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),fx,dfdx,xo,x,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
      x = (v0**2*sin(2*angrad))/g
      do 100 it=1,itmax
         xo=x
         dfdx =  c(1) + 2*c(2)*x + 3*c(3)*x**2
     &          + g*x/(v0*cos(angrad))**2-tan(angrad)

         fx=  c(1)*x + c(2)*x**2 + c(3)*x**3
     &       + .5*g*(x/(v0*cos(angrad)))**2-tan(angrad)*x
         dx = -fx/dfdx
         x=xo+dx
         write (*,2000)it,x,fx,dx 
         if(abs(x-xo).lt.eps*abs(x)) go to 200
  100    continue
      print *, 'Iteration failed to converge'
  200 dist = x
 2000 format(' it=',i2,', x= ',f7.1, ', fx = ',1p,e12.5,
     &       ', dx = ',e12.5)
      return
      end