program cannon
      implicit none
      real angle,angrad,range,dist,c(3)
      real :: vel=300.0
      data c /0., 6.76e-5,-6.76e-9/
c
c     Willie Stevenson   1/96
c     Modified by J. Mahaffy 2/97
c
c     angle - cannon inclination angle from horizontal in degrees
c     angrad - angle converted to radians
c     range -  range of the cannon shell in meters
c     vel   -  muzzle velocity of the cannon  (m/s)
c     c     -  array containing coefficients for the cubic fit
c              to relative ground eleveation:
c              h(x) = c(1)*x + c(2)*x**2 + c(3)*x**3
c
c     dist  -  function that calculates distance traveled by the
c              cannon shell given the inclination angle in radians
c              and muzzle velocity in meters per second as arguments
c              Assumes flat terrain and no air resistance.
c
c     input - subroutine to obtain an inclination angle in degrees
c             from the user and return both this angle and its value
c             in radians through the argument list
c
c     output - subroutine to print the shell range for the input
c              angle of inclination
c
c
      call input ( angle,angrad)
c
      range=dist(angrad,vel,c)
c
      call output (angle,range)
c
      stop
      end
c
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c
      subroutine input ( angle,angrad)
c     
c     Obtain inclination angle from user, Read it and change to radians
c                                          
      implicit none 
      real angle,angrad,pi
c
c     Willie Stevenson   1/96
c     
c     angle - cannon inclination angle in degrees
c     angrad - inclination angle changed to radians
c     pi - constant; 3.1415926536
c
      parameter (pi=3.1415926536) 
      print *, '    '
      print *, '    '
      print *, ' Provide the cannon inclination angle in degrees:'
      read *, angle
      angrad = angle*(pi/180)
      return
      end
c
      subroutine output (angle,range)
c
c     Print shell range for the input angle of inclination
c
c     Willie Stevenson   1/96
c
      implicit none 
      real angle,range
      print *, ' For inclination of ',angle,' degrees'
      print *, ' The shell range is ', range, ' meters'
      return
      end
c
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c   Insert your function (or functions) below this line
c
      function dist(angrad,v0,c)
      implicit none
      real angrad,range,v0,g,dist,eps,c(3),fx0,x0,x1,fx1,x,fx,x2,dx
      integer it, itmax
      parameter (g=9.807,eps=1.e-5,itmax=10)
c
c     J. Mahaffy 2/97
c
c     angrad - angle of inclination in radians
c     v0 - initial velocity of cannonball in meters per second
c     g - gravitational constant in meters per second squared
c     c     -  array containing coefficients for the cubic fit
c              to relative ground eleveation:
c              h(x) = c(1)*x + c(2)*x**2 + c(3)*x**3
c
c     dist - name of function, distance traveled in meters
c
c     Note that I defined g as a parameter, since its value will
c     never be changed.
c
c     NOW CALCULATE THE DISTANCE of TRAVEL
c
      x0 = (v0**2*sin(2*angrad))/g
      x1=0.999*x0
      do 100 it=1,itmax
         fx0=  c(1)*x0 + c(2)*x0**2 + c(3)*x0**3
     &       + .5*g*(x0/(v0*cos(angrad)))**2-tan(angrad)*x0
         fx1=  c(1)*x1 + c(2)*x1**2 + c(3)*x1**3
     &       + .5*g*(x1/(v0*cos(angrad)))**2-tan(angrad)*x1
         dx=fx1*(x1-x0)/(fx1-fx0)
         x0=x1
         x1=x1-dx
         write (*,2000)it,x1,fx
         if(abs(x1-x0).lt.eps*abs(x1)) go to 200
  100    continue
      print *, 'Iteration failed to converge'
  200 dist = x1
 2000 format(' it=',i2,', x= ',f7.1, ', fx = ',1p,e12.5)
    
      return
      end