      program rep08
c
c     Gravitational Fall Problem with Air Friction
c     subroutines:
c        onestep(x,y,vx,vy,t,dt,ga)
c        force(vx,vy,fx,fy,ga)
c        solvexf(v0,theta,ga,dt,xf,tf)
c        solvey0(x,y,vx,vy,t,ga)
c
      implicit real*8 (a-h,o-z)
c Default Parameters
      dt    = 0.02d0
      a     = 0.2d0
      v0    = 30.0d0
      gam   = 1.0d0
c
      open(16,file="rep08.dat")
c
c     write(*,800) "# v0, gam, dt = ?"
c     read(*,*) v0, gam, dt
      write( *,800) "# v0, gam, dt = ", v0, gam, dt
      write(16,800) "# v0, gam, dt = ", v0, gam, dt
 800  format(a, 3(1x,f10.2))
c
      ga = gam*a*a
      tfx = 0.0d0
      xfx = 0.0d0
      thx = 45.0d0
c
      do i = 1, 89
        theta = dble(i)
        call solvexf(v0,theta,ga,dt,xf,tf)
        write(16,100) theta, xf, tf
        if(xf .gt. xfx) then
           xfx  = xf
           tfx  = tf
           thx = theta
        endif
      enddo
c
        write(16,*)
        thx1=thx
      do i = -10, 10
        theta = thx1 + i * 0.1d0
        call solvexf(v0,theta,ga,dt,xf,tf)
        write(16,100) theta, xf, tf
        if(xf .gt. xfx) then
           xfx  = xf
           tfx  = tf
           thx = theta
        endif
      enddo
c
        write(16,*)
        thx2=thx
      do i = -10, 10
        theta = thx2 + i * 0.01d0
        call solvexf(v0,theta,ga,dt,xf,tf)
        write(16,100) theta, xf, tf
        if(xf .gt. xfx) then
           xfx  = xf
           tfx  = tf
           thx = theta
        endif
      enddo
c
        write(16,*)
      write(16,100) thx, xfx, tfx
      write(16,800) "# Theta, Xf = ", thx, xfx
      write( *,800) "# Theta, Xf = ", thx, xfx
c
      theta = nint(10*thx)/10.0d0
      write( *,801) "# Final Answer: theta = ", theta
      write(16,801) "# Final Answer: theta = ", theta
c
 100  format(4(1x,f15.7))
 801  format(a, 1x,f8.1)
      end
c **********************************************************************
      subroutine onestep(x,y,vx,vy,t,dt,ga)
c **********************************************************************
      implicit real*8 (a-h,o-z)
c
c Force at x,y
      call force(vx,vy,fx,fy,ga)

c Trial Propagation and Force at x,y
c     x1  = x  + vx*dt/2
c     y1  = y  + vy*dt/2
      vx1 = vx + fx*dt/2
      vy1 = vy + fy*dt/2
      call force(vx1,vy1,fx1,fy1,ga)
c
      dx =  vx1*dt   ! Modified Euler
      dy =  vy1*dt   ! Modified Euler
      dvx = fx1*dt   ! Modified Euler
      dvy = fy1*dt   ! Modified Euler

c Real Propagation in Modified Euler method
       x =  x + dx
       y =  y + dy
      vx = vx + dvx
      vy = vy + dvy
       t =  t + dt
c
      end
c **********************************************************************
      subroutine force(vx,vy,fx,fy,ga)
c **********************************************************************
      implicit real*8 (a-h,o-z)
      g=9.8d0
      v = sqrt(vx*vx+vy*vy)
      fx = - ga*vx*v
      fy = - ga*vy*v - g
      end
c **********************************************************************
      subroutine solvexf(v0,theta,ga,dt,xf,tf)
c **********************************************************************
      implicit real*8 (a-h,o-z)
c Constant
      pi    = 4.0d0*atan(1.0d0)
c
      t  = 0.0d0
      x  = 0.0d0
      y  = 0.0d0
      vx = v0*cos(theta*pi/180)
      vy = v0*sin(theta*pi/180)
c
 10   continue
        call onestep(x,y,vx,vy,t,dt,ga)
      if( y .ge. 0.0d0 ) goto 10
c
      call solvey0(x,y,vx,vy,t,ga)
      xf = x
      tf = t
c
      end
c **********************************************************************
      subroutine solvey0(x,y,vx,vy,t,ga)
c **********************************************************************
c   This subroutine gives the position where y = 0.
      implicit real*8 (a-h,o-z)
c
      eps = 1.0d-12
      do n = 1, 1000
        if(abs(y).lt.eps) goto 20
        dt = - y / vy
        call onestep(x,y,vx,vy,t,dt,ga)
      enddo
        write(*,*) " I cannot find the position where y = 0."
20    continue
c
      end
