Hypersonic Aero-Thermal Coupled (dim=13)

def _atmosphere_density(altitude):
    """Exponential atmosphere model."""
    h = max(altitude, 0.0)
    return _RHO0 * np.exp(-h / _H_SCALE)

def _gravity_accel(pos):
    """Spherical gravity: returns acceleration vector (toward center)."""
    r = np.linalg.norm(pos)
    r = max(r, _R_EARTH * 0.5)
    g_mag = _MU_EARTH / (r * r)
    return -g_mag * pos / r

def _mach_number(speed, altitude):
    """Approximate Mach from speed; sound speed decreases with altitude."""
    # Simplified: assume isothermal atmosphere for sound speed
    h = max(altitude, 0.0)
    a = _SPEED_OF_SOUND_SL * np.exp(-h / (2.0 * 42000.0))
    return speed / max(a, 1.0)

def hypersonic_aero_thermal_coupled_rhs(t, y):
    pos = y[0:3]
    vel = y[3:6]
    T = y[6:10].copy()
    T = np.maximum(T, 200.0)
    mass = max(y[12], 100.0)

    speed = np.linalg.norm(vel)
    speed = max(speed, 1e-6)
    v_hat = vel / speed

    altitude = np.linalg.norm(pos) - _R_EARTH
    rho_inf = _atmosphere_density(altitude)
    mach = _mach_number(speed, altitude)

    cd = np.interp(mach, _MACH_TABLE, _CD_TABLE)
    q_dyn = 0.5 * rho_inf * speed * speed
    drag = q_dyn * cd * _S_REF
    f_drag = -drag * v_hat

    g_accel = _gravity_accel(pos)

    if t < _T_BURN and mass > _M0 - _MDOT * _T_BURN:
        f_thrust = _THRUST * v_hat
        dm_dt = -_MDOT
    else:
        f_thrust = np.zeros(3)
        dm_dt = 0.0

    accel = (f_drag + f_thrust) / mass + g_accel

    # Sutton-Graves stagnation-point heating
    rho_inf_safe = max(rho_inf, 1e-20)
    q_stag = _SUTTON_GRAVES_COEFF * np.sqrt(rho_inf_safe / _R_NOSE) * speed**3

    # TPS thermal response: 4-layer lumped conduction
    # Surface BC: net heat flux = stagnation heating - radiative cooling
    q_rad = _TPS_EPSILON * _SIGMA_SB * T[0]**4
    q_net_surface = q_stag - q_rad

    # Inter-layer conduction fluxes (positive = toward substrate)
    q_01 = _TPS_KEFF_01 * (T[0] - T[1]) / _TPS_DX_01
    q_12 = _TPS_KEFF_12 * (T[1] - T[2]) / _TPS_DX_12
    q_23 = _TPS_KEFF_23 * (T[2] - T[3]) / _TPS_DX_23

    dTdt = np.empty(4)
    dTdt[0] = (q_net_surface - q_01) / _TPS_RHOC[0]
    dTdt[1] = (q_01 - q_12) / _TPS_RHOC[1]
    dTdt[2] = (q_12 - q_23) / _TPS_RHOC[2]
    dTdt[3] = q_23 / _TPS_RHOC[3]          # insulated back face

    # q_stag_norm tracks normalized stagnation heat flux (for diagnostics)
    q_peak_ref = 5e6
    dq_norm_dt = (q_stag / max(q_peak_ref, 1.0) - y[10]) / 10.0  # relaxation

    dy = np.empty(13)
    dy[0:3] = vel
    dy[3:6] = accel
    dy[6:10] = dTdt
    dy[10] = dq_norm_dt
    dy[11] = 0.0            # alpha_deg held constant
    dy[12] = dm_dt
    return dy