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 comparable 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:
%
%    04 November 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ORBITAL_DECAY\n' );
  fprintf ( 1, '  MATLAB version\n' );
  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' );

  x0 = [ 1.0, 0.0, 0.0, 0.8 ];
  tmin = 0.0;
  tmax = 20.0 * 3.895;

  [ t, xvec ] = ode45 ( @two_body_ode, [ tmin, tmax ], x0 );
%
%  Plot the result.
%
  figure ( );

  hold ( 'on' );

% tvec = linspace ( tmin, tmax, 1000 );
% xvec = deval ( sol, tvec );
  plot ( xvec(:,1), xvec(:,3), 'b-' );
  plot ( xvec(1,1), xvec(1,3), 'g.', 'MarkerSize', 40 )
  plot ( xvec(end,1), xvec(end,3), 'r.', 'MarkerSize', 40 )
  plot ( 0.0, 0.0, 'bo', 'MarkerSize', 40 );

  grid ( 'on' );

  title ( 'After twenty orbits, we are almost elliptical' );
  axis ( 'equal' );
  hold ( 'off' );

  filename = 'orbital_decay';
  print ( '-dpng', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ORBITAL_DECAY:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
