Reference: FM04.pdf, Sec. 2 “Mathematical Model”, Fig. 1.
Given a rocket’s mass properties, its aerodynamic coefficients, its propulsion profile, and the atmosphere it flies through, predict its trajectory: position, velocity, attitude and angular rates as functions of time, from launch to ground impact.
This is an initial value problem (IVP): a system of coupled nonlinear ordinary differential equations (ODEs), integrated forward in time from a known state at t=0 (the “Initial Conditions” box in Fig. 1) until a stopping condition is reached (altitude ≤ 0, the “Stop Simulation” box).
The paper’s 6-DOF model rests on five explicit assumptions. Each has a teaching point:
| # | Assumption | Why it matters |
|---|---|---|
| a | The flying body is rigid. | No structural bending/flexing modes — 6 DOF (3 translation + 3 rotation) fully describe the motion. A flexible-body model would need many more states. |
| b | All equations are referred to a body-fixed frame. | Moments of inertia are constant in this frame (they’d be time-varying in an inertial frame as the body rotates) — this is why body axes are the natural choice for rotational dynamics. |
| c | Aerodynamic coefficients are computed in the body-fixed frame. | Forces/moments from tables (Table 1) apply directly to body-axis states without extra rotation. |
| d | The Earth model includes ellipsoidal shape, rotation, gravity variation. | For short-range artillery trajectories (tens of km, under two minutes of flight) Earth curvature/rotation effects are small but not always negligible — see equations_of_motion.py’s optional include_earth_rotation flag (which now implements both the Eq. (3) Coriolis term on body rates AND a full rotating/curved-Earth Navigation Equation with extra V_N,V_E,V_D states — spherical-Earth approximation, see docs/coordinate-systems.md) and the flat-Earth vs rotating-Earth assignment in assignments.md. |
| e | The atmosphere varies with altitude (temperature, sonic speed, density). | Aerodynamic forces scale with dynamic pressure q̄ = ½ρV² and Mach number M = V/a; both change substantially over a multi-km-altitude trajectory, so a constant-atmosphere assumption would be a poor approximation. |
This lab’s implementation (src/simulator/equations_of_motion.py) uses a
12-state vector, matching every quantity that flows through Fig. 1’s
signal-flow diagram:
x = [u, v, w, body-axis velocity components [m/s]
p, q, r, body-axis angular rates [rad/s]
φ, θ, ψ, Euler angles (roll, pitch, yaw) [rad]
N, E, D] geodetic position (North/East/Down) [m]
When include_earth_rotation=True, three more states are appended
(V_N, V_E, V_D, geodetic-frame velocity) for the rotating/curved-Earth
Navigation Equation — see docs/coordinate-systems.md. This is off by
default and the 12-state form above is unchanged in that case.
See state-variables page in the GUI and
equations.md for what produces the rate of each of
these 12 quantities.
u̇, v̇, ẇ from forces (thrust +
aerodynamic + gravity) and Coriolis-like body-rate coupling terms.ṗ, q̇, ṙ from aerodynamic moments and gyroscopic
coupling through the inertia tensor, using the paper’s general
(non-axisymmetric) form with cross-inertia term Izx. RocketParams
defaults Izx=0, which makes this collapse exactly to the simpler
axisymmetric-body case-study formulas.φ̇, θ̇, ψ̇ from body rates — how attitude
angles evolve given the angular velocity.Ṅ, Ė, Ḋ from body velocity rotated into the
geodetic frame by the direction cosine matrix L_BE.These four groups are coupled: forces/moments depend on velocity,
attitude and altitude; velocity and attitude evolve from those same
forces/moments. That coupling is exactly why this must be solved
numerically (Section numerical-methods.md) rather
than in closed form.
state_derivative(t, x, rocket, atmo, aero, ...) -> dx/dt
implements all four groups in one function
(src/simulator/equations_of_motion.py), calling out to:
rocket.py — time-varying mass/inertia (boost vs. free-flight phase),atmosphere.py — density, sonic speed vs. altitude,aerodynamics.py — Mach-interpolated coefficients → forces/moments,frames.py — the direction cosine matrix and Euler kinematics.An integrator (integrators.py) then repeatedly evaluates this function to
advance the state in time — see numerical-methods.md.
The block-diagram style of Fig. 1 is itself worth teaching: it is exactly how you would structure this code in any language — a small “physics kernel” function (forces/moments → accelerations) wrapped by an integration loop, with clearly separated sub-models (atmosphere, aerodynamics, mass properties) that could each be swapped out independently. This separation-of-concerns is why
src/simulator/has one file per sub-model instead of one monolithic script.
See assignments.md Exercises 1, 6.