program cannon_drift
      implicit none
      real x,vx,ax, y,vy,ay, vx_ini,vy_ini, x_range,y_range,delta_t
      real mass,k
      integer pgopen,isymbol,i_again

      call input_n_init(vx_ini,vy_ini,delta_t,mass,k)

      if (pgopen('/xwin').le.0) stop
      call pgpap(5.0, 1.0)

      isymbol = 20
c     isymbol = 21

 100  write(*,*) 'How far (x-dir) and how high (y-dir) is the plot ?'
      read (*,*) x_range, y_range
      call pgeras
      call pgenv(0.0,x_range,0.0,y_range,0,0)

      x = 0.0
      y = 1.0
      vx = vx_ini
      vy = vy_ini

      do while (y.gt.0.0)
        call fly_drift (x,vx,ax,y,vy,ay,delta_t,mass,k)
        call pgpt (1,x,y,isymbol)
      end do

 102  write(*,*) 'Again ? (1=Yes;0=No)'
      read(*,*) i_again
      if (i_again.eq.1) goto 100
      if (i_again.eq.0) stop
      goto 102

      call pgclos

      end



      subroutine input_n_init(vx,vy,delta_t,mass,k)
      implicit none
      real vx,vy,v,angle,ke,delta_t,mass,k
      integer i_gun_powder
      write(*,*) 'How maqny packs of gun powder (integer) ?'
      read (*,*) i_gun_powder
      write (*,*) 'What is the angle of shooting ?'
      read (*,*) angle
      write(*,*) 'What is the small time-interval delta_t ?'
      read (*,*) delta_t
      write(*,*) 'What is the mass of the particle ?'
      read (*,*) mass
      write(*,*) 'What is the coefficient of air resistence ?'
      read (*,*) k
      angle = angle * 3.14159268/180
      ke = 1.0*i_gun_powder
      v = sqrt(2.0*ke/mass)
      vx = v*cos(angle)
      vy = v*sin(angle)
      end



      subroutine fly_drift(x,vx,ax,y,vy,ay,delta_t,mass,k)
      implicit none
      real x,vx,ax,y,vy,ay,delta_t,theda,g
      real v_old,v_old_x,v_new_x,v_old_y,v_new_y,mass,v_mid_x
      real a_mid_x,a_mid_y,a_drift_x,a_drift_y,a_drift,k,v_mid_y

      g = -9.8

      v_old_x = vx
      v_old_y = vy

      v_old = sqrt(v_old_x**2 + v_old_y**2)
      a_drift = -1*k*v_old/mass
      theda = atan(v_old_y/v_old_x)
      a_drift_x = a_drift*cos(theda)
      a_drift_y = a_drift*sin(theda)
      a_mid_x = a_drift_x
      a_mid_y = a_drift_y + g
      v_mid_x = v_old_x + a_mid_x*0.5*delta_t
      v_mid_y = v_old_y + a_mid_y*0.5*delta_t

      v_new_x = v_old_x + a_mid_x*delta_t
      x = x + v_mid_x*delta_t

      v_new_y = v_old_y + a_mid_y*delta_t
      y = y + v_mid_y*delta_t

      vx = v_new_x
      vy = v_new_y

      end