function grav_example_euler1(initial_distance, timestep)
% Program to illustrate the danger of close approaches 
%    in gravitational computation.  We start with a ball
%    held a distance 'initial_distance' above a very 
%    massive object of zero size.  The ball falls freely
%    toward the massive object, eventually flying past it
%    with erroneously large velocity.
%
% User supplies the initial position of the ball and 
%    the value of the timestep for numerical integration.
%    We use simple Euler's method to compute the motion.
%
% Print out lines which include mechanical energy during each timestep.
%
% Usage: grav_example_euler1 (initial_distance, timestep)
%
%    where
%             initial_distance       is initial position of ball
%                                    above the massive object (m)
%
%             timestep               is size of steps in time (s)
%
% Returns: nothing
%
% MWR 12/12/2012


% initialize variables
%
%    mass of massive object (kg)
m_obj = 1.0E15;
%    mass of ball (kg)
m_ball = 1.0;
%    current distance from massive object (m)
y = initial_distance;
%    new distance from massive object (m)
y_new = 0;
%    current velocity (positive upward) (m/s)
vy = 0;
%    new velocity (m/s)
vy_new = 0;
%    time since start of simulation (s)
t = 0;
%    universal constant of gravitation (N*m^2/kg^2)
G = 6.67E-11;



% compute initial mechanical energy
ke = 0.5*m_ball*vy*vy;
gpe = -G*m_obj*m_ball/abs(y);
initial_energy = ke + gpe;
% fprintf(1, ' ke %9.3f  gpe %9.3f  tot %9.3f \n', ke, gpe, energy);

% loop over time; continue as long as ball is less than the 
%    initial distance from object; in other words, stop when
%    the ball has gone too far, indicating errors in calculations
%while (abs(y) <= initial_distance) 

% this works better for testing with specific input values
while (t < 100)
    
     %  compute mechanical energy
     ke = 0.5*m_ball*vy*vy;
     gpe = -G*m_obj*m_ball/abs(y);
     energy = ke + gpe;
     frac_error = (energy - initial_energy) / abs(initial_energy);
   
     %  compute accel at current position; take care with the sign
     accel = -G*m_obj*m_ball/(y*y);
     if (y < 0)
         accel = 0.0 - accel;
     end
     
     %  print current time, position velocity, accel, energy, frac error in energy
     fprintf(1, ' t %9.3f  y %9.3f  vy %9.3f  ay %9.3f  E %12.4e  f %11.4e \n', t, y, vy, accel, energy, frac_error);
     
     
     %  use const accel to compute new velocity
     vy_new = vy + accel*timestep;
     
     %  use current velocity to compute new position
     y_new = y + vy*timestep;

     %  update position and velocity for next time through loop
     vy = vy_new;
     y = y_new;
    
     %  update time
     t = t + timestep;

% end of loop over time
end

%  compute mechanical energy
ke = 0.5*m_ball*vy*vy;
gpe = -G*m_obj*m_ball/abs(y);
energy = ke + gpe;
frac_error = (energy - initial_energy) / abs(initial_energy);

%  compute accel at current position; take care with the sign
accel = -G*m_obj*m_ball/(y*y);
if (y < 0)     
    accel = 0.0 - accel;
end

%  print current time, position velocity, accel, energy, frac error in energy
fprintf(1, ' t %9.3f  y %9.3f  vy %9.3f  ay %9.3f  E %12.4e  f %11.4e \n', t, y, vy, accel, energy, frac_error);



% end of program
end
