is a two-dimensional flight simulator. It simulates the longitudinal flight of an aircraft using three degrees-of-freedom (3DOF) equations of motion. This means the aircraft can fly forward, pitch, ascend and descend, but does not turn, yaw or roll. You will also see that it is not possible to invert the aircraft because I made some math approximations. But that's what one gets in less than 10 lines of BASIC code that can run realtime on the right machine! This project became a wee bit of an obsession and I am glad it is finally finished. The code and instructions for both Atari 800XL and BBC Micro Model B are
, I learned what game developers have known since the beginning: simulating physics with little bit of math can yield complex play action in only a few lines of code. The first digital video game exploited real physics:
I wanted to apply this technique to a game for 2017. Listening to two Atari podcasts also helped lead to me to choose to write FLTSIM2D. First, Kevin Savetz on
described how he wanted to write a text adventure for the 2016 10-liner and produced
. After hearing this, I started to think about what kind of game could be done for the contest that was completely different from typical arcade or puzzle style games. Second, Rob McMullen made an episode on Atari flight simulators for his
. It was fun to hear about
, which I played for endless hours on my Atari 800XL as a teenager. These thoughts came together (physics, text, and flight simulator), which led to FLTSIM2D.
Why 2-D? My college professor always said, "when you face a problem, simplify it down to something you know how to solve, and then start adding back in the details a little at a time." I had no idea how to do a full motion simulator and it would surely be too slow in BASIC. So 2-D should be easy, right? Lift, drag, thrust, and gravity - should be pretty simple. Well, it was in the end, but it took a while to get there. I learned enough to be dangerous from
; however, my airplane cartwheeled and the floating-point overflowed in the first version I wrote. After another try in BASIC, I resorted to a spreadsheet to try to first get the math right. I drew a free-body diagram (after recalling my statics and dynamics professors saying, "always start with a free-body diagram.") Ah! I had a sign wrong. Switch a positive to a negative, and voila, I had a stable system. I didn't know when I started this project that an aircraft is a second-order system, which means it can oscillate. At this point, I was not completely satisfied my math was correct. I'm pretty sure I had two different coordinate systems mixed together. I also tried not to get too discouraged when I found someone had made this the topic of their
. After all, I did have something kind of working. As I said ... obsession.
on 6-DOF differential equations of motion. This lecture was a turning point for me. Here's the code with some description of each command and equation with links to online resources as best as I can remember.
'-----------------------------------------------------------
'
' 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
$options +optimize, optimize=-convert_percent-const_replace, optimize=+const_folding
'Tested on Altirra
'-----------------------------------------------------------
First, a bit about platforms. I wrote this in TurboBASIC XL using the parser tool noted above. I started on Atari because, well, it's Atari. Unfortunately, my beloved machine was just too slow. I was able to run it in real time using an 65816 accelerator chip clocked at 7 MHz. This is easy to try out using the Altirra emualtor but is not legal for the contest. After the recommendation from an online community to check out the BBC Micro, I found it could run at real time with a co-processor. The BBC Micro text mode is several times faster than Atari, the 6502 is clocked faster, and a second processor off loads the display and I/O functions. That made it work. I put together a line-by-line comparison of the TBXL source and Atari and BBC optimized codes. The rest of this description is made using the TBXL source because it keeps the original variable names.
DIM K$(1)
Need a string to hold the keyboard input. As seen below, I used the original controls from Flight Simulator II.
DT = .1 : REM (S) TIME STEP
The time step value is a critical parameter in numerical solutions to nonlinear differential equations. I had to adjust it to keep the aircraft stable. Unfortunately, 0.1 seconds is the largest I could go. At 0.4 seconds, the stock Atari would have run in realtime. Alas...
CLS : REM CLEAR THE SCREEN
?
?,"FLIGHT SIMULATOR 2D"
POKE 752,1 : REM TURN OFF CURSOR
This sets up the screen. It's in glorious 40-column text mode! I hope you are an instrument-rated simulator pilot because I can't afford the clock cycles to draw the horizon.
REPEAT
Here begins the game loop.
BASETIME = TIME
I record the start time of each loop iteration so I can display the loop rate relative to real time below.
'ATMOSPHERE
SIGMA=(1-Z*8.0E-05) : REM LINEAR APPROX FOR RELATIVE AIR DENSITY
Every flight simulator needs an atmosphere model. Mine is a linear approximation of the relative density of the US Standard Atmosphere. Note, all my calculations are in SI units, so Z is in meters.
STALL = ((U+2*FLAPS)>29) ! (Z<1) : REM DETERMINE IF NOT STALLED
FLTSIM2D is based on a Cessna 172 because I assumed I could find the most information online about this general aviation aircraft. Fortunately, a consulting firm published a
free analysis with many specifications and parameters on the Cessna 172. Among the specs is the stall speed with flaps up and down. This equation is close enough for the simulator. The flag
STALL is in negative logic: 0 if stalled; 1 if not stalled. My aircraft's stall speed is 29 m/s (56 knots) with flaps up and 23 m/s ( 51 knots) with flaps fully extended. (The variable
FLAPS is in units of 10 degrees.) Also, if the aircraft is very close to the ground, it is not stalled.
QSW = 8.1*(U*U+W*W)*1.225*SIGMA : REM DYNAMIC PRESSURE * WING AREA
The engineering of aircraft depends upon many coefficients (lift, drag, moment) that are multiplied by dynamic pressure and wing area to compute forces, which is computed by this equation. The velocity computes U and W are squared, but it is faster to self-multiply them. (I just realized while writing this that I did not combine the 8.1*1.225 into 9.923. Some wasted compute cycles!)
'ANGLE OF ATTACK(S)
UNOSING = 1/(U + (U=0)) : REM U != 0
SLOPE = W*UNOSING : REM SLOPE OF WIND
SLOPE = -(SLOPE<=-1) + (SLOPE>-1)*SLOPE : REM LIMIT SLOPE TO <+/-1
SLOPE = (SLOPE>=1) + (SLOPE<1)*SLOPE : REM LIMIT SLOPE TO <+/-1
ALPHA = SLOPE - .22 * SLOPE^3 : REM WING ANGLE OF ATTACK WRT WIND
This group of statements computes the angle-of-attack (AOA) of the wing. As seen below, the lift and drag coefficients vary with AOA. The AOA is computed by taking the arc-tangent of the slope of the wind velocity with respect to the wing. (This is done in the aircraft reference coordinate frame, where U and W are the longitudinal and transverse velocity components, respectively). The first line avoids a singularity (divide by zero) when U=0. I just force it to be 1, if necessary. This works fine on the ground while waiting to take off. It really wreaks havoc in the air if you find yourself is a severe stall. The next three lines find the slope and limit it to +/- 1. The limiting is done because I use a cubic polynomial approximation of arc-tangent instead of ATN() in the last line for speed. The ability to loop the loop was something I had to give up early on to make a fast simulator possible.
ALPHAT = ALPHA + OMEGA*UNOSING*4.3 + 8.33E-3*DLTA + .0863*CL - .0873 : REM TAILPLANE AOA. INCLUDES ELEVATOR AND DOWNWASH TERMS
This statement computes the AOA of the tailplane and packs in a lot of physics. There are four adjustments made to the wing AOA. In reverse order: the tailplane on the Cessna 172 is pitched 5 degrees (0.0873 radians) relative to the aircraft. Second, the wing creates a downwash on the tailplane, which acts to change the AOA. It turns out this term is incredibly important for aircraft stability. I don't know if this stability is real-life or simply numerical, but I tested the simulator with and without it and had to have it in here. I found the equation in another lecture and collapsed all the constants into the 0.0863 coefficient. The variable CL is the coefficient of lift computed on the previous loop iteration. The next term is the effect of the elevator position on tailplane AOA. This term lets the elevator control change the pitch of the aircraft and I found the equation in the Elevator Design chapter of an online book on aircraft design. This simulator is brought to you by the online generosity of university professors everywhere. The last term takes into account the relative motion of the tailplane with respect to the wind due to aircraft pitch rate. The tailplane is a stabilizing element of the aircraft - the wind is always trying to push it back to a neutral position.
'LIFT & DRAG COEFFICIENTS
CLT = .4*ALPHAT -.24*ALPHAT^3 : REM TAILPLANE LIFT COEFFICIENT WITH AREA RATIO.
This section begins computation of the all important lift and drag coefficients. First, the tailplane lift coefficient is modeled using a cubic polynomial valid for a wide range of AOA. The Cessna 172 uses a NACA 0012 airfoil, which is also used in wind turbines. The turbine engineers are concerned with the entire 360 degree range of AOA's. I noticed a sine function is a nice fit, although it does not include details of how the airfoil stalls at about 20 degree. But, hey, it's only a BASIC aircraft. I also factor in the wing-to-tailplane area ratio, which allows me to simplify calculation of the total lift and drag and pitching moment later on in the program.
CL = 0.3+0.16*FLAPS+4.8*ALPHA+12*ALPHA*ABS(ALPHA)-46*ALPHA^3 : REM WING LIFT COEFFICIENT. CUBIC FIT TO GET "DOUBLE HOOKS"
CL = CL * STALL : REM IF STALLED NO WING LIFT
Here, the wing lift coefficient is computed. It is a bit more complicated than the tailplane because the Cessna 172 wing is a NACA 2412 asymmetric airfoil. A completed third-order polynomial models this nicely. The sign of the squared term is flipped with negative AOA's using the ABS() function. The effect of flap position is also included with an additional constant term. The lift is finally nullified if the wing is below stall speed. I broke this up into two statements at one point to fit it in 10 lines and never removed it. It works.
CLL = CL+CLT : REM TOTAL A/C LIFT COEFFICIENT
CD = .025 + .0575*CLL*CLL : REM DRAG COEFFICIENT USING OSWALD EFFICIENCY
The total lift is the sum of the wing and tailplane lift. The drag coefficient is computed using the Oswald efficiency concept. Fortunately, I found this theory when I was needing to compact the code. Otherwise, I would need two more polynomial equations for the wing and tail drag coefficients. Oswald lets me use a single quadratic for both.
'ENGINE WITH DROPOFF DUE TO ALTITUDE
THRUST=(SIGMA-.05)*THROTTLE*UNOSING*1100 : REM INVERT PROPULSIVE POWER EQUATION
THRUST=THRUST*(THRUST<=2000)+2000*(THRUST>2000) : REM LIMIT THRUST TO 2000 N
The propulsion system was almost more difficult to figure out than the aerodynamic forces. At the start of the development, the thrust was just a scaled value of the throttle setting. This model, it turns out, is like having a jet turbine powered aircraft. It was convenient to have something that just worked, although I knew it was not a good model when the aircraft was at 25,000 ft and still climbing! (The Cessna 172 has a ceiling of about 14,000 ft.) I read and read about propeller and combustion engine propulsion. Most of the treatments are thermodynamic in nature. Work = Force x Distance; differentiating with respect to time yields Power = Force x Speed. To determine the thrust from applied power, just divide by the speed. Of course this equation blows up when the aircraft is stopped at the end of the runway waiting to take off. I learned that real flight simulators have dynamic models instead of thermodynamic models for propulsion . Engine throttle creates torque, which accelerates the propeller, which moves air and creates thrust. I found a great page on propeller efficiency vs. rotation rate (RPM), but was never able to arrive at a model that would determine RPM. My disposable time was running short and the contest deadline was approaching. Fortunately, the aircraft simulator book gave me "permission" to make power proportional to throttle (not entirely accurate), to invoke the thermodynamic model, and to limit the thrust arbitrarily. That is what the above two statements do. Oh, and the aircraft ceiling? That is made possible by engine power dropoff vs. altitude, the equation for which was gleaned from equation 7 on the propeller performance page: the engine torque (but I approximate with power) is proportional to relative air density. In the end, the parameters here could use some adjustment, but we'll just go with it. This model has 2000 N (450 lbs) of static thrust for acceleration down the runway and an altitude ceiling of a about 14,000 ft. There is probably too much power available because at low altitudes the aircraft can climb at 1000 fpm, which is too fast. Oh well. It works. Whew.
'TIME STEP EQUATIONS OF MOTION
MY = -QSW*(.0308 + CL*(.28-0.1*FLAPS) + CLT*4.3) : REM COMPUTE PITCHING MOMENT, + IS NOSE UP
Here begins the group of statements that compute the motion of the aircraft. This statement computes and sums the moments. The most difficult part was keeping track of the signs of the three terms - get one wrong and the aircraft is unstable and tumbles to a numerical death. The pitching moment calculation was made clear by equation 3.1 in
this lecture on static longitudinal stability and control. I hope I'm applying it appropriately. The three terms are the pitching moment caused by air flowing around the wing, the wing lift acting at the end of a short moment arm (whose length is changed by the flap position) and the lift of the tailplane acting at the end of the tail. These terms are multiplied by the dynamic-pressure wing-area product to compute the total moment in units of Newton meters.
OMEGA = OMEGA + MY*5.48e-5 : REM UPDATE PITCHING RATE WITH PITCHING MOMENT AND IYY
OMEGA = OMEGA * ((Z>.01) ! (OMEGA>0)) : REM NO ROTATION UNLESS POSITIVE WHEN ON GROUND
The moment causes an acceleration in the pitch direction and changes the angular velocity OMEGA about the pitch axis. The first statement combines DT and the moment of inertia IYY into a single constant. The IYY value was provided by these guys, which was really nice of them because I couldn't find it anywhere else. The second statement limits the rotation direction to nose up when the aircraft is on the ground. The next group of statements implement the force calculations and equations of motion for linear translation.
U = U + ( (THRUST-QSW*CD)*1E-3 - SINTH*9.81 - OMEGA*W ) * DT : REM UPDATE LONGITUDINAL VELOCITY
W = W + ( -QSW*CLL*1E-3 + COSTH*9.81 + OMEGA*U ) * DT : REM UPDATE TRANSVERSE VELOCITY
W = W * ((Z>.01) ! (W<0)) : REM NO VERTICAL VELOCITY IF ON GROUND, UNLESS ITS UP
The calculations are performed in the aircraft reference frame: U is the longitudinal velocity (along the axis of the aircraft) and W is the transverse (vertical with respect to the aircraft) velocity. There are four forces effecting U: thrust, drag, gravity and Coriolis. There are three forces accelerating the aircraft along W: lift, gravity and Coriolis. Gravity is projected onto U and W by the sine and cosine of the pitch (computed below). The Coriolis forces are present because the forces and velocities are computed in the rotating aircraft reference frame. The transverse velocity is nearly vertical and limited to the up direction when the aircraft is on the ground.
'INERTIAL FRAME UPDATE
T=T+DT : REM STEP TIME FORWARD
THETA = THETA + OMEGA*DT : REM UPDATE PITCH ANGLE
COSTH = 1 - 0.49*THETA*THETA : REM QUADRATIC APPROX TO COSINE
SINTH = THETA -.15*THETA^3 : REM CUBIC APPROX TO SINE
VX = (U*COSTH + W*SINTH) : REM UPDATE VELOCITY OVER GROUND
VZ = (U*SINTH - W*COSTH) : REM UPDATE VERTICAL RATE
X = X + VX*DT : REM UPDATE GROUND TRACK POSITION
Z = Z + VZ*DT : REMO UPDATE ALTITUDE
Z = Z * (Z>0) : REM RESTRICT Z TO ABOVE GROUND
This group of statements computes states in the inertial (Earth) reference frame. First time is stepped forward. The pitch angle THETA is updated. The cosine and sine values of pitch are computed using polynomial approximations. This code executes in about half the speed of calling the COS() and SIN() functions. With the sine and cosine values, the U and W velocities are rotated into X and Z coordinates, which are then updated with the time step. Finally, the altitude is limited to positive values only.
OKTOLAND = OKTOLAND ! (Z>9) : REM BE A LITTLE FORGIVING ON THE TAKEOFF
This flag is used to end the game; however, it is not set until the aircraft is at least 9 meters off the ground. The plane has a tendency to bounce on takeoff. To avoid the bounce being counted as a landing, you must get surely off the ground.
'PILOT INPUT. USE KEY INPUTS SAME AS SUBLOGIC FLIGHT SIM II
K$=INKEY$
THROTTLE = THROTTLE + 5*((K$="\1F")*(THROTTLE<100)-(K$="\1E")*(THROTTLE>0))
DLTA = DLTA + 0.5*((K$="T")*(DLTA<23)-(K$="B")*(DLTA>-28)) : REM ELEVATOR UP/DOWN
FLPIN = FLPIN + ((K$="N")*(FLPIN<3)-(K$="Y")*(FLPIN>0)) : REM FLAPS IN/OUT
FLAPS = 0.05*FLPIN+0.95*FLAPS : REM A LITTLE IIR EXPONENTIAL RESPONSE TO MODEL THE FLAP DRIVE MOTOR
The above is the input logic. FLTSIM2D uses a subset of the control keys used in Flight Simulator II. Fortunately, the manual was preserved on Atari Mania (among others), so I didn't have to go by a fading memory. The throttle uses left and right arrow keys. (The BBC Micro version uses ,< and .> keys). The elevator is adjusted using T and B and the flaps using Y and N. To slow the flap actuation and not cause a huge step function in lift and drag, the flap value is smoothed by an exponential filter. I tried smoothing other inputs, but it only seemed to exasperate the phugoid.
POSITION 2,3
?"KTAS",INT(U*1.94)," " : REM KNOTS
?"PITCH",INT(THETA*573)*.1," " : REM DEGREES WITH 1 DECIMAL PLACE
?"ALT.",INT(Z*3.28)," " : REM FEET
?"VRATE",INT(VZ*197)," " : REM FEET/MINUTE
?
?"ELEV.", DLTA," " : REM ELEVATOR POSITION
?"POWER", THROTTLE," " : REM % POWER
?"FLAPS ", INT(FLAPS*10+0.5)," " : REM FLAP POSITION
?
?"STALL ", CHR$(161-116*STALL); : REM STALL INDICATOR
?CHR$(253-221*STALL) : REM BEEP IF STALLING!
?
?"DIST.", INT(X*.9)," " : REM DISTANCE TRAVELED OVER GROUND
?"TIME",T,,INT(600/(TIME-BASETIME));"% " : REM TIME AND SIMULATOR RATE RELATIVE TO REAL TIME
Obviously, this is the output section with all the print statements. I grouped the outputs to try to follow the dashboard from Flight Sim. II: speed, pitch, altitude, vertical rate. he internal calculations are done in SI units, but the output are in customary imperial units. So meters per second are turned into knots or feet per minute. Meters are turned into feet or yards. Radians are turned into degrees. The next group are the control surface position readouts: elevator, throttle and flaps. Then there's a stall indicator, which uses the BELL character to sound the alarm. Finally, distance (GPS?) and time are displayed. The execution rate of the simulator with respect to real time is displayed as a percentage. The TIME function returns 1/60 second ticks (for NTSC video) on the Atari. The BBC Micro returns centisecond (0.01 second) ticks. Values >100% are faster than real time.
UNTIL OKTOLAND & (Z=0) : REM END GAME IF BACK ON GROUND
?
IF VZ >-2 : REM CHECK TO SEE IF VERTICAL RATE IS SLOW ENOUGH
? "TOUCHDOWN"
ELSE
? "YOU CRASHED!\FD\FD\FD" : REM OH NO! VERTICAL TOO FAST!
ENDIF
