function orbital_decay ( )

%*****************************************************************************80
%
%% orbital_decay() simulates the decay of an orbit to an equilibrium.
%
%  Discussion:
%
%    Given two massive bodies subject to gravity, it is possible to write down
%    differential equations describing their motion.  These equations are
%    simpler to formulate in the frame of reference in which the center of 
%    mass of the two bodies does not move.  If one body is much more massive 
%    than the other, then our calculations in this new frame are essentially
%    the same as in the original geometry.  This is the case when one body
%    is the sun, and another a planet.  
%
%    By giving an arbitrary initial condition for the system, we are likely
%    to have a planet whose orbit is not the classical ellipse we would
%    expect.  However, if we follow the orbit over time, we will see it
%    gradually decay to an ellipse.
%
%    This simulation would need to be modified if we wanted to consider
%    the behavior of two bodies of com-parable mass, and expected to see
%    them both moving, or, even in the sun-planet case, if we wanted to
%    allow the sun to have a velocity while we stayed in a fixed frame
%    of reference.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 April 2021
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'orbital_decay():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  This simulation follows a small body for 20 orbits\n' );
  fprintf ( 1, '  around a relatively massive body - such as Mercury around\n' );
  fprintf ( 1, '  the sun.\n' );
  fprintf ( 1, '  Kepler''s equations for a two body system are used.\n' );
  fprintf ( 1, '  Initially, the orbit is NOT an ellipse, but as time passes,\n' );
  fprintf ( 1, '  the orbit decays into an elliptical shape.\n' );

  [ t0, y0, tstop ] = two_body_parameters ( );
  tstop = 20.0 * 3.895;
  [ t0, y0, tstop ] = two_body_parameters ( t0, y0, tstop );

  f = @two_body_deriv;
  tspan = [ t0, tstop ];

  [ t, y ] = ode45 ( f, tspan, y0 );
%
%  Plot the result.
%
  figure ( );
  clf ( );
  hold ( 'on' );
    plot ( y(:,1), y(:,3), 'b-' );
    plot ( 0.0, 0.0, 'r.', 'markersize', 40 );
  hold ( 'off' );

  grid ( 'on' );
  title ( 'After twenty orbits, we are almost elliptical' );
  axis ( 'equal' );

  filename = 'orbital_decay';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'orbital_decay():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end

