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