function solar_system(ending_time, time_step)
%   Compute the positions and velocities of bodies in the solar system
%   over a period of time.  We use initial values for Jan 1, 2011,
%   as provided by JPL's Horizons ephemeris system.
%
%   The input arguments are
%
%          ending_time              is the duration of the simulation
%                                        (days)
%
%          timestep                   is the timestep used in the numerical
%                                        integration (seconds)
%
%   The program prints lines to the standard output (screen) with the 
%   following format (the example below is shown on two lines for
%   clarity, but it should be on a single line in the real output):
%
%     ID   5    t   0.000E00  x  0.000E00  y  9.000E07   z  9.023E08   
%                            vx  1.250E03  vy  0.000E00  vz  0.000E00
%
%   where 
%            ID         is followed by an integer identifying the body
%                           (0 = Sun, 1 = Mercury, etc.)
%
%             t         is followed by the time (seconds)
%
%             x         is followed by the X position of the body (meters)
%
%             y         is followed by the Y position of the body (meters)
%
%             z         is followed by the Z position of the body (meters)
%
%            vx         is followed by the X velocity of the body (meters/sec)
%
%            vy         is followed by the Y velocity of the body (meters/sec)
%
%            vz         is followed by the Z velocity of the body (meters/sec)
%

%
% initialize positions, velocities, masses


% main loop over time

   % loop over each body i
   for i = 1 : 10 

      % set total force on body i to zero in x, y, z

      % loop over each other body j to calc forces it exerts on body i
      for j = 1 : 10

         % compute forces exerted by body j on body i in x, y, z

         % add those forces to running total of forces on body i in x, y z

      end

      % now, given the forces on body i, compute acceleration in x, y, z

      % compute the NEW velocity of body i in x, y, z

      % compute the NEW position of body i in x, y, z

   end


   % now, after NEW positions and velocities have been computed for 
   %    all objects, using their current positions,
   %    we can safely replace the current positions and velocities 
   %    with the new values
   for i = 1 : 10 

      % update positions and velocities for body i in x, y, z

      % print out values for body i 
   end


   % compute energies for all objects, if desired
  
   % update time

end   


