! simulation of velocity vs. time for a large cannon
! include the effect of air resistance 
! Program to accompany "Computational Physics" by N. Giordano
! Copyright Prentice Hall 1997
program cannon
   option nolet
   library "sglib*","sgfunc*"
   dim x(5000), y(5000)
   call initialize(dt,v_init,theta,A_m)
   call calculate(x,y,dt,v_init,theta,A_m)
   call display(x,y)
end
sub initialize(dt,v_init,theta,A_m)
   dt = 0.25	     ! seconds - all unit mks
   v_init = 700 ! m/s
   theta = 55   ! firing angle in degrees
   A_m = 4e-5   ! A2/m
end sub
! x() and y() are the position of the projectile
! dt = time step   v_init = initial speed   theta = launch angle
! A_m = proportional to drag force = A2/m
sub calculate(x(),y(),dt,v_init,theta,A_m)
   option angle degrees    ! use degrees rather than radians
   x(1) = 0
   y(1) = 0
   vx = v_init * cos(theta)
   vy = v_init * sin(theta)
   nmax = size(x)    ! this is the number of elements in the array x()
   for i = 2 to nmax
      x(i) = x(i-1) + vx * dt            ! Euler method equations
      y(i) = y(i-1) + vy * dt
      f = A_m * sqr(vx^2 + vy^2)         ! drag force from air resistance
      vy = vy - 9.8 * dt - f * vy * dt
      vx = vx - f * vx * dt
      if y(i) <= 0 then exit for         ! shell has hit the ground
   next i
   a = -y(i) / y(i-1)              ! interpolate to find landing point
   x(i) = (x(i) + a*x(i-1)) / (1+a)
   y(i) = 0
   mat redim x(i),y(i)
end sub
sub display(x(),y())     ! display the results with datagraph
   call settitle("Cannon shell y vs. x")
   call sethlabel("x (m)")
   call setvlabel("y (m)")
   call datagraph(x,y,0,1,"black")
end sub