function inclass_grav(timestep, duration)
%
% A very simple illustration of the dangers of numerical integration
% of gravity. We follow the motion of a very massive object
% (fixed in place) and a tiny object (free to move). They start
% at rest. The tiny object falls towards the big one.
% We follow the motion for "duration" seconds with the given "timestep"
%
% Arguments: timestep (input) size of timestep (seconds)
%
% duration (input) number of seconds to
% run the simulation
%
% MWR 4/28/2003
% mass of each object (kg)
m1 = 1e9;
m2 = 1;
% initial position of each object (m). Note that
% the massive object is fixed in place and never moves
x1 = 100;
x2 = 0;
% initial velocity of each object (m/s). Again, note that the
% massive object is fixed in place and can never move
v1 = 0.0;
v2 = 0.0;
% gravitational constant, in MKS units
G = 6.67E-11;
% initial total energy: GPE + KE
dist = x2 - x1;
initial_energy = -G*m1*m2/abs(dist) + 0.5*m2*v2*v2;
fprintf(1, 'initial energy is %9.4e\n', initial_energy);
t = 0;
while (t < duration)
% calculate the distance between the object here ...
dist =
% calculate current total energy
current_energy =
% now find the gravitational force between the objects ...
force =
% and the acceleration of the tiny object
accel2 =
% use simple Euler's method to calculate new velocity and position
new_v2 =
new_x2 =
% print out information on the CURRENT (not new) quantities
fprintf(1, ' %9.4f %10.5f %10.5f %9.4e \n', t, x2, v2, current_energy);
% update variables for next time through the loop
v2 = new_v2;
x2 = new_x2;
t = t + timestep;
end
% end of program