API

γ_AIR = 1.4

Adiabatic index or ratio of specific heats (dry air at 20º C).

R_AIR = 287.05287  (J/(Kg·K))

Specific gas constant for dry air.

gD = 9.80665  (m/s²)

Down component of gravity acceleration at Earth surface at 45º geodetic latitude.

1. Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control and simulation: dynamics, controls design, and autonomous systems. John Wiley & Sons. Equation (page 33)

Ellipsoid(a, b, f, finv, e2, ϵ2, name)
Ellipsoid(a, finv, name)

Earth ellipsoid model.

Fields

• a::Real: semi-major axis (m).
• b::Real: semi-minor axis (m).
• f::Real: flattening ($f$).
• finv::Real: inverse of flattening ($f^{-1}$).
• e2::Real: eccentricity squared.
• ϵ2::Real: second eccentricity squared.
• name::String: ellipsoid name.

Notes

$a$ is the semi-major axis, $b$ is the semi-minor axis. Ellipsoids are normally defined by $a$ and $f^{-1}$. The rest of parameters are derived.

• Flattening:

\[f = 1 - \frac{b}{a}\]

• Eccentricity (first eccentricity):

\[e^2 = \frac{a^2 - b^2}{a^2} = f (2 - f)\]

• Second eccentricity:

\[ϵ^2 = \frac{a^2 - b^2}{b^2} = \frac{f (2 - f)}{(1 - f)^2}\]

References

1. Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation: dynamics, controls design, and autonomous systems. John Wiley & Sons. Section 1.6 (page 23)
2. Rogers, R. M. (2007). Applied mathematics in integrated navigation systems. American Institute of Aeronautics and Astronautics. Chapter 4 (Nomenclature differs from 1).
3. Bowring, B. R. (1976). Transformation from spatial to geographical coordinates. Survey review, 23(181), 323-327.

ecef2llh(xecef, yecef, zecef; ellipsoid=WGS84)

Transform ECEF coordinates to geodetic latitude, longitude (rad) and ellipsoidal height (m) for the given ellipsoid (default ellipsoid is WGS84).

Notes

• The transformation is direct without iterations as [1] introduced the need to iterate for near Earth positions.
• [2] is an update of increased accuracy of [1]. The former is used in this implementation although the latter implementation is commented in the code.
• Model becomes unstable if latitude is close to 90º. An alternative equation can be found in [2] equation (16) but has not been implemented.

References

1. Bowring, B. R. (1976). Transformation from spatial to geographical coordinates. Survey review, 23(181), 323-327.
2. Bowring, B. R. (1985). The accuracy of geodetic latitude and height equations. Survey Review, 28(218), 202-206.

llh2ecef(lat, lon, height; ellipsoid=WGS84)

Transform geodetic latitude, longitude (rad) and ellipsoidal height (m) to ECEF for the given ellipsoid (default ellipsoid is WGS84).

References

1. Rogers, R. M. (2007). Applied mathematics in integrated navigation systems. American Institute of Aeronautics and Astronautics. (Page 75, equations 4.20, 4.21, 4.22)

body2horizon(x, y, z, ψ, θ, ϕ)

Transform the vector coordintes (x, y, z) given in body axis to local horizon given the Euler angles (ψ, θ, ϕ) (rad).

body2horizon(x, y, z, q0, q1, q2, q3)

Transform the vector coordintes (x, y, z) given in body axis to local horizon given the quaternions $q_0, q_1, q_2, q_3$.

body2wind(x, y, z, α, β)

Transform the vector coordintes (x, y, z) given in body axis to wind given the angle of attack (α) and the angle of sideslip (β) (rad).

ecef2horizon(x, y, z, lat, lon)

Transform the vector coordintes (x, y, z) given in ECEF (Earth Fixed Earth Centered) coordinates to local horizon coordinates using geodetic latitude and longitude (rad).

euler_angles(q0, q1, q2, q3)

Calculate Euler angles (ψ, θ, ϕ) (rad) given the quaternions $q_0, q_1, q_2, q_3$.

horizon2body(x, y, z, ψ, θ, ϕ)

Transform the vector coordintes (x, y, z) given in local horizon axis to body given the Euler angles (ψ, θ, ϕ) (rad).

horizon2body(x, y, z, q0, q1, q2, q3)

Transform the vector coordintes (x, y, z) given in local horizon axis to body given the quaternions $q_0, q_1, q_2, q_3$.

horizon2ecef(x, y, z, lat, lon)

Transform the vector coordintes (x, y, z) given in local horizon axis to ECEF (Earth Fixed Earth Centered) using geodetic latitude and longitude (rad).

horizon2wind(x, y, z, χ, γ, μ)

Transform the vector coordintes (x, y, z) given in local horizon axis to wind given the velocity angles (χ, γ, μ) (rad).

quaternions(ψ, θ, ϕ)

Calculate quaternion representation given the Euler angles (ψ, θ, ϕ) (rad).

rotation_matrix_zyx(q0, q1, q2, q3)

Calculate the rotation matrix $R_{ab}$ that transforms a vector in b frame to a frame ($v_{a} = R_{ab} v_{b}$) given the quaternions $q_0, q_1, q_2, q_3$.

rotation_matrix_zyx(α1, α2, α3)

Calculate the rotation matrix $R_{ab}$ that transforms a vector in b frame to a frame ($v_{a} = R_{ab} v_{b}$) given the rotations (α1, α2, α3) (rad).

Frame b is obtained from frame a performing three intrinsic rotations of magnitude α1, α2 and α3 in ZYX order.

wind2body(x, y, z, α, β)

Transform the vector coordintes (x, y, z) given in wind axis to body given the angle of attack (α) and the angle of sideslip (β) (rad).

wind2horizon(x, y, z, χ, γ, μ)

Transform the vector coordintes (x, y, z) given in wind axis to local horizon given the velocity angles (χ, γ, μ) (rad).

atmosphere_isa(height)

Calculate temperature, pressure, density and sound velocity for the given geopotential height according to International Standard Atmosphere 1976.

References

1.U.S. Standard Atmosphere, 1976, U.S. Government Printing Office, Washington, D.C., 1976

Source: https://en.wikipedia.org/wiki/U.S.StandardAtmosphere

pqr_2_quat_dot(p, q, r, q0, q1, q2, q3)

Transform body angular velocity (p, q, r) [rad/s] to quaternion rates [1/s].

See also

pqr_2_ψθϕ_dot

References

1. Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control and simulation: dynamics, controls design, and autonomous systems. John Wiley & Sons. Equation (1.8-15) (page 51).

pqr_2_ψθϕ_dot(p, q, r, ψ, θ, ϕ)

Transform body angular velocity (p, q, r) [rad/s] to Euler angles rates (ψdot, θdot, ϕ_dot) [rad/s] given the euler angles (θ, ϕ) [rad] using kinematic angular equations.

See also

ψθϕ_dot_2_pqr, pqr_2_quat_dot

References

1. Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control and simulation: dynamics, controls design, and autonomous systems. John Wiley & Sons. Equation (1.4-4) (page 20)

rate_of_climb_constrain_no_wind(γ, α, β, ϕ)

Calculate pitch angle (θ rad) to obtain a flight path angle (γ rad) at certain angle of attack (α rad), angle of sideslip (β rad) and roll (ϕ rad).

Inertial velocity and aerodynamic velocity are equivalent in the absence of wind:

$v_{cg-e}^e = R_{eb} R_{bw} v_{cg-air}^w$

References

1. Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation: dynamics, controls design, and autonomous systems. John Wiley & Sons. (Section 3.6, equation 3.6-3, page 187)

rigid_body_acceleration(acc_P, ω, ω_dot, r_PQ)

Calcualte rigid body acceleration field.

Return the acceleration of a point Q of a rigid solid given the acceleration of a point P (accP), the rotational velocity of the solid (ω), the rotational acceleration of the solid (ωdot) and the relative position of Q with respect to P.

$a_{10}^{Q} = a_{10}^{P} + \omega_{10} \times (\omega_{10} \times r^{PQ}) + \dot{\omega}_{10} \times r^{PQ}$

being:

• $a_{10}^{Q}$ the acceleration of point Q, fixed to 1, wrt 0
• $\omega_{10}$ the angular velocity of the solid 1 wrt 0
• $\dot{\omega}_{10}$ the angular acceleration of the solid 1 wrt 0
• $r^{PQ}$ the position of Q wrt P ($r^{Q}-r^{P}$)

References

1. Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation: dynamics, controls design, and autonomous systems. John Wiley & Sons. (Section 1.3, Formaula 1.3-14c, page 26)

rigid_body_velocity(vel_P, ω, r_PQ)

Calculate rigid solid velocity field.

Return velocity of a point Q of a rigid solid given the velocity of a point P (vel_P), the rotational velocity of the solid (ω) and the relative position of Q with respect to P.

If the reference frame 1 is attached to the solid and the velocity is calculated with respect to reference frame 0:

$v_{10}^{Q} = v_{10}^{P} + \omega_{10} \times r^{PQ}$

being:

• $v_{10}^{Q}$ the velocity of point Q, fixed to 1, wrt 0
• $\omega_{10}$ the angular velocity of the solid 1 wrt 0
• $r^{PQ}$ the position of Q wrt P ($r^{Q}-r^{P}$)

Every vector needs to be expressed in the same coordinate system.

References

1. Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation: dynamics, controls design, and autonomous systems. John Wiley & Sons. (Section 1.3, page 26)

tasαβ_dot_to_uvw_dot(tas, α, β, tas_dot, α_dot, β_dot)

Obatain body velocity derivatives given velocity in wind axis and its derivatives.

See also

uvw_dot_to_tasαβ_dot

References

1. Morelli, Eugene A., and Vladislav Klein. Aircraft system identification: Theory and practice. Williamsburg, VA: Sunflyte Enterprises, 2016. Derived from equation 3.32 (page 44).

uvw_dot_to_tasαβ_dot(u, v, w, u_dot, v_dot, w_dot)

Calculate time derivatives of velocity expressed as TAS, AOA, AOS.

Notes

Note that tas here is not necessarily true air speed. Could also be inertial speed in the direction of airspeed. It will concide whith TAS for no wind.

See also

tasαβ_dot_to_uvw_dot

References

1. Morelli, Eugene A., and Vladislav Klein. Aircraft system identification: Theory and practice. Williamsburg, VA: Sunflyte Enterprises, 2016. Equation 3.33 (page 44).
2. Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control and simulation: dynamics, controls design, and autonomous systems. John Wiley & Sons. Equation (2.3-10) (page 81).

uvw_to_tasαβ(u, v, w)

Calculate true air speed (TAS), angle of attack (α) and angle of side slip (β) from velocity expressed in body axis.

Notes

This function assumes that u, v, w are the body components of the aerodynamic speed. This is not true in genreal (wind speed different from zero), as u, v, w represent velocity with respect to an inertial reference frame.

References

1. Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control and simulation: dynamics, controls design, and autonomous systems. John Wiley & Sons. Equation (2.3-6b) (page 78).

ψθϕ_dot_2_pqr(ψ_dot, θ_dot, ϕ_dot, ψ, θ, ϕ)

Transform Euler angles rates (ψdot, θdot, ϕ_dot) [rad/s] to body angular velocity (p, q, r) [rad/s] given the euler angles (θ, ϕ) [rad] using kinematic angular equations.

See also

pqr_2_ψθϕ_dot

References

1. Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control and simulation: dynamics, controls design, and autonomous systems. John Wiley & Sons. Equation (1.4-3) (page 20)

coordinated_turn_bank(ψ_dot, α, β, tas, γ[, g])

Calculate roll angle (ϕ) [rad] for a given turn rate, angle of attack, angle of sideslip, tas and flight path angle in the absence of wind for a coordinated turn bank.

Imposes sum of forces along y body axis equal to zero.

Arguments

• ψ_dot: turn rate (rad/s).
• α: angle of attack (rad).
• β: angle of sideslip (rad).
• tas: true air speed (m/s).
• γ: flight path angle (rad).
• g: gravity (m/s²). Optional. Default value is gD.

steiner_inertia(cg, inertia_g, mass, p2)

Calculate the inertia tensor of a rigid solid at point 2 (p2), given the inertia tensor (inertia_g) at the center of gravity (cg) and the mass of the system.

$r = p_2 - cg$

$I_{2} = I_{cg} + m (r^T · r I - r · r^T)$

translate_forces_moments(forces_1, moments_1, p1, p2)

Calculate equivalent moment at point 2 (p2) given the forces (forces1) and moments (moments1) at point 1 (p1).

$r = p_1 - p_2$

$M_{2} = M_{1} + r \times F_{1}$

cas2eas(cas, ρ, p)

Calculate equivalent airspeed from calibrated airspeed, density (ρ) and pressure (p) at the current altitude.

cas2tas(cas, ρ, p)

Calculate true airspeed from calibrated airspeed, density (ρ) and pressure (p) at the current altitude.

compressible_qinf(tas, p, a)

Calculate compressible dynamic pressure from Mach number and static pressure (p)

Two different models are used depending on the Mach number:

• Subsonic case: Bernouilli's equation compressible form.
• Supersonic case: to be implemented.

References

1. Ward, D. T. (1993). Introduction to flight test engineering. Elsevier Science Ltd. (page 12)
2. Fundamentals of Aerdynamics, 5th edition, J.D.Anderson Jr (page 550)

eas2cas(eas, ρ, p)

Calculate calibrated airspeed from equivalent airspeed, density (ρ) and pressure (p) at the current altitude.

eas2tas(qc, ρ)

Calculate true airspeed from equivalent airspeed and density at current altitude (ρ).

References

1. Ward, D. T. (1993). Introduction to flight test engineering. Elsevier Science Ltd. (page 13, formula 2.15)

incompressible_qinf(tas, ρ)

Calculate incompressible dynamic pressure from true airspeed (tas) and density (ρ) at current altitude.

References

1. Ward, D. T. (1993). Introduction to flight test engineering. Elsevier Science Ltd. (page 13, formula 2.14)

qc2cas(qc)

Calculate calibrated airspeed from ASI (Air Speed indicator), differential pressure between impact pressure and static pressure.

$qc = p_t - p_s$

References

1. Ward, D. T. (1993). Introduction to flight test engineering. Elsevier Science Ltd. (page 13, formula 2.13)

qc2eas(qc, p)

Calculate equivalent airspeed from ASI (Air Speed indicator), differential pressure between impact pressure and static pressure (qc = pt - ps) and p.

References

1. Ward, D. T. (1993). Introduction to flight test engineering. Elsevier Science Ltd.

qc2tas(qc, ρ, p)

Calculate true airspeed from ASI (Air Speed indicator), differential pressure between impact pressure and static pressure (qc = pt - ps), rho and p.

References

1. Ward, D. T. (1993). Introduction to flight test engineering. Elsevier Science Ltd. (page 12, based on formula 2.11)

tas2cas(cas, ρ, p)

Calculate true airspeed from calibrated airspeed, density (ρ) and pressure (p) at the current altitude.

tas2eas(tas, ρ)

Calculate equivalent airspeed from true airspeed and density at current altitude (ρ).

References

1. Ward, D. T. (1993). Introduction to flight test engineering. Elsevier Science Ltd. (page 13, formula 2.15)