PROGRAM drag
! assume drag force proportional to v^2
CALL initial(y,y0,v,t,g,vt2,dt,dt_2)
CALL print_table(y,y0,t,dt,nshow)
DO until t >= 0
   CALL Euler_Richardson(y,v,t,g,vt2,dt,dt_2)
LOOP
CALL print_table(y,y0,t,dt,nshow)      ! values at t = 0
LET counter = 0
DO while t <= 0.80
   CALL Euler_Richardson(y,v,t,g,vt2,dt,dt_2)
   LET counter = counter + 1
   IF mod(counter,nshow) = 0 then
      CALL print_table(y,y0,t,dt,nshow)
   END IF
LOOP
END

SUB initial(y,y0,v,t0,g,vt2,dt,dt_2)
    LET t0 = -0.132               ! initial time (see Table 3.1)
    LET y0 = 0
    LET y = y0
    LET v = 0                     ! initial velocity
    INPUT prompt "terminal velocity (m/s) = ": vt
    LET vt2 = vt*vt
    LET dt = 0.001
    LET dt_2 = 0.5*dt
    LET g = 9.8
END SUB

SUB Euler_Richardson(y,v,t,g,vt2,dt,dt_2)   ! Euler-Richardson method
    LET v2 = v*v
    LET a = -g*(1 - v2/vt2)
    LET vmid = v + a*dt_2         ! velocity at midpoint
    LET ymid = y + v*dt_2         ! position at midpoint
    LET vmid2 = vmid*vmid
    LET amid = -g*(1 - vmid2/vt2)      ! acceleration at midpoint
    LET v = v + amid*dt
    LET y = y + vmid*dt
    LET t = t + dt
END SUB

SUB print_table(y,y0,t,dt,nshow)
    IF t < 0 then
       LET show_time = 0.1        ! time interval between output
       ! choice of dt = 0.001 convenient
       LET nshow = int(show_time/dt)
       CLEAR
       PRINT "time","displacement"
       PRINT
    END IF
    LET distance_fallen = y0 - y
    PRINT t,distance_fallen
END SUB