ballistics_engine/

atmosphere.rs

//! Enhanced atmospheric calculations for ballistics.
//!
//! This module provides Rust-accelerated implementations of atmospheric calculations
//! with full ICAO Standard Atmosphere support for improved accuracy at all altitudes.

use std::cmp::Ordering;

/// ICAO Standard Atmosphere layer definitions
#[derive(Debug, Clone)]
struct AtmosphereLayer {
    /// Base geopotential height of this layer (m)
    base_altitude: f64,
    /// Base temperature at layer start (K)
    base_temperature: f64,
    /// Base pressure at layer start (Pa)
    base_pressure: f64,
    /// Temperature lapse rate (K/m)
    lapse_rate: f64,
}

/// ICAO Standard Atmosphere constants
const G_ACCEL_MPS2: f64 = 9.80665;
const R_AIR: f64 = 287.0531; // Specific gas constant for dry air (J/(kg·K))
const GAMMA: f64 = 1.4; // Heat capacity ratio for air
const GEOPOTENTIAL_EARTH_RADIUS_M: f64 = 6_356_766.0;
const MIN_GEOMETRIC_ALTITUDE_M: f64 = -5000.0;
const MAX_GEOMETRIC_ALTITUDE_M: f64 = 84000.0;
const MIN_STANDARD_GEOPOTENTIAL_HEIGHT_M: f64 = -5000.0;
const MAX_STANDARD_GEOPOTENTIAL_HEIGHT_M: f64 = 84000.0;

/// CIPM constants for precise air density calculation
const R: f64 = 8.314472; // Universal gas constant
const M_A: f64 = 28.96546e-3; // Molar mass of dry air (kg/mol)
const M_V: f64 = 18.01528e-3; // Molar mass of water vapor (kg/mol)

/// ICAO Standard Atmosphere layer data up to 84 km
/// Pressures calculated using barometric formula between layers
const ICAO_LAYERS: &[AtmosphereLayer] = &[
    // Troposphere (-5 - 11 km; the sea-level base extrapolates smoothly below 0 m)
    AtmosphereLayer {
        base_altitude: 0.0,
        base_temperature: 288.15, // 15°C
        base_pressure: 101325.0,  // 1013.25 hPa
        lapse_rate: -0.0065,      // -6.5 K/km
    },
    // Tropopause (11 - 20 km)
    AtmosphereLayer {
        base_altitude: 11000.0,
        base_temperature: 216.65, // -56.5°C
        base_pressure: 22632.1,   // 226.32 hPa
        lapse_rate: 0.0,          // Isothermal
    },
    // Stratosphere 1 (20 - 32 km)
    AtmosphereLayer {
        base_altitude: 20000.0,
        base_temperature: 216.65, // -56.5°C
        base_pressure: 5474.89,   // 54.75 hPa
        lapse_rate: 0.001,        // +1 K/km
    },
    // Stratosphere 2 (32 - 47 km)
    AtmosphereLayer {
        base_altitude: 32000.0,
        base_temperature: 228.65, // -44.5°C
        base_pressure: 868.02,    // 8.68 hPa
        lapse_rate: 0.0028,       // +2.8 K/km
    },
    // Stratopause (47 - 51 km)
    AtmosphereLayer {
        base_altitude: 47000.0,
        base_temperature: 270.65, // -2.5°C
        base_pressure: 110.91,    // 1.11 hPa
        lapse_rate: 0.0,          // Isothermal
    },
    // Mesosphere 1 (51 - 71 km)
    AtmosphereLayer {
        base_altitude: 51000.0,
        base_temperature: 270.65, // -2.5°C
        base_pressure: 66.94,     // 0.67 hPa
        lapse_rate: -0.0028,      // -2.8 K/km
    },
    // Mesosphere 2 (71 - 84 km)
    AtmosphereLayer {
        base_altitude: 71000.0,
        base_temperature: 214.65, // -58.5°C
        base_pressure: 3.96,      // 0.04 hPa
        lapse_rate: -0.002,       // -2.0 K/km
    },
];

/// Calculate ICAO Standard Atmosphere conditions at any altitude.
///
/// This function implements the full ICAO Standard Atmosphere model with all
/// atmospheric layers up to 84 km altitude.
///
/// # Arguments
/// * `altitude_m` - Geometric altitude above mean sea level in meters (-5000 to 84000; values
///   outside are clamped). It is converted to geopotential height before layer evaluation.
///
/// # Returns
/// Tuple of (temperature_k, pressure_pa)
pub(crate) fn calculate_icao_standard_atmosphere(altitude_m: f64) -> (f64, f64) {
    let geometric_altitude = altitude_m.clamp(MIN_GEOMETRIC_ALTITUDE_M, MAX_GEOMETRIC_ALTITUDE_M);
    let geopotential_height = geometric_to_geopotential_height_m(geometric_altitude).clamp(
        MIN_STANDARD_GEOPOTENTIAL_HEIGHT_M,
        MAX_STANDARD_GEOPOTENTIAL_HEIGHT_M,
    );

    // Find the appropriate atmospheric layer
    let layer = ICAO_LAYERS
        .iter()
        .rev()
        .find(|layer| geopotential_height >= layer.base_altitude)
        .unwrap_or(&ICAO_LAYERS[0]);

    let height_diff = geopotential_height - layer.base_altitude;
    let temperature = layer.base_temperature + layer.lapse_rate * height_diff;

    let pressure = if layer.lapse_rate.abs() < 1e-10 {
        // Isothermal layer
        layer.base_pressure * (-G_ACCEL_MPS2 * height_diff / (R_AIR * layer.base_temperature)).exp()
    } else {
        // Non-isothermal layer
        let temp_ratio = temperature / layer.base_temperature;
        layer.base_pressure * temp_ratio.powf(-G_ACCEL_MPS2 / (layer.lapse_rate * R_AIR))
    };

    (temperature, pressure)
}

/// Convert physical geometric altitude to the geopotential height used by ICAO layer tables.
#[inline]
fn geometric_to_geopotential_height_m(geometric_altitude_m: f64) -> f64 {
    GEOPOTENTIAL_EARTH_RADIUS_M * geometric_altitude_m
        / (GEOPOTENTIAL_EARTH_RADIUS_M + geometric_altitude_m)
}

/// Resolve the station-pressure override for an air-density calculation.
///
/// Altitude and pressure are redundant inputs for density. The rule:
/// * An explicitly-supplied pressure is the authoritative STATION pressure (already
///   altitude-reduced); it is returned as `Some` and used directly, so altitude is NOT
///   double-counted.
/// * When pressure is left at the sea-level standard default (≈1013.25 hPa) while a real
///   altitude is given, the caller meant "standard atmosphere at this altitude": return
///   `None` so [`calculate_atmosphere`] derives the station pressure from altitude (ICAO
///   standard) instead of silently using sea-level density.
///
/// Without this, `--altitude` with the default pressure produced sea-level density (altitude
/// had no effect on drag). The ±0.5 hPa tolerance covers the `29.92 inHg ≈ 1013.21 hPa`
/// conversion, and `>1 m` avoids triggering at sea level. (Mirrors the existing
/// `pressure != 29.92` "user override" sentinel used elsewhere in the CLI.)
pub fn resolve_station_pressure(pressure_hpa: f64, altitude_m: f64) -> Option<f64> {
    const SEA_LEVEL_HPA: f64 = 1013.25;
    if (pressure_hpa - SEA_LEVEL_HPA).abs() < 0.5 && altitude_m.abs() > 1.0 {
        None // pressure left at default + real altitude → derive station pressure from altitude
    } else {
        Some(pressure_hpa) // explicit station pressure is authoritative
    }
}

/// Resolve the temperature override for an air-density calculation, mirroring
/// [`resolve_station_pressure`].
///
/// * An explicitly-supplied temperature is authoritative (returned as `Some`).
/// * When temperature is left at the sea-level standard default (15 °C) while a real altitude
///   is given, the caller meant "standard atmosphere at this altitude": return `None` so
///   [`calculate_atmosphere`] applies the ICAO lapse-rate temperature for that altitude
///   (≈ −6.5 °C/km).
///
/// Without this, `--altitude` with the default temperature held the air at 15 °C, which
/// under-estimates density (warm air is thinner) by ~2.4% at 1 km up to ~7% at 3 km versus the
/// standard atmosphere — validated against py_ballisticcalc, which derives both temperature and
/// pressure from altitude. The 0.1 °C tolerance matches the `59 °F = 15.0 °C` default exactly,
/// and `>1 m` avoids triggering at sea level. A shooter at a genuinely non-standard temperature
/// at altitude should pass an explicit temperature (same contract as station pressure).
pub fn resolve_station_temperature(temperature_c: f64, altitude_m: f64) -> Option<f64> {
    const SEA_LEVEL_TEMP_C: f64 = 15.0;
    if (temperature_c - SEA_LEVEL_TEMP_C).abs() < 0.1 && altitude_m.abs() > 1.0 {
        None // temperature left at default + real altitude → derive ICAO lapse temperature
    } else {
        Some(temperature_c) // explicit temperature is authoritative
    }
}

/// Return the station temperature and pressure that [`calculate_atmosphere`] will use after
/// applying the default-at-altitude resolution rules.
///
/// A thin, byte-identical delegate to
/// [`resolve_station_conditions_with_pressure_mode`] at [`PressureReferenceMode::Absolute`]
/// (MBA-1397) — kept as a separate name/signature so no existing caller needs to change.
pub fn resolve_station_conditions(
    temperature_c: f64,
    pressure_hpa: f64,
    altitude_m: f64,
) -> (f64, f64) {
    resolve_station_conditions_with_pressure_mode(
        temperature_c,
        pressure_hpa,
        altitude_m,
        PressureReferenceMode::Absolute,
    )
}

/// Whether a station-pressure input is already absolute (station) pressure, or a sea-level-
/// corrected altimeter setting (QNH / "barometer" reading) that must be reduced to station
/// pressure at the shooter's altitude before use (MBA-1397).
///
/// Kestrel, AB Mobile, Shooter, JBM, and Hornady 4DOF all let the user declare which one they
/// are entering. Previously this engine only supported [`Absolute`](Self::Absolute) — an
/// entered sea-level-corrected barometer/METAR value at a real altitude was silently treated
/// as station pressure, over-stating air density and under-stating drop.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
#[cfg_attr(feature = "cli", derive(clap::ValueEnum))]
pub enum PressureReferenceMode {
    /// The input is already absolute station pressure — today's meaning, and the default, so
    /// every existing caller that never sets this is byte-identical to pre-MBA-1397 behavior.
    #[default]
    Absolute,
    /// The input is a sea-level-corrected altimeter setting (QNH) and must be reduced to
    /// station pressure at the target altitude via [`reduce_qnh_to_station_pressure`].
    Qnh,
}

/// Reduce a sea-level-corrected altimeter setting (QNH) to the station pressure at
/// `altitude_m` (MBA-1397), using the inverse of the ICAO troposphere barometric formula:
///
/// `station = QNH * (1 + L*h/T0)^(-g/(L*R))`
///
/// which, for the troposphere's published lapse rate `L` and sea-level temperature `T0`, is
/// exactly the standard ICAO Doc 7488 / NOAA altimeter-setting reduction
/// `station = QNH * (1 - 0.0065*h/288.15)^5.25588`, where `h` is GEOPOTENTIAL height, not
/// the geometric altitude the caller supplies — this function converts the one to the other
/// first (`geometric_to_geopotential_height_m`), as ISA requires. Substituting geometric
/// altitude directly into that expression reproduces the reduction only to about 0.04 hPa
/// at 1500 m, which is why the worked values below are not what the bare formula gives.
///
/// Reuses the SAME troposphere layer (`ICAO_LAYERS[0]`) and physical constants
/// (`G_ACCEL_MPS2`, `R_AIR`) that `calculate_icao_standard_atmosphere` (private to this crate)
/// uses for its own non-isothermal-layer pressure formula, rather than re-hardcoding the
/// lapse rate / sea-level temperature / exponent a second time — this is that same formula,
/// solved for a caller-given sea-level pressure instead of the fixed 1013.25 hPa standard.
/// `altitude_m` is clamped and converted to geopotential height exactly like
/// `calculate_icao_standard_atmosphere`, so a QNH of exactly 1013.25 hPa reduces to precisely
/// the ICAO standard station pressure at every altitude (this is what makes the
/// omitted-pressure default byte-identical under either [`PressureReferenceMode`]).
///
/// The reduction is only physically meaningful within the troposphere (the altitude range
/// real shooting takes place in); it is not re-derived per-layer for higher altitudes.
pub fn reduce_qnh_to_station_pressure(qnh_hpa: f64, altitude_m: f64) -> f64 {
    let geometric_altitude = altitude_m.clamp(MIN_GEOMETRIC_ALTITUDE_M, MAX_GEOMETRIC_ALTITUDE_M);
    let geopotential_height = geometric_to_geopotential_height_m(geometric_altitude).clamp(
        MIN_STANDARD_GEOPOTENTIAL_HEIGHT_M,
        MAX_STANDARD_GEOPOTENTIAL_HEIGHT_M,
    );
    let layer = &ICAO_LAYERS[0]; // troposphere: the only layer the altimeter-setting reduction covers
    let temp_ratio = 1.0 + layer.lapse_rate * geopotential_height / layer.base_temperature;
    let exponent = -G_ACCEL_MPS2 / (layer.lapse_rate * R_AIR);
    qnh_hpa * temp_ratio.powf(exponent)
}

/// Mode-aware counterpart of [`resolve_station_pressure`] (MBA-1397). Does NOT change that
/// function's signature or behavior — [`PressureReferenceMode::Absolute`] is a pure
/// passthrough to it — so this is a new function alongside the original, not a replacement.
///
/// [`PressureReferenceMode::Qnh`] reduces `pressure_hpa` (interpreted as a QNH) to station
/// pressure via [`reduce_qnh_to_station_pressure`] and returns it as `Some(..)`
/// UNCONDITIONALLY, bypassing the default-sentinel ambiguity `resolve_station_pressure` must
/// use to infer omission from a bare `f64`. This matters: a QNH the caller explicitly declared
/// should never be silently re-interpreted as an omitted default merely because the REDUCED
/// number happens to land within the sentinel's tolerance of the 1013.25 hPa sea-level
/// constant (a real, reachable case — e.g. a high-pressure day's QNH reducing to ~1013 hPa at
/// a few hundred meters of elevation).
pub fn resolve_station_pressure_with_mode(
    pressure_hpa: f64,
    altitude_m: f64,
    mode: PressureReferenceMode,
) -> Option<f64> {
    match mode {
        PressureReferenceMode::Absolute => resolve_station_pressure(pressure_hpa, altitude_m),
        PressureReferenceMode::Qnh => {
            Some(reduce_qnh_to_station_pressure(pressure_hpa, altitude_m))
        }
    }
}

/// Mode-aware counterpart of [`resolve_station_conditions`] (MBA-1397): identical for
/// [`PressureReferenceMode::Absolute`] (byte-identical delegate — [`resolve_station_conditions`]
/// itself now forwards here), and reduces an explicit QNH pressure to station pressure via
/// [`resolve_station_pressure_with_mode`] for [`PressureReferenceMode::Qnh`] before applying
/// the same standard-atmosphere fallback used when either input is left at its sea-level
/// default.
pub fn resolve_station_conditions_with_pressure_mode(
    temperature_c: f64,
    pressure_hpa: f64,
    altitude_m: f64,
    pressure_mode: PressureReferenceMode,
) -> (f64, f64) {
    let temp_override = resolve_station_temperature(temperature_c, altitude_m);
    let press_override =
        resolve_station_pressure_with_mode(pressure_hpa, altitude_m, pressure_mode);
    let (std_temp_k, std_pressure_pa) = calculate_icao_standard_atmosphere(altitude_m);
    let temp_c = temp_override.unwrap_or(std_temp_k - 273.15);
    let pressure_hpa = press_override.unwrap_or(std_pressure_pa / 100.0);
    (temp_c, pressure_hpa)
}

// ---------------------------------------------------------------------------------------
// MBA-1366: density altitude as a direct atmosphere input
// ---------------------------------------------------------------------------------------

/// NWS/FAA published pressure-altitude constants. Duplicated here (not imported) because the
/// forward density-altitude formula this inverts,
/// `pdf_dope_card::calculate_density_altitude(_altitude_ft, pressure_inhg, temp_f)`, lives in
/// the `ballistics` BINARY crate's CLI-only `pdf_dope_card` module (`src/pdf_dope_card.rs`) —
/// this LIBRARY crate has no dependency on it and cannot call it. `src/main.rs`'s
/// `density_altitude_round_trips_through_the_dope_card_formula` test is the cross-crate proof
/// that the two stay numerically consistent; these constants must match that function's
/// literals exactly (145_366.45, 0.190_284, 3.57, 59.0, and the "120 ft/degC" correction) for
/// that round trip to hold.
const DA_NWS_SEA_LEVEL_HPA: f64 = 1013.25;
const DA_NWS_PRESSURE_ALT_FT: f64 = 145_366.45;
const DA_NWS_EXPONENT: f64 = 0.190_284;
const DA_NWS_SEA_LEVEL_TEMP_F: f64 = 59.0;
const DA_NWS_LAPSE_F_PER_1000FT: f64 = 3.57;
const DA_FT_TO_M: f64 = 0.3048;

/// Back-solve an ISA-equivalent station altitude/temperature/pressure from a directly-declared
/// density altitude (MBA-1366) — the single-value atmosphere entry mode Shooter, Nosler, AB
/// Analytics, Ballistic AE, and TRASOL all support as an alternative to altitude + pressure +
/// temperature.
///
/// This is the exact algebraic inverse of the published NWS/FAA density-altitude model this
/// engine already uses to REPORT density altitude
/// (`pdf_dope_card::calculate_density_altitude`):
///
/// ```text
/// pressure_alt_ft = 145366.45 * (1 - (station_hpa / 1013.25)^0.190284)
/// isa_temp_f      = 59.0 - pressure_alt_ft / 1000.0 * 3.57
/// density_alt_ft  = pressure_alt_ft + (120*5/9) * (station_temp_f - isa_temp_f)
/// ```
///
/// Solved for `station_hpa` given a target `density_alt_ft` and a station temperature `K = 120
/// * 5/9` is the "120 ft per °C" density-altitude correction rule, expressed in ft per °F):
///
/// ```text
/// pressure_alt_ft = (density_alt_ft - K*(station_temp_f - 59.0)) / (1 + K*3.57/1000.0)
/// station_hpa     = 1013.25 * (1 - pressure_alt_ft/145366.45)^(1/0.190284)
/// ```
///
/// # Temperature precedence (MBA-1366)
/// * `explicit_temperature_c == None`: ISA-at-density-altitude is the default. Algebraically
///   this collapses `pressure_alt_ft` to EXACTLY `density_alt_ft`: the correction term above is
///   zero by construction whenever the station temperature equals the ISA temperature at that
///   pressure altitude, which is self-consistently true when temperature is left at ISA — i.e.
///   omitting an explicit temperature is equivalent to typing `--altitude <density_altitude>`
///   with no `--temperature`/`--pressure` override.
/// * `explicit_temperature_c == Some(t)`: `t` wins outright (returned unchanged as
///   `temperature_c`, never re-derived) and the station pressure is re-solved so the resulting
///   density altitude still reproduces the input exactly — density is honored either way; only
///   the implied pressure (and therefore station altitude) differs from the ISA-default case.
///
/// # Height convention
/// The NWS/FAA formula being inverted performs no geopotential correction — it treats its
/// altitude as a plain height, exactly like `pdf_dope_card::calculate_density_altitude` does
/// (whose own `_altitude_ft` parameter is unused and undocumented as geometric-vs-geopotential
/// for that reason). Consistent with that, the `altitude_m` this returns is GEOMETRIC altitude
/// — the same convention every other engine altitude input uses (`--altitude`,
/// [`AtmosphericConditions::altitude`](crate::cli_api::AtmosphericConditions::altitude)) —
/// rather than being run through `geometric_to_geopotential_height_m` a second time here; that
/// conversion is `calculate_icao_standard_atmosphere`'s job (both private to this crate — see
/// their doc comments elsewhere in this file), applied once this altitude re-enters the normal
/// altitude-lapse pipeline (this is what "back-solve an ISA-equivalent atmosphere" means —
/// unlike [`get_direct_atmosphere`], which would freeze speed of sound and bypass that pipeline
/// entirely).
///
/// # Interaction with [`PressureReferenceMode`] (QNH)
/// Density altitude supersedes pressure outright — callers must not also run a declared
/// `--pressure`/QNH value through [`resolve_station_pressure_with_mode`] when a density
/// altitude is present; see the CLI/WASM call sites, which emit a note to that effect.
///
/// # Returns
/// `(altitude_m, temperature_c, pressure_hpa)`, meant to be written directly into
/// [`crate::cli_api::AtmosphericConditions`] (and the paired `BallisticInputs` environment
/// fields, so powder-temperature sensitivity and the moist speed of sound both see the real
/// resolved temperature) via the Authoritative
/// (`TrajectorySolver::new_with_resolved_station_atmosphere`) constructor.
pub fn resolve_atmosphere_for_density_altitude(
    density_altitude_m: f64,
    explicit_temperature_c: Option<f64>,
) -> (f64, f64, f64) {
    let density_altitude_ft = density_altitude_m / DA_FT_TO_M;
    const K: f64 = 120.0 * 5.0 / 9.0; // ft of density-altitude correction per °F station-vs-ISA delta

    let (pressure_alt_ft, temperature_c) = match explicit_temperature_c {
        Some(temp_c) => {
            let temp_f = temp_c * 9.0 / 5.0 + 32.0;
            let denom = 1.0 + K * DA_NWS_LAPSE_F_PER_1000FT / 1000.0;
            let pressure_alt_ft =
                (density_altitude_ft - K * (temp_f - DA_NWS_SEA_LEVEL_TEMP_F)) / denom;
            (pressure_alt_ft, temp_c)
        }
        None => {
            // ISA-at-DA default: the correction term vanishes by construction (see doc
            // comment above), so the pressure altitude IS the density altitude exactly.
            let pressure_alt_ft = density_altitude_ft;
            let isa_temp_f =
                DA_NWS_SEA_LEVEL_TEMP_F - pressure_alt_ft / 1000.0 * DA_NWS_LAPSE_F_PER_1000FT;
            let isa_temp_c = (isa_temp_f - 32.0) * 5.0 / 9.0;
            (pressure_alt_ft, isa_temp_c)
        }
    };

    let pressure_hpa = DA_NWS_SEA_LEVEL_HPA
        * (1.0 - pressure_alt_ft / DA_NWS_PRESSURE_ALT_FT).powf(1.0 / DA_NWS_EXPONENT);
    let altitude_m = pressure_alt_ft * DA_FT_TO_M;

    (altitude_m, temperature_c, pressure_hpa)
}

/// Enhanced atmospheric calculation with ICAO Standard Atmosphere.
///
/// # Arguments
/// * `altitude_m` - Altitude in meters
/// * `temp_override_c` - Temperature override in Celsius (None for standard)
/// * `press_override_hpa` - Pressure override in hPa (None for standard)
/// * `humidity_percent` - Humidity percentage (0-100)
///
/// # Returns
/// Tuple of (air_density_kg_m3, speed_of_sound_mps)
pub fn calculate_atmosphere(
    altitude_m: f64,
    temp_override_c: Option<f64>,
    press_override_hpa: Option<f64>,
    humidity_percent: f64,
) -> (f64, f64) {
    // Get standard atmosphere conditions or use overrides
    let combined_overrides = temp_override_c.zip(press_override_hpa);
    let (temp_k, pressure_pa) = if let Some((temp_c, pressure_hpa)) = combined_overrides {
        // Both overrides provided
        (temp_c + 273.15, pressure_hpa * 100.0)
    } else {
        // Get ICAO standard conditions
        let (std_temp_k, std_pressure_pa) = calculate_icao_standard_atmosphere(altitude_m);

        let final_temp_k = if let Some(temp_c) = temp_override_c {
            temp_c + 273.15
        } else {
            std_temp_k
        };

        let final_pressure_pa = if let Some(press_hpa) = press_override_hpa {
            press_hpa * 100.0
        } else {
            std_pressure_pa
        };

        (final_temp_k, final_pressure_pa)
    };

    // Humidity clamp shared by the CIPM density and the moist speed of sound.
    let humidity_clamped = humidity_percent.clamp(0.0, 100.0);
    let temp_c = temp_k - 273.15;

    // Density: CIPM-2007 is the single canonical humid-air density model. Every solver
    // (cli_api / monte_carlo / ffi / fast_trajectory) reaches this one formula through
    // calculate_atmosphere, so there is no second (Arden-Buck ideal-gas) density path to drift.
    let density = calculate_air_density_cimp(temp_c, pressure_pa / 100.0, humidity_clamped);

    // Speed of sound uses an ideal-gas moist-air mixture approximation, first order in the
    // water-vapor mole fraction. Extracted into `moist_speed_of_sound` so the integrators can
    // share it; its vapor pressure comes from the SAME IAPWS saturation formula
    // (`enhanced_saturation_vapor_pressure`) + CIPM enhancement factor that the density above
    // uses. Only the vapor fraction is shared; the acoustic relation remains ideal-gas.
    let speed_of_sound = moist_speed_of_sound(temp_k, pressure_pa, humidity_clamped);

    (density, speed_of_sound)
}

/// Speed of sound using an ideal-gas moist-air mixture approximation, first order in the
/// water-vapor mole fraction.
///
/// The water-vapor mole fraction is derived from the SAME IAPWS saturation vapor pressure
/// (`enhanced_saturation_vapor_pressure`) and CIPM-2007 enhancement factor used by
/// [`calculate_air_density_cimp`], so a single vapor formula feeds both density and c. Only the
/// vapor fraction is shared with the CIPM path; the acoustic relation itself remains ideal-gas.
/// The mixture applies first-order humidity corrections to the dry-air heat-capacity ratio and
/// gas constant, then evaluates `sqrt(gamma * R * T)`. It is not Cramer's full real-gas
/// polynomial and does not include Cramer's pressure or carbon-dioxide terms.
///
/// # Arguments
/// * `temp_k` - Temperature in Kelvin
/// * `pressure_pa` - Total (station) pressure in Pa
/// * `humidity_percent` - Relative humidity percentage (0-100)
///
/// # Returns
/// Speed of sound in m/s
pub fn moist_speed_of_sound(temp_k: f64, pressure_pa: f64, humidity_percent: f64) -> f64 {
    let humidity_clamped = humidity_percent.clamp(0.0, 100.0);
    let temp_c = temp_k - 273.15;

    // Water-vapor partial pressure p_v = RH * f * p_sv, matching CIPM's x_v exactly. p_sv is in
    // hPa (enhanced_saturation_vapor_pressure returns hPa), so convert to Pa before forming the
    // mole fraction against the Pa total pressure.
    let p_sv_hpa = enhanced_saturation_vapor_pressure(temp_k);
    let f = enhanced_enhancement_factor(pressure_pa, temp_c);
    let vapor_pressure_pa = humidity_clamped / 100.0 * f * p_sv_hpa * 100.0;

    // Cap the mole fraction at the physical maximum of 1 and guard pressure_pa == 0 (a 0 hPa
    // override would otherwise give +Inf -> NaN speed of sound).
    let mole_fraction_vapor = (vapor_pressure_pa / pressure_pa.max(f64::MIN_POSITIVE)).min(1.0);

    // Heat-capacity ratio and gas constant for moist air (mole-fraction coefficients). 0.378 is
    // the dry-air molecular-weight ratio (0.6078 would belong to specific humidity, not mole
    // fraction).
    let gamma_moist = GAMMA * (1.0 - mole_fraction_vapor * 0.062);
    let r_moist = R_AIR * (1.0 + 0.378 * mole_fraction_vapor);

    (gamma_moist * r_moist * temp_k).sqrt()
}

/// Enhanced air density calculation using CIPM formula with ICAO atmosphere.
///
/// # Arguments
/// * `temp_c` - Temperature in Celsius
/// * `pressure_hpa` - Pressure in hPa
/// * `humidity_percent` - Humidity percentage (0-100)
///
/// # Returns
/// Air density in kg/m³
pub fn calculate_air_density_cimp(temp_c: f64, pressure_hpa: f64, humidity_percent: f64) -> f64 {
    let t_k = temp_c + 273.15;

    // Enhanced saturation vapor pressure calculation
    let p_sv = enhanced_saturation_vapor_pressure(t_k);

    let pressure_pa = pressure_hpa * 100.0;

    // Enhanced enhancement factor with temperature dependence. CIPM constants use Pa.
    let f = enhanced_enhancement_factor(pressure_pa, temp_c);

    // Vapor pressure with clamping. p_sv is in hPa (enhanced_saturation_vapor_pressure
    // returns hPa — its critical-pressure constant is 220640 hPa), so p_v is in hPa too.
    let p_v = humidity_percent.clamp(0.0, 100.0) / 100.0 * f * p_sv;

    // Convert the vapor pressure to Pa BEFORE forming the mole fraction: the divisor below
    // is in Pa. Dividing the hPa p_v by the Pa total made x_v 100x too small, which erased
    // the humidity term and returned essentially dry-air density (e.g. 15 C / 1013.25 hPa /
    // 50% RH gave ~1.2254 instead of the CIPM-2007 moist value ~1.2211 — moist air is
    // LIGHTER than dry air).
    let p_v_pa = p_v * 100.0;

    // Floor the pressure divisor (mirrors calculate_atmosphere): a 0 hPa pressure would
    // otherwise make x_v = +Inf -> NaN density. No-op for all valid (>0) pressures.
    let p_pa = pressure_pa.max(f64::MIN_POSITIVE);

    // Mole fraction of water vapor (capped at the physical maximum of 1)
    let x_v = (p_v_pa / p_pa).min(1.0);

    // Enhanced compressibility factor. CIPM virial constants use Pa.
    let z = enhanced_compressibility_factor(p_pa, t_k, x_v);

    // Calculate density with enhanced precision
    // Note: parentheses are important here for correct operator precedence
    ((p_pa * M_A) / (z * R * t_k)) * (1.0 - x_v * (1.0 - M_V / M_A))
}

/// Enhanced saturation vapor pressure calculation.
/// Uses the IAPWS-IF97 formulation for high precision.
#[inline(always)]
fn enhanced_saturation_vapor_pressure(t_k: f64) -> f64 {
    // IAPWS-IF97 coefficients for better accuracy
    const A: [f64; 6] = [
        -7.85951783,
        1.84408259,
        -11.7866497,
        22.6807411,
        -15.9618719,
        1.80122502,
    ];

    // Ensure temperature is positive and reasonable
    let t_k_safe = t_k.max(173.15); // -100°C minimum

    let tau = 1.0 - t_k_safe / 647.096; // Critical temperature of water
    let ln_p_ratio = (647.096 / t_k_safe)
        * (A[0] * tau
            + A[1] * tau.powf(1.5)
            + A[2] * tau.powf(3.0)
            + A[3] * tau.powf(3.5)
            + A[4] * tau.powf(4.0)
            + A[5] * tau.powf(7.5));

    220640.0 * ln_p_ratio.exp() // Critical pressure in hPa (22.064 MPa)
}

/// CIPM-2007 enhancement factor `f = alpha + beta*p + gamma*t^2` (p in Pa, t in Celsius).
#[inline(always)]
fn enhanced_enhancement_factor(p: f64, t: f64) -> f64 {
    const ALPHA: f64 = 1.00062;
    const BETA: f64 = 3.14e-8;
    const GAMMA: f64 = 5.6e-7;

    ALPHA + BETA * p + GAMMA * t * t
}

/// CIPM-2007 compressibility factor `Z` (virial expansion, second order in `p/T`).
#[inline(always)]
fn enhanced_compressibility_factor(p: f64, t_k: f64, x_v: f64) -> f64 {
    // CIPM-2007 molar virial coefficients (p in Pa, t in Celsius).
    const A0: f64 = 1.58123e-6;
    const A1: f64 = -2.9331e-8;
    const A2: f64 = 1.1043e-10;
    const B0: f64 = 5.707e-6;
    const B1: f64 = -2.051e-8;
    const C0: f64 = 1.9898e-4;
    const C1: f64 = -2.376e-6;
    const D: f64 = 1.83e-11;
    const E: f64 = -0.765e-8;

    // Ensure temperature is positive
    let t_k_safe = t_k.max(173.15); // -100°C minimum
    let t = t_k_safe - 273.15;
    let p_t = p / t_k_safe;

    let z_second_order =
        1.0 - p_t * (A0 + A1 * t + A2 * t * t + (B0 + B1 * t) * x_v + (C0 + C1 * t) * x_v * x_v);

    let z_third_order = p_t * p_t * (D + E * x_v * x_v);

    z_second_order + z_third_order
}

/// Convert an `(x, y)` position in the shot-aligned frame to true world altitude.
///
/// The engine rotates gravity by `shooting_angle_rad`, so shot-frame X follows the inclined
/// line of fire and Y is perpendicular to it in the vertical plane. Atmosphere lookup needs the
/// world-vertical projection of that position, added to the station altitude.
#[inline]
pub(crate) fn shot_frame_altitude(
    base_altitude_m: f64,
    downrange_m: f64,
    shot_y_m: f64,
    shooting_angle_rad: f64,
) -> f64 {
    base_altitude_m
        + downrange_m * shooting_angle_rad.sin()
        + shot_y_m * shooting_angle_rad.cos()
}

/// Enhanced local atmospheric calculation with variable lapse rates.
///
/// # Arguments
/// * `altitude_m` - Query geometric altitude above mean sea level in meters
/// * `base_alt` - Base geometric altitude above mean sea level in meters
/// * `base_temp_c` - Base temperature in Celsius
/// * `base_press_hpa` - Base pressure in hPa
/// * `base_ratio` - Base density ratio
///
/// # Returns
/// Tuple of (air_density_kg_m3, speed_of_sound_mps)
pub fn get_local_atmosphere(
    altitude_m: f64,
    base_alt: f64,
    base_temp_c: f64,
    base_press_hpa: f64,
    base_ratio: f64,
) -> (f64, f64) {
    let (temp_k, _pressure_pa, density) =
        local_temp_pressure_density(altitude_m, base_alt, base_temp_c, base_press_hpa, base_ratio);

    // Dry speed of sound. 401.874 ~ gamma * R_air; kept exactly for back-compat with existing
    // callers (get_local_atmosphere_humid uses the precise moist formula instead).
    let speed_of_sound = (temp_k * 401.874).sqrt();

    (density, speed_of_sound)
}

/// Humidity-aware companion to [`get_local_atmosphere`]: identical local density, but the speed
/// of sound is the moist-air value ([`moist_speed_of_sound`]) evaluated at the LOCAL temperature
/// and pressure.
///
/// [`get_local_atmosphere`] is intentionally left unchanged (dry speed of sound) for
/// API/back-compat; call this variant only where a real relative humidity is available.
///
/// # Arguments
/// * `altitude_m` - Query geometric altitude above mean sea level in meters
/// * `base_alt` - Base (station) geometric altitude above mean sea level in meters
/// * `base_temp_c` - Base temperature in Celsius
/// * `base_press_hpa` - Base pressure in hPa
/// * `base_ratio` - Base density ratio (density / 1.225)
/// * `humidity_percent` - Relative humidity percentage (0-100)
///
/// # Returns
/// Tuple of (air_density_kg_m3, moist_speed_of_sound_mps)
pub fn get_local_atmosphere_humid(
    altitude_m: f64,
    base_alt: f64,
    base_temp_c: f64,
    base_press_hpa: f64,
    base_ratio: f64,
    humidity_percent: f64,
) -> (f64, f64) {
    let (temp_k, pressure_pa, density) =
        local_temp_pressure_density(altitude_m, base_alt, base_temp_c, base_press_hpa, base_ratio);
    (density, moist_speed_of_sound(temp_k, pressure_pa, humidity_percent))
}

/// Shared local temperature / pressure / density computation for [`get_local_atmosphere`] and
/// [`get_local_atmosphere_humid`]. Returns `(local_temp_k, local_pressure_pa, density_kg_m3)`.
#[inline]
fn local_temp_pressure_density(
    altitude_m: f64,
    base_alt: f64,
    base_temp_c: f64,
    base_press_hpa: f64,
    base_ratio: f64,
) -> (f64, f64, f64) {
    let base_temp_k = base_temp_c + 273.15;

    // A non-finite endpoint would make the boundary walk fail to advance. The
    // previous single-column formula also produced non-finite outputs here.
    if !altitude_m.is_finite() || !base_alt.is_finite() {
        return (f64::NAN, f64::NAN, f64::NAN);
    }

    let base_geopotential_m = geometric_to_geopotential_height_m(base_alt);
    let target_geopotential_m = geometric_to_geopotential_height_m(altitude_m);
    let (temp_k, pressure_hpa) = integrate_local_atmosphere_layers(
        base_geopotential_m,
        target_geopotential_m,
        base_temp_k,
        base_press_hpa,
    );

    // Enhanced density calculation
    let density_ratio = base_ratio * (base_temp_k * pressure_hpa) / (base_press_hpa * temp_k);
    let density = density_ratio * 1.225;

    (temp_k, pressure_hpa * 100.0, density)
}

/// Carry arbitrary station temperature and pressure through each crossed ICAO
/// layer in geopotential-height coordinates. The layer table supplies lapse rates and boundaries
/// only; station deviations from the standard atmosphere remain anchored at the converted base.
fn integrate_local_atmosphere_layers(
    base_geopotential_m: f64,
    target_geopotential_m: f64,
    mut temp_k: f64,
    mut pressure_hpa: f64,
) -> (f64, f64) {
    let mut current_alt = base_geopotential_m;

    if target_geopotential_m > current_alt {
        while current_alt < target_geopotential_m {
            // At an exact boundary, ascent starts in the higher layer.
            let layer_index = ICAO_LAYERS
                .iter()
                .rposition(|layer| current_alt >= layer.base_altitude)
                .unwrap_or(0);
            let segment_end = ICAO_LAYERS
                .get(layer_index + 1)
                .map_or(target_geopotential_m, |next| {
                    target_geopotential_m.min(next.base_altitude)
                });
            (temp_k, pressure_hpa) = integrate_local_atmosphere_segment(
                temp_k,
                pressure_hpa,
                segment_end - current_alt,
                ICAO_LAYERS[layer_index].lapse_rate,
            );
            current_alt = segment_end;
        }
    } else {
        while current_alt > target_geopotential_m {
            // At an exact boundary, descent starts in the lower layer. The
            // strict comparison is what makes a cross-layer round trip reversible.
            let layer_index = ICAO_LAYERS
                .iter()
                .rposition(|layer| current_alt > layer.base_altitude)
                .unwrap_or(0);
            let segment_end = if layer_index == 0 {
                target_geopotential_m
            } else {
                target_geopotential_m.max(ICAO_LAYERS[layer_index].base_altitude)
            };
            (temp_k, pressure_hpa) = integrate_local_atmosphere_segment(
                temp_k,
                pressure_hpa,
                segment_end - current_alt,
                ICAO_LAYERS[layer_index].lapse_rate,
            );
            current_alt = segment_end;
        }
    }

    (temp_k, pressure_hpa)
}

#[inline]
fn integrate_local_atmosphere_segment(
    base_temp_k: f64,
    base_pressure_hpa: f64,
    height_diff: f64,
    lapse_rate: f64,
) -> (f64, f64) {
    let temp_k = base_temp_k + lapse_rate * height_diff;
    let pressure_hpa = if lapse_rate.abs() < 1e-10 {
        base_pressure_hpa * (-G_ACCEL_MPS2 * height_diff / (R_AIR * base_temp_k)).exp()
    } else {
        let temp_ratio = temp_k / base_temp_k;
        base_pressure_hpa * temp_ratio.powf(-G_ACCEL_MPS2 / (lapse_rate * R_AIR))
    };

    (temp_k, pressure_hpa)
}

/// Determine local lapse rate based on altitude and atmospheric layer.
#[cfg(test)]
#[inline(always)]
fn determine_local_lapse_rate(altitude_m: f64) -> f64 {
    // Find the current atmospheric layer to get appropriate lapse rate
    let layer = ICAO_LAYERS
        .iter()
        .rev()
        .find(|layer| altitude_m >= layer.base_altitude)
        .unwrap_or(&ICAO_LAYERS[0]);

    layer.lapse_rate
}

/// Direct atmosphere calculation for simple cases.
///
/// # Arguments
/// * `density` - Pre-computed air density
/// * `speed_of_sound` - Pre-computed speed of sound
///
/// # Returns
/// Tuple of (air_density, speed_of_sound) - just passes through the values
#[inline(always)]
pub fn get_direct_atmosphere(density: f64, speed_of_sound: f64) -> (f64, f64) {
    (density, speed_of_sound)
}

/// Legacy function name for backwards compatibility
pub fn calculate_air_density_cipm(temp_c: f64, pressure_hpa: f64, humidity_percent: f64) -> f64 {
    calculate_air_density_cimp(temp_c, pressure_hpa, humidity_percent)
}

/// A single downrange-referenced atmosphere zone:
/// `(temp_c, pressure_hpa, humidity_percent, until_distance_m)`.
///
/// The T/P/H are the STATION-REFERENCED conditions (defined at the shooter base altitude) that
/// apply from the previous segment's threshold out to `until_distance_m`. This mirrors
/// [`crate::wind::WindSegment`]'s `(speed, angle, until_distance)` shape so the two segmented
/// models compose the same way (wind by X, atmosphere by X, altitude lapse by Y).
pub type AtmoSegment = (f64, f64, f64, f64);

/// Downrange-segmented atmosphere handler (MBA-1137), the density analogue of
/// [`crate::wind::WindSock`].
///
/// Holds a set of station-referenced atmosphere zones ordered by their `until_distance_m`
/// threshold and answers a stateless downrange lookup ([`AtmoSock::atmo_for_range`]).
///
/// The zone T/P/H are the base (shooter-altitude) conditions for that stretch of range; the
/// solver swaps them into the same local-atmosphere altitude-lapse pipeline that a single-station
/// solve uses, so the downrange (X) zone and the vertical (Y) altitude lapse compose orthogonally
/// without double-counting (the zone sets the base tuple, the lapse multiplies on top of it).
#[derive(Debug, Clone)]
pub struct AtmoSock {
    /// Zones sorted ascending by `until_distance_m` (segment slot 3).
    segments: Vec<AtmoSegment>,
}

impl AtmoSock {
    /// Create a new `AtmoSock` from station-referenced atmosphere zones.
    ///
    /// Each segment is `(temp_c, pressure_hpa, humidity_percent, until_distance_m)`. Segments are
    /// sorted by `until_distance_m`, with NaN thresholds ordered last.
    pub fn new(mut segments: Vec<AtmoSegment>) -> Self {
        segments.sort_by(|a, b| match (a.3.is_nan(), b.3.is_nan()) {
            (true, true) => Ordering::Equal,
            (true, false) => Ordering::Greater,
            (false, true) => Ordering::Less,
            (false, false) => a.3.partial_cmp(&b.3).unwrap(),
        });
        AtmoSock { segments }
    }

    /// True when this sock carries no zones (a lookup falls back to sea-level ISA).
    pub fn is_empty(&self) -> bool {
        self.segments.is_empty()
    }

    /// Stateless downrange lookup of the active zone's `(temp_c, pressure_hpa, humidity_percent)`.
    ///
    /// Selection matches [`crate::wind::WindSock::vector_for_range_stateless`]: the first segment
    /// whose `until_distance_m` STRICTLY exceeds `downrange_m` wins (thresholds are upper-exclusive).
    /// Unlike wind — which returns zero past the last threshold — the LAST zone is used for any
    /// distance at or beyond the final threshold (there is no "zero atmosphere"). An empty sock
    /// returns the sea-level ISA reference `(15 C, 1013.25 hPa, 0% RH)`.
    ///
    /// This is stateless and safe for numerical integration (the same X may be queried repeatedly
    /// or out of order across RK substeps).
    pub fn atmo_for_range(&self, downrange_m: f64) -> (f64, f64, f64) {
        if self.segments.is_empty() {
            return (15.0, 1013.25, 0.0); // sea-level ISA fallback
        }
        // NaN X can't be ordered; use the first (nearest) zone deterministically.
        if downrange_m.is_nan() {
            let s = self.segments[0];
            return (s.0, s.1, s.2);
        }
        for seg in &self.segments {
            if downrange_m < seg.3 {
                return (seg.0, seg.1, seg.2);
            }
        }
        // Beyond the final threshold: hold the last zone (no zeroing).
        let last = self.segments[self.segments.len() - 1];
        (last.0, last.1, last.2)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn inclined_shot_frame_position_maps_to_world_altitude() {
        let base_altitude = 100.0;
        let downrange = 1_000.0;
        let shot_y = 10.0;
        let angle = std::f64::consts::FRAC_PI_6;
        let expected = base_altitude + downrange * angle.sin() + shot_y * angle.cos();

        let actual = shot_frame_altitude(base_altitude, downrange, shot_y, angle);
        assert!(
            (actual - expected).abs() < 1e-12,
            "30-degree shot at x=1000/y=10 should be at {expected} m, got {actual} m"
        );
        assert_eq!(
            shot_frame_altitude(base_altitude, downrange, shot_y, 0.0),
            base_altitude + shot_y,
            "flat-fire altitude must remain byte-identical"
        );
        let downhill = shot_frame_altitude(base_altitude, downrange, shot_y, -angle);
        let expected_downhill = base_altitude - downrange * angle.sin() + shot_y * angle.cos();
        assert!((downhill - expected_downhill).abs() < 1e-12);
    }

    // ---- MBA-1136: CIPM-2007 as the single canonical humid-air density ----

    /// Gate 1: dry sea-level (15 C, 1013.25 hPa, 0% RH) must stay at the ISA reference — density
    /// 1.225 +- 0.002 kg/m^3 and speed of sound 340.3 +- 0.6 m/s. CIPM-2007 at 0% RH reduces to
    /// dry-air ideal gas to within rounding, so this is essentially unchanged from the pre-CIPM
    /// baseline (baseline was 1.225012 / 340.294; now 1.225521 / 340.294 — the tiny density bump
    /// is CIPM compressibility + exact molar mass, the speed of sound is bit-identical).
    #[test]
    fn test_mba1136_dry_sea_level_reference() {
        let (density, sos) = calculate_atmosphere(0.0, Some(15.0), Some(1013.25), 0.0);
        assert!(
            (density - 1.225).abs() < 0.002,
            "dry sea-level density {density} not within 1.225 +- 0.002"
        );
        assert!(
            (sos - 340.3).abs() < 0.6,
            "dry sea-level speed of sound {sos} not within 340.3 +- 0.6"
        );
    }

    /// Gate 2: humid air (15 C, 1013.25 hPa, 50% RH) is the CIPM-2007 value (~1.2211 +- 0.002),
    /// STRICTLY lighter than dry air at the same T/P, with a speed of sound slightly ABOVE dry.
    #[test]
    fn test_mba1136_humid_lighter_than_dry() {
        let (dry_rho, dry_sos) = calculate_atmosphere(0.0, Some(15.0), Some(1013.25), 0.0);
        let (moist_rho, moist_sos) = calculate_atmosphere(0.0, Some(15.0), Some(1013.25), 50.0);

        assert!(
            (moist_rho - 1.2211).abs() < 0.002,
            "50% RH density {moist_rho} not within CIPM 1.2211 +- 0.002"
        );
        assert!(
            moist_rho < dry_rho,
            "moist air ({moist_rho}) must be lighter than dry ({dry_rho})"
        );
        assert!(
            moist_sos > dry_sos,
            "moist speed of sound ({moist_sos}) must exceed dry ({dry_sos})"
        );
    }

    /// Gate 3: density is monotone-decreasing in humidity (100% RH < 50% RH < 0% RH).
    #[test]
    fn test_mba1136_density_monotone_in_humidity() {
        let (rho_0, _) = calculate_atmosphere(0.0, Some(15.0), Some(1013.25), 0.0);
        let (rho_50, _) = calculate_atmosphere(0.0, Some(15.0), Some(1013.25), 50.0);
        let (rho_100, _) = calculate_atmosphere(0.0, Some(15.0), Some(1013.25), 100.0);
        assert!(
            rho_100 < rho_50 && rho_50 < rho_0,
            "humidity monotonicity violated: 100%={rho_100}, 50%={rho_50}, 0%={rho_0}"
        );
    }

    /// rank 28: `calculate_atmosphere`'s density is exactly `calculate_air_density_cimp` — there
    /// is a single canonical humid-air density path (no separate Arden-Buck ideal-gas density).
    #[test]
    fn test_mba1136_atmosphere_density_is_cipm() {
        for (t, p, h) in [
            (15.0, 1013.25, 0.0),
            (15.0, 1013.25, 50.0),
            (30.0, 1000.0, 80.0),
            (-10.0, 1020.0, 20.0),
        ] {
            let (density, _) = calculate_atmosphere(0.0, Some(t), Some(p), h);
            let cipm = calculate_air_density_cimp(t, p, h);
            assert_eq!(
                density, cipm,
                "calculate_atmosphere density must equal CIPM at {t}C/{p}hPa/{h}%"
            );
        }
    }

    /// rank 9: the extracted `moist_speed_of_sound` is exactly what `calculate_atmosphere`
    /// returns (behavior-identical extraction), across dry and humid conditions.
    #[test]
    fn test_mba1136_moist_speed_of_sound_extraction() {
        for (t, p, h) in [
            (15.0, 1013.25, 0.0),
            (15.0, 1013.25, 50.0),
            (25.0, 900.0, 100.0),
        ] {
            let (_, sos) = calculate_atmosphere(0.0, Some(t), Some(p), h);
            let extracted = moist_speed_of_sound(t + 273.15, p * 100.0, h);
            assert_eq!(
                sos, extracted,
                "extracted moist_speed_of_sound must match calculate_atmosphere at {t}C/{p}hPa/{h}%"
            );
        }
    }

    /// rank 9: local-atmosphere reference values (base: 500 m, 10 C, 950 hPa, ratio 1.05).
    /// The 1500 m values include the geometric-to-geopotential height conversion.
    #[test]
    fn test_mba1136_get_local_atmosphere_reference() {
        let (d0, c0) = get_local_atmosphere(500.0, 500.0, 10.0, 950.0, 1.05);
        assert!(
            (d0 - 1.286250000000).abs() < 1e-9,
            "local density@500m drifted: {d0}"
        );
        assert!(
            (c0 - 337.328657395129).abs() < 1e-9,
            "local sos@500m drifted: {c0}"
        );
        let (d1, c1) = get_local_atmosphere(1500.0, 500.0, 10.0, 950.0, 1.05);
        assert!(
            (d1 - 1.165238292559).abs() < 1e-9,
            "local density@1500m drifted: {d1}"
        );
        assert!(
            (c1 - 333.435546617978).abs() < 1e-9,
            "local sos@1500m drifted: {c1}"
        );
    }

    #[test]
    fn fractional_station_altitude_is_not_quantized() {
        let station_altitude_m = 500.25;
        let station_temp_c = 10.0;
        let station_pressure_hpa = 950.0;
        let station_density_ratio = 1.05;

        let (density, speed_of_sound) = get_local_atmosphere(
            station_altitude_m,
            station_altitude_m,
            station_temp_c,
            station_pressure_hpa,
            station_density_ratio,
        );

        let expected_density = station_density_ratio * 1.225;
        let expected_speed_of_sound = ((station_temp_c + 273.15) * 401.874).sqrt();
        assert!((density - expected_density).abs() < 1e-12);
        assert!((speed_of_sound - expected_speed_of_sound).abs() < 1e-12);
    }

    /// rank 9: `get_local_atmosphere_humid` returns the SAME density as `get_local_atmosphere`,
    /// and at 0% RH its speed of sound reduces to the dry value (within the 401.874-vs-gamma*R
    /// constant rounding). At real humidity the speed of sound exceeds the dry value.
    #[test]
    fn test_mba1136_get_local_atmosphere_humid() {
        let (d_dry, c_dry) = get_local_atmosphere(1500.0, 500.0, 10.0, 950.0, 1.05);
        let (d_h0, c_h0) = get_local_atmosphere_humid(1500.0, 500.0, 10.0, 950.0, 1.05, 0.0);
        assert_eq!(d_dry, d_h0, "humid variant must not change density");
        assert!(
            (c_h0 - c_dry).abs() < 1e-3,
            "0% RH humid sos {c_h0} should match dry sos {c_dry}"
        );
        let (_, c_h80) = get_local_atmosphere_humid(1500.0, 500.0, 10.0, 950.0, 1.05, 80.0);
        assert!(c_h80 > c_dry, "humid sos {c_h80} should exceed dry {c_dry}");
    }

    #[test]
    fn test_icao_standard_atmosphere() {
        // Test sea level
        let (temp, press) = calculate_icao_standard_atmosphere(0.0);
        assert!((temp - 288.15).abs() < 0.01);
        assert!((press - 101325.0).abs() < 1.0);

        // The table's 11 km tropopause base is geopotential height; convert it back to the
        // geometric altitude accepted by the public atmosphere contract.
        let geometric_tropopause_m =
            GEOPOTENTIAL_EARTH_RADIUS_M * 11000.0 / (GEOPOTENTIAL_EARTH_RADIUS_M - 11000.0);
        let (temp_11km, press_11km) = calculate_icao_standard_atmosphere(geometric_tropopause_m);
        assert!((temp_11km - 216.65).abs() < 0.01);
        assert!((press_11km - 22632.1).abs() < 1.0);

        // Test stratosphere
        let (temp_25km, _) = calculate_icao_standard_atmosphere(25000.0);
        assert!(temp_25km > 216.65); // Temperature increases in stratosphere
    }

    #[test]
    fn standard_atmosphere_extends_below_sea_level() {
        let altitude_m = -430.0;
        let (temp_k, pressure_pa) = calculate_icao_standard_atmosphere(altitude_m);

        assert!((temp_k - 290.945_189_079_054).abs() < 1e-9);
        assert!((pressure_pa - 106_598.763_552_437).abs() < 0.1);

        let (station_temp_c, station_pressure_hpa) =
            resolve_station_conditions(15.0, 1013.25, altitude_m);
        assert!((station_temp_c - 17.795_189_079_054).abs() < 1e-9);
        assert!((station_pressure_hpa - 1_065.987_635_524).abs() < 1e-6);

        let (sea_density, _) = calculate_atmosphere(0.0, None, None, 0.0);
        let (below_sea_density, _) = calculate_atmosphere(altitude_m, None, None, 0.0);
        assert!((below_sea_density - 1.276_908_642_79).abs() < 1e-6);
        assert!(below_sea_density > sea_density * 1.04);

        assert_eq!(
            calculate_icao_standard_atmosphere(-6000.0),
            calculate_icao_standard_atmosphere(MIN_GEOMETRIC_ALTITUDE_M)
        );
    }

    #[test]
    fn standard_atmosphere_converts_geometric_to_geopotential_height() {
        let cases: [(f64, f64, f64); 3] = [
            (10_000.0, 223.252_092_647_979, 26_499.901_600_244),
            (30_000.0, 226.509_083_611_330, 1_197.032_108_466),
            (84_000.0, 190.841_043_736_102, 0.531_525_514_935),
        ];
        for (geometric_m, expected_temp_k, expected_pressure_pa) in cases {
            let (temp_k, pressure_pa) = calculate_icao_standard_atmosphere(geometric_m);
            let pressure_tolerance = expected_pressure_pa.max(1.0_f64) * 1e-6;

            assert!((temp_k - expected_temp_k).abs() < 1e-9);
            assert!((pressure_pa - expected_pressure_pa).abs() < pressure_tolerance);
        }
    }

    #[test]
    fn local_atmosphere_walks_icao_layers_continuously() {
        for altitude_m in [
            10999.0, 11000.0, 11001.0, 11050.0, 19999.0, 20000.0, 20001.0, 25000.0,
        ] {
            let (local_temp_k, local_pressure_pa, _) =
                local_temp_pressure_density(altitude_m, 0.0, 15.0, 1013.25, 1.0);
            let (standard_temp_k, standard_pressure_pa) =
                calculate_icao_standard_atmosphere(altitude_m);

            assert!(
                (local_temp_k - standard_temp_k).abs() < 1e-9,
                "local temperature diverged from ICAO at {altitude_m} m: local={local_temp_k}, standard={standard_temp_k}"
            );
            assert!(
                ((local_pressure_pa - standard_pressure_pa) / standard_pressure_pa).abs() < 5e-5,
                "local pressure diverged from ICAO at {altitude_m} m: local={local_pressure_pa}, standard={standard_pressure_pa}"
            );
        }

        let (density_below, sound_below) =
            get_local_atmosphere(10999.0, 0.0, 15.0, 1013.25, 1.0);
        let (density_above, sound_above) =
            get_local_atmosphere(11001.0, 0.0, 15.0, 1013.25, 1.0);
        assert!(
            (density_below - density_above).abs() < 0.001,
            "density jumped across 11 km: below={density_below}, above={density_above}"
        );
        assert!(
            (sound_below - sound_above).abs() < 0.1,
            "speed of sound jumped across 11 km: below={sound_below}, above={sound_above}"
        );
    }

    #[test]
    fn local_atmosphere_preserves_nonstandard_station_offset_and_round_trips() {
        let base_alt = 7500.0;
        let base_temp_c = 5.0;
        let base_pressure_hpa = 410.0;
        let base_ratio = 0.72;

        let (high_temp_k, high_pressure_pa, high_density) = local_temp_pressure_density(
            25000.0,
            base_alt,
            base_temp_c,
            base_pressure_hpa,
            base_ratio,
        );
        assert!((high_temp_k - 260.244_615_053_376).abs() < 1e-9);
        assert!((high_pressure_pa / 100.0 - 40.964_358_485_456).abs() < 1e-9);
        assert!((high_density - 0.094_186_400_274).abs() < 1e-9);

        let (back_temp_k, back_pressure_pa, back_density) = local_temp_pressure_density(
            base_alt,
            25000.0,
            high_temp_k - 273.15,
            high_pressure_pa / 100.0,
            high_density / 1.225,
        );
        assert!((back_temp_k - (base_temp_c + 273.15)).abs() < 1e-9);
        assert!((back_pressure_pa / 100.0 - base_pressure_hpa).abs() < 1e-8);
        assert!((back_density - base_ratio * 1.225).abs() < 1e-9);
    }

    #[test]
    fn test_enhanced_atmosphere_sea_level() {
        let (density, speed) = calculate_atmosphere(0.0, None, None, 0.0);
        assert!((density - 1.225).abs() < 0.01);
        assert!((speed - 340.0).abs() < 1.0);
    }

    #[test]
    fn test_resolve_station_pressure_contract() {
        // Default sea-level pressure + real altitude => derive from altitude (None).
        assert_eq!(resolve_station_pressure(1013.25, 2000.0), None);
        // 29.92 inHg ≈ 1013.21 hPa is also treated as the default (within tolerance).
        assert_eq!(resolve_station_pressure(1013.21, 2000.0), None);
        // An explicit, non-default station pressure is authoritative (Some, used directly).
        assert_eq!(resolve_station_pressure(850.0, 2000.0), Some(850.0));
        // At/near sea level the default is used directly (no derivation needed).
        assert_eq!(resolve_station_pressure(1013.25, 0.0), Some(1013.25));
    }

    #[test]
    fn qnh_reduction_matches_hand_computed_value_at_nonzero_altitude() {
        // Hand-computed (Python, double precision) from the ICAO inverse-barometric formula.
        // NOTE the geopotential step — 1500 m geometric is 1499.646 m geopotential, and it is
        // the geopotential height that enters the exponent:
        //   h_geo   = 6356766*1500/(6356766+1500)                    = 1499.6461... m
        //   station = 1030.0 * (1 - 0.0065*h_geo/288.15)^5.255875601466713
        //           = 859.5753123926447 hPa
        // Feeding the raw 1500 m in instead yields 859.5378567 hPa, so a reader checking this
        // pin against the bare formula would wrongly conclude the implementation is off.
        let reduced = reduce_qnh_to_station_pressure(1030.0, 1500.0);
        assert!(
            (reduced - 859.575_312_392_644_7).abs() < 1e-9,
            "reduced={reduced}"
        );
        // The reduction must strictly lower the pressure versus the raw QNH input.
        assert!(reduced < 1030.0);
    }

    #[test]
    fn qnh_reduction_is_identity_at_sea_level() {
        // At h=0 the geopotential height is 0, so the ratio is exactly 1.0 and QNH passes
        // through unchanged (bit-exact: 1.0f64.powf(x) == 1.0).
        assert_eq!(reduce_qnh_to_station_pressure(1030.0, 0.0), 1030.0);
        assert_eq!(reduce_qnh_to_station_pressure(950.5, 0.0), 950.5);
    }

    #[test]
    fn qnh_sea_level_standard_matches_icao_standard_atmosphere_at_every_altitude() {
        // A QNH of exactly 1013.25 hPa (the sea-level standard) must reduce to precisely the
        // ICAO standard station pressure at any altitude -- this is what keeps the omitted-
        // pressure default byte-identical whether the caller declares Absolute or Qnh.
        for altitude_m in [0.0, 500.0, 2000.0, 4500.0] {
            let (_, std_pressure_pa) = calculate_icao_standard_atmosphere(altitude_m);
            let reduced = reduce_qnh_to_station_pressure(1013.25, altitude_m);
            assert!(
                (reduced - std_pressure_pa / 100.0).abs() < 1e-9,
                "altitude={altitude_m} reduced={reduced} std={}",
                std_pressure_pa / 100.0
            );
        }
    }

    #[test]
    fn resolve_station_pressure_with_mode_absolute_matches_resolve_station_pressure() {
        // PressureReferenceMode::Absolute must be byte-identical to resolve_station_pressure
        // for every case that function's own contract test exercises.
        for (pressure_hpa, altitude_m) in
            [(1013.25, 2000.0), (1013.21, 2000.0), (850.0, 2000.0), (1013.25, 0.0)]
        {
            assert_eq!(
                resolve_station_pressure_with_mode(
                    pressure_hpa,
                    altitude_m,
                    PressureReferenceMode::Absolute
                ),
                resolve_station_pressure(pressure_hpa, altitude_m)
            );
        }
    }

    #[test]
    fn resolve_station_pressure_with_mode_qnh_bypasses_the_default_sentinel() {
        // Constructed so the REDUCED value coincidentally lands within the absolute-mode
        // sentinel's +/-0.5 hPa tolerance of 1013.25 hPa at a real (>1 m) altitude -- exactly
        // the case that would silently discard an explicit reading under the old,
        // mode-blind resolve_station_pressure.
        let qnh = 1030.0;
        let altitude_m = 138.078_300_203_223_02;
        let reduced = reduce_qnh_to_station_pressure(qnh, altitude_m);
        assert!(
            (reduced - 1013.25).abs() < 0.5,
            "test fixture must land inside the sentinel band; reduced={reduced}"
        );

        // Qnh mode must still return Some(reduced) -- never None -- regardless of the
        // coincidence, and must equal the reduction, not a derived-from-altitude value.
        assert_eq!(
            resolve_station_pressure_with_mode(qnh, altitude_m, PressureReferenceMode::Qnh),
            Some(reduced)
        );

        // Sanity: the OLD absolute-mode function, given that same reduced number directly,
        // WOULD have discarded it (returned None, deriving ICAO-standard pressure instead) --
        // demonstrating why Qnh mode must never be routed back through the plain sentinel
        // check on its output.
        assert_eq!(resolve_station_pressure(reduced, altitude_m), None);
        let (_, std_pressure_pa) = calculate_icao_standard_atmosphere(altitude_m);
        assert!(
            (std_pressure_pa / 100.0 - reduced).abs() > 10.0,
            "fixture must show the sentinel misfire actually changes the answer materially"
        );
    }

    #[test]
    fn resolve_station_conditions_with_pressure_mode_absolute_matches_resolve_station_conditions()
    {
        for (temp_c, pressure_hpa, altitude_m) in
            [(15.0, 1013.25, 2000.0), (-5.0, 850.0, 2000.0), (15.0, 1013.25, 0.0)]
        {
            assert_eq!(
                resolve_station_conditions_with_pressure_mode(
                    temp_c,
                    pressure_hpa,
                    altitude_m,
                    PressureReferenceMode::Absolute
                ),
                resolve_station_conditions(temp_c, pressure_hpa, altitude_m)
            );
        }
    }

    #[test]
    fn resolve_station_conditions_with_pressure_mode_qnh_reduces_pressure_only() {
        // -5.0 C is not the sea-level default sentinel (15.0 C), so it stays authoritative --
        // isolating that only the pressure resolution differs between the two modes.
        let (temp_c, pressure_hpa) = resolve_station_conditions_with_pressure_mode(
            -5.0,
            1030.0,
            1500.0,
            PressureReferenceMode::Qnh,
        );
        // Temperature resolution is untouched by pressure mode: an explicit, non-default
        // temperature stays authoritative exactly as resolve_station_conditions would give.
        assert_eq!(temp_c, -5.0);
        assert!((pressure_hpa - 859.575_312_392_644_7).abs() < 1e-9);
    }

    #[test]
    fn test_altitude_affects_density_with_default_pressure() {
        // Regression: with the default pressure, altitude MUST lower density (previously the
        // air-density path ignored altitude whenever pressure was the sea-level default).
        let press = resolve_station_pressure(1013.25, 0.0);
        let (rho_sea, _) = calculate_atmosphere(0.0, Some(15.0), press, 50.0);
        let press_alt = resolve_station_pressure(1013.25, 2000.0);
        let (rho_2km, _) = calculate_atmosphere(2000.0, Some(15.0), press_alt, 50.0);
        assert!(
            rho_2km < rho_sea * 0.9,
            "density at 2000 m ({rho_2km}) should be well below sea level ({rho_sea})"
        );

        // But an explicit station pressure stays authoritative (no altitude double-count):
        // density with an explicit pressure is independent of the altitude field.
        let p = resolve_station_pressure(900.0, 2000.0);
        let (rho_a, _) = calculate_atmosphere(2000.0, Some(15.0), p, 50.0);
        let (rho_b, _) = calculate_atmosphere(0.0, Some(15.0), p, 50.0);
        assert!(
            (rho_a - rho_b).abs() < 1e-9,
            "explicit pressure must ignore altitude"
        );
    }

    #[test]
    fn test_resolve_station_temperature_contract() {
        // Default 15 C + real altitude => derive ICAO lapse temperature (None).
        assert_eq!(resolve_station_temperature(15.0, 2000.0), None);
        // An explicit, non-default temperature is authoritative (Some, used directly).
        assert_eq!(resolve_station_temperature(-5.0, 2000.0), Some(-5.0));
        assert_eq!(resolve_station_temperature(30.0, 2000.0), Some(30.0));
        // At/near sea level the default is used directly (no derivation needed).
        assert_eq!(resolve_station_temperature(15.0, 0.0), Some(15.0));
    }

    #[test]
    fn test_altitude_only_default_matches_full_icao_standard() {
        // Regression: resolving BOTH temperature and pressure for an altitude-only query (defaults
        // left in place) must equal the fully-standard atmosphere at that altitude — i.e. altitude
        // now drives temperature (ICAO lapse) AND pressure, not just pressure. Validated against
        // py_ballisticcalc to ~0.04%. Previously the air held 15 C, leaving density ~7% too thin
        // (warm) at 3 km.
        for alt in [1000.0, 2000.0, 2500.0, 3000.0] {
            let t = resolve_station_temperature(15.0, alt);
            let p = resolve_station_pressure(1013.25, alt);
            let (rho_resolved, _) = calculate_atmosphere(alt, t, p, 0.0);
            let (rho_std, _) = calculate_atmosphere(alt, None, None, 0.0);
            assert!(
                (rho_resolved - rho_std).abs() < 1e-9,
                "alt {alt}: altitude-only default density {rho_resolved} should equal the full \
                 ICAO standard {rho_std}"
            );
            // And it must be denser than the old temperature-held-at-15C behavior (colder = denser).
            let (rho_warm, _) = calculate_atmosphere(alt, Some(15.0), p, 0.0);
            assert!(
                rho_resolved > rho_warm,
                "alt {alt}: lapse-temperature density {rho_resolved} should exceed 15 C-held {rho_warm}"
            );
        }
    }

    #[test]
    fn test_enhanced_atmosphere_with_humidity() {
        let (density_dry, speed_dry) = calculate_atmosphere(0.0, None, None, 0.0);
        let (density_humid, speed_humid) = calculate_atmosphere(0.0, None, None, 80.0);

        // Humid air should be less dense
        assert!(density_humid < density_dry);
        // Humid air should have slightly higher speed of sound
        assert!(speed_humid > speed_dry);
    }

    #[test]
    fn test_enhanced_atmosphere_stratosphere() {
        // Test in stratosphere where temperature increases
        let (density_20km, speed_20km) = calculate_atmosphere(20000.0, None, None, 0.0);
        let (density_30km, speed_30km) = calculate_atmosphere(30000.0, None, None, 0.0);

        // Density should decrease with altitude
        assert!(density_30km < density_20km);
        // Speed of sound should increase due to temperature increase
        assert!(speed_30km > speed_20km);
    }

    #[test]
    fn test_enhanced_cimp_density() {
        let density = calculate_air_density_cimp(15.0, 1013.25, 0.0);
        assert!((density - 1.225).abs() < 0.01);

        // Test with humidity
        let density_humid = calculate_air_density_cimp(15.0, 1013.25, 50.0);
        assert!(density_humid < density);
    }

    #[test]
    fn test_cipm_moist_air_matches_python_reference() {
        // Regression for the hPa/Pa mole-fraction slip: p_v (hPa) was divided by the total
        // pressure in Pa, making x_v 100x too small and erasing the humidity effect entirely.
        // Reference values computed with the validated Python implementation
        // (ballistics.physics.atmosphere_icao.calculate_air_density_cipm_icao), same cases as
        // the Flask suite's tests/test_atmosphere.py::TestCalculateAirDensityCIPM. Tolerance
        // 0.1% (matches that suite's Rust-vs-Python assertion).
        let cases = [
            (15.0, 1013.25, 50.0, 1.221125867723075),
            (30.0, 1000.0, 80.0, 1.1344071877123691),
            (-10.0, 1020.0, 20.0, 1.3500610713710515),
        ];
        for (temp_c, pressure_hpa, humidity_pct, expected) in cases {
            let density = calculate_air_density_cipm(temp_c, pressure_hpa, humidity_pct);
            let rel_err = ((density - expected) / expected).abs();
            assert!(
                rel_err < 1e-3,
                "CIPM density at {temp_c} C / {pressure_hpa} hPa / {humidity_pct}% RH: \
                 got {density}, expected {expected} (rel err {rel_err:.2e} >= 1e-3)"
            );
        }

        // Moist air must be materially lighter than dry air at the same temp/pressure
        // (the broken version returned a difference of only ~4e-5 kg/m^3).
        let dry = calculate_air_density_cipm(15.0, 1013.25, 0.0);
        let moist = calculate_air_density_cipm(15.0, 1013.25, 50.0);
        assert!(
            dry - moist > 3e-3,
            "humidity effect too small: dry {dry} vs 50% RH {moist}"
        );
    }

    #[test]
    fn test_variable_lapse_rates() {
        // Test that lapse rates change appropriately with altitude
        let lapse_tropo = determine_local_lapse_rate(5000.0);
        let lapse_strato = determine_local_lapse_rate(25000.0);

        assert!((lapse_tropo - (-0.0065)).abs() < 0.0001);
        assert!(lapse_strato > 0.0); // Positive lapse rate in stratosphere
    }

    // ---- MBA-1137: AtmoSock stateless downrange lookup (mirrors the WindSock tests) ----

    #[test]
    fn test_atmo_sock_empty_falls_back_to_isa() {
        let sock = AtmoSock::new(vec![]);
        assert!(sock.is_empty());
        // Empty sock returns the sea-level ISA reference regardless of distance.
        assert_eq!(sock.atmo_for_range(0.0), (15.0, 1013.25, 0.0));
        assert_eq!(sock.atmo_for_range(500.0), (15.0, 1013.25, 0.0));
    }

    #[test]
    fn test_atmo_sock_single_segment_holds_beyond_last() {
        // One zone until 100 m; it must apply BOTH before and beyond the threshold (unlike wind,
        // which zeroes past the last segment — atmosphere holds the last zone).
        let sock = AtmoSock::new(vec![(25.0, 1000.0, 30.0, 100.0)]);
        assert_eq!(sock.atmo_for_range(50.0), (25.0, 1000.0, 30.0));
        assert_eq!(sock.atmo_for_range(100.0), (25.0, 1000.0, 30.0)); // beyond last -> hold
        assert_eq!(sock.atmo_for_range(5000.0), (25.0, 1000.0, 30.0));
    }

    #[test]
    fn test_atmo_sock_boundary_is_upper_exclusive() {
        // A zone's until_distance_m is exclusive: a query exactly at the boundary rolls to the
        // next zone (mirrors WindSock::test_wind_sock_boundary_is_upper_exclusive).
        let sock = AtmoSock::new(vec![
            (30.0, 1010.0, 80.0, 100.0), // hot/humid near zone
            (-5.0, 900.0, 10.0, 200.0),  // cold/thin far zone
        ]);
        // Just below 100 m -> first zone.
        assert_eq!(sock.atmo_for_range(99.999), (30.0, 1010.0, 80.0));
        // Exactly 100 m -> second zone.
        assert_eq!(sock.atmo_for_range(100.0), (-5.0, 900.0, 10.0));
        // Beyond the last boundary -> hold the last zone (NOT zeroed).
        assert_eq!(sock.atmo_for_range(200.0), (-5.0, 900.0, 10.0));
        assert_eq!(sock.atmo_for_range(1e6), (-5.0, 900.0, 10.0));
    }

    #[test]
    fn test_atmo_sock_sorts_unordered_segments() {
        // Segments supplied out of order are sorted by until_distance so the lookup is monotone.
        let sock = AtmoSock::new(vec![
            (-5.0, 900.0, 10.0, 200.0),
            (30.0, 1010.0, 80.0, 100.0),
        ]);
        assert_eq!(sock.atmo_for_range(50.0), (30.0, 1010.0, 80.0));
        assert_eq!(sock.atmo_for_range(150.0), (-5.0, 900.0, 10.0));
    }

    #[test]
    fn test_atmo_sock_orders_nan_thresholds_last() {
        let positive_nan = f64::from_bits(0x7ff8_0000_0000_0001);
        let negative_nan = f64::from_bits(0xfff8_0000_0000_0002);
        let sock = AtmoSock::new(vec![
            (97.0, 997.0, 97.0, positive_nan),
            (30.0, 1030.0, 30.0, 300.0),
            (98.0, 998.0, 98.0, negative_nan),
            (10.0, 1010.0, 10.0, 100.0),
            (99.0, 999.0, 99.0, f64::NAN),
            (20.0, 1020.0, 20.0, 200.0),
        ]);

        let thresholds: Vec<_> = sock.segments.iter().map(|segment| segment.3).collect();
        assert_eq!(&thresholds[..3], &[100.0, 200.0, 300.0]);
        assert!(thresholds[3..].iter().all(|threshold| threshold.is_nan()));
        assert_eq!(sock.atmo_for_range(f64::NAN), (10.0, 1010.0, 10.0));
        assert_eq!(sock.atmo_for_range(350.0), (99.0, 999.0, 99.0));
    }

    #[test]
    fn test_atmo_sock_nan_uses_first_zone() {
        let sock = AtmoSock::new(vec![
            (30.0, 1010.0, 80.0, 100.0),
            (-5.0, 900.0, 10.0, 200.0),
        ]);
        // NaN can't be ordered; deterministically use the nearest (first) zone rather than panic.
        assert_eq!(sock.atmo_for_range(f64::NAN), (30.0, 1010.0, 80.0));
    }

    // ---- MBA-1366: density altitude as a direct atmosphere input ----

    /// With no explicit temperature, the resolved altitude must equal the input density
    /// altitude EXACTLY (this is the algebraic identity documented on
    /// `resolve_atmosphere_for_density_altitude`: the correction term vanishes when station
    /// temperature is left at ISA).
    #[test]
    fn density_altitude_default_temp_pressure_alt_equals_density_alt() {
        for da_ft in [0.0_f64, 1000.0, 3000.0, 7500.0, -500.0] {
            let da_m = da_ft * DA_FT_TO_M;
            let (altitude_m, temp_c, pressure_hpa) =
                resolve_atmosphere_for_density_altitude(da_m, None);
            assert!(
                (altitude_m - da_m).abs() < 1e-6,
                "da_ft={da_ft}: altitude_m {altitude_m} should equal da_m {da_m}"
            );
            // Sanity: pressure must be a plausible station pressure (strictly decreasing with DA).
            assert!(pressure_hpa > 0.0 && pressure_hpa < 1100.0);
            assert!(temp_c.is_finite());
        }
        // At DA=0 the result must be exactly the sea-level ISA reference.
        let (alt0, temp0, press0) = resolve_atmosphere_for_density_altitude(0.0, None);
        assert!(alt0.abs() < 1e-9);
        assert!((temp0 - 15.0).abs() < 1e-6);
        assert!((press0 - 1013.25).abs() < 1e-6);
    }

    /// Higher density altitude at the ISA default must yield a lower station pressure (thinner
    /// air), matching the physical meaning of density altitude.
    #[test]
    fn density_altitude_pressure_decreases_monotonically() {
        let (_, _, p0) = resolve_atmosphere_for_density_altitude(0.0, None);
        let (_, _, p1) = resolve_atmosphere_for_density_altitude(1000.0, None);
        let (_, _, p2) = resolve_atmosphere_for_density_altitude(3000.0, None);
        assert!(p1 < p0);
        assert!(p2 < p1);
    }

    /// An explicit temperature must be honored EXACTLY (never re-derived), while the implied
    /// pressure altitude (and therefore the resolved altitude/pressure) differs from the
    /// ISA-default case -- proving density is still honored (a different pressure/altitude
    /// combination that reproduces the SAME density altitude at the warmer/colder temperature).
    #[test]
    fn density_altitude_explicit_temperature_overrides_isa_default_but_honors_density() {
        let da_m = 1000.0 * DA_FT_TO_M;
        let (default_alt_m, default_temp_c, default_pressure_hpa) =
            resolve_atmosphere_for_density_altitude(da_m, None);

        // A HOTTER-than-ISA explicit station temperature at the same density altitude implies a
        // LOWER pressure altitude: hot (less dense) air already accounts for some of the
        // "thinness" that defines the density altitude, so less physical elevation is needed to
        // reach the rest of it (equivalently: the correction term subtracts from DA before the
        // pressure inversion, since a warmer station reads MORE dense-feeling per unit of actual
        // elevation than ISA does).
        let hot_temp_c = default_temp_c + 20.0;
        let (hot_alt_m, hot_temp_out, hot_pressure_hpa) =
            resolve_atmosphere_for_density_altitude(da_m, Some(hot_temp_c));
        assert_eq!(hot_temp_out, hot_temp_c, "explicit temperature must be honored exactly");
        assert!(
            hot_alt_m < default_alt_m,
            "hotter station temp at same DA should imply a lower pressure altitude: \
             hot={hot_alt_m} default={default_alt_m}"
        );
        assert!(hot_pressure_hpa > default_pressure_hpa);

        // A COLDER-than-ISA explicit temperature implies a HIGHER pressure altitude.
        let cold_temp_c = default_temp_c - 20.0;
        let (cold_alt_m, cold_temp_out, cold_pressure_hpa) =
            resolve_atmosphere_for_density_altitude(da_m, Some(cold_temp_c));
        assert_eq!(cold_temp_out, cold_temp_c);
        assert!(cold_alt_m > default_alt_m);
        assert!(cold_pressure_hpa < default_pressure_hpa);
    }

    /// Round trip using ONLY this crate's own math (the forward NWS pressure-altitude formula,
    /// re-derived independently in the test rather than reusing the production function) -- an
    /// implementation-independent algebraic sanity check. The REAL cross-crate round trip
    /// against the production `pdf_dope_card::calculate_density_altitude` forward function
    /// lives in `src/main.rs` (that function is CLI-binary-only and unreachable from this
    /// library crate).
    #[test]
    fn density_altitude_round_trips_through_hand_derived_forward_formula() {
        fn forward_density_altitude_ft(pressure_hpa: f64, temp_c: f64) -> f64 {
            let pressure_alt_ft =
                145_366.45 * (1.0 - (pressure_hpa / 1013.25).powf(0.190_284));
            let isa_temp_f = 59.0 - (pressure_alt_ft / 1000.0) * 3.57;
            let temp_f = temp_c * 9.0 / 5.0 + 32.0;
            pressure_alt_ft + (120.0 * 5.0 / 9.0) * (temp_f - isa_temp_f)
        }

        for da_ft in [0.0_f64, 500.0, 2500.0, 6000.0] {
            let da_m = da_ft * DA_FT_TO_M;
            // Default (ISA-at-DA) branch.
            let (_, temp_c, pressure_hpa) = resolve_atmosphere_for_density_altitude(da_m, None);
            let round_tripped_ft = forward_density_altitude_ft(pressure_hpa, temp_c);
            assert!(
                (round_tripped_ft - da_ft).abs() < 1e-6,
                "da_ft={da_ft}: round-tripped {round_tripped_ft} (temp_c={temp_c}, pressure_hpa={pressure_hpa})"
            );

            // Explicit-temperature branch.
            let explicit_temp_c = 25.0;
            let (_, temp_out, pressure_hpa_explicit) =
                resolve_atmosphere_for_density_altitude(da_m, Some(explicit_temp_c));
            assert_eq!(temp_out, explicit_temp_c);
            let round_tripped_explicit_ft =
                forward_density_altitude_ft(pressure_hpa_explicit, temp_out);
            assert!(
                (round_tripped_explicit_ft - da_ft).abs() < 1e-6,
                "da_ft={da_ft} (explicit {explicit_temp_c}C): round-tripped {round_tripped_explicit_ft}"
            );
        }
    }
}