*FLIGHTv2.m* provides a six-degree-of-freedom rigid-body simulation of an aircraft, as well as trimming calculations and the generation of a linearized model at any flight condition. It is a tutorial program, heavily commented to aid interpretation. Aircraft control histories, initial conditions, flag settings, and other program control actions are entered by changing variables contained in the code; there is no separate user interface. Version 2 is a major update to Version 1, including a revised business jet model based on [B-1] and a supersonic jet trainer model.

The code has been designed for simplicity and clarity, leaving a challenge to find ways to make the program run faster or more efficiently. Possible additions include simulation of "frozen" random turbulence or microburst wind shear, flight in aggressive maneuvering, and reproduction of accident scenarios. The early phase of rocket launch trajectories or the late phase of atmospheric entry can be simulated by declaring the proper initial conditions and vehicle model.

V/STOL aircraft flight simulation can be accommodated within the program framework. Aircraft mass is assumed to be constant during the simulation. Modeling mass variation with thrust, atmospheric state, and control setting requires the addition of at least one element to the state variable. With additional programming, the suite of control effectors can be expanded, and closed-loop control can be incorporated. The flat-Earth dynamic equations can be replaced by round, rotating earth equations (Section 3.3 of *Flight Dynamics*) for 6-DOF launch, orbit, and reentry simulation.

The MATLAB code is written in chapter format as follows:

**TABLE OF CONTENTS**

- CONVERSION FACTORS
- USER INPUTS
- ANALYSIS FLAGS
- MODEL and FLIGHT CONDITION SELECTION FLAGS
- NOMINAL ALTITUDE, AIRSPEED, or MACH NUMBER
- INITIAL and FINAL TIMES for Simulation
- ADDITIONAL INITIAL CONDITIONS
- TEST INPUTS at INITIAL CONDITION
- CONTROL PERTURBATION HISTORY (Test Inputs: rad or 100%)

- MODEL SELECTION, STATE VECTOR and CONTROL INITIALIZATION
- TRIM CALCULATION (for Steady Level Flight at Initial V and h)
- LTI STABILITY-AND-CONTROL-DERIVATIVE Calculation at Trimmed Flight Condition
- FLIGHT PATH SIMULATION, with Euler or Quaternion Option
- Euler Angles
- Quaternion

- Control and State History PLOTS
- NOMENCLATURE

The program provides the option of computing rotation (direction cosine) matrices using either Euler angles or quaternions. *FLIGHTv2.m* and its associated functions could be translated to similar open-source languages, such as Python, Scilab, or GNU Octave. No explicit or implicit warranties are made regarding the accuracy or correctness of the computer code. The code is subject to continuing change, and user suggestions for improving the code are welcome.

The twelfth-order mathematical models are motivated by actual aircraft; however, they should not be construed as accurate representations of any particular aircraft. Aerodynamic coefficients are either angle-of-attack-dependent or Mach-number-dependent functions. Inertial properties (total mass and moments of inertia) are constant during a simulation run. The aircraft dynamic model is specified by *ACtype* (aircraft type, see *FLIGHTv2.m*) and *MODEL* (angle-of-attack-dependent or Mach-dependent math model). The user determines which aircraft to simulate by setting these flags in the program. Aircraft characteristics are described in a function such as *BusinJetA.m* or *BusinJetM.m*. *ModelBuild.m* is an adjunct program for creating mathematical models using handbook methods, as described in Section B.4. Other aircraft types can be simulated in user-defined m-files.

*FLIGHTv2.m* is the executive file (or script) that calls program functions. Initial conditions are defined, the three primary features (trim, linearization, and simulation) are enabled, and computed output is generated. The nonlinear equations of motion take the vector form, d**x**/dt = **f**[(**x**(t), **u**(t),**w**(t)], **x**(0) given; **x**, **u**, and **w** are the state, control, and disturbance vectors. Primary features are described as follows:

- Trimming for steady, level flight (i.e., finding the equilibrium condition) is accomplished by minimizing a quadratic cost function
*J*of axial acceleration, normal acceleration, and pitch acceleration. The longitudinal trimming parameters (stabilator angle, throttle setting, and pitch angle) are found using the downhill simplex algorithm implemented in*fminsearch.m*. The optimal value of*J*is zero, and the search typically terminates when the error in*J*is less than 10^{-12}. - Small perturbations from the trimmed condition are described by the linear, time-invariant (LTI) model, dΔ
**x**/dt =**F**Δ**x**(t) +**G**Δ**u**(t) +**L**Δ**w**(t). The matrices are generated by*numjac.m*, a numerical evaluation of the Jacobian matrices associated with the nonlinear equations of motion at the equilibrium condition. The twelfth-order LTI model matrices are saved to the disk files*Fmodel.mat*and*Gmodel.mat*. - MATLAB's
*ode45.m*is the default integrator for the nonlinear equations of motion. The integrator calculates the state history,**x**(t) in [t(0), t(f)];*ode23.m*,*ode15s.m*, or any other MATLAB integration function can be selected as an alternative. The state history is displayed in time plots, with angles converted from the radians used in calculation to degrees.

The user can change the units of plotted quantities or add additional plots. Any result (e.g., numerical values of the state history) can be displayed in the MATLAB Command Window simply by removing the semi-colon at the end of the line of code. The *ACtype* flag selects either a low-angle-of-attack, Mach-dependent model or a high-angle-of-attack, low-subsonic model. The *RM* flag selects either Euler-angle or quaternion calculation of the rotation matrix.

The aircraft model must be determined by running a script such as *BusinessJetA.m*, *Business JetM.m*, or *ModelBuildJT.m* prior to running *FLIGHTv2.m* (Sections B.2.1, B.2.2, B.4.2).

The inertial, geometric, and aerodynamic characteristics of the aircraft type are specified by the parameter *ACtype*, which selects a function such as *BusinJetA.m* or *BusinJetM.m* to download the tabulated parameters. The specified parameters are saved in a mat-file that is loaded by the run-time function *AlphaModel.m* or *MachModel.m*, as selected by *MODEL*. The run-time function is called by the equations-of-motion function *EoMv2.m* or *EomQv2.m* for use in trimming and simulating aircraft motions and for calculating a linearized model at the trim condition.

Trim conditions for the angle-of-attack-dependent and Mach-dependent models of an aircraft type are approximately the same at one flight condition. For the business jet example, this condition is an altitude of 3,050 m and Mach number of 0.3. For all other flight conditions, trim values, LTI models, and simulated transient responses are different. Large-amplitude responses for Mach-independent aerodynamic coefficients are computed using *AlphaModel.m*. Small-amplitude responses for aerodynamic coefficients that vary with Mach number are computed using *MachModel.m*. If aerodynamic data for concurrent angle of attack and Mach number variations are available, an alternative function that follows the *AlphaModel.m* and *MachModel.m* structures could be scripted. (The scalar arguments for table interpolation in each function would be replaced by vector arguments of angle of attack and Mach number.)

The *BusinJetA.m* script expresses the inertial, geometric, aerodynamic, and thrust properties of the aircraft in incompressible flow (approximating *M* < 0.5). The angle-of-attack ranges from -10 to +90 deg, and the sideslip-angle range is about +/-8 degrees. Low-angle-of-attack static coefficients up to stall (angle of attack = 14 deg) are based on the wind tunnel data given in [B-1]. High-angle-of-attack coefficients (i.e., beyond 35 deg) are Newtonian estimates based on the aircraft planform. The low- and high-angle data are connected by cubic curve fits to the values and slopes of the low- and high-alpha coefficients at the matching conditions. Control surface effects are linear in their respective deflections and are tapered to zero at 90-deg angle of attack. Static lateral-directional effects are linear in sideslip angle. Dynamic coefficients are assumed to be functions of the static coefficients, as presented in Section 3.4 of the book. Thrust available is proportional to the air density ratio raised to a power that fits both static thrust at sea level and the thrust required for flight at maximum Mach number (0.81) and altitude of 13,720 m (45,000 ft). It is independent of angle of attack. *BusinJetA.m* is called prior to *FLIGHT.m* initialization. Tabulated parameters are saved in *AlphaModelData.mat*, which is loaded by *AlphaModel.m*.

The *BusinJetM.m* script uses aerodynamic and dimensional data contained in [B-1] with estimates of inertial properties of a generic business jet. Mass, inertial, geometric, aerodynamic, and thrust properties are tabulated in *BusinJetM.m*. Details of the configuration are extrapolations of the subsonic wind-tunnel data based on sweep and aspect ratio of the wing and tail. Estimates of Mach effects are based on the Prandtl factor or the modified Helmbold equation (Section 2.4 of *Flight Dynamics*). Thrust is modeled as in the alpha-dependent case. *BusinJetM.m* is called during *FLIGHT.m* initialization. Tabulated parameters are saved in the mat-file *MachModelData.mat*, which is loaded by *MachModel.m*.

The equations of motion for functions *EoM.m* and *EoMQ.m* are written using the flat-earth assumption, as described in Section 3.2 of *Flight Dynamics*. Linear and translational rates are expressed in body axes, linear position is expressed in earth-relative axes, and Euler angles (*EoM.m*) or quaternions (*EoMQ.m*) describe the orientation of the body frame with respect to the earth frame. Input and output orientations are expressed as Euler angles. If quaternion calculation of the rotation matrix is selected, transformations to and from quaternions are internal to the program. The lift and drag coefficients produced by *AlphaModel.m* or *MachModel.m* are transformed to body-axis coefficients; the remaining coefficients and thrust produced by the function are expressed in body axes. Effects of the wind are accounted for in the aerodynamic calculations. With the Euler-angle option, an ad hoc limit on the cosine of the pitch angle is imposed to prevent singular calculations in near-vertical flight. This expediency introduces a small error when the pitch angle is near +/-90 deg. The quaternion option is non-singular throughout the flight envelope and does not require the approximation.

The function *event.m* specifies a stopping condition that terminates the simulation before the final time if the altitude goes below zero.

Trim control settings are calculated by minimizing a quadratic function of longitudinal accelerations (i.e., rates of change of axial velocity, normal velocity, and pitch rate) contained in *TrimCost.m*. Equal weight is given to each component in the cost function, *J*; the defining matrix, **R**, is an identity matrix. The stabilator angle, throttle setting, and pitch angle (which equals the angle of attack in steady, level flight) are adjusted to minimize *J*.

*DCM.m* computes the direction-cosine (or rotation) matrix as a function of Euler angles. *RMQ.m* computes the rotation matrix with quaternions. The rotation matrices transform vectors from the earth-relative frame of reference to the body-axis frame.

*LinModelv2.m* calls function *EoMv2.m* to generate **f**(**x**,**u**) and appends seven dummy elements to account for the controls. Thus, **xdotj** is a function of [**x**^{T} **u**^{T}], and *numjac* calculates the Jacobian matrix evaluated at the nominal values of state and control. The stability matrix, **F**, is contained in the upper-left (12 x 12) partition, and the control-effect matrix, **G**, is contained in the upper-right (12 x 7) partition of the Jacobian.

*NatFreq* presents complex-valued eigenvalues, natural frequencies, damping ratios, and time constants that identify oscillations, convergence, and divergence in the system's modes of motion.

*WindField.m* produces a three-component wind vector as a function of altitude, with linear interpolation between tabulated points. The first point is tabulated at negative altitude to assure that no computational problems occur at zero altitude. The wind vector is rotated to body axes for application in *FLIGHT.m* and *EoM.m*.

1976 U.S. Standard Atmosphere air density, air pressure, air temperature, and sound speed are generated as functions of altitude by *Atmos.m*. The first point is tabulated at negative altitude (-1000 m) to assure that no computational problems occur at zero altitude.

The *ModelBuild.m* script uses handbook methods to estimate inertial and aerodynamic properties of a subject aircraft from available aircraft specifications. The specific constants and equations contained in the script depend on the configuration being modeled; for example, tailed and tailless aircraft have different control surfaces, thrust models are different for reciprocating and jet engines, and thrust offset effects are different for single- and multi-engine models. The script first computes *AlphaModelData.mat* and then *MachModelData.mat*, expanding on and combining the calculations of Sections B.2.1 and B.2.2. *ModelBuild.m* is run prior to executing *FLIGHTv2.m*. A prototype is presented in *ModelBuildJT.m*.

**TABLE OF CONTENTS**

- CONVERSION FACTORS and CONSTANTS
- USER INPUTS and DERIVED PARAMETERS
- DIMENSIONS of AIRCRAFT
- AIRCRAFT POWERPLANT
- COMPONENTS: Masses, Ratios, Centers of Mass, and Aerodynamic Pressure
- Masses
- Centers of Mass and Pressure
- Geometric Properties, Simulation Symbols

- AIRCRAFT: Mass and Moments/Products of Inertia
- Component Dimensions
- Moments and Products of Inertia

- INCOMPRESSIBLE AERODYNAMICS
- ANGLE OF ATTACK Table
- LONGITUDINAL AERODYNAMICS
- LIFT COEFFICIENT
- DRAG COEFFICIENT
- PITCHING MOMENT COEFFICIENT
- CONTROL EFFECT COEFFICIENTS

- LATERAL-DIRECTIONAL AERODYNAMICS
- SIDE FORCE COEFFICIENT
- ROLLING MOMENT COEFFICIENT
- YAWING MOMENT COEFFICIENT
- CONTROL EFFECT COEFFICIENTS

- SAVE AlphaModelData.mat

- COMPRESSIBLE AERODYNAMICS
- MACH NUMBER TABLE
- LONGITUDINAL AERODYNAMICS
- LIFT COEFFICIENT
- DRAG COEFFICIENT
- PITCHING MOMENT COEFFICIENT
- CONTROL EFFECT COEFFICIENTS

- LATERAL-DIRECTIONAL AERODYNAMICS
- SIDE FORCE COEFFICIENT
- ROLLING MOMENT COEFFICIENT
- YAWING MOMENT COEFFICIENT
- CONTROL EFFECT COEFFICIENTS

- SAVE MachModelData.mat

- NOMENCLATURE

The prototype *ModelBuild.m* script, presents a model for a generic jet trainer with supersonic capability. The script *ModelBuildJT.m* generates both alpha-dependent and Mach-dependent models, with output stored as *AlphaModelData.mat* and *MachModelData.mat*. The incompressible-flow model, identified as *ACType* = *JetTrainA*, is the first to be computed, and its output provides input for *AlphaModel.m*; hence, *MODEL* = *Alph*. The compressible flow model is specified by *JetTrainM* and *Mach*, with Mach-dependent coefficients extrapolated from *M* = 0 via Prandtl-Glauert-Goethert-Ackeret and Helmbold functions contained in *MachSimFcnJT.m*.

The generic jet trainer is similar to the USN *Super Hornet*, although it is smaller and has a single engine. Aerodynamic and inertial properties of the trainer are estimated from three-view drawings, which are scaled according to the wingspan ratio of the two aircraft. Scaling allows areas, chord and fuselage lengths, taper ratios, and sweep angles to be estimated. Masses and moments of inertia of components are estimated using conventional structural analysis of simple geometric shapes and the parallel-axis theorem, as in [R-1]. For this conceptual design, static and dynamic aerodynamic coefficients are estimated with handbook methods, following Sections 2.4 and 3.4 of *Flight Dynamics, Second Edition* and [F-1]. Thrust for the trainer model is based on estimated characteristics of a current powerplant.

The aerodynamic coefficients are tabulated separately as functions of angle of attack (-10 deg to 90 deg) and Mach number (0 to 1.5). The alpha-dependent data are stored in *AlphaModelData.mat*, which serves as input to the simulation function *AlphaModel.m*. The Mach-dependent data are stored in *MachModelData.mat*, which serves as input to the simulation function *MachModel.m*.

This function provides a smooth approximation in a region between scalar aerodynamic coefficients that are tabulated for low and high values of the independent variable. Angle-of-attack-dependent coefficients may be tabulated for partially attached or fully separated upper-surface flows at low and high angles. For the alpha-dependent model, the cubic fit provides estimates in the intermediate range from stall to post-stall that match the values and slopes of the coefficients at the interior points. For the Mach-dependent case, the cubic fit provides estimates in the transonic regime, matching the values and slopes of the coefficients at the interior points.

F-1 Fink, R., and Hoak, D., "USAF Stability and Control DATCOM," USAF AFFDL, Wright Patterson AFB, Dayton, 1978.

R-1 Raymer, D. P., "Aircraft Design: A Conceptual Approach," American Institute of Aeronautics and Astronautics, Reston, 1989 (sixth edition, 2018).

**Flight Dynamics, Second Edition**.

Last updated on June 21, 2023.

© Copyright 2023 by Robert F. Stengel. All rights reserved.

MATLAB is a registered trademark of The MathWorks, Inc.