F-16 Subsonic aircraft aero inertial, and autopilot models

Goal

Provide flight dynamic aero, propulsion, mass properties, guidance, stability and control subsystems for a maneuvering F-16 subsonic vehicle in order to generate flight simulation verification time history data.

Acronyms

ANL: Aircraft nose left
ANR: Aircraft nose right
ANU: Aircraft nose up
AP: Autopilot
CCW: Counter-clockwise
CM: Center-of-mass
CWFN: Clockwise from North
DOF: Degrees of freedom
FORTRAN: FORMula TRANslator
KEAS: Knots equivalent airspeed
IDL: International Date Line
LWD: Left wing down
nm: nautical mile
MAC: Mean aerodynamic chord
RWD: Right wing down
SAS: Stability Augmentation System
TED: Trailing edge down
TEL: Trailing edge left

Assumptions and Limitations

The various models are specified in accordance with AIAA S-119 [S119], also known as DAVE-ML. See http://daveml.org for more information.

Flight envelope is limited to the vicinity of 10,000 ft MSL and 287.8 knots equivalent airspeed (around Mach 0.5).

Scope of effort

A representative but unclassified and open-source flight vehicle model is provided for use in generating or comparing subsonic maneuvering in atmospheric flight, in a non-proprietary format that is both human- and text-readable.

The S-119 format is designed for unambiguous representation of flight dynamic models in an archivable format.

To this end, separate system models of a single-engine fighter plane are provided in the following files:

F16_aero.dml, mod O, dated 2013-09-11:: Subsonic aerodynamics model
F16_prop.dml, dated 2012-08-07:: Simple propulsion model
F16_inertia.dml, dated 2012-08-07:: Mass properties model
F16_control.dml, rev E, dated 2013-09-12:: SAS and three-axis autopilot model
F16_gnc.dml, rev B, dated 2013-09-12:: Maneuver generating autopilot

It is intended that either the F16_control.dml or F16_gnc.dml models be used, depending on the requirements for trajectory generation. The aerodynamic, propulsion and inertia models should always be included.

Components

Each component model is described in detail below.

F16_aero.dml

This model is based on [TM212145] by Garza and Morelli, which implements aero data published in [TP1538] by Nguyen. The Garza model was written in Matlab®, a proprietary model format from The Mathworks, Inc. that has been re-written to avoid proprietary formats. It encodes a non-linear subsonic flight vehicle model with data table aero coefficients and is representative of a moderate-fidelity engineering or flight training simulation model.

Table 1. Input Signals - *F16_aero.dml*
Name	Units	Sign	Description [nominal FTV value]
`trueAirspeed`	ft/s	>0	True airspeed, ft per second
`angleOfAttack`	deg	+ANU	Angle of attack
`angleOfSideslip`	deg	+ANL	Angle of sideslip
`bodyAngularRate_Roll`	rad/s	+RWD	Body roll rate
`bodyAngularRate_Pitch`	rad/s	+ANU	Body pitch rate
`bodyAngularRate_Yaw`	rad/s	+ANR	Body yaw rate
`elevatorDeflection`	deg	+TED	Elevator deflection
`aileronDeflection`	deg	+LWD	Aileron deflection (right - left)/2
`rudderDeflection`	deg	+TEL	Rudder deflection

Table 2. Output Signals - *F16_aero.dml*
Name	Units	Sign	Description
`referenceWingChord`	ft	const.	Longitudinal reference length
`referenceWingSpan`	ft	const.	Lateral reference length
`referenceWingArea`	ft^2	const.	Reference area
`aeroBodyForceCoefficient_X`	n.d.	+fwd	Total coefficient of force along the body X-axis acting about the moment reference center
`aeroBodyForceCoefficient_Y`	n.d.	+right	Total coefficient of force along the body Y-axis acting about the moment reference center
`aeroBodyForceCoefficient_Z`	n.d.	+down	Total coefficient of force along the body Z-axis acting about the moment reference center
`aeroBodyMomentCoefficient_Roll`	n.d.	+RWD	Total coefficient of moment around the moment reference center in the body X-axis (rolling moment)
`aeroBodyMomentCoefficient_Pitch`	n.d.	+ANU	Total coefficient of moment around the moment reference center in body Y-axis (pitching moment)
`aeroBodyMomentCoefficient_Yaw`	n.d.	+ANR	Total coefficient of moment around the moment reference center in the body Z-axis (yawing moment)

F16_prop.dml

This model contains data tables and equations found in [Stevens03] that were originally realized in FORTRAN. It models thrust as a function of Mach, altitude, and power lever angle; if the commanded power lever angle is greater than 50, augmented (afterburner) thrust is included.

Table 3. Input Signals - *F16_prop.dml*
Name	Units	Sign	Description [nominal FTV value]
`powerLeverAngle`	pct	0-100	Throttle (power lever angle), 0 - 100. 50 is MIL (max dry) thrust; 100 is MAX (burner)
`altitudeMSL`	ft	+up	Geometric altitude above reference geoid
`mach`	n.d.	>0	Flight Mach number

Table 4. Output Signals - *F16_prop.dml*
Name	Units	Sign	Description
`thrustBodyForce_X`	lbf	+fwd	Steady-state thrust of engine for given inputs in body X axis
`thrustBodyForce_Y`	lbf	+right	Steady-state thrust of engine for given inputs in body Y axis (always 0)
`thrustBodyForce_Z`	lbf	+down	Steady-state thrust of engine for given inputs in body Z axis (always 0)
`thrustBodyMoment_Roll`	ft-lbf	+RWD	Steady-state moment of engine for given inputs about the body X axis (always 0)
`thrustBodyMoment_Pitch`	ft-lbf	+RWD	Steady-state moment of engine for given inputs about the body Y axis (always 0)
`thrustBodyMoment_Roll`	ft-lbf	+RWD	Steady-state moment of engine for given inputs about the body Z axis (always 0)

F16_inertia.dml

This models a simple constant-mass condition with fixed values for products and cross-products of inertia, from [Stevens03]. For this test, the input value for the CG position should be 25% (a floating point value of 25.0, as the input is in percent).

Table 5. Input Signals - *F16_inertia.dml*
Name	Units	Sign	Description [nominal FTV value]
`vrsPositionOfCM`	pct	+aft	Longitudinal location of the center of mass in percent of MAC. Input value (CG location) for these tests should be 25%.

Table 6. Output Signals - *F16_inertia.dml*
Name	Units	Sign	Description
`bodyMomentOfInertia_Roll`	slug-ft^2	>0	Rolling moment of inertia about the body X axis
`bodyMomentOfInertia_Pitch`	slug-ft^2	>0	Pitching moment of inertia about the body Y axis
`bodyMomentOfInertia_Yaw`	slug-ft^2	>0	Yawing moment of inertia about the body Z axis
`bodyProductOfInertia_ZX`	slug-ft^2	normal	Cross-product of inertia in the body X-Z plane (no sign reversal)
`bodyProductOfInertia_XY`	slug-ft^2	normal	Cross-product of inertia in the body X-Y plane (no sign reversal)
`bodyProductOfInertia_YZ`	slug-ft^2	normal	Cross-product of inertia in the body Y-Z plane (no sign reversal)
`totalMass`	slug	>0	Total mass of the vehicle (20,500 lbm)
`bodyPositionOfCmWrtMrc_X`	ft	+fwd	Longitudinal location of the center of mass relative to moment reference point
`bodyPositionOfCmWrtMrc_Y`	ft	+right	Lateral location of the center of mass relative to moment reference point
`bodyPositionOfCmWrtMrc_Z`	ft	+down	Vertical location of the center of mass relative to moment reference point

F16_control.dml

This model contains optional stability augmentation for piloted flight as well as a one, two, or three-axis autopilot with course capture-and-track capability. It is the basis for the more specialized F16_gnc.dml that follows, but can be used as-is for arbitrary maneuvering.

Stability augmentation is provided by two separate linear quadratic regulator, full-state-feedback gain matrices developed by the author that are optimized for flight at 10,000 ft and 287 KEAS. Flight at other speed/altitude combinations may be sub-optimal or even unstable.

Pilot inputs are converted to surface and power lever angle commands. If stabilityAugmentationOn_disc is true, the vehicle’s dynamics will be stabilized. If autopilotOn_disc is true, stability augmentation will be engaged and the vehicle will respond to airspeed, altitude, and heading commands (equivalentAirspeedCommand, altitudeMslCommand and trueBaseCourseCommand which acts as a desired heading input). If lateralDeviationError is non-zero, the vehicle will attempt to intercept and track a desired course along the direction supplied in trueBaseCourseCommand.

While trimming the F-16 model to an equilibrium state, stabilityAugmentationOn_disc and autopilotOn_disc should both be set to false. It is recommended that the simulation’s trim feature manipulate trimmedPilotControl_throttle and trimmedPilotControl_long to trim the longitudinal state of the vehicle. If the simulation’s trim feature is designed to manipulate pilotControl_throttle and pilotControl_long, then the simulation should set trimmedPilotControl_throttle and trimmedPilotControl_long to zero. (The F16_control.dml file specifies a non-zero initial condition for these two trim variables that is near their trim values for straight and level flight at 10,000 ft MSL and Mach 0.5.)

Table 7. Input Signals - *F16_control.dml*
Name	Units	Sign	Description [nominal FTV value]
`pilotControl_throttle`	0 → +1	+incr	Pilot throttle control position
`pilotControl_long`	-1 → +1	+aft	Pilot longitudinal control position
`pilotControl_lat`	-1 → +1	+right	Pilot lateral control position
`pilotControl_yaw`	-1 → +1	+right	Pilot rudder pedal position
`trimmedPilotControl_throttle`	0 → +1	+incr	Trimmed position of throttle
`trimmedPilotControl_long`	-1 → +1	+aft	Trimmed position of pitch control
`stabilityAugmentationOn_disc`	0, 1	+true	Stability augmentation engage flag (discrete)
`autopilotOn_disc`	0, 1	+true	Autopilot engage flag (discrete)
`equivalentAirspeedCommand`	KEAS	>0	Desired equivalent airspeed (autopilot input)
`altitudeMslCommand`	ft	>0	Desired absolute altitude (autopilot input)
`lateralDeviationError`	ft	+right	Lateral deviation error from desired course for autopilot
`trueBaseCourseCommand`	deg	+CWFN	True heading of desired ground track (autopilot input)
`altitudeMsl`	ft	+up	Geometric altitude above mean sea level (SAS and autopilot feedback)
`equivalentAirspeed`	KEAS	+up	Equivalent airspeed (SAS and autopilot feedback)
`angleOfAttack`	deg	+ANU	Angle of attack (SAS feedback)
`angleOfSideslip`	deg	+ANL	Angle of sideslip (SAS feedback)
`eulerAngle_Roll`	deg	+RWD	Roll angle (SAS feedback)
`eulerAngle_Pitch`	deg	+ANU	Pitch angle (SAS feedback)
`eulerAngle_Yaw`	deg	+CWFN	Heading angle (autopilot feedback)
`bodyAngularRate_Roll`	rad/s	+RWD	Body roll rate (SAS feedback)
`bodyAngularRate_Pitch`	rad/s	+ANU	Body pitch rate (SAS feedback)
`bodyAngularRate_Yaw`	rad/s	+ANR	Body yaw rate (SAS feedback)

Table 8. Output Signals - *F16_control.dml*
Name	Units	Sign	Description
`elevatorDeflection`	deg	+TED	Elevator command
`aileronDeflection`	deg	+LWD	Aileron deflection (right - left)/2
`rudderDeflection`	deg	+TEL	Rudder deflection
`powerLeverAngle`	pct	0-100	Throttle (power lever angle), 0 - 100. 50 is MIL (max dry) thrust; 100 is MAX (burner)

F16_gnc.dml

This model can be used as a replacement for the control law described in F16_control.dml and is designed to perform a 3-nm CCW circular orbit around one of two spots on the Earth: the North Pole or around the equator/international data line intersection, for the purposes of simulation checkout and verification.

The primary difference between F16_gnc.dml and F16_control.dml is the GNC version has no trueBaseCourseCommand; the behavior of the vehicle when autopilotOn_disc is true depends on the value of selectCircumnavigator_disc: if 0, the vehicle will fly to and circle around the equator/international date line intersection; if 1, the vehicle will fly to and circle the North Pole at a distance of 3 nautical miles at the specified altitude and speed.

These circles take approximately 3.5 minutes at 10,000 ft and 287 KEAS.

See the description of F16_control.dml for notes on manipulating the control system to trim the F-16 model.

Table 9. Input Signals - *F16_gnc.dml*
Name	Units	Sign	Description [nominal FTV value]
`pilotControl_throttle`	0 → +1	+incr	Pilot throttle control position
`pilotControl_long`	-1 → +1	+aft	Pilot longitudinal control position
`pilotControl_lat`	-1 → +1	+right	Pilot lateral control position
`pilotControl_yaw`	-1 → +1	+right	Pilot rudder pedal position
`trimmedPilotControl_throttle`	0 → +1	+incr	Trimmed position of throttle
`trimmedPilotControl_long`	-1 → +1	+aft	Trimmed position of pitch control
`stabilityAugmentationOn_disc`	0, 1	+true	Stability augmentation engage flag (discrete)
`autopilotOn_disc`	0, 1	+true	Autopilot engage flag (discrete)
`selectCircumnavigator_disc`	0, 1	+Npole	Selects point to circle; 0 = equator/IDL, 1 = North Pole
`geLatitude`	deg	+north	Geodetic latitude of vehicle’s center of mass (nav feedback)
`geLongitude`	deg	+east	Geodetic longitude of vehicle’s center of mass east of PM (nav feedback)
`equivalentAirspeedCommand`	KEAS	>0	Autopilot commanded equivalent airspeed (autopilot input)
`altitudeMslCommand`	ft	>0	Autopilot commanded absolute altitude (autopilot input)
`altitudeMsl`	ft	+up	Geometric altitude above mean sea level (SAS and autopilot feedback)
`equivalentAirspeed`	KEAS	+up	Equivalent airspeed (SAS and autopilot feedback)
`angleOfAttack`	deg	+ANU	Angle of attack (SAS feedback)
`angleOfSideslip`	deg	+ANL	Angle of sideslip (SAS feedback)
`eulerAngle_Roll`	deg	+RWD	Roll angle (SAS feedback)
`eulerAngle_Pitch`	deg	+ANU	Pitch angle (SAS feedback)
`eulerAngle_Yaw`	deg	+CWFN	Heading angle (autopilot feedback)
`bodyAngularRate_Roll`	rad/s	+RWD	Body roll rate (SAS feedback)
`bodyAngularRate_Pitch`	rad/s	+ANU	Body pitch rate (SAS feedback)
`bodyAngularRate_Yaw`	rad/s	+ANR	Body yaw rate (SAS feedback)

Table 10. Output Signals - *F16_gnc.dml*
Name	Units	Sign	Description
`elevatorDeflection`	deg	+TED	Elevator command
`aileronDeflection`	deg	+LWD	Aileron deflection (right - left)/2
`rudderDeflection`	deg	+TEL	Rudder deflection
`powerLeverAngle`	pct	0-100	Throttle (power lever angle), 0 - 100. 50 is MIL (max dry) thrust; 100 is MAX (burner)

Implementation

The AIAA S-119 standard is an XML-based encoding scheme that provides both human- and machine-readable algorithm and data table specifications in an unambiguous way, suitable for exchange over distance and, for archival purposes, time. While specifying models in this fashion may initially present a challenge to the implementer, investment in automated scripts or parsers should greatly simplify future exchange of models. The nature of XML should make these models understandable even without reference to the formal grammar definition document (available from http://daveml.org).

The grammar is simple enough that equivalent source code can be implemented by hand in most cases. If preferred, several parsers are available, some open and some by request. See the Tools section of http://daveml.org for a current list.

These parsers are intended either for use at run-time (an interpreting parser) or at compile-time (a translating parser that emits C, Simulink®, FORTRAN or other programming source code). A parser that can do both is also available.

S-119 models may also contain their own self-verification check-cases that will help ensure proper implementation of the model by the chosen parser. At present, the aero (F16_aero.dml) and propulsion (F16_prop.dml) models contain verification checkcase definitions and data.

Variable names used by these various models, where appropriate, reflect the suggested naming convention described in [S119].

Results

Using the models cited above, a six-DOF non-linear aircraft simulation was realized in a popular dynamic simulation and control design tool with an oblate rotating spheroidal Earth with a constant 32.174 ft/s^2 gravity field and a tabular US 1976 Standard Atmosphere model. The center of mass was set at 25% (the single input to the inertia model) and the aircraft was trimmed for level flight, and then various modes of autopilot were engaged and recorded; these results are given below for the convenience of other implementers.

Trimmed flight conditions

It is important to note that for proper trim, the summed aero and propulsion forces and moments obtained from the aero (F16_aero.dml) and engine (F16_prop.dml) models must be transferred from the moment reference center (35% of MAC behind the wing leading edge) to the center of mass obtained from the mass properties model (F16_inertia.dml).

The trim scheme employed was for the values given below were obtained by trimming for wings-level, horizontal, un-accelerated flight using elevator, throttle, and pitch attitude/angle of attack (since flight path angle is zero).

Table 11. Trimmed level flight conditions
Name	Value	Units
Altitude	10,013	ft MSL
True airspeed	565.6854	ft/s
C.M. position	25.0	percent
Pitch attitude	2.6538	deg +ANU
Longitudinal pilot stick position	12.96	percent +AFT
Horizontal tail position	-3.2410	deg +TED
Pilot throttle setting	13.9019	percent max
Power Lever Angle	13.9019	percent max

The F-16 Trim Conditions - Case 11 drawing shows the pitch trim solution graphically.

Graphical depiction of subsonic F-16 trim solution

(click on the image for a higher-resolution one)

A note on initial closed-loop responses

Several comments should preface the results.

In the examples that follow, the basic F-16 is under control of two Linear, Quadratic Regulator (LQR) autopilots, one for pitch and one for roll/yaw. These two LQR autopilots are at the heart of both the simple altitude, velocity, heading and course autopilot (F16_control.dml) and the navigating autopilot (F16_gnc.dml). To simplify this autopilot implementation, perfect state feedback is assumed (no state estimator model is used), and the control system has no integrated states. This results in:

large excursions in control inputs to changes in sharp commands that tend to saturate the control surfaces
Type 0 response, which implies no integral error feedback. The commanded altitude, airspeed, course, etc. are not tracked exactly; there is usually a small steady offset once steady-state behavior is achieved a few seconds after the sim is initiated. This is exacerbated due to small differences in trim conditions from truly steady-state. This is most apparent in the shorter duration step inputs where these small biases in, say, heading angle take a few seconds to damp out.

Standard autopilot - F16_control

The F16_control.dml model can be used to follow simple altitude, velocity, heading and ground track commands. The F16_control model has two inputs that serve to engage successive levels to augmentation: if stabilityAugmentationOn_disc is on (i.e. greater than 0.5), the control law will provide damping through the LQR controllers; if this is off (i.e. less than 0.5) changes in the pilot inputs (pilotControl_throttle, pilotControl_long, pilotControl_lat, and pilotControl_yaw) result is bare F-16 dynamics.

If autopilotOn_disc is on, the pilot inputs are ignored, stability augmentation is forced on, and the vehicle tracks commands from altitudeMslCommand, equivalentAirspeedCommand, trueBaseCourseCommand and lateralDeviationError inputs. The first two inputs should be self-explanatory, but the last two need some explanation.

The autopilot will point the aircraft heading to match trueBaseCourseCommand as long as lateralDeviationError is zero. This is similar behavior to a heading-hold function on a small single-engine airplane. If however a course deviation is fed back from some sensor into the lateralDeviationError input, the autopilot’s behavior is modified so that it modifies the vehicle’s heading from the initial trueBaseCourseCommand in order to reduce lateralDeviationError. This is similar to a NAV mode on an aircraft equipped with radio or satellite navigation receivers.

The four examples given below show the behavior of the F16_control model in an example 6DOF simulation model of the F-16 subsonic vehicle. These maneuvers all begin from the trimmed condition given previously.

Altitude command step change response

The Altitude change plot below shows the response of the vehicle with the F16_control autopilot in operation, for the case with autopilotOn_disc on, and a step change of 100 ft in the commanded altitude is inserted at 5 seconds after the simulation starts.

Altitude change time history

(click on the image for a higher-resolution one)

Velocity command step change response

The Velocity change plot below shows the response of the vehicle with the F16_control autopilot in operation, for the case with autopilotOn_disc on, and with a speed reduction of 10 knots equivalent airspeed inserted at 5 seconds after the simulation starts.

Velocity change time history

(click on the image for a higher-resolution one)

Heading command step change response

The Heading change strip chart and Heading change track plots below show the response of the vehicle with the F16_control autopilot in operation, for the case with autopilotOn_disc on, and with a heading change of 15 degrees to the right of the original course, at 15 seconds after the simulation starts. In this case, the lateralDeviationError signal is kept at zero throughout the maneuver, so this is just a change in heading but is not a tracking task.

Heading change time history

(click on the image for a higher-resolution one)

Heading change ground track

(click on the image for a higher-resolution one)

Lateral offset error step response

The Lateral Course Offset change strip chart and Lateral Course Offset change track plots below show the response of the vehicle with the F16_control autopilot in operation, for the case with autopilotOn_disc on, and with a course offset change of 2,000 ft to the left of the original course, at 20 seconds after the simulation starts. In this case, the lateralDeviationError signal is driven by the step (initially in the amount -2,000 ft or left of course) while trueBaseCourseCommand is kept steady at 90 deg (true East).

It can be seen that the vehicle banks right to a new intercept heading, then gradually returns to the original heading to follow the new ground track, 2,000 to the right of and parallel with the original one.

Lateral Course Offset change time history

(click on the image for a higher-resolution one)

Lateral Course Offset change ground track

(click on the image for a higher-resolution one)

Advanced autopilot - F16_gnc

The autopilot logic provided in F16_gnc.dml provides two pre-planned operational modes: one is to fly to and circle the North pole; the other is to fly to and circle the equator/International Date Line (IDL) intersection. These maneuvers were chosen to exercise the difficulty in dealing with flight near the poles, as longitude nears a singularity, and in circling the IDL when both latitude and longitude change sign.

The F16_gnc model has three inputs that serve as switches to affect this behavior. For either circling navigators to be engaged, both the stabilityAugmentationOn_disc and autopilotOn_disc inputs need to be turned on (by setting them to floating-point values of 1.0). [Off is implied with a value of 0.0 provided].

Circling the equator - International Date Line

With the F16_gnc autopilot input circumnavigatorSelectPolar set to 0.0, the autopilot will attempt to intercept and intercept a 3-nm counter-clockwise circle around the equator/IDL intersection.

Using the example simulation, the F-16 model was initialized at the trimmed flight conditions from a point located at 0 latitude, -179.95 degrees (west) longitude, heading northeast (45 deg true) A plot of the resulting ground track Equator-IDL circumnavigation ground track of a 360 second run, plotted in latitude/longitude below, indicates the vehicle makes a left turn at the 30 degree bank limit until on a heading to intercept the circle in the northwest quadrant, then intercepts and turns left again to track the circle around the intersection. It continues to maintain approximately 10,013 feet and 585.6 ft/s while doing so as shown in the strip chart image Equator-IDL circumnavigation time history.

Equator-IDL circumnavigation ground track

(click on the image for a higher-resolution one)

Equator-IDL circumnavigation time history

(click on the image for a higher-resolution one)

Circling the North Pole

Similarly, with the F16_gnc autopilot input circumnavigatorSelectPolar set to 1.0, the autopilot will attempt to intercept and intercept a 3-nm counter-clockwise circle around the North Pole.

The F-16 example simulation was initialized at the same trimmed flight conditions from a point located at 89,95 geodetic latitude, 0 degrees longitude, heading ESE (100 deg true). A plot of the resulting ground track of a 360 second run, plotted in latitude/longitude North pole circumnavigation ground track, indicates the vehicle makes a left turn at the 30 degree bank limit until on a heading to intercept the circle in the northeast quadrant, then intercepts and turns left again to track the circle around the intersection. It continues to maintain approximately 10,013 feet and 585.6 ft/s while doing so North pole circumnavigation time history.

North pole circumnavigation ground track

(click on the image for a higher-resolution one)

North pole circumnavigation time history

(click on the image for a higher-resolution one)

Change history

Table 12. Revisions to this file
2013-05-03	v1	Original version
2013-05-29	v2	Corrected units on angular rate inputs
2013-06-06	v3	Added trimmed control values as inputs per Michael Madden
2013-09-19	v4	Updated rev numbers; updated I/O lists; added trimmed values; added plots
2013-10-22	v5	Added notes in a couple places that the CM position input to inertial model should be 25% for all tests.
2013-10-22	v6	Emphasized that stab aug is forced on if any autopilot function is requested in either controller

References

[S119] Anon.: "Flight Dynamics Model Exchange Standard." American National Standard ANSI/AIAA S-119-2011, American Institute of Aeronautics and Astronautic, Washington, DC, March 2011. Available from http://arc.aiaa.org/doi/book/10.2514/4.867965 (as of May 2013)
[TM212145] Garza, Fredrico R. and Morelli, Eugene A.: "A Collection of Nonlinear Aircraft Simulations in MATLAB." NASA TM-2003-212145, National Aeronautics and Space Administration, Hampton, VA, 2003
[TP1538] Nguyen, Luat T., et al.: "Simulator Study of Stall/Post-Stall Characteristics of a Fighter Airplane with Relaxed Longitudinal Static Stability." NASA TP-1538, National Aeronautics and Space Administration, Hampton, VA, 1979
[Stevens03] Stevens, Brian L. and Lewis, Frank L.: "Aircraft Control and Simulation, second edition." John Wiley & Sons, Hoboken, NJ, 2003