The Apogee simulator uses the US Standard Atmosphere 1976 (USSA76) for atmospheric properties.

Overview¶

Atmospheric drag significantly affects rocket trajectories during the lower portion of ascent. We need two key properties:

1. Density \(\rho(h)\) - for drag force calculation
2. Speed of Sound \(a_s(h)\) - for Mach number

USSA76 Model¶

The US Standard Atmosphere 1976 defines temperature, pressure, and density as functions of geometric altitude.

Temperature Profile¶

The atmosphere is divided into layers with different lapse rates:

Density Computation¶

Within each layer, density follows from the hydrostatic equation and ideal gas law.

Isothermal Layer: $$ \rho(h) = \rho_b \exp\left(-\frac{g_0 M (h - h_b)}{R^* T_b}\right) $$

Gradient Layer: $$ \rho(h) = \rho_b \left(\frac{T_b}{T_b + L(h - h_b)}\right)^{1 + g_0 M / (R^* L)} $$

Where: - \(\rho_b\) = density at layer base [kg/m³] - \(T_b\) = temperature at layer base [K] - \(L\) = lapse rate [K/m] - \(M\) = molar mass of air (0.0289644 kg/mol) - \(R^*\) = universal gas constant (8.31447 J/(mol·K))

Speed of Sound¶

For an ideal gas:

Where: - \(\gamma = 1.4\) (ratio of specific heats for air) - \(R = R^*/M = 287.05\) J/(kg·K)

Implementation¶

Lookup Table Approach¶

Rather than evaluating the piecewise functions repeatedly, we build a precomputed table.

Implementation:

def build_atmosphere_table(*, z_max_m: float, dz_m: float) -> AtmosphereTable:
    """Build discrete USSA76 table for interpolation [Section 4.1]."""
    z = np.arange(0.0, z_max_m + dz_m, dz_m, dtype=float)
    ds = ussa1976.compute(z=z, variables=["rho", "cs"])

    rho = np.asarray(ds["rho"].values, dtype=float)
    cs = np.asarray(ds["cs"].values, dtype=float)

    return AtmosphereTable(z_m=jnp.asarray(z), rho=jnp.asarray(rho), cs=jnp.asarray(cs))

Interpolation¶

Linear interpolation between table points:

def rho_at(self, h_m: Array) -> Array:
    """Interpolate density rho(h) [Section 4.1, Eq. 15]."""
    return jnp.interp(h_m, self.z_m, self.rho, left=self.rho[0], right=self.rho[-1])

def cs_at(self, h_m: Array) -> Array:
    """Interpolate speed of sound cs(h) [Section 4.1, Eq. 15]."""
    return jnp.interp(h_m, self.z_m, self.cs, left=self.cs[0], right=self.cs[-1])

Default Parameters¶

Aerodynamic Drag¶

Mach Number¶

Physical meaning: - \(M < 1\): Subsonic - \(M = 1\): Sonic (transonic region) - \(M > 1\): Supersonic - \(M > 5\): Hypersonic

Drag Force¶

Where: - \(\rho(h)\) = atmospheric density [kg/m³] - \(v\) = velocity magnitude [m/s] - \(C_D(M)\) = drag coefficient (function of Mach) - \(A_{\text{ref}}\) = reference area [m²]

Dynamic Pressure¶

A key derived quantity:

Max-Q occurs when \(q\) is maximum (typically around 80-90 seconds into flight).

Implementation:

def compute_derived(*, y: Array, earth: EarthParams, stage: StageParams, atmos: AtmosphereTable) -> Derived:
    """Compute derived scalars per LaTeX Report (Section 4)."""
    r, v = y[_R], y[_V]

    h = r - earth.r_e
    h_clip = jnp.maximum(h, 0.0)

    rho = atmos.rho_at(h_clip)
    cs = atmos.cs_at(h_clip)
    mach = v / cs

    cd = stage.cd(mach)
    drag = 0.5 * rho * cd * stage.a_ref * v * v
    q = 0.5 * rho * v * v

    return Derived(h=h, rho=rho, cs=cs, mach=mach, drag=drag, q=q)

Drag Coefficient Model¶

The drag coefficient varies with Mach number:

Constant \(C_D\)¶

Simplest model (used in Apogee):

Implementation:

class ConstantCd:
    def __init__(self, cd: float = 0.3):
        self.cd = cd

    def __call__(self, mach: Array) -> Array:
        return jnp.full_like(mach, self.cd)

Transonic Variation (Future)¶

A more realistic model includes transonic drag rise:

Altitude Effects¶

Density vs. Altitude¶

Above ~100 km (Kármán line), atmospheric effects become negligible.

High-Altitude Handling¶

For altitudes above the table maximum:

# Add an extra point just above z_max so extrapolation above the table
# approaches vacuum (rho=0) while remaining C0-continuous at z_max.
z = np.concatenate([z, np.array([z[-1] + dz_m], dtype=float)])
rho = np.concatenate([rho, np.array([0.0], dtype=float)])

Falcon 9 Aerodynamic Parameters¶

Implementation:

diameter: float = 3.7  # Vehicle diameter [m]

a_ref = math.pi * (params.diameter / 2) ** 2  # Reference area

References¶

1. NOAA/NASA/USAF (1976). U.S. Standard Atmosphere, 1976. U.S. Government Printing Office. NOAA-S/T 76-1562.

2. Anderson, J.D. (2016). Introduction to Flight. McGraw-Hill. Chapter 3 - "The Standard Atmosphere".

3. Python Package: ussa1976

Next Steps¶

• Equations of Motion - Where drag is used
• Orbital Mechanics - Beyond atmosphere