You guys have the supplementary notes in PDF form for the aircraft modes of motion, and the last bit is best taught by way of an example.

# Final Examples#

## Aircraft Modes from Stability Derivatives#

Consider the following aircraft data for the Pipe Cherokee PA-28-180 flying at 50m/s.

$$m$$ = 1090kg

$$S$$ = 15

$$I_{yy}$$ = 1700

$$I_{xx}$$ = 3100

$$I_{zz}$$ = 1400

$$I_{xz}$$ = 0

$$U_e$$ = 50

$$\rho$$=1.06

$$\bar{c}$$ = 1.6

$$C_{L_0}$$ = 0.543

$$C_{D_0}$$ = 0.0615

$$b$$ = 9.11

$$X_u$$ = -0.06728

$$Z_u$$ = -0.396

$$M_u$$ = 0.0

$$X_w$$ = 0.02323

$$Z_w$$ = -1.729

$$M_w$$ = -0.2772

$$X_q$$ = 0.0

$$Z_q$$ = -1.6804

$$M_q$$ = -2.207

$$M_{\dot{w}}$$ = -0.0197

$$Z_{\delta_e}$$ = -17.01

$$M_{\delta_e}$$ = -44.71

$$Y_v$$ = -0.1444

$$L_v$$ = -0.1166

$$N_v$$ = 0.174

$$L_p$$ = -2.283

$$N_p$$ = -1.732

$$L_r$$ = 1.053

$$N_r$$ = -1.029

$$Y_{\delta_r}$$ = 2.113

$$L_{\delta_r}$$ = 0.6133

$$N_{\delta_r}$$ = -6.583

$$L_{\delta_a}$$ = 3.101

$$N_{\delta_a}$$ = 0.0

The following example has been completed using Python, and I’m a better coder than when I first wrote the examples in MATLAB nearly three years ago, so my codes look more like programs than a list of calculator instructions.

I know that this can make things seem less accessible to some, but it’s smarter to keep data in a class like below.

# The means that I'm using the compute this is to store the Cherokee stability derivatives in a class,
# and then make a function to produce the stability deritatives. This is the easiest way to make readable code
# that's easily extensible to other aircraft without copying and pasting, or having all the derivatives live
# in the general namespace

import numpy as np
import sympy as sp
from IPython.display import display, Math, Latex, Markdown
import plotly.graph_objects as go

# All data for Piper Cherokee PA-28-180
class Cherokee():
h = 1200
V = 50
S = 15
m = 1090
Ixx = 1300
Izz = 1400
Ixz = 0
Iyy = 1700
CL0 = 0.543
CD0 = 0.0615
b = 9.11
cbar = 1.3
U0 = 50
theta_0=0

Ue = V

g = 9.80665

Xu = -0.06728
Zu = -0.396
Mu = 0.0
Xw = 0.02323
Zw = -1.729
Mw = -0.2772
Xq = 0.0
Zq = -1.6804
Mq = -2.207
Mdw = -0.0197
Zde = -17.01
Mde = -44.71

Yv = -0.1444
Lv = -0.1166
Nv = 0.174
Lp = -2.283
Np = -1.732
Lr = 1.053
Nr = -1.029
Ydr = 2.113
Ldr = 0.6133
Ndr = -6.583
Lda = 3.101
Nda = 0.0

# Make A and B matrices for longitudinal motion
def MakeLongitudinal(ac):
# Determine the M* derivatives (assuming Theta_e = 0)
Mu_star = ac.Mu + ac.Mdw*ac.Zu;
Mw_star = ac.Mw + ac.Mdw*ac.Zw;
Mq_star = ac.Mq + ac.Mdw*ac.Ue;
Mth_star = 0;
Mde_star = ac.Mde + ac.Mdw*ac.Zde;

if not hasattr(ac, 'theta0'):
ac.theta0 = 0

A = np.array([[ac.Xu, ac.Xw, 0, -ac.g*np.cos(ac.theta0)],
[ac.Zu, ac.Zw, ac.U0, -ac.g*np.sin(ac.theta0)],
[Mu_star, Mw_star, Mq_star, Mth_star],
[0, 0, 1, 0]])

# Make the control/input matrix, for good measure
B = np.array([, [ac.Zde], [Mde_star], ])
return A, B

def MakeLateral(ac):
# Making the starred terms
Imess = ac.Ixx * ac.Izz / (ac.Ixx * ac.Izz - ac.Ixz**2)
I2 = ac.Ixz / ac.Ixx

# Lstarred terms
Lvstar = Imess * (ac.Lv + I2 * ac.Nv)
Lpstar = Imess * (ac.Lp + I2 * ac.Np)
Lrstar = Imess * (ac.Lr + I2 * ac.Nr)
Ldrstar = Imess * (ac.Ldr + I2 * ac.Ndr)
Ldastar = Imess * (ac.Lda + I2 * ac.Nda)

# Nstarred terms
I2 = ac.Ixz / ac.Izz
Nvstar = Imess * (ac.Nv + I2 * ac.Lv)
Npstar = Imess * (ac.Np + I2 * ac.Lp)
Nrstar = Imess * (ac.Nr + I2 * ac.Lr)
Ndrstar = Imess * (ac.Ndr + I2 * ac.Ldr)
Ndastar = Imess * (ac.Nda + I2 * ac.Lda)

# Make the lateral matrix
A = np.matrix([[ac.Yv, 0, -ac.U0, ac.g*np.cos(ac.theta_0), 0],
[Lvstar, Lpstar, Lrstar, 0, 0],
[Nvstar, Npstar, Nrstar, 0, 0],
[0, 1, np.tan(ac.theta_0), 0, 0],
[0, 0, 1/np.cos(ac.theta_0), 0, 0]])

# And the control matrix
if not hasattr(ac, 'Yda'):
ac.Yda = 0

B = np.matrix([[ac.Ydr, ac.Yda], [Ldrstar, Ldastar], [Ndrstar, Ndastar], [0, 0], [0, 0]])

return A, B

# Make the two sets of matrices
Alon, Blon = MakeLongitudinal(Cherokee)
Alat, Blat = MakeLateral(Cherokee)

print("The longitudinal A and B matrices for this aircraft are")
display(Math('[A_{lon}] = ' + sp.latex(sp.Matrix(Alon).evalf(5))))
display(Math('[B_{lon}] = ' + sp.latex(sp.Matrix(Blon).evalf(5))))

print("The lateral/directional A and B matrices for this aircraft are")
display(Math('[A_{lat}] = ' + sp.latex(sp.Matrix(Alat).evalf(5))))
display(Math('[B_{lat}] = ' + sp.latex(sp.Matrix(Blat).evalf(5))))

The longitudinal A and B matrices for this aircraft are

$\begin{split}\displaystyle [A_{lon}] = \left[\begin{matrix}-0.06728 & 0.02323 & 0 & -9.8067\\-0.396 & -1.729 & 50.0 & 0\\0.0078012 & -0.24314 & -3.192 & 0\\0 & 0 & 1.0 & 0\end{matrix}\right]\end{split}$
$\begin{split}\displaystyle [B_{lon}] = \left[\begin{matrix}0\\-17.01\\-44.375\\0\end{matrix}\right]\end{split}$
The lateral/directional A and B matrices for this aircraft are

$\begin{split}\displaystyle [A_{lat}] = \left[\begin{matrix}-0.1444 & 0 & -50.0 & 9.8067 & 0\\-0.1166 & -2.283 & 1.053 & 0 & 0\\0.174 & -1.732 & -1.029 & 0 & 0\\0 & 1.0 & 0 & 0 & 0\\0 & 0 & 1.0 & 0 & 0\end{matrix}\right]\end{split}$
$\begin{split}\displaystyle [B_{lat}] = \left[\begin{matrix}2.113 & 0\\0.6133 & 3.101\\-6.583 & 0\\0 & 0\\0 & 0\end{matrix}\right]\end{split}$
# Get the eigenvalues
eigs_lon, eigv_lon = np.linalg.eig(Alon)

# We know these will always be a complex pair, so can separate. You could get the second one
# by fixing the index, but I like to do boolean indexing
eig1 = eigs_lon
eig2 = eigs_lon[eigs_lon.real != eig1.real]

# A quick function to print the single eigenvalue as a complex pair
def makecomplexpair(eigin, roundto=4):
return str(eigin.round(roundto)).replace('+', '±').strip('(').strip(')')

print(f"The eigenvalues associated with longtudinal motion are the complex pair {makecomplexpair(eig1)}, and {makecomplexpair(eig2)}")

The eigenvalues associated with longtudinal motion are the complex pair -2.4663±3.4056j, and -0.0279±0.2452j

# do the same for the lateral-directional modes
# Get the eigenvalues
eigs_lat, eigv_lat = np.linalg.eig(Alat)

# Find the dutch roll - the only part with a non-zero imaginary component
eig_DR = eigs_lat[eigs_lat.imag > 0]
other_lat_eigs = eigs_lat[eigs_lat.imag == 0]

print(f"The eigenvalues associated with lateral-direction motion are:\nThe complex pair {makecomplexpair(eig_DR)}")
# This next line is complicated but elegant because I like to enjoy how wonderfully easy it is to manipulate strings in Python by contrast to MATLAB
print(f"and {len(other_lat_eigs)} other real eigenvalues:  {' '.join(str(other_lat_eigs.real.round(4)).split())[2:].replace(' ', ' and ').strip('[').strip(']')}")

The eigenvalues associated with lateral-direction motion are:
The complex pair -0.3468±3.3718j
and 3 other real eigenvalues:  0. and -2.7823 and 0.0194


### Longitudinal Modes#

The two eigenvalues are -2.4663±3.4056j and -0.0279±0.2452j, and we can easily discriminate between the two modes because of the respective size. Recall that the damped natural frequency is given by the imaginary part of the eigenvalue, so we can see:

# A bit of logic to determine which mode is which
if eig1.imag > eig2.imag:
SP = eig1
PH = eig2
else:
SP = eig2
PH = eig1

print(f"From inspection of the imaginary part of the eigenvalues, {makecomplexpair(SP)} is clearly the short period mode, with a period of {2*np.pi/SP.imag:1.4}s\n")
print(f"From inspection of the imaginary part of the eigenvalues, {makecomplexpair(PH)} is clearly the Phugoid mode, with a period of {2*np.pi/PH.imag:1.4}s\n")

From inspection of the imaginary part of the eigenvalues, -2.4663±3.4056j is clearly the short period mode, with a period of 1.845s

From inspection of the imaginary part of the eigenvalues, -0.0279±0.2452j is clearly the Phugoid mode, with a period of 25.63s


Since both modes are stable, the time to half amplitude can be found from the real part of the eigenvalue

def timetohalf(eig):
return -0.69/eig.real

for e, n in zip([SP, PH], ['Short Period', "Phugoid"]):
print(f"The time to half amplitude of the {n} mode is {timetohalf(e):1.2f}s")

The time to half amplitude of the Short Period mode is 0.28s
The time to half amplitude of the Phugoid mode is 24.75s


The undamped natural frequency, $$\omega_n$$, is given by the absolute value of the eigenvalue,

$\omega_n=\sqrt{\Re\left({\lambda}\right)^2 + \Im\left({\lambda}\right)^2}$

and hence the damping ratio, $$\zeta$$ can be determine either from comparison with the standard form characteristic equation and

$$\zeta = -\frac{\Re(\lambda)}{\omega_n}$$

or from the definition of the damped natural frequency

$$\zeta = \sqrt{1-\frac{\omega_d^2}{\omega_n^2}}$$

for e, n in zip([SP, PH], ['Short Period', "Phugoid"]):
print(f"The undamped natural frequency of the {n} mode is {np.abs(e):1.4f}rad/s")

zeta = -e.real/np.abs(e)
print(f"which gives a damping ratio, from method 1 of {zeta:1.4f}")

zeta = np.sqrt(1-e.imag**2/np.abs(e)**2)
print(f"and a damping ratio, from method 2 of {zeta:1.4f}\n")

The undamped natural frequency of the Short Period mode is 4.2048rad/s
which gives a damping ratio, from method 1 of 0.5865
and a damping ratio, from method 2 of 0.5865

The undamped natural frequency of the Phugoid mode is 0.2468rad/s
which gives a damping ratio, from method 1 of 0.1130
and a damping ratio, from method 2 of 0.1130


the two methodologies obviously are, and were always going to be equivalent, but it’s nice to show redundancy in methods.

### Lateral Directional Modes#

From the A matrix for lateral directional motion the eigenvalues

-0.3468±3.3718j, 0.0,
-2.7823, 0.0194

are found. The first value is the only complex pair, and therefore has to be the Dutch Roll mode (the only oscillatory lateral/directional mode). We can see that is stable (negative real part) and the following can be yielded using the same as before:

print(f"From inspection of the imaginary part of the eigenvalue, {makecomplexpair(eig_DR)}, the Dutch Roll mode, has a period of {2*np.pi/eig_DR.imag:1.4}s\n\
and a time to half amplitude of {-0.69/eig_DR.real:1.4}s.")

From inspection of the imaginary part of the eigenvalue, -0.3468±3.3718j, the Dutch Roll mode, has a period of 1.863s
and a time to half amplitude of 1.99s.


For the remaining eigenvalues, one simply denotes the aircraft’s ability to yaw - the zero value shows that the aircraft has directional freedom and no real concept of heading stability (you could at most say the aircraft possesses neutral heading stability). The other two eigenvalues are real-only; -2.7823 and 0.0194.

We know that one denotes the spiral mode and one denotes the roll mode. The absolute magnitude of the damping indicates that -2.7823 is the roll mode, whilst it can also be seen that 0.0194 is unstable so cannot be the roll mode.

In summary:

print(f"For this aircraft the spiral mode is unstable with a time to double amplitude of {0.69/other_lat_eigs.real:1.2f}s")

print(f"For this aircraft the roll mode is stable (which is *has* to be) with a time to half amplitude of {0-0.69/other_lat_eigs.real:1.2f}s")

For this aircraft the spiral mode is unstable with a time to double amplitude of 35.56s
For this aircraft the roll mode is stable (which is *has* to be) with a time to half amplitude of 0.25s


## Extracting Aircraft Modes from Data#

Say you’ve been given some data for an aircraft subjected to a roll disturbance in a datafile - the information you have is roll attitude vs time in a text file.

The data is the free-response of the aircraft, just as seen in the examples following the impulsive forcing of the aircraft.

Plotting these data vs time is the first step, and we can see: