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