Aerodynamic and Wind Policies
All external forces and moments on the rigid body are expressed as spatial wrenches in the body frame:
where \(\boldsymbol{\tau}_B \in \mathbb{R}^3\) are moments (indices 0–2) and \(\mathbf{F}_B \in \mathbb{R}^3\) are forces (indices 3–5), consistent with the Featherstone spatial-vector convention used throughout Aetherion.
Aerodynamic wrenches are computed by AeroPolicy objects passed as
template parameters to RigidBody::VectorField.
Wind effects are factored out into a separate WindModel layer so that
any combination of drag model and wind profile can be composed without
modifying the core integrator.
Atmosphere-Relative Velocity
The key quantity for all aerodynamic force computations is the atmosphere-relative velocity in the body frame. The atmosphere co-rotates with Earth; a body orbiting without atmosphere-relative motion experiences zero drag:
where
\(\mathbf{v}_B = \nu_B^{3:5}\) is the body-frame linear velocity (lower half of the body twist),
\(R \in SO(3)\) is the rotation from body to ECI frame,
\(\boldsymbol{\omega}_\oplus = [0,\,0,\,\omega_\oplus]^\top\) is Earth’s spin vector in the ECI frame (\(\omega_\oplus = 7.2921150 \times 10^{-5}\) rad/s),
\(\mathbf{r}_\mathrm{ECI}\) is the vehicle’s ECI position, and
\(\mathbf{v}_\mathrm{wind,ECEF}\) is the ambient wind in the ECEF frame, rotated to ECI at the current Earth Rotation Angle \(\theta_\mathrm{ERA}(t) = \omega_\oplus\,t\).
When there is no wind, \(\mathbf{v}_\mathrm{wind,ECEF} = \mathbf{0}\) and
the expression reduces to the Earth-surface velocity subtraction used in
DragOnlyAeroPolicy.
The true airspeed (atmosphere-relative speed) reported in the output CSV is
All drag and damping formulas below use \(\mathbf{v}_\mathrm{rel}\) and \(V\).
Drag-Only Policy (DragOnlyAeroPolicy)
Models a sphere or any body with purely resistive (no-lift) aerodynamics. The drag force in the body frame is
where \(\rho(h)\) is the US Standard Atmosphere 1976 density at the geocentric altitude \(h \approx \|\mathbf{r}_\mathrm{ECI}\| - R_e\).
The corresponding body-frame wrench has zero moment components:
Parameters: \(C_D\) (drag coefficient), \(S_\mathrm{ref}\) (reference area [m²]).
AD-safe: yes — the square-root singularity at \(V = 0\) is avoided by adding a small floor \(\varepsilon = 10^{-15}\) before taking \(\sqrt{V^2 + \varepsilon}\).
Used in: Scenario 6 (sphere with drag), Scenarios 9 & 10 (cannonball).
Drag and Rotary-Damping Policy (BrickDampingAeroPolicy)
Extends Drag-Only Policy (DragOnlyAeroPolicy) with aerodynamic moment damping proportional to the body angular rates. This is the standard non-dimensional rate formulation used for the tumbling-brick scenarios.
Damping moments use the standard dimensionless rate convention:
where \(b\) is the span reference length and \(\bar{c}\) is the chord reference length.
The rolling, pitching, and yawing damping moments are
The full wrench combines drag force and damping moments:
Parameters: \(C_D,\,S\) (drag); \(b,\,\bar{c}\) (reference lengths); \(C_{l_p},\,C_{m_q},\,C_{n_r}\) (rotary derivatives).
Note: Cross-coupling derivatives (\(C_{l_r},\,C_{n_p}\), etc.) are
stored in AerodynamicParameters and default to zero.
Used in: Scenario 3 (tumbling brick with damping).
Wind Models
All wind-aware policies receive the ambient wind via a WindModel — a lightweight struct providing a single templated method:
template<class S>
Eigen::Matrix<S,3,1> velocity_ecef(
const Eigen::Matrix<S,3,1>& r_eci, // ECI position [m]
S t_s) const; // simulation time [s]
The method returns the wind velocity in the ECEF frame [m/s] and must
be callable for both S = double (evaluation) and S = CppAD::AD<double>
(Jacobian recording).
Three built-in models are provided.
ZeroWind
Calm atmosphere. Reduces Wind-Aware Drag Policy (WindAwareDragPolicy<Wind>) to Drag-Only Policy (DragOnlyAeroPolicy).
ConstantECEFWind
The helper ConstantECEFWind::from_ned() converts a NED wind
specification at the launch geodetic position into the ECEF frame:
where \(\phi_0,\,\lambda_0\) are the launch geodetic latitude and longitude.
Used in: Scenario 7 (steady 20 ft/s eastward wind).
LinearWindShear
Wind varies linearly with geocentric altitude in each NED component:
where \(h = \|\mathbf{r}_\mathrm{ECI}\| - R_e\) is the geocentric altitude, \(\alpha\) is the altitude gradient [m/s per m], and \(\beta\) is the sea-level intercept [m/s]. The NED wind is converted to ECEF once at construction from the launch latitude/longitude.
For NASA TM-2015-218675 Scenario 8 the formula is given explicitly:
“Vwind = (0.003h − 20) ft/s from west; h is height MSL in ft.”
In SI units: \(v_E(h) = 0.003\,h_\text{m} - 6.096\) m/s (\(\alpha_E = 0.003\), \(\beta_E = -6.096\), \(\alpha_N = \beta_N = 0\)). The wind is eastward above 2 032 m (6 667 ft), zero there, and westward below.
Used in: Scenario 8 (2D wind shear).
GeodesicCallbackWind
Wraps an arbitrary user-provided callable
std::function<Eigen::Vector3d(double lat_rad,
double lon_rad,
double alt_m,
double t_s)>
that returns the NED wind at the current geodetic position. This enables integration with external weather services or pre-computed atmospheric data fields.
AD compatibility strategy — the implicit Newton solver in the Radau IIA
integrator records a CppAD tape that calls the wind model with
\(S = \mathtt{AD<double>}\) values. Since an external API cannot be
called with AD arguments, the implementation uses the
frozen-Jacobian approximation:
During the double path (residual evaluation), the callback is invoked with plain
doublecoordinates and the resulting ECEF wind vector is cached.During the AD path (Jacobian recording), the cached value is returned as a constant.
This is equivalent to treating the wind as piecewise-constant over each Newton step — a standard and physically reasonable approximation since wind fields change on timescales much longer than an integration step.
To add a custom wind profile, implement velocity_ecef(r_eci, t) and
register the type:
// In namespace Aetherion::FlightDynamics:
template<> struct is_wind_model<MyWind> : std::true_type {};
Wind-Aware Drag Policy (WindAwareDragPolicy<Wind>)
General drag policy that accepts any registered WindModel as a template
parameter. It subsumes Drag-Only Policy (DragOnlyAeroPolicy) (using ZeroWind) and
extends it to steady or shear winds.
The atmosphere-relative velocity is computed as
where \(R_z(\theta) = R_\mathrm{ECEF\to ECI}(\theta)\) rotates the ECEF wind into the ECI frame at the current ERA \(\theta_\mathrm{ERA}(t) = \omega_\oplus\,t\).
The drag force is then identical to Drag-Only Policy (DragOnlyAeroPolicy):
The policy is instantiated as, for example:
using ShearPolicy =
WindAwareDragPolicy<LinearWindShear>;
ShearPolicy p{ CD, S_ref,
LinearWindShear::from_ned(
0.0, 0.003, // gradient (N, E) [m/s per m]
0.0, -6.096, // intercept (N, E) [m/s at h=0]
lat0, lon0) };
Snapshot Air-Data Corrections
MakeSnapshot1() computes the reported true airspeed (TAS),
Mach number, and dynamic pressure from the Earth-relative NED velocity
\(\mathbf{v}_\mathrm{NED}\):
For wind-aware simulators (Scenarios 7 and 8), the snapshot
snapshot() override subtracts the ambient wind before computing
air-data quantities:
so that the reported TAS, Mach, and \(\bar{q}\) match the NASA reference convention. For Scenario 8 (power-law shear), the wind-at-altitude correction uses the same exponent \(n\) as the policy.
Policy Composition Summary
Policy |
Drag force |
Damping moments |
Wind subtraction |
Scenario(s) |
|---|---|---|---|---|
|
✗ |
✗ |
✗ |
1, 2, 9, 10 (dragless) |
|
✓ |
✗ |
Earth rotation only |
6 |
|
✓ (opt.) |
✓ |
Earth rotation only |
3 |
|
✓ |
✗ |
Earth rotation only |
— |
|
✓ |
✗ |
Earth rotation + constant wind |
7 |
|
✓ |
✗ |
Earth rotation + linear altitude shear |
8 |
|
✓ |
✗ |
Earth rotation + API wind |
user-defined |