Program projectile
c         author Amit Lath 
c         date 9/19/07
cccccccccccccccccccccccccccc
c  This program will numerically solve the
c  coupled differential equation for a
c  projectile with quadratic and linear 
c  drag coefficients.
c  The variables are "b" for linear drag
c  and "c" for quad drag.  These are calculated
c  from "beta" and "gamma" as given in Taylor
c  "Classical Mechanics" 2nd printing, eq. 2.4
c   - 2.6.
ccccccccccccccccccccccccccccccccccc
           implicit none
c
c  Define variables for x, y, xdot, xdoubledot
c   as well as ydot, ydoubledot
cc

c                  The original x,y,vx,vy
           real*8  x0, y0, v0, vx0, vy0
c
c                  The angle of the gun (taken from x axis)
           real*8  theta
c
c                  The iteration variable (x, x+dx, etc)        
           real*8  x_i, y_i, x_j, y_j
           real*8  vx_i, vy_i, vx_j, vy_j
c
c                  The various drag variables from Taylor
           real*8  c, b, beta, gamma, diameter
c
c                  time derivatives of x, y 
           real*8  xd, xdd, yd, ydd           
c
c                  Time:  final, tau 
           real*8  t, dt, t_final, tau_b, tau_c
c
c                  Gravitational constant
           real*8  g
c
c                  Mass of projectile
           real*8  mass
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c   Now we need to set the various constants
c    to do the problem.  We will define the x0, y0 = 0
c    for now.  Change them if needed
c    We will define v0, and calculate vx0, vy0
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc

cc    INITIAL POSITION,  CHANGE HERE (in meters)
       x0=0.0
       y0=0.0
     
cc     INITIAL VELOCITY, CHANGE HERE (in meters/sec)
       v0 = 30.0
c
cc     ANGLE OF LAUNCH,  CHANGE HERE (in degrees)  
       theta = 50.0

c       DIAMETER, MASS OF PROJECTILE, CHANGE HERE (in meters, kg)
        diameter = 0.07
        mass     = 0.15
c
c       Gravitational constant, CHANGE HERE (in m/s^2)
        g=9.81
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

       vx0 = v0*cos(theta*3.1415927/180.0)
       vy0 = v0*sin(theta*3.1415927/180.0)
c
c      Now we need to calculate the coefficients of drag
c       From Taylor...  (beta and gamma in SI units)

        beta = 1.6e-4
        gamma = 0.25

        b = beta*diameter
        c = gamma*diameter*diameter


c
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Now we have all the variables we need to get started
c     doing the calculation.   We need to figure out how long
c     (in time) to do the calculation.  We find the "tau" for
c     the linear and quad cases.

c       tau_b = mass/b
       tau_c = mass/(c*vx0)        

c      Take t_final to be "non drag" time + a little bit
c       t_final = 2.5*vy0/g
        t_final = 8.0
c
c      Take 5000 time steps from 0 to t_final
       dt = t_final/5000

       t=0.0
       x_i = x0
       y_i = y0
       vx_i = vx0
       vy_i = vy0

       do while (t.le.t_final)

c
c         Turn the differntial equations  (2.19, 2.25, 2.61)
c         into difference equations.  NOTE: this is not a
c         very efficient or accurate way to do this!  There are
c         better ways (Runge-Kutte, etc).
c

          vx_j = dt/mass*
     .           ((-c*sqrt(vx_i**2 + vy_i**2)*vx_i) +
     .            (-b*vx_i)) +
     .            vx_i       

          vy_j = -g*dt +
     .           dt/mass*
     .           ((-c*sqrt(vx_i**2 + vy_i**2)*vy_i) +
     .            (-b*vy_i)) +
     .            vy_i         
         
c         Now find x, y from the integral of vx, vy.

           x_j = x_i + vx_i*dt
           y_j = y_i + vy_i*dt

c         Now write out the variables

           write(6,100) t, x_i, y_i, vx_i, vy_i

cc  Comment out the above line, and uncomment the two lines
cc   below if you want to output the "no drag" values
cc

c           write(6,100) t, x0+vx0*t, y0+vy0*t-0.5*g*t*t, 
c     .                                     vx0, vy0-g*t


 100       format(f9.3,1x,f9.3,1x,f9.3,1x,f9.3,1x,f9.3)
c          Now go on to next step
           t = t+dt
           x_i = x_j
           y_i = y_j
           vx_i = vx_j
           vy_i = vy_j

       end do  ! t .le. t_final

       end