c vertical fall with friction - eq of motion: mdv/dt=-mg-bv
c SI units m~kg v ~m/s a~ m/s^2 b~kg/s
real m
data b, m, g/6.7E-2 ,0.2,9.8/
a(y1,y0)= -(g+(b/m)*(y1-y0)/dt)
vscale=((m*g)/b)
tscale=vscale/g
print*,'vterminal,tscale=',vscale,tscale
print*,' '
dt=tscale/100.
nstep=int((3.*tscale)/dt)
kprint=int(float(nstep)/60.)
kount=kprint
c initial conditions
y0=0.
v0=0.
y1=y0+v0*dt
t=0.
write(6,100)t, y0, v0
do 10 i=2,nstep
t=dt*float(i)
y2=2.*y1-y0+dt**2*a(y1,y0)
v=(y2-y1)/dt
if(i.eq.kount)then
kount=kount+kprint
print 100,t,y2,v
endif
y0=y1
y1=y2
10 continue
100 format(2x,'t,y,v=',3(4x,e11.4))
stop
end
vterminal,tscale= 29.2537308 2.98507452
t,y,v= 0.0000E+00 0.0000E+00 0.0000E+00
t,y,v= 0.1194E+00 -0.5205E-01 -0.8689E+00
t,y,v= 0.2388E+00 -0.2397E+00 -0.1987E+01
t,y,v= 0.3582E+00 -0.5576E+00 -0.3062E+01
t,y,v= 0.4776E+00 -0.1001E+01 -0.4094E+01
t,y,v= 0.5970E+00 -0.1564E+01 -0.5085E+01
t,y,v= 0.7164E+00 -0.2242E+01 -0.6038E+01
t,y,v= 0.8358E+00 -0.3032E+01 -0.6952E+01
t,y,v= 0.9552E+00 -0.3928E+01 -0.7831E+01
t,y,v= 0.1075E+01 -0.4926E+01 -0.8675E+01
t,y,v= 0.1194E+01 -0.6023E+01 -0.9486E+01
t,y,v= 0.1313E+01 -0.7214E+01 -0.1027E+02
t,y,v= 0.1433E+01 -0.8496E+01 -0.1101E+02
t,y,v= 0.1552E+01 -0.9865E+01 -0.1173E+02
t,y,v= 0.1672E+01 -0.1132E+02 -0.1242E+02
t,y,v= 0.1791E+01 -0.1285E+02 -0.1309E+02
t,y,v= 0.1910E+01 -0.1446E+02 -0.1372E+02
t,y,v= 0.2030E+01 -0.1615E+02 -0.1433E+02
t,y,v= 0.2149E+01 -0.1790E+02 -0.1492E+02
t,y,v= 0.2269E+01 -0.1972E+02 -0.1549E+02
t,y,v= 0.2388E+01 -0.2161E+02 -0.1603E+02
t,y,v= 0.2507E+01 -0.2357E+02 -0.1655E+02
t,y,v= 0.2627E+01 -0.2558E+02 -0.1705E+02
t,y,v= 0.2746E+01 -0.2765E+02 -0.1753E+02
t,y,v= 0.2866E+01 -0.2978E+02 -0.1799E+02
t,y,v= 0.2985E+01 -0.3196E+02 -0.1844E+02
t,y,v= 0.3104E+01 -0.3420E+02 -0.1886E+02
t,y,v= 0.3224E+01 -0.3648E+02 -0.1927E+02
t,y,v= 0.3343E+01 -0.3881E+02 -0.1967E+02
t,y,v= 0.3463E+01 -0.4119E+02 -0.2004E+02
t,y,v= 0.3582E+01 -0.4361E+02 -0.2041E+02
t,y,v= 0.3701E+01 -0.4607E+02 -0.2076E+02
t,y,v= 0.3821E+01 -0.4857E+02 -0.2109E+02
t,y,v= 0.3940E+01 -0.5112E+02 -0.2141E+02
t,y,v= 0.4060E+01 -0.5370E+02 -0.2172E+02
t,y,v= 0.4179E+01 -0.5631E+02 -0.2202E+02
t,y,v= 0.4299E+01 -0.5896E+02 -0.2230E+02
t,y,v= 0.4418E+01 -0.6165E+02 -0.2258E+02
t,y,v= 0.4537E+01 -0.6436E+02 -0.2284E+02
t,y,v= 0.4657E+01 -0.6711E+02 -0.2309E+02
t,y,v= 0.4776E+01 -0.6988E+02 -0.2334E+02
t,y,v= 0.4896E+01 -0.7269E+02 -0.2357E+02
t,y,v= 0.5015E+01 -0.7552E+02 -0.2379E+02
t,y,v= 0.5134E+01 -0.7838E+02 -0.2401E+02
t,y,v= 0.5254E+01 -0.8126E+02 -0.2421E+02
t,y,v= 0.5373E+01 -0.8416E+02 -0.2441E+02
t,y,v= 0.5493E+01 -0.8709E+02 -0.2460E+02
t,y,v= 0.5612E+01 -0.9004E+02 -0.2479E+02
t,y,v= 0.5731E+01 -0.9302E+02 -0.2496E+02
t,y,v= 0.5851E+01 -0.9601E+02 -0.2513E+02
t,y,v= 0.5970E+01 -0.9902E+02 -0.2529E+02
t,y,v= 0.6090E+01 -0.1021E+03 -0.2545E+02
t,y,v= 0.6209E+01 -0.1051E+03 -0.2560E+02
t,y,v= 0.6328E+01 -0.1082E+03 -0.2574E+02
t,y,v= 0.6448E+01 -0.1113E+03 -0.2588E+02
t,y,v= 0.6567E+01 -0.1144E+03 -0.2602E+02
t,y,v= 0.6687E+01 -0.1175E+03 -0.2614E+02
t,y,v= 0.6806E+01 -0.1206E+03 -0.2627E+02
t,y,v= 0.6925E+01 -0.1237E+03 -0.2638E+02
t,y,v= 0.7045E+01 -0.1269E+03 -0.2650E+02
t,y,v= 0.7164E+01 -0.1301E+03 -0.2660E+02
t,y,v= 0.7284E+01 -0.1333E+03 -0.2671E+02
t,y,v= 0.7403E+01 -0.1365E+03 -0.2681E+02
t,y,v= 0.7522E+01 -0.1397E+03 -0.2691E+02
t,y,v= 0.7642E+01 -0.1429E+03 -0.2700E+02
t,y,v= 0.7761E+01 -0.1461E+03 -0.2709E+02
t,y,v= 0.7881E+01 -0.1494E+03 -0.2717E+02
t,y,v= 0.8000E+01 -0.1526E+03 -0.2725E+02
t,y,v= 0.8119E+01 -0.1559E+03 -0.2733E+02
t,y,v= 0.8239E+01 -0.1591E+03 -0.2741E+02
t,y,v= 0.8358E+01 -0.1624E+03 -0.2748E+02
t,y,v= 0.8478E+01 -0.1657E+03 -0.2755E+02
t,y,v= 0.8597E+01 -0.1690E+03 -0.2762E+02
t,y,v= 0.8716E+01 -0.1723E+03 -0.2768E+02
t,y,v= 0.8836E+01 -0.1756E+03 -0.2774E+02
