F-16 High-Fidelity 6DOF Flight Dynamics Simulation

Framework
Achraf khelifi touhami
November 11, 2025

Contents

1 List of Abbreviations 3

2 Abstract 4

3 Introduction 4
3.1 Equations of Motion Block . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.2 Translational Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.3 Rotational Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.4 Euler and Quaternion Kinematics . . . . . . . . . . . . . . . . . . . . . . . . 8
3.5 NED Position Rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.6 Body-Fixed Velocities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

4 Actuator Models 8
4.1 Control Surface Actuators . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
4.2 Thrust Vector Control System . . . . . . . . . . . . . . . . . . . . . . . . . . 10
4.3 Control Input Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

5 Atmospheric Model (ISA) 11

5.1 Reference Constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
5.2 Temperature and Pressure Calculation . . . . . . . . . . . . . . . . . . . . . 12
5.3 Derived Quantities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

6 Atmospheric Turbulence Modeling 13

6.1 Dryden Turbulence Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
6.2 Discrete-Time Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . 13
6.3 Turbulence-Induced Aerodynamic Perturbations . . . . . . . . . . . . . . . . 14
6.3.1 Linear Velocity Perturbations . . . . . . . . . . . . . . . . . . . . . . 14
6.3.2 Angular Rate Perturbations . . . . . . . . . . . . . . . . . . . . . . . 14

1
7 Aerodynamic Force and Moment 15
7.1 Aerodynamic Force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
7.2 Aerodynamic Moment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

8 Thrust Vectoring 16
8.1 Force Components . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
8.2 Induced Moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

9 Sensors Model 17

10 Sensor Model Architecture 17

10.1 Bias dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
10.2 Measurement noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
10.3 Inertial Measurement Unit (IMU) . . . . . . . . . . . . . . . . . . . . . . . . 18
10.3.1 Accelerometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
10.3.2 Gyroscope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

11 Simulation Setup 20

12 Results 21

13 Conclusion 23

14 References 24

2
1 List of Abbreviations
• V : Airspeed (m/s)
• α: Angle of attack (rad)
• β : Sideslip angle (rad)
• γ : Flight-path angle (rad)
• µ : Bank angle (rad)
• p, q, r : Body angular rates (rad/s)
• ϕ, θ, ψ : Euler angles (rad)
• χ : Heading angle (rad)
• u, v, w : Body-axis velocities (m/s)
• xE , yE , zE : Earth-fixed NED position (m)
• q0 , q1 , q2 , q3 : Quaternion components
• D, L, Y : Drag, lift, side forces (N)
• La , Ma , Na : Aerodynamic moments (N·m)
• Tx , Ty , Tz : Thrust force components (N).
• LT , MT , NT : Thrust-induced moments (N·m).
• uturb , vturb , wturb : Dryden turbulence velocity components (m/s).
• pturb , qturb , rturb : Angular turbulence velocity components (rad/s).
• σu , σv , σw : Turbulence intensity parameters (m/s).
• Lu , Lv , Lw : Turbulence length scales (m).
• m : Mass (kg)
• g : Gravity (9.81 m/s2 )
• Ixx , Iyy , Izz , Ixz : Inertia tensor elements (kg·m2 )
• ∆ : Ixx Izz − Ixz
2

• IMU: Inertial Measurement Unit

• bias(t) : bias process
• noise(t) : measurement noise
• H(s) : Sensor Transfer Function

3
2 Abstract
This project presents a comprehensive high-fidelity open-loop simulation framework for the
F-16 Fighter jet, integrating six degrees of freedom (6DOF) flight dynamics. The simulation
captures complex nonlinear flight behaviour through a 22-state dynamic model while incor-
porating realistic sensor measurements from inertial measurement units (IMU),The sensor
models provide realistic measurement data incorporating systematic errors and stochastic
noise. The framework features translational and rotational dynamics, quaternion-based at-
titude representation, aerodynamic force and moment modelling, actuator dynamics with
physical constraints, thrust vector control (TVC) effects, and atmospheric turbulence mod-
elling using the Dryden turbulence model. The simulation environment, developed in MAT-
LAB/Simulink, employs fourth-order Runge-Kutta (RK4) numerical integration. The frame-
work enables a comprehensive evaluation of flight control systems under realistic operational
conditions.

3 Introduction
The F-16 stands as one of the most versatile and widely deployed fighter aircraft in modern
aviation history. Its exceptional manoeuvrability, coupled with sophisticated flight control
systems, has established it as a benchmark for high-performance military aircraft design.
Understanding and accurately modelling the complex flight dynamics of such aircraft, along
with their associated systems, is paramount for advancing aerospace engineering, developing
next-generation control systems, and ensuring mission success in demanding operational
environments.
Modern fighter aircraft rely heavily on integrated sensor systems to provide accurate
state estimation for flight control, navigation, and mission systems. These sensors, including
inertial measurement units (IMU), GNSS reicievers, barometric altimeter and others must
operate reliably across diverse flight conditions while maintaining precision in the presence
of various error sources.
Flight dynamics simulation of high-performance aircraft presents unique challenges due
to the inherently nonlinear, coupled, and multivariable nature of aerodynamic phenomena,
particularly at extreme flight conditions. The F-16’s capability to operate across a wide
flight envelope, including high angles of attack and aggressive manoeuvring scenarios, ne-
cessitates sophisticated mathematical models that capture the intricate interplay between
aerodynamic forces, inertial effects, control system dynamics, sensor measurements and ac-
tuator dynamics.
This report presents a comprehensive open-loop simulation framework for the F-16 air-
craft, employing a high-fidelity 22-state nonlinear mathematical model that encompasses
complete six-degree-of-freedom dynamics. The simulation architecture integrates transla-
tional and rotational motion equations, quaternion-based attitude representation to avoid
singularities, detailed actuator dynamics with realistic physical constraints, advanced thrust
vector control capabilities, comprehensive atmospheric turbulence modelling, and sophisti-
cated sensor models incorporating both systematic bias and stochastic noise components.

4
 The mathematical foundation is based on the first principles of Newtonian mechan-
ics, incorporating the Newton-Euler equations of motion along with appropriate coordinate
transformations between principal reference frames. The sensor models follow established
s-domain transfer function representations with a dual-noise architecture that captures both
low-frequency bias drift and high-frequency measurement noise characteristic of real sensor
systems.
The simulation environment, developed within MATLAB’s Simulink framework, provides
a robust computational platform for numerical integration using different Numerical Solvers
to solve a set of ordinary differential equations (ODEs), including the Euler Forward and
fourth-order Runge-Kutta schemes. In this simulation, the 4th-order Runge-Kutta solver
has been chosen as it is more accurate than the Euler forward solver, but at a higher com-
putational cost.
The model architecture is designed with modularity in mind, facilitating future exten-
sions for closed-loop flight control system development and validation. Instead of using
built-in Simulink blocks to model the 6-DOF nonlinear dynamics, a user-defined function
(MATLAB Function Block) was used to provide design freedom for implementing various
design concepts. Figure 1 illustrates the input-output structure of the developed simulation
framework, demonstrating the interconnected nature of the various subsystems.

Figure 1: Fighter Aircraft Simulation Model Architecture

3.1 Equations of Motion Block

The equations of motion governing aircraft dynamics are formulated through the application
of Newton-Euler principles, which synthesize Newton’s second law of motion with coordinate
transformation matrices operating between principal reference frames. The mathematical
description of aircraft motion necessitates the establishment of three fundamental coordinate
systems, each serving a distinct analytical purpose as illustrated in Figure 2:

• Inertial Frame (North-East-Down, NED): The inertial reference frame consti-

tutes an Earth-fixed coordinate system wherein the positive Xe axis is oriented north-
ward, the positive Ye axis eastward, and the positive Ze axis directs downward in
alignment with the local gravity vector. This frame provides the absolute spatial ref-
erence for aircraft position and trajectory determination.

5
 • Body-Fixed Frame: The body axis system is rigidly affixed to the aircraft’s center
of gravity (CG) and adheres to conventional aerospace definitions. The positive Xb
axis extends forward through the aircraft nose, the positive Yb axis projects laterally
toward the starboard wingtip, and the positive Zb axis points downward perpendicular
to the aircraft’s longitudinal-lateral plane.

• Wind Frame: The wind axis system is aligned with the aircraft’s velocity vector rel-
ative to the surrounding air mass. The positive Xw axis is collinear with the velocity
vector (pointing in the direction of flight), the positive Zw axis lies in the plane of sym-
metry perpendicular to Xw (pointing downward), and the positive Yw axis completes
the right-handed orthogonal system (pointing to the right). This frame is particularly
advantageous for aerodynamic force representation, as the drag force acts along the
negative Xw axis and the lift force along the negative Zw axis.

These coordinate systems, along with their associated transformation matrices, form
the mathematical foundation for expressing the complete six-degree-of-freedom motion of
the aircraft, encompassing both translational and rotational dynamics in three-dimensional
space.

Figure 2: Definition of Axes and Aerodynamical Angles.

6
3.2 Translational Dynamics

Y sin β − mg sin γ − D Tx cos β cos α + Ty sin β + Tz sin β cos α

V̇ = + (1)
m m
mg cos γ cos µ − L −Tx sin α + Tz cos α
α̇ = q − (r sin α + p cos α) tan β + + (2)
mV cos β + 0.1 mV cos β
−Tx sin β cos α + Ty cos β − Tz sin β sin α Y cos β + mg cos γ sin µ
β̇ = −r cos α + p sin α + +
mV mV
(3)
−Y cos β sin µ + L cos µ Ty Tx

γ̇ = − cos β + cos α sin β sin µ + sin α cos µ
mV mV mV
Tz

+ sin α sin β sin µ − cos α cos µ (4)
mV
p cos α + r sin α L
Y + Ty
µ̇ = + tan β + sin µ tan γ + cos β cos µ tan γ
cos β mV mV
mg Tx sin α − Tz cos α

− tan β cos µ cos γ + tan β + sin µ tan γ
V mV
Tx cos α + Tz sin α
− sin β cos µ tan γ (5)
mV
Y cos β cos µ + L sin µ Tx

χ̇ = + sin α sin µ − cos α sin β cos µ
mV cos γ mV cos γ
Ty Tz

+ cos µ cos β − cos µ sin α sin β + cos α sin µ (6)
mV cos γ mV cos γ

3.3 Rotational Dynamics

2
(Ixx − Iyy )Ixz + Ixz Izz pq (Iyy Izz − Izz 2
Izz (La + LT ) + Ixz (Na + NT ) − Ixz )qr
ṗ = + +
∆ ∆ ∆
(7)
Ma + MT (Ixx − Izz )pr Ixz (p2 − r2 )
q̇ = − − (8)
Iyy Iyy Iyy
2 2
Ixz (La + LT ) + Ixx (Na + NT ) (Ixx − Ixx Iyy + Ixz )pq Ixz (Ixx − Iyy + Izz )qr
ṙ = + − (9)
∆ ∆ ∆
2
where ∆ = Ixx · Izz − Ixz and Ixx ,Ixx ,Ixx are the moments of inertia.

7
3.4 Euler and Quaternion Kinematics

ϕ̇ = p + q sin ϕ tan θ + r cos ϕ tan θ (10)

θ̇ = q cos ϕ − r sin ϕ (11)
q sin ϕ + r cos ϕ
ψ̇ = (12)
cos θ
q̇ = 21 q ⊗ [0, p, q, r]⊤ (13)

3.5 NED Position Rates

ẋE = V cos χ cos γ, ẏE = V sin χ cos γ, żE = −V sin γ (14)

3.6 Body-Fixed Velocities

Tx − D cos α cos β + L sin α

u̇ = rv − qw + g sin θ + (15)
m
Ty + Y
v̇ = −ru + pw − g cos θ sin ϕ + , (16)
m
Tz − D sin α cos β − L cos α
ẇ = qu − pv − g cos θ cos ϕ + . (17)
m
where V is flight speed, m is the mass of the aircraft fighter, and g is the gravity accel-
eration. α and β represent the angle of attack and sideslip angle, respectively. γ, χ and µ
denote flight path angle, velocity heading angle and roll angle around the velocity vector,
respectively. p, q and r represent roll, pitch and yaw angular rate in the body coordinate
system, respectively. xE , yE and zE are the flight position coordinates.

4 Actuator Models
The actuators are critical components of the flight control system, providing the motive
power needed to move the flight control surfaces and manage thrust vectoring. The F-16
flight control system encompasses both conventional aerodynamic control surfaces (aileron,
elevator, rudder) and advanced thrust vector control (TVC) capabilities for enhanced ma-
noeuvrability.

4.1 Control Surface Actuators

The primary flight control surfaces, namely aileron (δa ), elevator (δe ), and rudder (δr ), are
shown in Figure 3. These conventional control surfaces provide the primary means of aircraft
attitude and trajectory control through aerodynamic moment generation.

8
 Figure 3: F-16 aerodynamic control surfaces

The actuator mathematical models for the control surfaces are represented as second-
order transfer functions with pure time delay

δi (s) −sτ ωn2

=Ke (18)
δi,com (s) s2 + 2ζωn s + ωn2

where:

• K: steady-state gain,

• ωn : natural frequency [rad/s],

• ζ: damping ratio.

• τ : time delay .

• i ∈ {a, e, r} (aileron, elevator, rudder)

For implementation, the delay term e−sτ can be approximated using 1-order Padé expansions
as:
1 − τ2 s
e−sτ ≈ (19)
1 + τ2 s
for K = 1 , ωn = 75 rad/s, ζ = 0.7 and τ = 0.008 s. We will have

752 1 − 0.008


δi (s) 2
ss
=

(20)
δi,com (s) 1 + 0.004 s s2 + 105 s + 752

9
4.2 Thrust Vector Control System
In addition to conventional control surfaces, the F-16 model incorporates thrust vector con-
trol capabilities through one movable engine nozzle. The TVC system provides two degrees
of freedom: pitch vectoring (δz ) and yaw vectoring (δy ), enabling direct thrust-based control
moments that are particularly effective at high angles of attack where aerodynamic control
authority may be reduced.
The magnitude of engine thrust TE and the deflection angles of the nozzle are controlled
through dedicated actuator systems with similar dynamic characteristics but with a delay
in time.
The enhanced actuator model framework with first-order Padé approximation:

752 1 − 0.005


δy (s) 2
s
= 1.0 · 0.005

(21)
δy,com (s) 1+ 2
s (s2 + 2(0.7)(75)s + 752 )
752 1 − 0.005


δz (s) 2
s
= 1.0 · (22)
1 + 0.005


δz,com (s) 2
s (s2 + 2(0.7)(75)s + 752 )
502 1 − 0.012


TE (s) 2
s
= 1.0 · (23)
1 + 0.012


TE,com (s) 2
s (s2 + 2(0.8)(50)s + 502 )
The TVC actuator parameters are:
• Nozzle deflection actuators (δy , δz ): ωn = 75 rad/s, ζ = 0.7, K = 1, τ = 8 ms

• Thrust magnitude actuator (TE ): ωn = 50 rad/s, ζ = 0.8,K = 1, τ = 12 ms

4.3 Control Input Constraints

All actuators operate within physical limitations that must be enforced to ensure realistic
simulation behaviour. The control input vector is defined as:

U = [δa , δe , δr , δy , δz , TE ]T (24)

The operational limits for each control input are:

−25 ≤ δe ≤ 25 (elevator deflection) (25)

−25 ≤ δa ≤ 25 (aileron deflection) (26)
−25 ≤ δr ≤ 25 (rudder deflection) (27)
−15 ≤ δy ≤ 15 (TVC pitch vectoring) (28)
−15 ≤ δz ≤ 15 (TVC yaw vectoring) (29)
1 N ≤ TE ≤ 130000 N (engine thrust) (30)

These saturation limits are enforced in the simulation to prevent unrealistic control com-
mands and ensure the actuator responses remain within their physical capabilities. The
combination of conventional control surfaces and thrust vectoring provides the F-16 with
exceptional control authority across its entire flight envelope.

10
 Figure 4 illustrates the Simulink implementation of the actuation system.

Figure 4: Simulink block implementation of the Actuation System

The complete actuator system parameters are summarised in Table 1.

Table 1: Actuator Model Parameters

Control Position Rate ωn ζ τ

Surface Limit Limit (rad/s) (ms)
δa,e,r ±25 ±80/s 75 0.6 8
δy,z ±15 ±60/s 75 0.7 5
TE 1 − 130 kN ±50 kN/s 50 0.8 12

This implementation of the actuator models provide a realistic representation of modern

flight control actuator systems, including the effects of digital processing delays that are char-
acteristic of contemporary fly-by-wire flight control architectures. The time delay parameters
can be adjusted based on specific hardware configurations and system requirements.

5 Atmospheric Model (ISA)

As an aircraft operates at various altitudes, the aerodynamic coefficients need to account
for changes in density as represented by a standard atmosphere model. The International
Standard Atmosphere (ISA) model used in this simulation provides standard atmospheric
values for air temperature T , pressure P , density ρ, speed of sound a, and dynamic pressure
qbar , as a function of the altitude h and true airspeed VT .

11
5.1 Reference Constants
The ISA model uses the following sea-level reference conditions:
T0 = 288.15 K (sea-level temperature) (31)
P0 = 101325 Pa (sea-level pressure) (32)
R = 287.05 J/(kg·K) (specific gas constant for air) (33)
g = 9.80665 m/s2 (standard gravity) (34)
γ = 1.4 (heat capacity ratio) (35)

5.2 Temperature and Pressure Calculation

The atmospheric properties are computed using the hydrostatic equation and the ideal gas
law based on altitude layers.
Troposphere (h ≤ h1 = 11000 m):
T = T0 + L1 h (36)
 − RLg
T 1
P = P0 (37)
T0
where L1 = −0.0065 K/m is the temperature lapse rate in the troposphere.
Above Troposphere (h > 11000 m):
For altitudes above 11 km, an isothermal approximation is used:
T1 = T0 + L1 h1 (38)
 − RLg
T1 1
P1 = P0 (39)
T0
T = T1 (40)
 
g(h − h1 )
P = P1 exp − (41)
RT1

5.3 Derived Quantities

Once temperature and pressure are determined, the remaining atmospheric properties are
calculated using fundamental thermodynamic relationships:
Air density:
P
ρ= (42)
RT
Dynamic pressure:
1
qbar = ρVT2 (43)
2
Speed of sound: p
a = γRT (44)
Mach number:
VT
M= (45)
a
12
6 Atmospheric Turbulence Modeling
Real-world flight operations are significantly influenced by atmospheric turbulence, which
originates from atmospheric pressure instabilities and non-uniform temperature distribu-
tions, introducing stochastic disturbances that affect aircraft stability, handling qualities,
and structural loads. The incorporation of accurate turbulence models is essential for com-
prehensive flight simulation and robust control system validation under realistic atmospheric
conditions.
The inherently stochastic nature of atmospheric turbulence cannot be modelled exactly,
necessitating statistical approaches that capture the essential spectral characteristics of tur-
bulent velocity fields. Several well-established turbulence models have been used in the liter-
ature to implement turbulence, among which the Dryden and von Kármán spectral models
have achieved widespread acceptance due to their physical realism, mathematical simplicity,
and computational efficiency.

6.1 Dryden Turbulence Model

The Dryden turbulence model is widely adopted in aerospace engineering due to its math-
ematical tractability and good correlation with measured atmospheric turbulence spectra.
The model generates correlated turbulence velocity components through shaping filters that
transform white noise inputs into realistic atmospheric disturbances with appropriate spec-
tral characteristics.
The Dryden model defines three turbulence velocity components corresponding to longi-
tudinal (uturb ), lateral (vturb ), and vertical (wturb ) atmospheric disturbances. These compo-
nents are generated using the following continuous-domain transfer functions:

r
2Lu 1
Gug (s) = σu · (46)
πV 1 + LVu s
r √
2Lv 1 + 2 V3Lv s
Gvg (s) = σv ·
2 (47)
πV 1 + 2L v
s
V
r √
2Lw 1 + 2 V3Lw s
Gwg (s) = σw ·
2 (48)
πV 1 + 2Lw sV

where σu , σv , σw are the turbulence intensity parameters, Lu , Lv , Lw are the turbulence

length scales, and V is the aircraft airspeed. The longitudinal turbulence exhibits first-order
dynamics, while the lateral and vertical components follow second-order characteristics to
capture the more complex cross-coupling effects observed in atmospheric measurements.

6.2 Discrete-Time Implementation

The continuous-domain Dryden filters are converted to discrete-time equivalents using the
bilinear transformation (Tustin method) to enable numerical integration within the simu-

13
lation framework. This transformation preserves the essential spectral characteristics while
ensuring numerical stability at the chosen integration timestep.
The bilinear transformation s = dt2 z+1
z−1
is applied to each transfer function, resulting in
discrete-time difference equations suitable for real-time implementation. The transformed
filters maintain appropriate filter memory states to ensure continuity and realistic temporal
correlation in the generated turbulence sequences.

6.3 Turbulence-Induced Aerodynamic Perturbations

Atmospheric turbulence modifies the relative wind velocity vector experienced by the aircraft.
The turbulence components in the body-fixed reference frame are denoted as uturb , vturb ,
and wturb , representing longitudinal, lateral, and vertical wind disturbances, respectively.
These components are generated through spectral filtering of white noise to produce realistic
atmospheric disturbances with appropriate temporal and spatial correlation characteristics.

6.3.1 Linear Velocity Perturbations

The relative wind velocity experienced by the aircraft is computed by subtracting the tur-
bulence velocity components from the aircraft’s body-frame velocity components:

urel = u − uturb (49)

vrel = v − vturb (50)
wrel = w − wturb (51)

The magnitude of the relative velocity and the aerodynamic angles are computed from
the perturbed velocity components:
q
Vrel = u2rel + vrel
2 2
+ wrel (52)
 
wrel
α = arctan (53)
urel
 
vrel
β = arcsin (54)
Vrel

where α is the angle of attack and β is the sideslip angle, both modified by the presence
of atmospheric turbulence.

6.3.2 Angular Rate Perturbations

In addition to linear velocity turbulence, atmospheric disturbances induce rotational per-
turbations on the aircraft. These angular rate disturbances arise from spatial gradients in
the turbulent velocity field and can be approximated using the following relations based on

14
MIL-STD-1797A:
π
pgust = (wturb − vturb ) (55)
4b
π
qgust = wturb (56)
3b
π
rgust = (vturb + wturb ) (57)
4b
where b is the aircraft wingspan and pgust , qgust , rgust represent roll, pitch, and yaw rate
disturbances, respectively. These angular perturbations modify the effective body angular
rates experienced by the aerodynamic surfaces, thereby affecting the damping and rate-
dependent moment contributions.
The relative angular rates become:
prel = p − pgust (58)
qrel = q − qgust (59)
rrel = r − rgust (60)

7 Aerodynamic Force and Moment

7.1 Aerodynamic Force
The aerodynamic forces are computed using the disturbed values Vrel , α, and β in place of
the undisturbed values.
The aerodynamic forces are initially computed in the wind-axis, where:
L = qbar S CL (Lift, perpendicular to Vrel ) (61)
D = qbar S CD (Drag, parallel to Vrel ) (62)
Y = qbar S CY (Side force) (63)
where S is the reference wing area, and qbar is the dynamic pressure calculated using the
relative velocity:
1 2
qbar = ρVrel (64)
2
CL ,CD , CY are the aerodynamic force coefficients which are evaluated using perturbed
aerodynamic angles.
For a typical high-performance fighter aircraft, such as the F-16, can be represented as:
qrel c̄
CL = CL0 + CLα α + CLq + CLδe δe (65)
2Vrel

CD = CD0 + CDα α2 + CDδe |δe | (66)

prel b rrel b
CY = CYβ β + CYp + C Yr + CYδa δa + CYδr δr (67)
2Vrel 2Vrel

15
7.2 Aerodynamic Moment
The aerodynamic moments about the aircraft center of gravity are similarly computed using
the relative velocity and angular rates.
The aerodynamic moments about the body-fixed axes are:
Laero = qbar S b Cℓ (68)
Maero = qbar S c̄ Cm (69)
Naero = qbar S b Cn (70)
where b is the wingspan and c̄ is the mean aerodynamic chord.
Cl ,Cm , Cn are the moment coefficients represented as:
prel b rrel b
Cℓ = Cℓβ β + Cℓp + Cℓr + Cℓδa δa + Cℓδr δr (71)
2Vrel 2Vrel

qrel c̄
Cm = Cm0 + Cmα α + Cmq + Cmδe δe (72)
2Vrel

prel b rrel b
Cn = Cnβ β + Cnp + Cnr + Cnδa δa + Cnδr δr (73)
2Vrel 2Vrel

Table 2: Aerodynamic Stability Derivatives

Derivative Value Derivative Value

CLα 5.8 rad−1 Cℓβ −0.12 rad−1
CDα 0.15 – Cℓp −0.40 –
C Yβ −0.75 rad−1 Cℓr 0.15 –
CLq 8.0 – Cmα −1.2 rad−1
C Yp 0.05 – Cmq −12.0 –
C Yr 0.30 – Cnβ 0.10 rad−1
CLδe 0.40 rad−1 Cnp −0.05 –
C Yδ r 0.15 rad−1 Cnr −0.20 –
Cℓδa 0.15 rad−1 Cmδe −1.5 rad−1
Cℓδr 0.01 rad−1 Cnδa 0.01 rad−1
– – Cnδr −0.12 rad−1

8 Thrust Vectoring
Modern fighter aircraft like the F-16 employ thrust vectoring systems to enhance manoeu-
vrability and control authority. The thrust vectoring system redirects the engine exhaust
flow through movable nozzles, creating additional control forces and moments which pro-
vide additional control authority, particularly beneficial during high-angle-of-attack flight or
when aerodynamic control surfaces become less effective.

16
8.1 Force Components
The thrust vector components in the body frame are calculated based on the nozzle deflection
angles. The thrust force decomposed into body-axis components according to the nozzle
orientation:

Tx = ζf TE cos(δz ) cos(δy ) (74)

Ty = ζf TE cos(δz ) sin(δy ) (75)
Tz = ζf TE sin(δz ) (76)

where δy and δz are the nozzle deflection angles in the yaw and pitch planes, respectively,TE
is the engine thrust magnitude, and ζf is the combined nozzle loss coefficient accounting for
efficiency losses in the vectoring system.

8.2 Induced Moments

The thrust forces applied at the engine location create additional moments about the air-
craft’s centre of gravity. These thrust-induced moments are calculated using the cross prod-
uct of the moment arm and thrust force vector:

LT = yE Tz − zE Ty (77)
MT = zE Tx − xE Tz (78)
NT = xE Ty − yE Tx (79)

where (xE , yE , zE ) represents the engine location relative to the centre of gravity in body
coordinates.

9 Sensors Model
Modern fighter aircraft are critically dependent on integrated sensor systems to provide
accurate state estimation for flight control, navigation, and mission systems. The F-16
simulation framework incorporates a comprehensive multisensor model that includes inertial
measurement units (IMU) with realistic error characteristics and dynamic responses.

10 Sensor Model Architecture

The multi-sensor suite is represented by a unified, five-stage processing chain that captures
the principal characteristics of practical sensing systems:
1. Dynamic filtering: sensor dynamics described by a transfer function.

2. Scale transformation: application of sensor-specific scale (gain) factors.

3. Bias injection: addition of slowly varying systematic offsets.

4. Noise corruption: addition of measurement noise (high-frequency).

17
 5. Time delay: latency due to digital processing and communication.

A general time-domain model for a single measurement channel is

y(t) = K H{u(t − T )} + b(t − T ) + n(t − T ) (80)

where H represents the sensor transfer function, K is the scale factor, u(t) the true input,
T the pure time delay, b(t) the bias process, and n(t) the additive measurement noise.
In the Laplace domain:

Y (s) = K H(s) e−sT U (s) + e−sT B(s) + N (s)

(81)

i.e. the delay operator e−sT applies to the complete processed signal.

10.1 Bias dynamics

Sensor bias denotes a slowly varying systematic error caused by temperature, ageing, calibra-
tion drift and similar effects. A common and practical model is the first-order Gauss-Markov
process
1
ḃ(t) = − b(t) + w(t) (82)
τ
where τ > 0 is the bias correlation time and w(t) is zero-mean white driving noise with
variance σw2 . The corresponding transfer function from the driving noise to the bias is

B(s) 1
= . (83)
W (s) τs + 1

where the bias variance is

σw2 τ 2σb2
σb2 = =⇒ σw2 = (84)
2 τ

10.2 Measurement noise

High-frequency measurement disturbances (e.g. ADC quantisation, electronic noise, small-
scale turbulence) are modelled as additive, zero-mean Gaussian white noise:

n(t) ∼ N (0, σn2 ), (85)

where σn2 denotes the noise variance.

10.3 Inertial Measurement Unit (IMU)

The IMU provides three-axis acceleration and angular rate measurements through accelerom-
eter and gyroscope sensor clusters. These sensors are fundamental to the aircraft’s inertial
navigation system and flight control loops. The sensor families follow the general architecture
above.

18
10.3.1 Accelerometer
The accelerometer measurement vector is modelled as

yacc (t) = Kacc Hacc {ftrue (t − T )} + bacc (t − T ) + nacc (t − T ), (86)

with a second-order dynamics approximation

ωn2
Hacc (s) = , (87)
s2 + 2ζωn s + ωn2

and bias dynamics governed by the Gauss-Markov model in which each axis follows the same
form as (3). Typical parameter values (representative) are:

• natural frequency: ωn = 100 rad/s (≈15.9 Hz)

• damping ratio: ζ = 0.7 (underdamped)

• scale factor: Kacc = 1.0 ± 0.005 (0.5% error)

• bias correlation time: τacc = 300 s

• bias std: σb,acc = 0.01 m/s2

√
• noise variance: σn,acc
2
= 0.001 m/s2 / Hz.

10.3.2 Gyroscope
The gyroscope measurement is modelled as

ω meas (t) = Kgyro Hgyro {ω true (t − T )} + bgyro (t − T ) + ngyro (t − T ) (88)

with a first-order dynamics approximation

1
Hgyro (s) = (89)
τgyro,d s + 1

and Gauss–Markov bias dynamics per axis. Representative parameters are:

• dynamics time constant: τgyro,d = 0.01 s (≈100 Hz bandwidth)

• scale factor: Kgyro = 1.0 ± 0.002 (0.2% error)

• bias correlation time: τgyro,b = 300 s

• bias std: σb,gyro = 0.001 rad/s

√
• noise variance: σn,gyro
2
= 1 × 10−4 rad/s/ Hz.

19
 Figure 5 illustrates the Simulink implementation of this architecture for the gyroscope,
demonstrating the cascaded processing stages from true altitude input to corrupted mea-
surement output.

Figure 5: Simulink block implementation of the gyroscope sensor model

11 Simulation Setup
The simulation has been implemented in MATLAB R2011a with the following initial condi-
tions and parametres:

• Integrator: RK4, ∆t = 0.001 s, T = 50 s.

• Initial: V = 200 m/s, α = 0.05 rad, β = 0, p = q = r = 0.

• u = 200, v = w = 0.

• Euler (0, 0.05, 0), xE , yE , zE = (0, 0, 10), quaternion from Euler.

• Aircraft: m = 9000 kg; Ixx = 12875, Iyy = 15673, Izz = 20552, Ixz = 1031, g = 9.80.

• the nozzle loss coefficient ζf = 1.

• S = 30 Wing area (m²).

• b = 9.5 Wing span (m).

• c = 3.5 Mean aerodynamic chord (m).

• Engine offset: (−0.5, 0, 0) m.

• Turbulence: Dryden model with light intensity, filter initialisation at zero states.

20
12 Results
The following results showcase both the flight dynamic response.

Figure 6: Angle of attack α.

Figure 7: Yaw rate r.

21
Figure 9: Body linear vertical velocity w.

Figure 8: Pitch rate q.

22
 Figure 10: Altitude h.

13 Conclusion
This work presents a comprehensive high-fidelity simulation framework for the F-16 fighter
aircraft, successfully integrating six-degree-of-freedom nonlinear flight dynamics with realis-
tic sensor models and atmospheric disturbances. The implementation demonstrates several
key achievements that establish a robust foundation for advanced flight control research.
The numerical architecture, employing fourth-order Runge-Kutta integration with 1 ms
timestep, maintained stable behaviour throughout extended simulation periods without nu-
merical instabilities. The adoption of quaternion-based attitude representation effectively
eliminated the gimbal lock singularities inherent to Euler angle formulations, ensuring stable
flight across the entire attitude envelope, including extreme pitch conditions near ±90.
The actuator subsystem accurately captures the physical constraints and dynamic char-
acteristics of modern flight control systems. Position and rate saturation limits, combined
with second-order transfer function dynamics and realistic time delays, provide authentic
representations of control surface and thrust vectoring behaviour during aggressive manoeu-
vring scenarios. The observed oscillatory responses in the angular rates (p, q, r), angle of
attack (α), and side angle (β) reflect the expected behavior of the open-loop system in the
absence of stabilising feedback control.
Integration of atmospheric turbulence through the discrete-time Dryden model success-
fully introduces realistic disturbances while maintaining numerical stability via bilinear trans-
formation. Turbulence forces and moments provide essential environmental excitation for a
complete validation of the control system under operational conditions.
The multi-sensor architecture implements a sophisticated dual-noise error model captur-
ing both systematic bias drift through Gauss-Markov processes and high-frequency mea-
surement noise. The integration of the IMU model with realistic transfer function dynamics,
scale factor errors, and time delays provides authentic measurement data.
This framework provides a robust testbed for advanced research in flight control design,

23
sensor fusion, enabling rigorous evaluation of control architectures under realistic operational
conditions such as high-fidelity dynamics, realistic actuator constraints and environmental
disturbances, etc..

14 References
1. Liu, J., Sun, M., Chen, Z., & Sun, Q. High AOA decoupling control for aircraft based
on ADRC.

2. Öznurlu, C. Data-Driven Model Discovery and Control of Lateral-Directional Fighter

Aircraft Dynamics.

3. Liu, J., Sun, M., Chen, Z., & Sun, Q. High AOA decoupling control for aircraft based
on ADRC.

4. Liu, J., Sun, M., Chen, Z., & Sun, Q. Output feedback control for aircraft at high
angle of attack based upon fixed-time extended state observer.

5. Ijaz, S., Fuyang, C., Hamayun, M. T., & Anwaar, H. Adaptive integral-sliding-mode
control strategy for manoeuvring control of the F16 aircraft subject to aerodynamic
uncertainty.

6. MUHARREM ULU. SLIDING MODE GUIDANCE OF AN AIR-TO-AIR MISSILE.

7. Groves, P. D. Principles of GNSS, Inertial, and Multisensor Integrated Navigation

Systems. Artech House, 2013.

8. Teuku Mohd Ichwanul Hakim, Ony Arifianto. Implementation of Dryden Continuous

Turbulence Model into Simulink for LSA-02 Flight Test Simulation, 2018.

9. Eduardo Pedó Gutkoski. Noise modelling in low-cost MEMS IMU sensors using ex-
perimental data, 2022.

24