      program prog101
      implicit real*8(a-h,o-z)
c
      ti = 0.0d0
      dt = 0.01d0
      xi = 0.0d0 
      yi = 0.0d0
      theta = atan(1.0d0)
      vxi = 20.0d0 * sin(theta) 
      vyi = 20.0d0 * cos(theta)   
      gan = 1.0d0
      a = 0.1d0
      ga = gan * a**2  
c
      t = ti
      x = xi
      y = yi
      vx = vxi
      vy = vyi
c
      open(16,file='prog101.dat')
      write(*,800)x,y,t
      write(16,800)x,y,y
c      write(16,800)x,y,t
c      write(16,810)x,y
c      
 100  continue
c
      call onestep(x,y,vx,vy,dt,ga)
      t = t + dt
      write(*,800)x,y,t
      write(16,800)x,y,y
c      write(16,800)x,y,t
c      write(16,810)x,y
      
      if( y.ge.0.0d0 ) goto 100      
c
      stop
 800  format(3(1x,f15.7))      
 810  format(2(1x,f15.7))      
      end
c---------------------------------------------
      subroutine onestep(x,y,vx,vy,dt,ga) 
c---------------------------------------------
      implicit real*8(a-h,o-z)
c
c x1,y1での微係数vx1,vy1を求める
      call force(vx,vy,fx,fy,ga) 
c      x1 = x + vx * dt /2
c      y1 = y + vy * dt /2
      vx1 = vx + fx * dt /2
      vy1 = vy + fy * dt /2
c
c x1,y1での微係数を使いx,y,vx,vyを進める
      call force(vx1,vy1,fx1,fy1,ga)
      x = x + vx1*dt
      y = y + vy1*dt
      vx = vx + fx1*dt
      vy = vy + fy1*dt
c
      return
      end
c---------------------------------------------
      subroutine force(vx,vy,fx,fy,ga)
c---------------------------------------------
      implicit real*8(a-h,o-z)
      m = 1.0d0
      g = 9.8d0
      v =sqrt(vx**2 + vy**2) 
      fx = - ga * v * vx /m
      fy = - ga * v * vy /m - g
      return
      end

