'-----------------------------------------------------------
'
'					FLTSIM2D
'
'Longitudinal (2-D: Z & X) flight simulator 
'for Atari 800XL 2017 8-bit BASIC 10-Liner
'
'Jeff Piepmeier
'March 4, 2017
'http://jeffpiepmeier.blogspot.com/
'http://github.com/jeffpiep/
'
'Parsed with TurboBASIC XL Parser Tool 
'http://github.com/dmsc/tbxl-parser
$options +optimize, optimize=-convert_percent-const_replace
'Tested on Altirra
'http://www.virtualdub.org/altirra.html
'-----------------------------------------------------------

DIM K$(1)

'SET UP CONSTANTS AND INITIAL CONDITIONS
PI = 3.14159265
DT = .1 : REM (S) TIME STEP (0.1 SEEMS STABLE, 0.2 MIGHT BE OK, 1.0 IS UNSTABLE)
'U = 0 : REM (M/S) REFERENCE FRAME X-AXIS SPEED
'W = 0 : REM (M/S) REFERENCE FRAME Z-AXIS SPEED (+ IS DOWN)
'X = 0 : REM (M) INERTIAL FRAME 
'Z = 0 : REM (M) INERTIAL FRAME
'OMEGA = 0 : REM (RAD/S) INITIAL THETA RATE
'THETA = 0*PI/180 : REM (RAD) INITIAL THETA
MASS = 900 : REM (KG)
G = 9.807 : REM (M/S^2) ACCEL OF GRAVITY
IYY = 1825 : REM (KG-M^2) MOMENT OF INERTIA
RHO = 1.225 : REM (KG/M^2) STP AIR DENSITY
SW = 16.2 : REM (M^2) WIND AREA
ST = 3.25 : REM (M^2) TAILPLANE AREA
ER = 0.415 : REM ELEVATOR CHORD TO TAILPLANE CHORD RATIO
CM = -0.044 : REM WING MOMENT COEFFICIENT
CHORD = 1.4 : REM (M) MEAN CHORD LENGTH
AR = 7.38 : REM ASPECT RATIO 
ARMW = - CHORD * 0.2 : rem 0.2 : REM (M) CG RELATIVE TO CENTER OF LIFT "+"==AFTWARD
'THROTTLE = 0 : REM (%) THROTTLE SEETING
ARMT = 4.3 : REM (M) TAILPLANE MOMENT ARM RELATIVE TO a.c.
'DLTA = 0 : REM (DEG) ELEVATOR SETTING - UP IS NEGATIVE 
'FLAPS = 0 : REM FLAPS UP
'OKTOLAND = 0 : REM STATE FLAG FOR ENDING GAME
C = 0.12 : REM ALTITUDE POWER DROP OFF PARAMETER

CLS 

REPEAT

'ATMOSPHERE 
SIGMA=(1-Z*8.0E-05) : REM LINEARY APPROX FOR RELATIVE AIR DENSITY
RHOZ = SIGMA*RHO : REM AIR DENSITY AT ALTITUDE
STALL = (U-2*FLAPS)>25  : REM DETERMINE IF NOT STALLED
Q = (U^2+W^2)*RHOZ*0.5*STALL : REM DYNAMIC PRESSURE

'PLANE
ALPHA = ATN(W/(U + (U=0))) : REM ANGLE OF ATTACH WRT WIND
CL = SIN(2*ALPHA*(EXP(-ABS(ALPHA)/PI/8)+1))*1.35+0.25+0.167*FLAPS : REM LIFT COEFFICIENT 
CD = SIN(ALPHA)^2+0.0223+9.2E-03+0.01 + FLAPS*0.02 : REM DRAG COEFFICIENT
LIFTW = Q*CL*SW : REM WING LIFT FORCE
DRAGW = Q*CD*SW : REM WING DRAG FORCE

'TAILPLANE
'baseline
ALPHAT = ATN((W+OMEGA*ARMT)/(U+(U=0))) + 1.15 * DLTA*PI/180 * ER + 2*CL/PI/AR - 5*PI/180
'no rotational velocity
'ALPHAT = ATN((W)/(U+(U=0))) + 1.15 * DLTA*PI/180 * ER + 2*CL/PI/AR 
'no downwash
'ALPHAT = ATN((W+OMEGA*ARMT)/(U+(U=0))) + 1.15 * DLTA*PI/180 * ER 
'no downwash no rotational velocity
'ALPHAT = ATN((W)/(U+(U=0))) + 1.15 * DLTA*PI/180 * ER 

CLT = SIN(2*ALPHAT) : REM TAILPLANE LIFT COEFFICIENT
CDT = 1.8*SIN(ALPHAT)^2 : REM TAILPLANE DRAG COEFFICIENT
LIFTT = Q*CLT*ST : REM TAILPLANE LIFT FORCE
DRAGT = Q*CDT*ST : REM TAILPLANE DRAG FORCE

'ENGINE
'POWER DROPOFF WITH ALTITUDE FROM C
THRUST=THROTTLE*950*(SIGMA-C)/(1-C)/(U+(U=0)) : REM INVERT PROPULSIVE POWER EQUATION, AND MAX POWER IS 95kW
THRUST=THRUST*(THRUST<=2000)+2000*(THRUST>2000) : REM LIMIT THRUST TO 2000 N

'FORCES AND MOMENTS
ZFORCE = -(LIFTW+LIFTT)
XFORCE = -(DRAGW+DRAGT)
MY = Q*CM*SW*CHORD + LIFTW*ARMW - LIFTT*ARMT : REM + IS NOSE UP
'THERE SHOULD BE TWO DIFFERENT MOMENT REGIMES, BUT SWITCHING BETWEEN
'REGIMES CAUSES A PITCH OSCILLATION AND A/C LOSES ALTITUDE RIGHT 
'AFTER TAKEOFF. JUST USE THE FLIGHT REGIME AND DON'T WORRY ABOUT IT

'TIME STEP EQUATIONS OF MOTION
OMEGA = OMEGA + (MY/IYY)*DT
OMEGA = OMEGA * ((Z>.01) ! (OMEGA>0)) : REM NO ROTATION UNLESS POSITIVE WHEN ON GROUND
U = U + ( (XFORCE+THRUST)/MASS - G*SINTH - OMEGA*W ) * DT  
W = W + (     ZFORCE/MASS + G*COSTH + OMEGA*U ) * DT
W = W * ((Z>.01) ! (W<0)) : REM NO VERTICAL VELOCITY IF ON GROUND, UNLESS ITS UP

'INERTIAL FRAME UPDATE
T=T+DT
THETA = THETA + OMEGA*DT
COSTH = COS(THETA)
SINTH = SIN(THETA)
VX = (U*COSTH + W*SINTH)
VZ = (U*SINTH - W*COSTH)
X = X + VX*DT
Z = Z + VZ*DT
Z = Z * (Z>0) : REM RESTRICT X TO ABOVE GROUND

OKTOLAND = OKTOLAND ! (Z>152) : REM MUST GET TO 500 FT BEFORE YOU CAN LAND

'PILOT INPUT
K$=INKEY$
THROTTLE = THROTTLE + 5*((K$="\1F")*(THROTTLE<100)-(K$="\1E")*(THROTTLE>0))
THROTTLE = THROTTLE + (100-THROTTLE)*((K$="/") & (THROTTLE=0)) - THROTTLE*((K$="/") & (THROTTLE>0))
DLTA = DLTA + 0.5*((K$="T")*(DLTA<23)-(K$="B")*(DLTA>-28))
FLPIN = FLPIN + ((K$="N")*(FLPIN<3)-(K$="Y")*(FLPIN>0))
FLAPS = 0.2*FLPIN+0.8*FLAPS

POSITION 2,2
?"TIME",T," "
?"PITCH",INT(THETA*180/PI*10)/10," "
?"ALT.",INT(Z*3.28084)," "
?"KTAS",INT(U*1.94384)," "
?"VRATE",INT(VZ*196.85)," "
?"DIST.",INT(X)," "
?"POWER", THROTTLE," "
?"ELEV.", DLTA," "
?"FLAPS ", INT(FLAPS*10+0.5)," "
?"STALL ", NOT STALL

'lprint T,THETA*180/PI


UNTIL OKTOLAND & (Z=0) 
IF VZ >-1
	? "TOUCHDOWN"
ELSE
	? "YOU CRASHED!"
ENDIF
