DECLARE FUNCTION CL! (VEL!)
DECLARE FUNCTION CD! (VEL!)
MASS = 20000 / 32.2 'Mass of orbiter in slugs
COEFTEMP = .5 * 214 / MASS 'Temperary variable for calculating drag and lift
HS = 26000
CLS
OPEN "J-231_I.txt" FOR INPUT AS #1
OPEN "J_231_M.txt" FOR OUTPUT AS #2
GMAXMAX = 0
'DO
H.OLD = H1
T.OLD = T1
'"INPUT #1, T1, H1, V1
H1 = 301022
V1 = 4652
Hdot = 56 * 3.3
H1 = H1 * 3.3
V1 = V1 * 3.3
IF H1 < 100000 THEN GOTO NEXT.INPUT
GMAX = 0!
DTIME = 1 'Step size used in the numerical integration
TIME = 0 'TIME is the simulation clock, it is
'not really needed
VEL = V1 'Velocity in ft/sec, scaler quantity
DELAZ = 0!
RANGE = 60000
RANGE = RANGE * 6080 'Range to the site converted to feet
'Hdot = (H1 - H.OLD) / (T1 - T.OLD)
H = H1
DTIME = 2
X = -RANGE 'The X coordinate of the orbiter, in feet
Y = 0 'The Y coordinate of the orbiter, in feet
'The origin of X,Y is the landing sight.
'The orbiter is initially placed on the negative Y axis
VELHORZ = SQR(VEL * VEL - Hdot * Hdot) 'Horizontal velocity, ft/sec
VX = VELHORZ * COS(-DELAZ / 57.29) 'The X velocity of the orbiter, ft/sec
VY = VELHORZ * SIN(-DELAZ / 57.29) 'The Y velocity of the orbiter, ft/sec
DT2 = DTIME * DTIME / 2! 'Constant to used in Taylor's series expansion
DT3 = DT2 * DTIME / 3! 'Constant to used in Taylor's series expansion
TSTART = TIMER 'TIMER is the real time returned from the
'computer's clock, used to test the speed
'of the software.
'Start of the BIG loop:
GOSUB FINDDERV
DO
GOSUB PROPTIME
GOSUB FINDDERV
WRITE #2, T, G
LOOP UNTIL VEL < 50 OR H < 100 OR D > 1500 OR ABS(Hdot) > 50000 OR X > 0 OR END.FLAG! = 1 '***
'Stop loop when:
' VEL < 2500 (TAEM interface, Mach 2.5)
' Drag > 150 (the sim probably blew-up)
' Hdot > 5000 (the sim probably blew-up)
' X > 0 (Flew past the site)
' END.FLAG! is set (time step has been adjusted to stop at Mach 2.5)
END2:
NEXT.INPUT:
IF GMAX > GMAXMAX THEN GMAXMAX = GMAX
PRINT T1, H1, V1
WRITE #2, T, G
PRINT "GMAX = ", GMAX
PRINT "GMAXMAX = ", GMAXMAX
'IF GMAX > 4 THEN END
LOOP
GOTO END1
FINDDERV:
' This routine finds the derivatives to use in the Taylor's
' Series expansions.
RHO = .0028964 * EXP(-H / HS) 'Density of the air in slugs/ft^3
RHODOT = -RHO * Hdot / HS! 'Derivative of RHO with time
CDTEMP = CD(VEL) 'Drag Coef as a function of VEL
CLTEMP = CL(VEL) 'Lift Coef as a function of VEL
'This is where Alpha comes into play
'CD and CL function were built to assume
'Alpha 31 stretch to Mach 8 then
'low alpha
DTEMP = COEFTEMP * VEL * VEL * CDTEMP 'temperary storage for drag
LTEMP = COEFTEMP * VEL * VEL * CLTEMP 'temperary storage for lift
D = DTEMP * RHO 'Drag acceleration, ft/sec/sec
L = LTEMP * RHO 'Lift acceleration, ft/sec/sec
'Acceleration = force/mass
G = SQR(D * D + L * L) / 32.2
IF G > GMAX THEN GMAX = G
VDOT = -D - Hdot * 32.2 / VEL 'Velocity dot, derivative
'with respect to time
DTEMPDOT = COEFTEMP * CDTEMP * 2 * VEL * VDOT 'Derivative of DTEMP with time
LTEMPDOT = COEFTEMP * CLTEMP * 2 * VEL * VDOT
DDOT = DTEMP * RHODOT + DTEMPDOT * RHO 'Deriv. of D with time
LDOT = LTEMP * RHODOT + LTEMPDOT * RHO 'Deriv. of L with time
BANK = -DELAZ * 2.2 'Bank angle in degrees, right is positive
IF BANK > 70 THEN BANK = 70 'Current built for 2.2 times Delaz
IF BANK < -70 THEN BANK = -70 'Limit to 70 degrees
COSB = COS(BANK / 57.29) 'Cosine of bank angle
SINB = SIN(BANK / 57.29) 'Sin of bank angle
' This section calculates the altitude derivatives:
HDDOT = L * COSB + (((VEL + 1359) / 25812) ^ 2! - 1!) * 32.2 - D * Hdot / VEL
'H Double dot -- the second derivative of altitude with time
HDDDOT = LDOT * COSB + 64.2 * (VEL + 1300!) * VDOT / (25500! * 25500!) - (DDOT * Hdot / VEL + D * HDDOT / VEL - D * Hdot * VDOT / (VEL * VEL))
'H Triple dot -- the third derivative of altitude with time
' This section calculates the horizontal position derivatives:
CVANGLE = VX / VELHORZ 'Cosine of the angle between the X axis
'and the horizontal velocity vector
SVANGLE = VY / VELHORZ 'Sin of the angle between the X axis
'and the horizontal velocity vector
SINFLIGHT = -Hdot / VEL 'Sin of the flight path angle
SINFLIGHTDOT = -DDOT / VEL + Hdot * VDOT / (VEL * VEL) '***
LV = L * SINFLIGHT * COSB 'The component of lift parellel to the
'the horizontal velocity vector
LP = -L * SINB 'The component of horizontal lift
'perpendicualar to the horizontal velocity
'vector
AX = -D * CVANGLE + LV * CVANGLE - LP * SVANGLE 'Acceleration in the
'X dirrection
AY = -D * SVANGLE + LV * SVANGLE + LP * CVANGLE 'Acceleration in the
'Y dirrection
LPDOT = -LDOT * SINB 'Derivative of LP with respect to time
LVDOT = LDOT * SINFLIGHT * COSB + L * SINFLIGHTDOT * COSB '****
CVANGLEDOT = AX / VELHORZ - VX * (-D + LV) / (VELHORZ * VELHORZ)
SVANGLEDOT = AY / VELHORZ - VY * (-D + LV) / (VELHORZ * VELHORZ)
AXDOT = -(DDOT * CVANGLE + D * CVANGLEDOT) + (LVDOT * CVANGLE + LV * CVANGLEDOT) + (LPDOT * SVANGLE + LP * SVANGLEDOT) 'The derivative of AX with time
AYDOT = -(DDOT * SVANGLE + D * SVANGLEDOT) + (LVDOT * SVANGLE + LV * SVANGLEDOT) + (LPDOT * CVANGLE + LP * CVANGLEDOT)'The derivative of AY with time
RETURN
PROPTIME:
' Find new values to use in the next time step
' Uses Taylor's Series expansion to propagate from the current time:
END.FLAG! = 0
TIME = TIME + DTIME 'Simulation time in seconds
H = H + Hdot * DTIME + HDDOT * DT2 + HDDDOT * DT3
Hdot = Hdot + HDDOT * DTIME + HDDDOT * DT2
X = X + VX * DTIME + AX * DT2 + AXDOT * DT3
Y = Y + VY * DTIME + AY * DT2 + AYDOT * DT3
RANGE = SQR(X * X + Y * Y)
VX = VX + AX * DTIME + AXDOT * DT2
VY = VY + AY * DTIME + AYDOT * DT2
VELHORZ = SQR(VX * VX + VY * VY) 'Horizontal velocity
VEL = SQR(Hdot * Hdot + VELHORZ * VELHORZ) 'Total velocity
DELAZTEMP = -(ATN(VY / VX) + ATN(-Y / X)) * 57.29
'Should really be the ATN2 function with the
'signs of the angle put in correctly.
IF ABS(DELAZTEMP) > 90 THEN
'Use old Delaz for printing.
'DELAZ should really use an
'ATAN2 function and if
'X becomes positive the
'DELAZ calculation is in
'error. But it means you
'have flown past the site
'and the program is stoped
'any way.
ELSE
DELAZ = DELAZTEMP
END IF
RETURN
PRINT.DATA:
' The following is a printout at each time step
' For the HP-48 you will need to store these parameter
' for future readout to the crew.
A$ = " #.#####^^^^ #.#####^^^^ #.#####^^^^ #.#####^^^^ #.#####^^^^ #.#####^^^^"
' PRINT #1, USING A$; TIME; VEL / 1000; DELAZ; RANGE / 6080; HDOT; D
RETURN
END1:
FUNCTION CD (VEL) STATIC
CD = .5
END FUNCTION
FUNCTION CL (VEL) STATIC
CL = .3 * CD(VEL)
'CL = 0!
END FUNCTION