Model

Model Overview

This page summarizes the physics model implemented in src/simulator/. For full derivations see Equations and the docs/ folder in the main repository.

The five modeling assumptions

From the source paper (Sec. 2), the 6-DOF model rests on:

#	Assumption	Consequence
a	Rigid body	6 DOF (3 translation + 3 rotation) fully describe motion; no structural modes
b	Body-fixed reference frame	Moments of inertia are constant in this frame
c	Aerodynamic coefficients computed in body axes	Table lookups apply directly to body-axis states
d	Earth model (ellipsoidal, rotating, gravity)	Optional fidelity toggle (include_earth_rotation)
e	Altitude-varying atmosphere	Density/Mach change substantially over a multi-km trajectory

State vector (12 states)

x = [u, v, w,       body-axis velocity           [m/s]
     p, q, r,       body-axis angular rates       [rad/s]
     φ, θ, ψ,       Euler angles (roll,pitch,yaw)  [rad]
     N, E, D]       geodetic position               [m]

Coordinate frames

• Body-fixed frame (F_B) — origin at C.G., x forward, y right, z down.
• Local geodetic / NED frame (F_E) — origin at launch, x North, y East, z Down.

Attitude uses a 3-2-1 (yaw→pitch→roll) Euler sequence. The direction cosine matrix L_BE (implemented in src/simulator/frames.py) transforms vectors from body axes to the geodetic frame and is used both for navigation (v_E = L_BE @ v_B) and for resolving Earth's rotation into body axes.

Gimbal lock occurs at pitch angle θ = ±90°, where the Euler-rate kinematic matrix is singular — see the Coordinate Frames lab page and Assignment Exercise 5 for a hands-on exploration.

Equation groups

Four coupled groups produce dx/dt, mirroring the source paper's Fig. 1 block diagram:

1. Translational dynamics (u̇,v̇,ẇ) — Newton's second law in a rotating body frame: thrust + aerodynamic force + gravity, plus Coriolis-like ω×v coupling terms.
2. Rotational dynamics (ṗ,q̇,ṙ) — Euler's equations for a rigid body, with gyroscopic coupling through the inertia tensor. For the paper's axisymmetric case (Iyy=Izz), roll decouples from pitch/yaw, but pitch and yaw remain gyroscopically coupled whenever spin p is nonzero ("coning"/epicyclic motion).
3. Kinematics equation (φ̇,θ̇,ψ̇) — converts body angular rates into Euler-angle rates.
4. Navigation equation (Ṅ,Ė,Ḋ) — rotates body velocity into the geodetic frame to update position.

Atmosphere

src/simulator/atmosphere.py implements the 1976 US Standard Atmosphere's troposphere (0–11 km, constant lapse rate) and lower stratosphere (11–20 km, isothermal) layers, giving temperature, pressure, density, and sonic speed as functions of altitude. This stands in for the source paper's unpublished atmosphere table (see FAQ).

Aerodynamics

src/simulator/aerodynamics.py interpolates five coefficients against Mach number: CA (axial force), CN_alpha (normal-force curve slope), Cl_p (roll damping), Cmq (pitch/yaw damping), Cm_alpha (pitching-moment curve slope — the key static-stability parameter for a fin-stabilized body). These feed the standard force/moment build-up equations using dynamic pressure q̄ = ½ρV²:

Tx_aero = -q̄·S·CA
Ty_aero =  q̄·S·CN_alpha·β
Tz_aero = -q̄·S·CN_alpha·α
L  = q̄·S·D·(Cl_p·p·D/2V)
M  = q̄·S·D·(Cm_alpha·α + Cmq·q·D/2V)
N  = q̄·S·D·(Cm_alpha·β + Cmq·r·D/2V)

⚠️ These specific coefficient values are a representative reconstruction, not the paper's exact Table 1 data — see FAQ.

Numerical integration

Three integrators are implemented in src/simulator/integrators.py: forward Euler (O(dt) global error), classical RK4 (O(dt⁴) global error, hand-written), and an adaptive scipy.integrate.solve_ivp (RK45) wrapper used as a ground-truth reference. See the Numerical Integrator lab page for a live convergence/stability demo.

Dispersion / sensitivity analysis

src/simulator/dispersion.py reproduces the paper's Table 2 one-parameter- at-a-time sensitivity sweeps (10 of its 12 published uncertainty parameters), returning range/drift/radial impact-point error vs. each parameter's perturbation — directly comparable to the paper's Figs. 10–21.