The maths behind the Spitfire's aerodynamic flight simulation
The sophisticated simulation of the Spitfire's flight model is one of the most lauded features of Aviator, and rightly so. The heart of the flight model is the ApplyFlightModel routine, where the aerodynamic forces on the plane are calculated and applied, all using 8-bit maths and 3D geometry. It's impressive stuff, so let's take a look at what's involved.
Before tackling the following, you might find it useful to read the deep dive on rotation matrices. For more information on the flight model calculations themselves, check out the deep dive on matching the code to the flight model, which takes the flight model's description from the manual and matches it up with the code.
To keep things as simple as possible, this deep dive takes a look at how the flight model calculations work for normal flight, i.e. when the plane is flying normally, and is not stalling, trying to land or experiencing alien-induced turbulence. For other aspects of the flight model, see the deep dives on stalling and recovery, ripping the wings off, take-offs and landings and on-ground calculations.
The aim of the flight model
---------------------------
The complex calculations of Aviator's flight model produce a relatively simple set of adjustments that we apply to the plane's velocity, turn rate, orientation and location. This updating process is run on each iteration of the main loop, resulting in a Spitfire that flies through the air in a mathematically accurate manner.
Specifically, the process takes all the relevant environmental conditions into consideration, and calculates the following:
- A vector (dxVelocity dyVelocity dzVelocity) that we add to the plane's velocity:
[ xVelocity ] [ xVelocity ] [ dxVelocity ] [ yVelocity ] = [ yVelocity ] + [ dyVelocity ] [ zVelocity ] [ zVelocity ] [ dzVelocity ]
- A vector (dxTurn dyTurn dzTurn) that we add to the plane's turn rate:
[ xTurn ] [ xTurn ] [ dxTurn ] [ yTurn ] = [ yTurn ] + [ dyTurn ] [ zTurn ] [ zTurn ] [ dzTurn ]
- A vector (dxRotation dyRotation dzRotation) that we add to the plane's orientation:
[ xRotation ] [ xRotation ] [ dxRotation ] [ yRotation ] = [ yRotation ] + [ dyRotation ] [ zRotation ] [ zRotation ] [ dzRotation ]
- The plane's new location in space, which is simply the current location plus the updated velocity vector:
[ xPlane ] [ xPlane ] [ xVelocity ] [ yPlane ] = [ yPlane ] + [ yVelocity ] [ zPlane ] [ zPlane ] [ zVelocity ]
The challenge, then, is to calculate the three (dx dy dz) vectors. Let's see what's involved.
The flight model's forces
-------------------------
The ApplyFlightModel routine and its subroutines calculate a number of different forces while calculating the three (dx dy dz) vectors. Some of these are linear forces, while others are angular forces; it's important to understand the difference.
- Linear forces are pretty easy to understand, as they act in a straight line. For example, the engine pushes the plane along in a straight line, so it's a linear force. Each axis in a linear force's vector represents the linear force along that axis.
- Angular forces are a bit more complicated, as they effectively rotate the object to which they apply around its centre of gravity. For example, the ailerons apply angular forces to the wings that rotate the plane around the z-axis, rolling the plane (as the z-axis points forward out of the plane's nose). These angular forces are also known as moments, or torque, and each axis in an angular force's vector represents the rotational force around that axis, with a positive force rotating clockwise, and a negative force rotating anti-clockwise.
Given this, here's a list of the forces that the flight model calculates, starting with the linear forces:
- (xLiftDrag yLiftDrag zLiftDrag)
- Feeds into the (dxVelocity dyVelocity dzVelocity) calculation
- Describes the lift, drag and side forces caused by the plane moving through the air
- xLiftDrag moves the plane right/left
- yLiftDrag moves the plane up/down
- zLiftDrag moves the plane forwards/backwards
- (0 yFlapsLift 0)
- Feeds into the dyVelocity calculation
- Describes the lift that's provided by the flaps when they are on
- The only non-zero component is in the y-axis as this represents lift
- The x- and z-axis components are always zero, so we don't bother to store them
- yFlapsLift lifts the plane depending on the forward velocity and the status of the flaps
- (xLinear yLinear zLinear)
- Feeds into the (dxVelocity dyVelocity dzVelocity) calculation
- Is used to calculate the combined effect of the above two forces, plus the effect of the engine
And here are the angular forces used in the calculation:
- (xMoments yMoments zMoments)
- Feeds into the (dxTurn dyTurn dzTurn) calculation
- Describes the angular forces caused by the plane moving through the air
- xMoments pitches the plane
- yMoments yaws the plane
- zMoments rolls the plane
- (xControls yControls zControls)
- Feeds into the (dxTurn dyTurn dzTurn) calculation
- Describes the effect of the flight controls
- xControls pitches the plane (elevator)
- yControls yaws the plane (rudder)
- zControls rolls the plane (ailerons)
- (0 0 zSlipMoment)
- Feeds into the dzTurn calculation
- Describes the rolling moment due to wing forces resulting from aircraft slip
- The only non-zero component is in the z-axis, as this represents a roll
- The x- and y-axis components are always zero, so we don't bother to store them
- zSlipMoment rolls the plane depending on side force/velocity
Let's now look at the various stages in calculating these forces, and how they get combined to work out the effect on the plane. The links will take you to the relevant parts of the source code, or if you want to match the following calculations to the description of the flight model in the game's manual, see the deep dive on matching the code to the flight model.
Calculation 1: Linear forces
----------------------------
See code: ApplyFlightModel (Part 2 of 7)
We start by rotating the plane's velocity vector from the outside world's frame of reference to the plane's frame of reference, using matrix 1:
[ xVelocityP ] [ xVelocity ] [ yVelocityP ] = matrix1 x [ yVelocity ] [ zVelocityP ] [ zVelocity ]
See code: ApplyAerodynamics (Part 1 of 3), ApplyAerodynamics (Part 3 of 3)
Next, we calculate the linear forces due to the lift, drag and side forces. These three forces are:
- x-axis: Side forces due to aircraft slip
- y-axis: Wing lift
- z-axis: Aircraft drag (parasitic)
and they are calculated as follows:
[ xLiftDrag ] = [ xVelocityP * 2 * 4 ] [ yLiftDrag ] = [ yVelocityP * 2 * 4 ] * maxv * airDensity [ zLiftDrag ] = [ zVelocityP * 2 ]
where:
airDensity = ~yPlaneHi * 2 maxv = max(|xVelocityP|, |yVelocityP|, |zVelocityP|)
The airDensity factor reduces the effect of all the aerodynamic-related forces at higher altitudes, due to the lower density of the air (which therefore has less effect when interacting with the plane).
Finally, if zLiftDrag is positive then air is passing over the wings in a front-to-back direction, so we also do the following calculation that gives us the lift due to the flaps (though if the flaps are not extended, then this value is reduced to zero by the scaling process in calculation 3, so we only get lift from the flaps when they are on):
[ xFlapsLift ] [ 0 ] [ yFlapsLift ] = [ zLiftDrag ] [ zFlapsLift ] [ 0 ]
As mentioned above, we only store yFlapsLift in this calculation because xFlapsLift and zFlapsLift are always zero.
We now have the basic linear forces that act on the plane.
Calculation 2: Angular forces
-----------------------------
See code: ApplyAerodynamics (Part 1 of 3), ApplyAerodynamics (Part 3 of 3)
Next we calculate the basic angular forces, starting with the moments:
[ xMoments ] = [ yVelocityP * 2 - xTurn * 250 / 256 ] [ yMoments ] = [ xVelocityP * 2 - yTurn * 250 / 256 ] * maxv * airDensity [ zMoments ] = [ -zTurn * 2 ]
The above calculates the following angular forces:
- x-axis: Pitching moment due to tail forces, due to both linear motion and pitch rate
- y-axis: Yawing moment due to aircraft side, fin and rudder forces, due to both aircraft slip and yaw rate
- z-axis: Rolling moment due to wing forces resulting from roll rate
Then we calculate the rolling moment due to wing forces resulting from aircraft slip (xLiftDrag being the linear side force that we calculated above):
[ xSlipMoment ] [ 0 ] [ ySlipMoment ] = [ 0 ] [ zSlipMoment ] [ xLiftDrag ]
As mentioned above, we only store zSlipMoment in the middle calculation because xSlipMoment and ySlipMoment are always zero.
See code: ApplyFlightControl
Next we calculate the angular forces from the plane's controls:
[ xControls ] = [ zLiftDrag * elevatorPosition ] [ yControls ] = [ zLiftDrag * rudderPosition ] [ zControls ] = [ zLiftDrag * aileronPosition ]
The above calculates the following angular forces:
- x-axis: Pitching moment due to elevator deflection
- y-axis: Yawing moment due to rudder demand
- z-axis: Rolling moment due to aileron deflection
See code: ApplyAerodynamics (Part 1 of 3), ApplyAerodynamics (Part 3 of 3)
Finally, if zLiftDrag is positive then air is passing the plane in a front-to-back direction, so we also do the following to calculate the pitching moment due to centre of pressure, which is affected by the position of the flaps:
xMoments = xMoments + (zLiftDrag * 8) when the flaps are off
xMoments - (zLiftDrag * 4) when the flaps are on
We now have the basic angular forces that act on the plane.
Calculation 3: Scale the forces
-------------------------------
See code: ScaleFlightForces
Next we scale all the forces we've calculated, first by the relevant force factor in forceFactor, and then by the relevant scale factor in scaleFactor (each component of each force has its own set of factors). The scaling process uses this formula:
scaledForce = unscaledForce * forceFactor * 2 ^ scaleFactor
In a sense, this is where the Spitfire's unique characteristics are encoded; the scale factors are the coefficients in the equation, and different planes would have different coefficients. When Geoff Crammond talks about writing Aviator, he talks about getting hold of a pilot's manual and other data for the Spitfire that he fed into the simulation, and I suspect that this is where a lot of that data is encapsulated.
Specifically, the forces are scaled as follows (in normal flight, with the undercarriage up, the flaps off and the engine on):
[ xLiftDragSc ] [ xLiftDrag ] 176 * 16 [ yLiftDragSc ] = [ yLiftDrag ] * 156 * 16 (39 * 16 when stalling) [ zLiftDragSc ] [ zLiftDrag ] 62 (+/- uc, flaps, engine) [ xFlapsLiftSc ] [ 0 ] n/a [ yFlapsLiftSc ] = [ yFlapsLift ] * 0 (152 * 4 when flaps on) [ zFlapsLiftSc ] [ 0 ] n/a [ xMomentsSc ] [ xMoments ] 212 [ yMomentsSc ] = [ yMoments ] * 201 / 4 [ zMomentsSc ] [ zMoments ] 204 / 2 [ xSlipMomentSc ] [ 0 ] n/a [ ySlipMomentSc ] = [ 0 ] * n/a [ zSlipMomentSc ] [ zSlipMoment ] 40 / 32 [ xControlsSc ] [ xControls ] 255 / 2 [ yControlsSc ] = [ yControls ] * 141 / 4 [ zControlsSc ] [ zControls ] 190 / 4
Note that these are all scalar multiplications, not matrix multiplications.
The force factors in the above list are for normal flight, but three of the factors for the linear forces change dynamically with the flying conditions. Specifically, this is how these three factors work:
- The force factor for yLiftDrag is:
- 156 * 16 in normal flight
- 39 * 16 when the plane is stalling
- The force factor for zLiftDrag is:
- 42 when everything that might cause drag is off, i.e. undercarriage up, flaps off and engine off
- +10 when the undercarriage goes down
- +200 when the flaps go on
- +20 when the engine is switched on
- So the factor is 52 when starting on the runway with the undercarriage down, flaps off and engine off
- And the factor is 62 in normal flight, with the undercarriage up, flaps off and engine on
- The force factor for yFlapsLift is:
- 0 when the flaps are off
- 152 * 4 when the flaps are on
All the other flight factors are constant.
If the plane is experiencing turbulence from an exploding alien, then additional calculations are implemented at this point - see the deep dive on explosions and turbulence for details.
Calculation 4: Effect on the turn rate
--------------------------------------
We now start the process of combining the various forces we've calculated, so we can apply them to the plane's turn rate, velocity and rotation.
See code: ApplyFlightModel (Part 1 of 7)
First, we set a gravity vector that pulls down along the y-axis, and rotate it from the outside world's frame of reference to the plane's frame of reference, using matrix 4:
[ xGravity ] [ 0 ] [ yGravity ] = matrix4 x [ -512 ] [ zGravity ] [ 0 ]
See code: ApplyTurnAndThrust (Part 1 of 2)
Next, we calculate the (dxTurn dyTurn dzTurn) vector, which contains the change in the plane's turn rate. We do this as follows:
[ dxTurn ] [ xMomentsSc ] [ xControlsSc ]
[ dyTurn ] = [ yMomentsSc ] + [ yControlsSc ]
[ dzTurn ] [ zMomentsSc ] [ zControlsSc ]
[ yGravity ] [ 0 ]
+ [ 0 ] - [ 0 ]
[ 0 ] [ zSlipMoment ]
This calculation adds the moments due to flight and the moments due to the controls, and subtracts the rolling moment due to wing forces resulting from aircraft slip, all of which we calculated above. It also adds in another angular force: a pitching moment due to the centre of gravity, which takes the gravitational force in yGravity and applies it as an angular force around the x-axis, giving us our pitching moment.
Calculation 5: Engine thrust
----------------------------
See code: ApplyTurnAndThrust (Part 2 of 2)
Next we create a vector in (xLinear yLinear zLinear) that combines all the various linear forces we calculated above into one vector, and then we add in the linear forces from the engine (but only if we are not going too fast). The calculation is as follows, starting with the high-speed calculation, and then the normal speed calculation.
If zVelocityPHi >= 48 (so forward speed >= 500 mph), then we calculate:
[ xLinear ] [ 0 ] [ xLiftDragSc ] [ yLinear ] = [ yFlapsLiftSc ] - [ yLiftDragSc ] [ zLinear ] [ 0 ] [ 5632 ]
Note that the 5632 number is actually (&EA zLinearLo), which is in the range -5632 and -5377 depending on the value of zLinearLo, though as it's such a large number, the low byte isn't that important. The effect is that going really fast puts a huge strain on the plane, resulting in a huge force pressing against the nose and the front edges of the wings and tail.
If zVelocityPHi < 48 (so forward speed < 500 mph), we calculate:
[ xLinear ] [ 0 ] [ xLiftDragSc ] [ 0 ] [ yLinear ] = [ yFlapsLiftSc ] - [ yLiftDragSc ] + [ 0 ] [ zLinear ] [ 0 ] [ zLiftDragSc ] [ zEngine ]
where zEngine is 0 if the engine is off, or the following if the engine is on:
zEngine = max(0, thrustScaled - (max(0, zVelocityP) / 16))
* airDensity / 512
where thrustScaled is the current thrust level in (thrustHi thrustLo), but:
- Doubled if thrust >= 1024 (i.e. the throttle is between 80% and 100%)
- Doubled if zVelocity is in the range 512 to 1023 (i.e. when the plane is rising vertically between 21.6 to 43.2 mph)
In effect, the above calculation takes the engine's thrust, reduces it with higher forward speeds, and reduces it again by the standard high-altitude effect. In the common case that the plane is travelling forwards and the engine is on, it simplifies to:
zEngine = (thrustScaled - (zVelocityP / 16)) * airDensity / 512
In this case it's a bit easier to see the effect of the current velocity and altitude on the engine force.
We now have the sum of all the linear forces on the plane, along with the engine forces. If the plane is on the ground, additional calculations are implemented at this point - see the deep dive on on-ground calculations for details.
Calculation 6: Effect on velocity
---------------------------------
See code: ApplyFlightModel (Part 6 of 7), AdjustVelocity
Next we calculate the (dxVelocity dyVelocity dzVelocity) vector, which contains the change in the plane's velocity. We do this by taking (xLinear yLinear zLinear), which contains all the linear forces on the plane, and rotating it from the plane's frame of reference to the outside world's frame using matrix 2.
We then subtract 16 from yLinearHi (which subtracts 4096 from yLinear) to apply the force of gravity, before doubling the result:
[ dxVelocity ] ( [ xLinear ] [ 0 ] ) [ dyVelocity ] = 2 * ( matrix2 x [ yLinear ] - [ 4096 ] ) [ dzVelocity ] ( [ zLinear ] [ 0 ] )
In all, this applies the linear forces to the plane, including gravity, which alter its velocity.
Calculation 7: Effect on orientation
------------------------------------
See code: AdjustTurn
We calculate the (dxRotation dyRotation dzRotation) vector by first updating the (xTurn yTurn zTurn) vector to its new value:
[ xTurn ] [ xTurn ] [ dxTurn ] [ yTurn ] = [ yTurn ] + [ dyTurn ] [ zTurn ] [ zTurn ] [ dzTurn ]
See code: ApplyFlightModel (Part 6 of 7), AdjustRotation
We then calculate the change in the plane's orientation by rotating the result by the reverse of the current roll angle using matrix 3:
[ dxRotation ] [ xTurn ] [ dyRotation ] = matrix3 x [ yTurn ] [ dzRotation ] [ zTurn ]
This alters the turn rate of the plane and calculates how we should be changing the orientation of the plane, based on the change in the turn rate that we already calculated in (dxTurn dyTurn dzTurn).
Calculation 8: Update the plane with the results
------------------------------------------------
See code: ApplyFlightModel (Part 6 of 7)
We now have all our (dx dy dz) vectors, so we can update the plane's location, velocity and orientation, giving us our final result (we updated the (xTurn yTurn zTurn) vector in the previous step):
[ xVelocity ] [ xVelocity ] [ dxVelocity ] [ yVelocity ] = [ yVelocity ] + [ dyVelocity ] [ zVelocity ] [ zVelocity ] [ dzVelocity ] [ xRotation ] [ xRotation ] [ dxRotation ] [ yRotation ] = [ yRotation ] + [ dyRotation ] [ zRotation ] [ zRotation ] [ dzRotation ] [ xPlane ] [ xPlane ] [ xVelocity ] [ yPlane ] = [ yPlane ] + [ yVelocity ] [ zPlane ] [ zPlane ] [ zVelocity ]
Note that the plane is flipped over if the rotation calculation makes it point backwards.
And that is Aviator's flight model - an impressive bit of coding, particularly given that this is all on an 8-bit computer. To see how the above calculations map to the description of the flight model in the manual, see the deep dive on matching the code to the flight model, and for other aspects of the flight model, see the deep dives on stalling and recovery, ripping the wings off, take-offs and landings and on-ground calculations.