The met_rocket simulation was written to estimate a launch angle which would deploy a parachute camera above the Overberg Test Range. It was written in extended FORTRAN 77 for the g77 compiler. Fortunately, the fpt engineering tool can process code in this style.
The entire program is only 200 lines long.
C ***************************************************************************** C met_rocket.for 13-Jul-06 Stephen Day Mandy et al C 18-Oct-06 John Collins C C ***************************************************************************** C 18-Oct-06 John Collins Removed the (very few) F90 features C ***************************************************************************** C PROGRAM met_rocket C IMPLICIT NONE C C ***************************************************************************** C REAL 1 g 1 ,p_burn_time 1 ,p_burn_thrust 1 ,p_initial_mass 1 ,p_stage_2_mass 1 ,p_fuel_mass 1 ,p_launch_vel 1 ,frametime PARAMETER ( 1 g = 32.2D0 ! fpsps 1 ,p_burn_time = 2.11D0 ! secs 1 ,p_burn_thrust = 4021.0D0*g ! lb-ft 1 ,p_initial_mass = 67.2D0 ! lb 1 ,p_stage_2_mass = 10.0D0 ! lb 1 ,p_fuel_mass = 37.6D0 ! lb 1 ,p_launch_vel = 190.0 ! fps 1 ,frametime = 0.001D0 ! secs 1 ) C INTEGER 1 frame 1 ,monitor_interval 1 ,stage ! Stage, before/after separation 1 ,u ! Logical unit number C REAL 1 dtr 1 ,time 1 ,xdd ! Down range acceleration, fpsps 1 ,xd ! Down range velocity, f/s 1 ,x ! Down range distance, ft 1 ,pxdd ! Last frame accelleration 1 ,pxd ! Last frame velocity 1 ,hdd ! Height acceleration, fpsps 1 ,hd ! Height velocity, fps 1 ,h ! Height, ft 1 ,phdd ! Last frame height acceleration 1 ,phd ! Last frame height velocity 1 ,theta ! Angle of flight to ground, radians 1 ,thetad ! Launch angle to ground, degrees 1 ,thrust 1 ,mass ! lb 1 ,drag 1 ,vel ! Direction of flight 1 ,atm_pressure ! Atmospheres 1 ,c_drag(2,2) ! Drag coeff, linear/quad, stage 1/2 C C ***************************************************************************** C OPEN (UNIT = 7,FILE = 'met_rocket.txt',STATUS = 'UNKNOWN') DO u = 6,7 WRITE (u,'(/,''Met Rocket Simulation'')') WRITE (u,'(''====================='',/)') ENDDO WRITE (6,'(''The program writes the output data to the file '' 1 ,''met_rocket.txt in the current'',/,''directory. '' 1 ,''Edit the program as marked to get interactice '' 1 ,''control of all of'',/,''the parameters'',/)') C C Initial conditions dtr = ATAN(1.0D0)/45.0D0 frame = 0 time = 0.0D0 hdd = 0.0D0 hd = 0.0D0 h = 0.0D0 xdd = 0.0D0 xd = 0.0D0 x = 0.0D0 phdd = 0.0D0 phd = 0.0D0 pxdd = 0.0D0 pxd = 0.0D0 C C Un-comment lines beginning CCC to control drag parameters C ========================================================= C WRITE (6,'(''Launch angle (degrees): '',$)') READ (5,*)thetad theta = thetad*dtr C c_drag(1,1) = 0.0D0 c_drag(2,1) = 0.0008D0 c_drag(1,2) = 0.0D0 c_drag(2,2) = 0.00021D0 monitor_interval = 1000 C CCC WRITE(6,'(''Drag stage 1 linear coeff: '',$)') CCC READ(5,*)c_drag(1,1) CCC WRITE(6,'(''Drag stage 1 quadratic coeff: '',$)') CCC READ(5,*)c_drag(2,1) CCC WRITE(6,'(''Drag stage 2 linear coeff: '',$)') CCC READ(5,*)c_drag(1,2) CCC WRITE(6,'(''Drag stage 2 quadratic coeff: '',$)') CCC READ(5,*)c_drag(2,2) CCC WRITE(6,'(''Monitor interval: '',$)') CCC READ(5,*)monitor_interval C WRITE (7,'(''Launch angle:'',T32,F12.2)')thetad WRITE (7,'(''Drag coefficients'')') WRITE (7,'('' Stage 1 linear:'',T32,F12.6)')c_drag(1,1) WRITE (7,'('' Stage 1 quadratic:'',T32,F12.6)')c_drag(2,1) WRITE (7,'('' Stage 2 linear:'',T32,F12.6)')c_drag(1,2) WRITE (7,'('' Stage 2 quadratic:'',T32,F12.6,/)')c_drag(2,2) C C ***************************************************************************** C DO WHILE (theta .GE. 0.0D0) frame = frame+1 C IF (time .LE. p_burn_time) THEN stage = 1 mass = p_initial_mass- 1 MAX(p_fuel_mass*(time/p_burn_time),0.0) thrust = p_burn_thrust ELSE stage = 2 mass = p_stage_2_mass thrust = 0.0D0 ENDIF C C Pressure in atmospheres atm_pressure = 10.0D0**(-h/50850.0D0) C vel = SQRT(xd*xd+hd*hd) C Assume drag is linear with pressure and quadratic with velocity drag = atm_pressure* 1 (c_drag(2,stage)*vel**2+c_drag(1,stage)*vel) C xdd = ((thrust-drag)/mass)*COS(theta) hdd = ((thrust-drag)/mass)*SIN(theta)-g C C Trapezoidal integration - AB2 xd = xd+frametime*(2*xdd-pxdd) hd = hd+frametime*(2*hdd-phdd) pxdd = xdd phdd = hdd x = x+frametime*(2*xd-pxd) h = h+frametime*(2*hd-phd) pxd = xd phd = hd IF (vel .GT. p_launch_vel) THEN theta = ATAN2(hd,xd) ENDIF time = time+frametime C IF (frame .EQ. monitor_interval*(frame/monitor_interval)) THEN DO u = 6,7 WRITE (u,'(/,''Frame:'',I8,'' Time:'',F10.2 1 ,'' Theta:'',F10.2,'' Pressure:'',F10.4)') 1 frame,time,theta/dtr,atm_pressure WRITE (u,'('' h '' 1 ,'' hd '' 1 ,'' hdd '' 1 ,'' x '' 1 ,'' xd '' 1 ,'' xdd '' 1 )') WRITE (u,'(6F13.5)')h,hd,hdd,x,xd,xdd ENDDO CCC WRITE (6,'(''<CR> to continue'')') CCC READ (5,*) ENDIF C ENDDO C C Final values DO u = 6,7 WRITE (u,'(//,''Maximum altitude'')') WRITE (u,'(/,''Frame:'',I8,'' Time:'',F10.2 1 ,'' Theta:'',F10.2,'' Pressure:'',F10.4)') 1 frame,time,theta/dtr,atm_pressure WRITE (u,'('' h '' 1 ,'' hd '' 1 ,'' hdd '' 1 ,'' x '' 1 ,'' xd '' 1 ,'' xdd '' 1 )') WRITE (u,'(6F13.5)')h,hd,hdd,x,xd,xdd WRITE (u,'(/,''End of run'')') ENDDO CLOSE (7) C C ***************************************************************************** C END ! *****************************************************************