The following computer listing (written in QuickBasic) provides a simplified aircraft simulation model for evaluating the effectiveness of various guidance schemes. This simulation runs on a personal computer, and the results obtained using it have been found to be comparable to those obtained on a full six-degrees-of-freedom simulator. This model was developed by J. Rene Barrios of the Honeywell Company.
The Wind Shear Simulation Model (WSSM) is a point mass three-degrees-of-freedom mathematical model which simulates the motion of an aeroplane in a vertical plane. The equations of motion, which are described in the wind axes, include the wind components of velocity and acceleration so that the aeroplane dynamics during a wind shear encounter are accurately modelled. This model has been used by several investigators to study the behaviour of an aeroplane during wind shear encounters.
Note: The Wind Shear Simulation Model provided at the end of this Appendix is an example written in Microsoft QuickBasic. Other programming languages such as Microsoft FORTRAN, C, or assembly language are also acceptable.
The Equations of Motion
The motion of a constant mass point in the vertical plane may be described by four equations of state and a control variable. For an aeroplane, it is convenient to use an orthogonal reference frame which is attached to the frame of the aeroplane, and its x-direction points in the direction of motion. Such a reference frame is the relative wind reference frame.
The following equations model the states of the aeroplane in the wind axes:
Vdt = g[(T.csalf) − D)/W − sngam] − Wxdt.csgam − Wzdt.sngam (1)
Gdt = {g[(T.snalf + L)/W − csgam] + Wxdt.sngam − Wzdt.csgam}/V (2)
Hdt = V.sngam + Wz (3)
Xdt = V.csgam + Wx (4) where:
| Vdt | Rate of change of true airspeed in knots/second |
|---|
| G | Gravitational constant in knots/second |
| T | Total engine thrust in lb |
| csalf | cos (alpha) |
| alpha | Angle of attack in radians |
| D | Total drag in lb. |
| W | Gross weight in lb |
| sngam | sin (gamma) |
| gamma | Flight path angle in radians |
| Wxdt | Inertial wind shear x-component in knots/second |
| Gdt | Rate of change of gamma in rad/second |
| snalf | sin (alpha) |
| L | Total lift in lb |
| V | True airspeed in knots |
| Hdt | Altitude rate in knots |
| Wz | Inertial wind z-component in knots |
| Xdt | Ground speed in knots |
In the above equations, the positive directions are upwards and forwards. This implies that tailwinds and updrafts are positive, while headwinds and downdrafts are negative. All the states can be determined from a given alpha; therefore, alpha is the control variable.
Since the model is that of a point mass, it is necessary to introduce the concept of alpha_command and actual alpha to account for the effect of the horizontal tail/elevator. This is done by introducing a lag between the alpha_command and the actual alpha. Therefore, any command that is given to the elevator or stabiliser can be interpreted as an alpha_command, which will cause a change in the angle of attack.
From equations 1, 2, 3 and 4, it can be seen that any change in alpha will produce a change in the longitudinal and normal accelerations, which in turn will change the states of the aeroplane.
The Path Control Function
The different segments of the trajectory flown by the WSSM are described by a series of alpha_commands, which are generated by the procedure explained below.
1.The aeroplane is trimmed for the initial conditions specified by the user. The initial conditions are usually specified as the altitude, gross weight, flaps, speed, flight path angle, and wind characteristics. The trimming operation consists in finding the angle of attack that satisfies the equations of state and will result in an unaccelerated motion.
2.After the initial trim, alpha_command must be specified for each segment of the trajectory, which usually consists of a climb or descent segment at constant speed or constant path angle, and guidance through wind disturbances. The wind disturbance is provided by wind models that can be selected at initialisation time.
3.In order to specify an alpha_command, the user must supply a subroutine where a quadratic function is defined in such a way that when minimised with respect to alpha, and constrained by the equations of state, the minimising alpha will produce the desired path in an optimal manner. For example, if we want to fly initially at a constant path angle, say 8 degrees, then the quadratic function may be defined by the expression:
(5)where:
| cst | Function to be minimised w.r.t alpha |
|---|
| dt
Gdt*dt | Time increment used in simulation in sec.
A predictive term which anticipates the change in gamma |
Other expressions follow:
The minimisation of the function cst is performed by a subroutine at each time frame and is totally transparent to the user, who has to supply only the objective function cst.
4.Each expression defining a different value of the objective function cst is called a ‘LAW’. The user selects the guidance law to be used during the wind shear encounter at the time the menu is displayed. This method allows the user to compare different guidance laws under the exact same conditions.
The Wind Models
The WSSM has two types of wind models: the Dallas-Ft Worth accident wind field, simulated by a quad_vortex model, and the constant shear model, which is user defined via the initial conditions menu.
Plotting Capabilities
The WSSM can plot up to 3 runs with 10 parameters per run. The length of each run should be kept under 60 seconds. This feature allows the user to compare different trajectories by overlaying the results.
The Programme
The WSSM is written in Microsoft QuickBasic, which is a highly structured language with a very friendly full-page editor. QuickBasic is very convenient for development, since it allows the user to stop execution, change the programme and continue executing. It also interfaces with Microsoft FORTRAN, C, or assembly language.
The procedure suggested for this application is that the WSSM be compiled without subroutines DETECT and GUIDE. DETECT and GUIDE can be separately compiled and put in a library called WNDSHR.QLB. These external subroutines may be written in Microsoft FORTRAN, C, or assembly language.
Listing of Programme
'************ AIRCRAFT FLIGHT PROFILE SIMULATION **********
DECLARE SUB PLOT ( )
DECLARE SUB TAKEOFF ( )
DECLARE SUB EULER ( )
DECLARE SUB MCRBRST ( )
DECLARE SUB WINDS ( )
DECLARE SUB OPT ( )
DECLARE SUB MIN (DM, M2, C1, C2, C3, M)
DECLARE SUB BEGIN ( )
DECLARE SUB VSHAKER ( )
DECLARE SUB COST ( )
DECLARE SUB LIMIT ( )
DECLARE SUB RATES ( )
DECLARE SUB THRUST ( )
DECLARE SUB ATMOS ( )
DECLARE SUB PRINTS ( )
DECLARE SUB DRAGS ( )
COMMON SHARED FLPS%, GEAR%, GEAR$, CL, CD, LIFT, DRAG, ALPHA
COMMON SHARED SEC, ALT, DST, HDOT, ALF, GAM, GAMREF, GREF, G
COMMON SHARED WSALERT%, WXO, WL, WX, WXDT, WZ, WZDT, DFW
COMMON SHARED WV, LC%, GM, GREFF, NOSAVE, GMO
COMMON SHARED DELTA, ISA, T0, SPDSND, VT, VC, MACH, A0, TAT, TAMF
COMMON SHARED THRST, EPR, TFCT, APPFLG%
COMMON SHARED SNGM, CSGM, CSAL, SNAL, VDOT, WG, GDOT, XDOT
COMMON SHARED AWX, AWZ, AU, AZ, VG, GRND, KF1, GMIN, KF2
COMMON SHARED ACMD, OLDALF, DT, HP, LP, ALFLIM
COMMON SHARED LAW%, GMR, ASS, CST, VTO, GCMD
COMMON SHARED OUTFILE$, DM, ALT1, PL$, TTT, WXDTO, TDX, TSH, WZO, TDZ, TSV
COMMON SHARED GM1, VTP, THETA
COMMON SHARED ALFRTE, PLMFLG%
'*******************************************
' MAIN PROGRAM *
'*******************************************
START: '< ----------------<< RE-RUNS START HERE
CLOSE : CLEAR
COLOR 15, 1: CLS : VIEW PRINT
LOCATE 8, 23: PRINT "WINDSHEAR SIMULATION"
LOCATE 10, 23: PRINT "FOR "
LOCATE 12, 23: PRINT "BOEING 737/200 "
LOCATE 23, 23: PRINT "TYPE " + CHR$(&H22) + "I" + CHR$(&H22) + " FOR INFORMATION" DO WHILE a$ = "" a$ = INKEY$
LOOP
IF a$ = "I" OR a$ = "i" THEN a$ = "": CLS
'---------------------- INFORMATION PAGE----------------------------------------- LOCATE 2, 2: PRINT "BOEING 737/200 INFORMATON"
LOCATE 3, 2: PRINT "JT8D-17 ENGINES"
LOCATE 5, 2: PRINT "-------------------------------------------------------------------------" LOCATE 7, 2: PRINT "ALLOWABLE WEIGHT RANGES..................: 75,000 TO 120,000 POUNDS" LOCATE 9, 2: PRINT "ALLOWABLE TAKEOFF FLAP SETTINGS.........: 1, 2, 5, 15, 20, 25 DEGREES" LOCATE 11, 2: PRINT "ALLOWABLE LANDING FLAP SETTINGS.........: 30, 40 DEGREES"
LOCATE 13, 2: PRINT "TAKEOFF EPR AT SEA LEVEL, STD. DAY...: 2.1 "
LOCATE 15, 2: PRINT "REFERENCE WING AREA.......................: 980 SQUARE FEET" LOCATE 17, 2: PRINT "REFERENCE TAKEOFF SPEED..................: V2 + 10"
LOCATE 19, 2: PRINT "REFERENCE LANDING SPEED..................: 1.3 Vs"
LOCATE 23, 2: PRINT "Press Any Key to Continue..."
DO: LOOP WHILE INKEY$ = ""
END IF
ANS$ = "2"
CLS
WHILE (ANS$ = "2")
LOCATE 10, 30: PRINT "Fly ..... 1"
LOCATE 12, 30: PRINT "Plot ..... 2"
LOCATE 14, 30: PRINT "Exit ..... 3"
LOCATE 18, 30: INPUT "Selection ....."; ANS$
IF ANS$ = "2" THEN
CALL PLOT
COLOR 15, 1
CLS
END IF
WEND
IF ANS$ = "3" THEN END
CALL BEGIN 'GET DATA/INITIALIZE VARIABLES
CALL THRUST 'INITALIZE THRUST
CALL TAKEOFF 'INITALIZE TAKEOFF
CALL COST 'SUBROUTINE COST
CALL PRINTS 'SUBROUTINE PRINT
FOR ICL% = 1 TO TTT ' TTT IS THE RUN TIME IN SECONDS
CALL THRUST ' SUBROUTINE EPR/THRUST
CALL WINDS ' SUBROUTINE WINDS
' CALL DETECT ' SUBROUTINE WINDSHEAR DETECTION
' SUPPLIED BY USER
' MUST RESIDE IN LIBRARY WNDSHR.QLB
CALL OPT ' SUBROUTINE OPTIMIZE
CALL LIMIT ' SUBROUTINE ALPHA RATE
CALL EULER ' SUBROUTINE INTEGRATE
CALL ATMOS ' SUBROUTINE ATMOSPHERE
CALL PRINTS ' SUBROUTINE PRINT
IF ALT < 0 THEN EXIT FOR
NEXT ICL%
PRINT "RUN IS COMPLETE"
PRINT "TYPE " + CHR$(&H22) + "D" + CHR$(&H22) + " FOR RUN DATA" a$ = ""
DO WHILE a$ = "" ' Wait for key to be pressed a$ = INKEY$
LOOP
VIEW PRINT: COLOR 15, 4: CLS
IF a$ = "D" OR a$ = "d" THEN a$ = ""
LOCATE 2, 2: PRINT "DATA FROM CURRENT RUN"
LOCATE 4, 2: PRINT "--------------------------------------------------------------" LOCATE 6, 2: PRINT "GROSS WEIGHT: "; WG; " POUNDS"
LOCATE 7, 2: PRINT "ISA DEVIATION: "; ISA; " DEG C"
LOCATE 8, 2: PRINT "FLAP POSITION: "; FLPS%; " DEGREES"
LOCATE 9, 2: PRINT "GEAR POSITION: "; GEAR$
LOCATE 11, 2: PRINT "CONTROL LAW: "; LAW%
LOCATE 12, 2: PRINT "GAMMA REFERENCE: "; GAMREF
LOCATE 13, 2: PRINT "PITCH LIMITING: "; PL$
IF PL$ = "YES" THEN
LOCATE 14, 2: PRINT "MAXIMUM PITCH: "; HP * 57.3; " DEGREES"
LOCATE 15, 2: PRINT "MINIMUM PITCH: "; LP * 57.3; " DEGREES"
END IF
LOCATE 16, 2: PRINT "TIME OF RUN: "; TTT * DT; " SECONDS"
IF DFW = 1 THEN
LOCATE 17, 2: PRINT "DALLAS/FW Wind Model"
ELSE
LOCATE 17, 2: PRINT "HORIZ. WIND MAGNITUDE "; WXO; " KNOTS"
LOCATE 18, 2: PRINT "HORIZ. SHEAR MAGNITUDE: "; WXDTO; " KNOTS/SECOND"
LOCATE 19, 2: PRINT "HORIZ. SHEAR DURATION: "; TDX; " SECONDS"
LOCATE 20, 2: PRINT "VERT. WIND MAGNITUDE: "; WZO * 1.689; " FT/SECOND"
LOCATE 21, 2: PRINT "VERT. WIND DURATION: "; TDZ; " SECONDS"
LOCATE 22, 2: PRINT "--------------------------------------------------------------"
END IF
IF LEN(OUTFILE$) = 0 THEN OUTFILE$ = "NONE"
LOCATE 23, 2: PRINT "OUTPUT FILE: "; OUTFILE$
LOCATE 24, 2: PRINT "Press Any Key to Continue...."
DO: LOOP WHILE INKEY$ = "" 'Wait for key to be pressed
END IF
GOTO START
END
SUB ATMOS STATIC
'******************************************************
' SUBROUTINE ATMOSPHERE *
'******************************************************
STATIC THETA
L% = ALT > 36089!
FISA = 1.8 * ISA
IF ALT > 36089 THEN
TMP = .7519 * T0
DELTA = .2234 * EXP((36089! - ALT) / 20806)
ELSE
TMP = T0 - .0035662 * ALT
DELTA = (TMP / T0) ^ 5.256
END IF
TAMB = TMP + FISA 'TAMBient in deg. R
TAMF = TAMB - 459.7 ' " " F
THETA = TAMB / T0
SQRTH = SQR(THETA)
SPDSND = A0 * SQRTH
IF VT > 0 THEN MACH = VT / SPDSND
VC = A0 * SQR(5 * (((1 + MACH * MACH / 5) ^ 3.5 - 1) * DELTA + 1) ^ .28571 - 5)
TAX = (TMP + FISA) * (1 + .2 * MACH * MACH) 'Deg. R
TAT = 5 * (TAX - 459.7 - 32) / 9 'Deg. C
IF INKEY$ < > "" THEN PRINT : INPUT "Press ENTER to continue...."; XXX
END SUB
SUB BEGIN STATIC
CLS : VIEW PRINT
'<--------------------------------- DATA_INPUT ------------------------------>
PRINT
INPUT "OUTPUT FILE (DEFAULT IS NO FILE) "; OUTFILE$
IF OUTFILE$ = " " THEN
NOSAVE = 1
ELSE
NOSAVE = 0
END IF
' CONSTANTS USED IN CALCULATIONS:
A0 = 661.478599# 'Speed of sound at sea level in knots
G = 19.07583 'Gravitational constant in knots/sec
T0 = 518.67 'Standard temperature at SL in deg Rankine
DT = .25 'Simulation time step in seconds
' ---------------------------------- INITIALIZATION OF VARIABLES-------------------------- GMIN = 0
VDOT = 0
ALT1 = 0
INPUT "TAKEOFF OR APPROACH (T/A) (Default is T)....."; ANS$
IF ANS$ = "a" OR ANS$ = "A" THEN
INPUT "ENTER ALTITUDE IN FEET (Default is 1000' . "; ALT1
IF ALT1 = 0 THEN ALT1 = 1000
APPFLG% = 1
TFCT = 1
END IF
ALT = ALT1
' ----------------------- CONFIGURTATION CONSTANTS------------------------------------------ ASS = 16.5 'Stick Shaker alpha in degrees
ASS = ASS / 57.3 ' " " " "radians
' ------------------------------------------- GROSS WEIGHT ENTRY------------------------------------------ PRINT : INPUT "ENTER GROSS WEIGHT IN POUNDS (Default is 110000) "; WG
IF WG = 0 THEN WG = 110000! ' DEFAULT SETTING
FL% = 0
WHILE (NOT FL%)
INPUT "ENTER FLAPS SETTING (Default is 0)............."; FLPS%
SELECT CASE FLPS%
CASE 0, 1, 2, 5, 15, 20, 25, 30, 40
FL% = -1
CASE ELSE
FL% = 0
PRINT "Invalid flaps setting"
PRINT "Only 0, 1, 2, 5, 15, 20, 25, 30, & 40 are supported"
PRINT
END SELECT
WEND
IF FLPS% < 15 THEN GEAR% = 1
IF FLPS% = 15 THEN INPUT "GEAR UP OR DOWN (1/0) (Default is Down)....."; GEAR%
IF GEAR% = 1 THEN
GEAR$ = " UP"
ELSE
GEAR$ = " DOWN"
END IF
INPUT "ENTER ISA DEV. IN DEGREES C (Default is 0)......"; ISA
PRINT
CALL VSHAKER ' COMPUTE V2+10 FOR FLAPS<33 OR 1.3Vs FOR FLAPS>32
PRINT " CONTROL LAW SELECTION:"
PRINT
PRINT " Speed = 1.1* V_stall = 1"
PRINT " Alpha = Stick Shaker Alpha = 2"
PRINT " Horizontal Acceleration = 0 = 3"
PRINT " 15_Degree Pitch = 4"
PRINT " Theoretical HONEYWELL/SPERRY = 5"
PRINT " User Defined = 6"
PRINT
INPUT " SELECT CONTROL LAW ....................... "; LAW%
IF LAW% = 0 THEN LAW% = 5
PRINT : PRINT
' ----------------------------------- GAMMA REFERENCE INPUT--------------------------------- IF LAW% > 4 THEN
INPUT "ENTER GAMMA REFERENCE IN DEGREES (Default is 0)........."; GMR
PRINT
GAMREF = GMR
GMR = GMR / 57.3: GMIN = GMR
END IF
' --------------------------------------- PITCH LIMITING SELECTION--------------------------------------------- INPUT "PITCH LIMITING DESIRED (y/n) (Default is NO)..............."; PL$
IF PL$ = "Y" OR PL$ = "y" THEN
PL$ = "YES"
INPUT " MAXIMUM PITCH ALLOWED IN DEGREES "; HP
INPUT " MINIMUM PITCH ALLOWED IN DEGREES "; LP
HP = HP / 57.3: LP = LP / 57.3: PL% = 1
ELSE
HP = 100
LP = -100
PL% = 0
PL$ = "NO"
END IF
CLS
' ----------------------------------------------- TIME FOR RUN----------------------------------------------------------------
PRINT
INPUT "ENTER TIME OF RUN IN SECONDS (Default is 45).............."; TTT
TTT = TTT / DT
IF TTT = 0 THEN TTT = 45 / DT ' DEFAULT SETTING
' ------------------------------------- WINDSHEAR SET UP---------------------------------------------------- INPUT "DALLAS/FW Wind Model (y/n)...(Default is constant Shear)....."; ANS$
IF ANS$ = "Y" OR ANS$ = "y" THEN
DFW = 1
ELSE
DFW = 0
PRINT
INPUT "MAGNITUED OF HORZ. WIND IN KNOTS......(Head wind < 0)......"; WXO
INPUT "MAGNITUED OF HORZ. SHEAR IN KT/SEC. (Dec. Perf. > 0)......"; WXDTO
INPUT "DURATION OF HORZ. SHEAR IN SEC.............(Default is 0)......"; TDX
INPUT "TIME FOR SHEAR TO START IN SEC..............(Default is 0)......"; TSH
PRINT
INPUT "MAGNITUED OF VERT. WIND IN FT/SEC. (Down Draft < 0)......"; WZO
WZO = WZO / 1.689 'Convert to knots
INPUT "DURATION OF VERT. WIND IN SEC...........(Default is 0)........"; TDZ
INPUT "TIME FOR SHEAR TO START IN SEC...........(Default is 0)......"; TSV
PRINT
END IF
' -------------------- OTHER SET UPS---------------------------------------------------- VT = VTO
WX = WXO
CALL ATMOS ' SUBROUTINE ATMOSPHERE
' ---------------------------------- HEADERS FOR SCREEN DISPLAY---------------------------------------------------- CLS : PRINT
PRINT "TIME ALT HDOT VT ALPHA GAMMA PITCH GREF WXDT WZ VDOT ALRT"
PRINT "(SEC) (FT) (FPM) (KTS) (DEG) (DEG) (DEG) (DEG) (KT/S) (FPS) (KT/S)"
PRINT STRING$(75, "-"): VIEW PRINT 5 TO 25
'**************************************
' SUBROUTINE INIT_OUTPUT FILE *
'**************************************
IF NOSAVE THEN ' CREATE OUTPUT FILE
ELSE
OPEN "O", 2, OUTFILE$
FMT$ = " ###.## ##### ##### #### #### ## ###.## ###.## ###.##"
FMT$ = FMT$ + " ###.## ###.# ###.# ##.## "
END IF
END SUB
SUB COST STATIC
'************************************************************
' SUBROUTINE COST *
'************************************************************
CALL DRAGS ' SUBROUTINE DRAG & LIFT
CALL RATES ' SUBROUTINE RATES
IF LC% = 0 THEN 'Constant gamma segment
FCT = (GM + GDOT * DT - GMO) ^ 2
GREFF = 57.3 * GMO
ELSE 'All guidance laws
SELECT CASE LAW%
CASE 1 ' ---------------------------- 1.1*Vstall----------------------------------- CST = (VT + VDOT * DT - 1.1 * 135) ^ 2
CASE 2 ' ---------------------------- Alpha = Ass----------------------------------- CST = (ALPHA - ASS) ^ 2
CASE 3 ' ---------------------------- Ax = 0-----------------------------------
CST = (VDOT - VT * GDOT * GM + WXDT) ^ 2
CASE 4 ' ---------------------------- 15 Degrees----------------------------------- CST = (GM + 3 * GDOT * DT + ALPHA - 15 / 57.3) ^ 2
CASE 5 ' ---------------------------- User Defined----------------------------------- PRINT "Not defined"
STOP
CASE 6 ' ---------------------------- User Supplied----------------------------------- 'User must supply a subroutine called GUIDE
'which must reside in the WNDSHR.QLB Library
'GUIDE can have a list of arguments
'As an example
'ALF = 57.3*ALPHA
'PTH = 57.3 * (ALPHA + GM)
' units : ft fpm kt deg g's *
'CALL GUIDE(ALT, HDOT, VC, ALF, PTH, AU, AZ, CST)
END SELECT
END IF
' CST is the Cost Function to be minimized
END SUB
SUB DRAGS STATIC
'*****************************************************
' SUBROUTINE DRAG FOR B737/200 *
'*****************************************************
X = 57.3 * ALPHA + 1
CF5 = 0: CF4 = 0: CF3 = 0: CF2 = 0
SELECT CASE FLPS%
CASE 0
CF1 = .091
CF0 = .0156
CASE 1
CF3 = -1.164058E-04
CF2 = 2.48561E-03
CF1 = .0905781
CF0 = .062114
CASE 2
CF0 = .101198
CF1 = .110993
CF2 = -.0015162
CF3 = 1.8931E-04
CF4 = -7.1427E-06
CF5 = -4.2776E-09
CASE 5
CF0 = .192638
CF1 = .123509
CF2 = -.0051477
CF3 = 6.4968E-04
CF4 = -3.0891E-05
CF5 = 4.1291E-07
CASE 10
CF0 = .249855
CF1 = .114005
CF2 = 7.1207E-04
CF3 = -9.9541E-05
CF4 = 7.0431E-06
CF5 = -2.3773E-07
CASE 15
CF0 = .40149
CF1 = .118723
CF2 = -6.4877E-04
CF3 = 6.6281E-05
CF4 = -1.6113E-07
CF5 = -1.4278E-07
CASE 25
CF0 = .592655
CF1 = .122433
CF2 = -.0026365
CF3 = 3.5963E-04
CF4 = -1.5579E-05
CF5 = 1.0894E-07
CASE 30
IF X < 4 THEN
CF1 = .12
CF0 = .72
ELSE
CF3 = -1.651192E-04
CF2 = 4.16461E-03
CF1 = 8.337061E-02
CF0 = .8350316
END IF
CASE 40
IF X < 4 THEN
CF1 = .12
CF0 = 1.08
ELSE
CF3 = -1.689903E-04
CF2 = 3.733285E-03
CF1 = 8.483822E-02
CF0 = 1.201596
END IF
CASE ELSE
PRINT "Flaps "; FLPS%; " not available....."
END
END SELECT 'For CL computation
CL = ((((CF5 * X + CF4) * X + CF3) * X + CF2) * X + CF1) * X + CF0
SELECT CASE FLPS% 'Low Speed Drag Polars
CASE 0
D0 = .013285: D1 = .052868: D2 = -.07182: D3 = .071561
CASE 1
D0 = .026143: D1 = .022358: D2 = -.00083: D3 = .016338
CASE 2
D0 = .070346: D1 = -.0852: D2 =.097453: D3 = -.01207
CASE 5
D0 = .045214: D1 = -.0178: D2 =.04373: D3 = .002101
CASE 10
D0 = -.04266: D1 = .19643: D2 = -.1152: D3 = .03966
CASE 15
IF GEAR% = 0 THEN
D0 = .034954: D1 = .098892: D2 = -.04187: D3 = .020496
ELSE
D0 = -.02822: D1 = .174631: D2 = -.0874: D3 = .029566
END IF
CASE 25
D0 = -.10416: D1 = .327506: D2 = -.17059: D3 = .043313
CASE 30
D0 = .124697: D1 = -.03348: D2 =.055295: D3 = -.00311
CASE 40
D0 = .124925: D1 = .052537: D2 =.006912: D3 = .0058
CASE ELSE
PRINT "Flaps "; FLPS% " not available...."
END
END SELECT
CD = ((D3 * CL + D2) * CL + D1) * CL + D0
Q = 1451770 * MACH * MACH * DELTA 'B737/200
LIFT = Q * CL
DRAG = Q * CD
END SUB
SUB EULER STATIC
'*****************************************************************
' SUBROUTINE EULER'S PREDICTOR/CORRECTOR *
' (INTEGRATION SUBROUTINE) *
'*****************************************************************
DTH = DT / 3600: DTM = DT / 60: SEC = SEC + DT: VTP = VT
CALL RATES ' SUBROUTINE RATES <<PREDICTOR>>
ALT1 = ALT: HDOT1 = HDOT: ALT = ALT + HDOT * DTM
GM1 = GM: GDOT1 = GDOT: GM = GM + GDOT * DT
DST1 = DST: XDOT1 = XDOT: DST = DST + XDOT * DTH
VT1 = VT: VDOT1 = VDOT: VT = VT + VDOT * DT
CALL RATES ' SUBROUTINE RATES <<CORRECTOR>>
ALT = ALT1 + (HDOT1 + HDOT) * DTM / 2
GM = GM1 + (GDOT1 + GDOT) * DT / 2
DST = DST1 + (XDOT1 + XDOT) * DTH /2
VT = VT1 + (VDOT1 + VDOT) * DT / 2
END SUB
SUB LIMIT STATIC
'*****************************************************************************
' SUBROUTINE ALPHA DOT AND PITCH LIMIT *
'*****************************************************************************
ALPHA = OLDALF + .25 * (ACMD - OLDALF) 'Pitch dynamics
CALL DRAGS ' SUBROUTINE DRAG (REQ'D FOR RATE SUB CALL)
IF PLMFLG% = 0 THEN EXIT SUB
OLDGM = GM
PLIM% = 0
DO WHILE (PLIM% = 0)
CALL RATES ' SUBROUTINE RATES
X = ALPHA + OLGM + GDOT * DT
IF X > HP THEN ALPHA = .9 * ALPHA
IF X < LP THEN ALPHA = 1.1 * ALPHA
IF ALPHA > ALFLIM THEN
ALPHA = ALFLIM
PLIM% = 1
END IF
LOOP
END SUB
SUB MCRBRST STATIC
IF MU1 = 0 THEN
MU1 = -37141!
AV = 5500: H1 = 2500: G3 = 3: J1 = -700: J2 = 800: J3 = 6.5
MU2 = -20000
BV = 12000: H2 = 2000: N1 = 200: N2 = 2500: N3 = 4
WX = 5
IF ALT > 1000 THEN
PRINT
PRINT " DFW data not available above 1000'"
PRINT " Please start at or below 1000'"
END
END IF
END IF
X = 6078 * DST: Y = ALT: A1 = AV: A2 = BV
NX1 = Y - H1: DENX1 = (Y - H1) ^ 2 + (X - A1) ^ 2
NY1 = X + J2 - A1: DENY1 = (Y + J1 - H1) ^ 2 + (X + J2 - A1) ^ 2
NX2 = Y - H2: DENX2 = (Y - H2) ^ 2 + (X - A2) ^ 2
NY2 = X + N2 - A2: DENY2 = (Y + N1 - H2) ^ 2 + (X + N2 - A2) ^ 2
NX3 = Y + H1: DENX3 = (Y + H1) ^ 2 + (X - A1) ^ 2
NY3 = X + J2 - A1: DENY3 = (Y + J1 + H1) ^ 2 + (X + J2 - A1) ^ 2
NX4 = Y + H2: DENX4 = (Y + H2) ^ 2 + (X - A2) ^ 2
NY4 = X + N2 - A2: DENY4 = (Y + N1 + H2) ^ 2 + (X + N2 - A2) ^ 2
XX = MU1 * (-NX1 / DENX1 + NX3 / DENX3) + MU2 * (NX2 / DENX2 - NX4 / DENX4)
WX = WX + .65 * (XX - WX) + 2 * G3
IF DST = 0 THEN WXP = WX
ZZ = MU1 * (NY1 / DENY1 - NY3 / DENY3) * J3 + MU2 * (-NY2 / DENY2 + NY4 / DENY4) * N3 WZ = WZ + .65 * (ZZ - WZ)
IF DST = 0 THEN WZP = WZ
WX5 = WX4: WX4 = WX3: WX3 = WX2: WX2 = WX1: WX1 = WX
WZ5 = WZ4: WZ4 = WZ3: WZ3 = WZ2: WZ2 = WZ1: WZ1 = WZ
IF WCNT% < 4 THEN WXDT = (WX - WXP) / DT: WXP = WX
IF WCNT% < 4 THEN WZDT = (WZ - WZP) / DT: WZP = WZ
IF WCNT% > 3 THEN WXDT = (26 * WX5 - 27 * WX4 - 40 * WX3 - 13 * WX2 + 54 * WX1) / (70 * DT) IF WCNT% > 3 THEN WZDT = (26 * WZ5 - 27 * WZ4 - 40 * WZ3 - 13 * WZ2 + 54 * WZ1) / (70 * DT) IF ABS(WXDT) > 15 THEN WXDT = 15 * SGN(WXDT)
IF ABS(WZDT) > 15 THEN WZDT = 15 * SGN(WZDT)
WCNT% = WCNT% + 1
END SUB
SUB MIN (DM, M2, C1, C2, C3, M) STATIC
'************************************************************
'SUBROUTINE MIN_CST BY LEAST SQUARES PARABOLA *
'************************************************************
ALPHA = M2 + DM 'INCREMENT ALPHA
CALL COST 'SUBROUTINE COST
IF DM < 0 THEN
C4 = CST
ELSE
SWAP C1, C3
C5 = CST
END IF
ALPHA = M2 - DM 'DECREMENT ALPHA
CALL COST 'SUBROUTINE COST
IF DM < 0 THEN
C5 = CST
ELSE
C4 = CST
END IF
M = ABS(DM) * (14 * C1 + 7 * C4 - 7 * C5 - 14 * C3) / (20 * C1 - 10 * C4 - 20 * C2 - 10 * C5 + 20 * C3)
END SUB
SUB OPT STATIC
'****************************************************************************
'SUBROUTINE OPTALF - DETERMINES THE ALPHA REQD FOR CMD GAMMA *
'****************************************************************************
OLDALF = ALPHA: GM1 = GM
CALL ATMOS ' SUBROUTINE ATMOSPHERE
CALL RATES ' SUBROUTINE RATES
DM = 1 / 57.3 ' SET ALPHA INCREMENT TO 1 DEGREE
C1 = 1E+20
C2 = 1E+20
C3 = 1E+20
OPTFLG% = 0
WHILE (OPTFLG% = 0)
CALL COST ' SUBROUTINE COST
C3 = C2: C2 = C1: C1 = CST
M3 = M2: M2 = M1: M1 = ALPHA
LGC% = C1 > C2 AND C3 = 1E+20
IF LGC% THEN
DM = -DM ' Reverse search direction
C1 = C2: C2 = CST: M1 = M2: M2 = ALPHA
ALPHA =ALPHA + 2 * DM
ELSE
IF C1 < C2 THEN
L% = ABS(OLDALF - ALPHA) / DT > ALFRTE OR ALPHA > ALFLIM OR ALPHA < -.08
IF L% THEN OPTFLG% = 1
ALPHA = ALPHA + DM
ELSE
DM = DM / 2
CALL MIN(DM, M2, C1, C2, C3, M)'Fit parabola & find minimum
ALPHA = M2 + M 'This is the optimum alpha
OPTFLG% = 1 'Set flag to terminate
END IF
END IF
WEND
ALFLIM = ASS 'SET ALPHA LIMIT TO ALPHA STICK SHAKER
SELECT CASE LAW%
CASE 4
ALFLIM = ASS - .035 'LIMIT TO SS MINUS 2 DEG
CASE 5, 6
ALFLIM = ASS - KF2
CASE ELSE
END SELECT
IF ALPHA < -.08 THEN ALPHA = -.08
IF ALPHA > ALFLIM THEN ALPHA = ALFLIM
ACMD = ALPHA 'SET ALPHA COMMAND TO COMPUTED ALPHA
END SUB
SUB PLOT
'************************************************************************************* '* PLOT ROUTINE *
'************************************************************************************* REM $DYNAMIC
' TWO DIMENSIONAL PLOTTER
DEFINT I-L, N
DIM F$(3) ' file name array
DIM DTA(3, 250, 15) ' data array
DIM TY$(14) ' title array (dependant variable)
TITLE$ = "HONEYWELL WINDSHEAR SIMULATION" ' main title
TX$ = "Time (s)" ' X title
TY$(1) = "Altitude ft "
TY$(2) = "Alt Rate fpm "
TY$(3) = "T A S kts "
TY$(4) = "Alpha deg "
TY$(5) = "Gamma deg "
TY$(6) = "Pitch deg "
TY$(7) = "G_ref deg "
TY$(8) = "Hz Shear kps "
TY$(9) = "Vt Wind fps "
TY$(10) = "Vt rate kps "
TY$(11) = "W/S Flag "
NV = 12
CLS
LOCATE 3, 15: PRINT "Enter the names of the data files you wish to plot."
FOR NC = 1 TO 3
LOCATE 6 + 2 * NC, 25 ' input
PRINT "FILENAME "; NC; " "; ' filenames
INPUT ; F$(NC) ' containing
IF F$(NC) = " " THEN EXIT FOR ' data
NEXT NC
NC = NC - 1 ' number of curves to plot
LOCATE 20, 15: PRINT "Reading from disk........."
FOR I = 1 TO NC
CLOSE
OPEN "I", #1, F$(I) ' open file for input
NP = 0
DO
NP = NP + 1 ' number of points
FOR J = 1 TO NV
INPUT #1, DTA(I, NP, J) ' read data
NEXT J
LOOP UNTIL EOF(1)
CLOSE
NEXT I
DO ' display all selected parameters
DO ' prompt user until a valid parameter is selected
100 CLS
LOCATE 3, 20: PRINT "Select the parameter you wish to plot."
FOR I = 1 TO NV - 1
LOCATE 4 + I, 30: PRINT TY$(I); " = "; I
NEXT I
LOCATE 21, 30: INPUT "parameter number (0 to exit)"; PARAM%
IF PARAM% = 0 THEN
CLS
EXIT SUB ' return to calling program
END IF
LOOP UNTIL 1 <= PARAM% AND PARAM% <= 14 'end of select loop
PARAM% = PARAM% + 1
DX = 5 ' x axis grid increment
GOSUB 400 ' find maximum x and y values
IF PLTFLG% = 1 THEN
PRINT "No information to plot...."
PRINT "Press any key to continue..."
DO: LOOP WHILE INKEY$ = ""
GOTO 100
END IF
GOSUB 600 ' grid and titles
FOR I = 1 TO NC
GOSUB 1110 ' plot graph
NEXT I
DO
LOOP WHILE INKEY$ = ""
CLS : SCREEN 0
LOOP
'************************************************************************************** 400 '* MAX SUBROUTINE *
'************************************************************************************** '
MAXX = DTA(1, 1, 1)
MAXY = DTA(1, 1, PARAM%)
MINY = DTA(1, 1, PARAM%)
FOR I = 1 TO NC
FOR J = 1 TO NP
IF DTA(I, J, 1) > MAXX THEN MAXX = DTA(I, J, 1)
IF DTA(I, J, PARAM%) > MAXY THEN MAXY = DTA(I, J, PARAM%)
IF DTA(I, J, PARAM%) < MINY THEN MINY = DTA(I, J, PARAM%)
NEXT J
NEXT I
PLTFLG% = 0
DY = (MAXY - MINY) / 15
IF DY = 0 THEN
PLTFLG% = 1
DY = 5
END IF
MAG = 10 ^ (INT(LOG(DY) / LOG(10))): DY = DY / MAG
IF DY <= 5 THEN
DY = 5
ELSE
DY = 10
END IF
DY = DY * MAG
IF INT(MAXX / DX) < > MAXX / DX THEN MAXX = INT(MAXX / DX + 1) * DX
IF INT(MAXY / DY) < > MAXY / DY THEN MAXY = INT(MAXY / DY + 1) * DY
IF INT(MINY / DY) < > MINY / DY THEN MINY = INT(MINY / DY) * DY
NUMX = MAXX / DX
NUMY = (MAXY - MINY) / DY
RETURN
600
'*************************************************************************************** '* GRID AND TITLES *
'*************************************************************************************** '
CLS
SCREEN 2 ' 640*200 monochrome graphics
KEY OFF
'
FOR J = 0 TO NUMX
Z = J * 580 / NUMX + 59
LINE (Z, 10) - (Z, 170) ' vertical grid line
Z = J * 71 / NUMX + 7 a = DX * J
IF a < > 0 THEN ' adjustment for D = INT(LOG(a) / LOG(10)) + 1 ' large numbers
IF D > 1 THEN Z = Z - D + 1
END IF
LOCATE 23, Z
PRINT a;
NEXT J
FOR J = 0 TO NUMY
Z = J * 160 / NUMY + 10
LINE (60, Z) - (640, Z) ' horizontal grid line
Z = 22 - J * 20 / NUMY
LOCATE Z, 2
Z = DY * J + MINY
AZ = ABS(Z)
IF INT(Z) = Z THEN
G$ = "######"
ELSEIF AZ < .1 THEN
G$ = "#.####"
ELSEIF AZ >= .1 AND AZ < 1 THEN
G$ = "##.###"
ELSEIF AZ >= 1 AND AZ < 10 THEN
G$ = "###.##"
ELSEIF AZ >= 10 AND AZ < 100 THEN
G$ = "####.#"
ELSE
G$ = "######"
END IF
PRINT USING G$; Z;
NEXT J
Z = (80 - LEN(TITLE$)) / 2 + 2
LOCATE 1, Z: PRINT TITLE$ ' print main title
LOCATE 24, 36: PRINT TX$; ' X axis title
LOCATE 8, 1 ' Y
FOR J = 1 TO LEN(TY$(PARAM% - 1)) ' axis
PRINT MID$(TY$(PARAM% - 1), J, 1) ' title
NEXT J
LOCATE 25, 10: PRINT "1"; ' curve
LINE (90, 195) - (130, 195)
LOCATE 25, 20: PRINT "2"; ' labels
FOR J = 0 TO 40 STEP 8
XX = 170 + J
PSET (XX, 195)
CIRCLE (XX + 80, 195) , 2
NEXT J
LOCATE 25, 30: PRINT "3";
RETURN
'***************************************************************************************** '* PLOTTING ROUTINE *
'***************************************************************************************** '
1110 FOR J = 1 TO NP
XX = 580 * DTA(I, J, 1) / MAXX + 60' calculate pixel X position
YY = 170 - 160 * (DTA(I, J, PARAM%) - MINY) / (MAXY - MINY)
IF J = 1 THEN GOTO 1170
IF I = 1 THEN LINE (XXOLD, YYOLD) - (XX, YY) ' line 1170
1170 XXOLD = XX: YYOLD = YY
IF I = 2 THEN PSET (XX, YY) ' point
IF I = 3 THEN CIRCLE (XX, YY), 2 ' circle
NEXT J
RETURN
END SUB
REM $STATIC
DEFSNG I-L, N
SUB PRINTS
'******************************************************************
' SUBROUTINE PRINT TO SCREEN AND FILE *
'******************************************************************
ACMDG = 57.3 * ACMD
ALF = 57.3 * ALPHA
GAM = 57.3 * GM
PITCH = ALF + GAM
WZX = 1.689 * WZ
IF NOSAVE = 0 THEN PRINT #2, SEC, ALT, HDOT, VT, ALF, GAM, PITCH, GREFF, WXDT, WZX, VDOT, WSALERT% FMT1$ = "###.## #### ##### ### ###.# ###.# ###.# ##.# ###.# ###.# ###.# #"
PRINT USING FMT1$; SEC, ALT, HDOT, VT, ALF, GAM, PITCH, GREFF, WXDT, WZX, VDOT, WSALERT%
END SUB
SUB RATES STATIC
'******************************************************************
' SUBROUTINE RATES *
'******************************************************************
SNGM = SIN(GM): CSGM = COS(GM): SNAL = SIN(ALPHA): CSAL = COS(ALPHA)
VDOT = G * ((THRST * CSAL - DRAG) / WG - SNGM) - WXDT * CSGM - WZDT * SNGM
GDOT = G * ((LIFT + THRST * SNAL) / WG - CSGM) + WXDT * SNGM - WZDT * CSGM
GDOT = GDOT / VT
HDOT = 101.28 * (VT * SNGM + WZ)
XDOT = VT * CSGM + WX
AWX = VDOT + WXDT * CSGM + WZDT * SNGM 'Inertial Acc. along Wind_x axis
AWZ = VT * GDOT - WXDT * SNGM + WZDT * CSGN 'Inertial Acc. along Wind_z axis
AU = (AWX * CSAL + AWZ * SNAL) / G 'LONG. ACCEL. ->=?
AZ = (AWX * SNGM + AWZ * CSGM) / G 'VERT. ACCEL. UP=?
VG = XDOT
GRND = (VT * GM + WZ) / (VT + WX) 'Gamma w/r ground
KF1 = 1
GHAT = GMIN * (1 + WX / VT)
IF WZ > -30 AND WZ < -20 THEN KF1 = 1 + .025 * (WZ + 20)
IF WZ <= -30 THEN KF1 = .75
DGAM = 57.3 * (20 * GDOT - (GHAT - GRND + (1 - KF1) * WZ / 152 + 20 * GDOT))
IF DGAM < 0 THEN
KF2 = (2 + .4 * DGAM)
ELSE
KF2 = 2
END IF
IF KF2 < 0 THEN KF2 = 0
KF2 = KF2 / 57.3
END SUB
SUB TAKEOFF STATIC
'********************************************************************
' SUBROUTINE INTIALIZE TAKEOFF *
'********************************************************************
IF APPFLG% = 0 THEN
ALPHA = .12
WHILE (LIFT <= WG)
CALL DRAGS
ALPHA = ALPHA + .01
WEND
GM = (THRST - DRAG) / WG 'COMPUTE POTENTIAL GAMMA
ELSE
GM = -3 / 57.3
ALPHA = 2 / 57.3
CALL DRAGS
TFCT = 1
CALL THRUST
T = DRAG - .052 * WG
IF T < 0 THEN T = .2 * THRST
TFCT = T / THRST
THRST = T
END IF
GMO = GM
CALL RATES
END SUB
SUB THRUST STATIC
'*******************************************************************
' SUBROUTINE EPR/THRUST *
'*******************************************************************
' TAKE-OFF THRUST FOR JT8D-17 ENGINES
VE = 1.668 * VT
R00 = 14688.74: R01 = -.65187546#: R02 = 6.7371E-05
R10 = -13.9295: R11 = .000751143#: R12 = -1.5405E-07
R20 = .014643: R21 = 5.3444E-07: R22 = -4.8907E-10
AA0 = (R02 * ALT + R01) * ALT + R00
AA1 = (R12 * ALT + R11) * ALT + R10
AA2 = (R22 * ALT + R21) * ALT + R20
THRST = 2 * ((AA2 * VT + AA1) * VT + AA0) 'Temp. = 100 F
IF APPFLG% = 1 THEN
IF LC% = 1 AND TFCT < 1 THEN
GMO = .136
TSPL = 5.5
'Engine Spool Up Time
TFCT = TFCT + DT / TSPL
END IF
IF TFCT > 1 THEN TFCT = 1
ELSE
TFCT = 1
END IF
THRST = TFCT * THRST
'' THRST = 2 * (((2.64159E-05 * VT + 5.110896E-03) * VT - 12.56476) * VT + 15550)
END SUB
SUB VSHAKER STATIC
'------------------------------------ COMPUTATION OF Vss AND V2------------------------- V2 = 145
VTO = V2 + 10' SETS INITAL SPEED EQUAL TO V2 + 10
SELECT CASE FLPS%
CASE 10
IF VTO < 150 THEN VTO = 150 ' TAKEOFF
CASE 18
IF VTO < 148 THEN VTO = 148 ' FLAP
CASE 22
IF VTO < 147 THEN VTO = 147 ' SETTINGS
CASE 33
VTO = 63.11225 + .222468 * WG / 1000 ' APPROACH
CASE 42 ' FLAP
VTO = 62.67386 + .21744 * WG / 1000 ' SETTINGS
CASE ELSE
END SELECT
END SUB
SUB WINDS STATIC
'******************************************************************
' SUBROUTINE WINDS *
'******************************************************************
'
IF TDX > 0 THEN
T1 = 4
T2 = TSH
T3 = T1 + T2
T4 = -4
T5 = T3 + TDX
T6 = T5 - T4
B1 = 3 * WXDTO / T1 ^ 2
A1 = -2 * B1 / (3 * T1)
B2 = 3 * WXDTO / T4 ^ 2
A2 = -2 * B2 / (3 * T4)
IF SEC > T2 AND SEC < = T3 THEN
X = SEC - T2
WXDT = (A1 * X + B1) * X * X
END IF
IF SEC > T5 AND SEC < = T6 THEN
X = SEC - T6
WXDT = (A2 * X + B2) * X * X
END IF
IF SEC > T6 THEN WXDT = 0
WX = WX + WXDT * DT
END IF
IF TDZ > 0 THEN
T1 = 4
T2 = TSV
T3 = T1 + T2
T4 = - 4
T5 = T3 + TDZ
T6 = T5 - T4
B1 = 3 * WZO / T1 ^ 2
A1 = -2 * B1 / (3 * T1)
B2 = 3 * WZO / T4 ^ 2
A2 = -2 * B2 / (3 * T4)
IF SEC > T2 AND SEC <= T3 THEN
X = SEC - T2
WZ = (A1 * X + B1) * X * X
WZC = WZ
END IF
IF SEC > T5 AND SEC <= T6 THEN
X = SEC - T6
WZ = (A2 * X + B2) * X * X
WZC = WZ
END IF
KALT = (-.0000011 * ALT + .00212) * ALT - .0251
IF KALT < 0 OR ALT <= 0 THEN KALT = 0
KALT = 1
WZ = KALT * WZC
IF SEC > T6 THEN WZ = 0
WZDT = (WZ - WZ1) / DT
WZ1 = WZ
END IF
IF DFW = 1 THEN CALL MCRBRST 'DALLAS Model
END SUB
[Amdt ETSO/16]