c c c
      program fall
      implicit none
c
c   Program to calculate the dynamics of a falling body
c
c   Arrays:
c     v   -   velocity at each time integration step
c     z   -   height at  each time integration step
c     t   -   time for each corresponding v and z
c     zreal - Actual height at time t(i) for comparison with computed z
c
c    John Mahaffy  4/15/95
c
      integer np
      parameter (np=2000)
c
c   In this program, I am using blank common just to save space in the
c   executable program file (a.out).  It only appears in the main program
c   Named common "const" communicates between subroutines.  By using the
c   "common/const/..." statement in the main routine, I insure that some
c   machines do not lose values in that common block when the first
c   subroutine containing "common/const/..." is exited.
c
c
      common v(np),z(np),t(np), zreal(np)
c
      real v,z,t,zreal
      common/const/dtmin,g
      real dtmin,g
      real dt
      integer nsteps
c
      call input(z,dt)
      call odesolve(v,z,t,dt,np,nsteps)
      call realans(t,z,nsteps,zreal)
      call output (t,z,zreal,v,nsteps)
      stop
      end
c
c
      block data cons
c
c
c    Initialize various constants
c
c    John Mahaffy  4/15/95
c
c     g  -  gravitational acceleration
c     dtmin - time step selected if 0.0 is entered
c
      common/const/dtmin,g
      real dtmin,g
      data dtmin,g/.001,9.807/
c
      end block data cons
c
      subroutine input (z,dt)
c
c   Obtain user input for initial height and solution time step
c
c    John Mahaffy  4/15/95
c
c  Output Arguments:
c     z(1)   -  initial height
c     dt     -  integration time step
c
      real z(*)
      common/const/dtmin,g
      real dtmin,g
c
      write(6,*) ' What is the initial height (meters)?'
      read *, z(1)
      write(6,*) ' What time step size do you want (seconds)?'
      read *, dt
      if(dt.le.0.) dt=dtmin
      return
      end
c
      subroutine odesolve(v,z,t,dt,np,nsteps)
c
c   Solve the Ordinary Differential Equation of motion for the body
c
c    John Mahaffy  4/15/95
c
c   Arguments:
c   Input
c     dt     -   timestep size
c     np     -   size of v,z, and t arrays
c   Output:
c     v    -   velocity
c     z    -   height
c     t    -   time
c     nsteps - last step in the integration
c
      real v(*),z(*),t(*),dt,rdt
      common/const/dtmin,g
      real dtmin,g
c
c   Solve the equation for a falling body
c
c      2
c     d z
c     ---2   =  - g
c     d t
c
c    NOTE:  When you try this you will see that solving this for of the
c           equation will result in the amplification of machine round-off
c           error.  Never use this approach for this type equation.
c
c   Set remaining initial conditions:
c
      t(1)=0.
      v(1)=0.
c
c    For the first step we need to fold in the initial velocity condition
c    v(1)=(z(2)-z(0))/(2*dt)  giving z(0)=z(1) - 2*dt*v(1)
c    Substituting this value of z(0) into the equation at t=0 gives:
c
      z(2) = z(1) - .5*g*dt**2
      t(2) = dt
c
c     Now loop through time steps until z goes negative or we run out of space
c
      rdt=.5/dt
      do 100 i=2,np-1
         z(i+1)= 2*z(i)-z(i-1)-dt**2*g
         v(i)=(z(i+1) - z(i-1)) * rdt
         t(i+1)=t(i)+dt
         if(z(i).lt.0.) go to 200
 100     continue
      write(6,*) 'Ran out of space to continue integration'
      write(6,*) ' Last height was ',z(np),' meters'
      i=np
 200  nsteps=i
c     return
      end
c
      subroutine realans(t,z,nsteps,zreal)
c
c     Get the values of the analytic solution to the differential equation
c     for each time point to check the numerical accuracy.
c
c    John Mahaffy  4/15/95
c
      common/const/dtmin,g
      real dtmin,g
      real t(*),z(*),zreal(*)
c
      do 10 i=1,nsteps
  10  zreal(i)=z(1)-.5*g*t(i)**2
      return
      end
c
      subroutine output(t,z,zreal,v,nsteps)
      implicit none
c
c   Outputs the full results of the time integration
c
c    John Mahaffy  4/15/95
c
      real v(*),z(*),t(*), zreal(*)
      integer nsteps,i
      print *, 'Look in the file fall.output for detailed results'
      open (11,file='fall.output')
      write(11,2000)
      do 300 i=1,nsteps
      write(11,2001) t(i),v(i),z(i),zreal(i)
 300  continue
 2000 format (33x,'Computed',8x,'Actual',/,
     &        6x,'Time',9x,'Velocity', 8x,'Height',8x,'Height')
 2001 format (1x,1p,4e15.7)
      return
c
      end subroutine output
c
c
c c