Appendix 6 to ETSO-C117b – Wind Shear Simulation Model

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:

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]