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 crate::gravity::Gravity;
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
// ICAO reference quantities are defined using standard Earth gravity, independently of the
// acceleration chosen for an internal simulation.
const G_ACCEL_MPS2: f64 = crate::constants::G_ACCEL_MPS2;
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)
}

/// Whether a pressure lies inside the legacy default sentinel's band, within 0.5 hPa of the
/// 1013.25 hPa sea-level standard. [`resolve_station_pressure`] reads an unmarked value there
/// as omitted at any altitude above 1 m.
pub(crate) fn in_legacy_default_pressure_band(pressure_hpa: f64) -> bool {
    const SEA_LEVEL_HPA: f64 = 1013.25;
    (pressure_hpa - SEA_LEVEL_HPA).abs() < 0.5
}

/// Whether a temperature lies inside the legacy default sentinel's band, within 0.1 °C of the
/// 15 °C sea-level standard. [`resolve_station_temperature`] reads an unmarked value there as
/// omitted at any altitude above 1 m.
pub(crate) fn in_legacy_default_temperature_band(temperature_c: f64) -> bool {
    const SEA_LEVEL_TEMP_C: f64 = 15.0;
    (temperature_c - SEA_LEVEL_TEMP_C).abs() < 0.1
}

/// Resolve the station-pressure override for an air-density calculation, by the legacy
/// default-sentinel convention.
///
/// Presence-blind legacy convention: a bare `f64` cannot say whether it was supplied, so a
/// sea-level default at a real altitude is read as omitted. The rule:
/// * A pressure within 0.5 hPa of the 1013.25 hPa sea-level standard at an altitude above 1 m
///   is read as "omitted, standard atmosphere at this altitude": `None`, so
///   [`calculate_atmosphere`] derives the station pressure from altitude (ICAO standard)
///   instead of flying sea-level density there.
/// * Any other value is returned as `Some` and used as the STATION pressure, so altitude is
///   not double-counted.
///
/// The ±0.5 hPa tolerance covers the `29.92 inHg ≈ 1013.21 hPa` conversion, and `>1 m` avoids
/// triggering at sea level. A real reading inside that band (QNH 1031 hPa at 150 m reduces to
/// 1012.86 hPa) is indistinguishable from the default here. A caller that knows presence uses
/// [`resolve_station_conditions_with_presence`], which flies a supplied value as given.
pub fn resolve_station_pressure(pressure_hpa: f64, altitude_m: f64) -> Option<f64> {
    if in_legacy_default_pressure_band(pressure_hpa) && 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, by the legacy
/// default-sentinel convention, mirroring [`resolve_station_pressure`].
///
/// Presence-blind legacy convention: a bare `f64` cannot say whether it was supplied, so a
/// sea-level default at a real altitude is read as omitted.
/// * A temperature within 0.1 °C of the 15 °C sea-level standard at an altitude above 1 m is
///   read as "omitted, standard atmosphere at this altitude": `None`, so
///   [`calculate_atmosphere`] applies the ICAO lapse-rate temperature for that altitude
///   (≈ −6.5 °C/km).
/// * Any other value is returned as `Some` and used as given.
///
/// Without the omitted case, `--altitude` with the default temperature held the air at 15 °C,
/// which under-estimates density by ~2.4% at 1 km up to ~7% at 3 km versus the standard
/// atmosphere (validated against py_ballisticcalc). The 0.1 °C tolerance matches the
/// `59 °F = 15.0 °C` default exactly, and `>1 m` avoids triggering at sea level. A real 15 °C
/// reading at altitude is indistinguishable from the default here; a caller that knows presence
/// uses [`resolve_station_conditions_with_presence`].
pub fn resolve_station_temperature(temperature_c: f64, altitude_m: f64) -> Option<f64> {
    if in_legacy_default_temperature_band(temperature_c) && 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.
///
/// Presence-blind: both values go through the legacy default sentinels
/// ([`resolve_station_temperature`], [`resolve_station_pressure`]). See
/// [`resolve_station_conditions_with_presence`].
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 altimeter-setting relation (MBA-1605):
///
/// `station = QNH * (1 + (L*h/T0) * (P0/QNH)^n)^(1/n)`, with `n = -L*R/g`
///
/// for the ICAO troposphere's lapse rate `L` (-0.0065 K/m), sea-level temperature `T0`
/// (288.15 K) and sea-level pressure `P0` (1013.25 hPa), so `n` = 0.190263. This is what an
/// altimeter setting means: an altimeter set to QNH reads field elevation on the ground, so the
/// station pressure's ISA pressure altitude is the QNH's own pressure altitude plus `h`. Written
/// that way it is the additive form `station = P0 * ((QNH/P0)^n + L*h/T0)^(1/n)`; the
/// algebraically equal form above is the one evaluated, because it returns QNH bit-exactly at
/// `h = 0`.
///
/// This is the altimeter-setting relation on ICAO constants with geopotential height. It has the
/// shape of the NWS/ASOS altimeter-setting algorithm but is not that algorithm, which uses its
/// own rounded constants (`n` = 0.190284 and `k` = 8.4228807e-5, `P0^n * 0.0065/T0` at 288 K),
/// adds the field elevation as entered rather than its geopotential height, and subtracts a
/// 0.3 hPa sensor offset. Even without that offset, NWS/ASOS gives a station pressure about
/// 0.1 hPa lower than this function at 1500 m and 0.3 to 0.45 hPa lower at 4000 to 5000 m (QNH
/// 950 to 1050 hPa).
///
/// Before 0.49.0 this function scaled the standard pressure ratio by QNH,
/// `station = QNH * (1 + L*h/T0)^(1/n)`, and called that the altimeter-setting reduction. It is
/// not: the two agree only at QNH = 1013.25 hPa or at `h = 0`, and the ratio form ran 0.49 hPa
/// low for a QNH of 1030 hPa at 1500 m.
///
/// `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. `altitude_m` is clamped and
/// converted to geopotential height exactly like `calculate_icao_standard_atmosphere`, and
/// `(P0/QNH)^n` is exactly 1 at QNH = 1013.25 hPa. At QNH = 1013.25 hPa the reduction therefore
/// equals the ICAO standard station pressure to rounding, but not bit for bit at every altitude
/// (for example 1012.0492647136854 against 1012.0492647136847 at 10 m; the two agree exactly at
/// 150 m). Callers that must stay bit-identical for an omitted pressure keep the mode they
/// already used for it.
///
/// 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. No
/// station pressure exists once the QNH's pressure altitude plus `h` reaches `T0/-L`
/// (44,331 m, where the troposphere's lapse rate reaches absolute zero), nor for a non-positive
/// QNH; there this returns NaN, which callers refuse as a non-finite pressure.
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 relation covers
    let height_term = layer.lapse_rate * geopotential_height / layer.base_temperature;
    let exponent = -G_ACCEL_MPS2 / (layer.lapse_rate * R_AIR); // 1/n
    let sea_level_hpa = layer.base_pressure / 100.0;
    let setting_term = (sea_level_hpa / qnh_hpa).powf(exponent.recip()); // (P0/QNH)^n
    qnh_hpa * (1.0 + height_term * setting_term).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.
///
/// Presence-blind: this is [`resolve_station_conditions_with_presence`] with
/// [`StationAtmospherePresence::NEITHER`], bit for bit. The temperature sentinel runs in BOTH
/// arms, so a typed 15 °C at altitude flies the ICAO lapse temperature here even with a
/// declared QNH; only the pressure sentinel is bypassed in the `Qnh` arm. A caller that knows
/// which values were supplied uses [`resolve_station_conditions_with_presence`].
pub fn resolve_station_conditions_with_pressure_mode(
    temperature_c: f64,
    pressure_hpa: f64,
    altitude_m: f64,
    pressure_mode: PressureReferenceMode,
) -> (f64, f64) {
    resolve_station_conditions_with_presence(
        temperature_c,
        pressure_hpa,
        altitude_m,
        pressure_mode,
        StationAtmospherePresence::NEITHER,
    )
}

/// Which of a station temperature / pressure pair the caller actually supplied (0.49.0).
///
/// A field that is not supplied is resolved by the pre-0.49.0 default-sentinel rule
/// ([`resolve_station_temperature`] / [`resolve_station_pressure_with_mode`]): the ICAO
/// standard at |altitude| > 1 m when the value is the sea-level default, the value itself
/// otherwise. That makes [`Self::NEITHER`] exactly the old behaviour of every presence-blind
/// API, which is why it is the `Default`, and makes a supplied field authoritative whatever
/// its value: a typed 15 °C, 59 °F, 1013.25 hPa or 29.92 inHg at altitude flies as typed.
///
/// Pass it to [`resolve_station_conditions_with_presence`], or to one of the
/// `*_with_station_presence` twins of the solver, zero, Monte Carlo, lead, hold-curve and
/// truing entry points.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default)]
pub struct StationAtmospherePresence {
    /// The station temperature was a reading the caller supplied.
    pub temperature_supplied: bool,
    /// The station pressure (or the QNH it is reduced from) was a reading the caller supplied.
    pub pressure_supplied: bool,
}

impl StationAtmospherePresence {
    /// Neither value was supplied: the pre-0.49.0 default-sentinel rule, bit for bit.
    pub const NEITHER: Self = Self::new(false, false);
    /// Both values were supplied: both fly as given (a declared QNH after its reduction).
    pub const BOTH: Self = Self::new(true, true);

    /// Presence from two flags.
    pub const fn new(temperature_supplied: bool, pressure_supplied: bool) -> Self {
        Self {
            temperature_supplied,
            pressure_supplied,
        }
    }

    /// Presence of a temperature reading and a pressure reading.
    pub const fn from_readings(temperature: StationReading, pressure: StationReading) -> Self {
        Self::new(temperature.supplied, pressure.supplied)
    }
}

/// A station temperature or pressure on a presence-carrying request (0.49.0): the value the
/// engine is handed and whether the caller supplied it.
///
/// On the wire it is a bare JSON number when supplied and ABSENT when omitted. The type's own
/// `Deserialize` rejects a JSON `null`; the truing DTOs (`TruingModelInputsV1`,
/// `WindTruingRequest`) read a `null` field as omitted, as the `card.*` requests do. An omitted
/// reading carries the surface's legacy default (59 °F / 29.92 inHg on
/// the truing DTOs, the metric CLI's 1013.25 hPa expressed in inHg, and so on), so it resolves
/// exactly as that default did before 0.49.0. Use [`Self::omitted`] for that, never a bare
/// `Option`: at |altitude| <= 1 m an omitted value flies as itself, so the default it carries
/// is part of the result.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct StationReading {
    /// The value handed to the engine.
    pub value: f64,
    /// Whether the caller supplied it. `false` means the legacy default-sentinel rule decides.
    pub supplied: bool,
}

impl StationReading {
    /// A reading the caller supplied: it flies as given.
    pub const fn supplied(value: f64) -> Self {
        Self {
            value,
            supplied: true,
        }
    }

    /// No reading: `default_value` is the surface's legacy default, resolved by the
    /// pre-0.49.0 rule (ICAO standard at a real altitude when it is the sea-level default).
    pub const fn omitted(default_value: f64) -> Self {
        Self {
            value: default_value,
            supplied: false,
        }
    }

    /// `true` when no reading was supplied. Used as `skip_serializing_if`.
    pub fn is_omitted(&self) -> bool {
        !self.supplied
    }
}

impl serde::Serialize for StationReading {
    fn serialize<S: serde::Serializer>(&self, serializer: S) -> Result<S::Ok, S::Error> {
        serializer.serialize_f64(self.value)
    }
}

impl<'de> serde::Deserialize<'de> for StationReading {
    fn deserialize<D: serde::Deserializer<'de>>(deserializer: D) -> Result<Self, D::Error> {
        <f64 as serde::Deserialize>::deserialize(deserializer).map(Self::supplied)
    }
}

/// Resolve the station temperature and pressure a solve flies, honouring which of them the
/// caller supplied (0.49.0). This is the presence-aware resolver every surface that knows
/// presence goes through.
///
/// For each field:
/// * supplied: returned as given; a pressure declared as QNH ([`PressureReferenceMode::Qnh`])
///   is first reduced to station pressure with [`reduce_qnh_to_station_pressure`];
/// * omitted: resolved exactly as before 0.49.0, through [`resolve_station_temperature`] /
///   [`resolve_station_pressure_with_mode`], then the ICAO standard at `altitude_m` when the
///   value is the sea-level default at a real altitude.
///
/// With [`StationAtmospherePresence::NEITHER`] the result is bit-identical to
/// [`resolve_station_conditions_with_pressure_mode`] (which now delegates here). With
/// [`StationAtmospherePresence::BOTH`] in [`PressureReferenceMode::Absolute`] it returns the
/// pair unchanged.
pub fn resolve_station_conditions_with_presence(
    temperature_c: f64,
    pressure_hpa: f64,
    altitude_m: f64,
    pressure_mode: PressureReferenceMode,
    presence: StationAtmospherePresence,
) -> (f64, f64) {
    let temp_override = if presence.temperature_supplied {
        Some(temperature_c)
    } else {
        resolve_station_temperature(temperature_c, altitude_m)
    };
    let press_override = if presence.pressure_supplied {
        Some(match pressure_mode {
            PressureReferenceMode::Absolute => pressure_hpa,
            PressureReferenceMode::Qnh => reduce_qnh_to_station_pressure(pressure_hpa, altitude_m),
        })
    } else {
        resolve_station_pressure_with_mode(pressure_hpa, altitude_m, pressure_mode)
    };
    if let (Some(temp_c), Some(pressure_hpa)) = (temp_override, press_override) {
        // Value-neutral fast path: the ICAO evaluation below would be discarded.
        return (temp_c, pressure_hpa);
    }
    let (std_temp_k, std_pressure_pa) = calculate_icao_standard_atmosphere(altitude_m);
    (
        temp_override.unwrap_or(std_temp_k - 273.15),
        press_override.unwrap_or(std_pressure_pa / 100.0),
    )
}

// ---------------------------------------------------------------------------------------
// Station-pressure plausibility
// ---------------------------------------------------------------------------------------

/// The lowest ratio of station pressure to the ICAO standard pressure at the station's
/// altitude that [`station_pressure_plausibility`] treats as weather. See that function for
/// where the band comes from.
pub const STATION_PRESSURE_MIN_PLAUSIBLE_RATIO: f64 = 0.85;
/// The highest ratio of station pressure to the ICAO standard pressure at the station's
/// altitude that [`station_pressure_plausibility`] treats as weather.
pub const STATION_PRESSURE_MAX_PLAUSIBLE_RATIO: f64 = 1.08;

/// How an explicit absolute station pressure compares with the ICAO standard pressure at the
/// altitude it was sent with. Returned by [`station_pressure_plausibility`].
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct StationPressurePlausibility {
    /// The ICAO standard pressure at the altitude, in pascals.
    pub icao_standard_pa: f64,
    /// The station pressure divided by `icao_standard_pa`.
    pub ratio: f64,
    /// Whether `ratio` lies outside [`STATION_PRESSURE_MIN_PLAUSIBLE_RATIO`] to
    /// [`STATION_PRESSURE_MAX_PLAUSIBLE_RATIO`] (both edges inside), or is not finite.
    pub implausible: bool,
}

/// Compare an explicit absolute station pressure, in pascals, with the ICAO standard pressure
/// at `altitude_m`, and say whether the two can describe the same station.
///
/// An explicit absolute pressure is authoritative wherever the engine takes one, so a pressure
/// that belongs to some other altitude is flown as it was sent. The two common cases are a
/// sea-level figure sent at a real elevation: a 101325 Pa placeholder, or a weather report's
/// altimeter setting (QNH) sent as station pressure instead of as a QNH. At Cusco's 3310.128 m
/// a 101325 Pa placeholder is 1.503 times the standard pressure there. The reverse, a true
/// station reading sent with the wrong altitude (omitted, or in the wrong unit), shows up as a
/// ratio far from 1 too. Callers use this to warn; nothing here changes a pressure.
///
/// The band, [`STATION_PRESSURE_MIN_PLAUSIBLE_RATIO`] to
/// [`STATION_PRESSURE_MAX_PLAUSIBLE_RATIO`], is wide enough that recorded weather never leaves
/// it. The WMO's table of weather and climate extremes (as of 31 July 2025) lists these
/// sea-level pressure records:
///
/// - 870 hPa in the eye of Typhoon Tip, 12 October 1979, the lowest (tornadoes excluded);
/// - 1083.8 hPa at Agata, Russia, 31 December 1968, the highest at a station below 750 m;
/// - 1089.1 hPa at Tosontsengel, Mongolia, 30 December 2004, the highest above 750 m.
///
/// Against the 1013.25 hPa standard those are 0.859, 1.070 and 1.075. They are sea-level
/// figures. At a station above sea level the high extreme moves toward standard, because the
/// record highs come from cold air: reduced to sea level through a column denser than the
/// standard one, a station pressure nearer standard becomes the record sea-level figure. Cold
/// air moves the low side the other way. Under a cold column an ordinary sea-level pressure
/// puts a high station further below standard than the sea-level figure is: a 3000 m station
/// under a 243 K column reads 0.948 times the ICAO standard there with 1013.25 hPa at sea
/// level, and 0.917 with 980 hPa. Both stay well above 0.85.
///
/// The band catches gross mismatches, not every placeholder: 101325 Pa stays inside it up to
/// about 645 m, where it is already 8% above standard.
///
/// This is for absolute station pressure only. A QNH is a sea-level figure by declaration, so
/// comparing it with the standard at the station's altitude would say nothing about it.
pub fn station_pressure_plausibility(
    station_pressure_pa: f64,
    altitude_m: f64,
) -> StationPressurePlausibility {
    let (_, icao_standard_pa) = calculate_icao_standard_atmosphere(altitude_m);
    let ratio = station_pressure_pa / icao_standard_pa;
    let implausible = !(STATION_PRESSURE_MIN_PLAUSIBLE_RATIO
        ..=STATION_PRESSURE_MAX_PLAUSIBLE_RATIO)
        .contains(&ratio);
    StationPressurePlausibility {
        icao_standard_pa,
        ratio,
        implausible,
    }
}

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

/// NWS pressure-altitude and FAA-rule constants — the ONE copy.
///
/// Both directions of the density-altitude model now live in this module: the forward
/// [`faa_rule_density_altitude_ft`] (the engine's reported density altitude) and its exact
/// inverse [`resolve_atmosphere_for_density_altitude`] (MBA-1366, the DA entry mode). They
/// used to live apart, with these constants duplicated here because the forward formula sat
/// in `pdf_dope_card` — first in the `ballistics` binary, then (from 0.33.0's Task 10) in
/// this library but behind `#[cfg(feature = "pdf")]`, which this unconditional module could
/// not import across without breaking the `--no-default-features` wasm32 build CI
/// type-checks.
///
/// The fix was to move the formula DOWN rather than import it UP: `atmosphere` is
/// unconditional, so the forward function lives here and
/// `pdf_dope_card::calculate_density_altitude` now delegates to it, bit for bit. That leaves
/// no second copy to drift, and it is what lets the `atmosphere.density_altitude` bridge
/// command exist in every build with the bridge — including wasm32, which has no PDF card.
///
/// `src/main.rs`'s `density_altitude_round_trips_through_the_dope_card_formula*` tests still
/// pin forward and inverse together, now through a single set of constants.
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;
/// The FAA rule of thumb, "120 ft of density altitude per °C away from standard", expressed
/// per °F because the NWS constants above are Fahrenheit.
const DA_FAA_RULE_FT_PER_DEG_F: f64 = 120.0 * 5.0 / 9.0;
const DA_FT_TO_M: f64 = 0.3048;

/// NWS pressure altitude, feet, from STATION pressure in hPa.
///
/// The published NWS station-pressure equation
/// (<https://www.weather.gov/media/epz/wxcalc/pressureAltitude.pdf>):
///
/// ```text
/// PA_ft = 145366.45 * (1 - (P_hPa / 1013.25)^0.190284)
/// ```
///
/// Station pressure, never an altimeter setting (MBA-643): a QNH already corrects to sea
/// level and would report the pressure altitude of the sea rather than of the firing point.
pub fn nws_pressure_altitude_ft(station_pressure_hpa: f64) -> f64 {
    DA_NWS_PRESSURE_ALT_FT
        * (1.0 - (station_pressure_hpa / DA_NWS_SEA_LEVEL_HPA).powf(DA_NWS_EXPONENT))
}

/// Density altitude by the FAA rule of thumb, feet: NWS pressure altitude plus 120 ft for
/// every °C the station is away from ISA temperature at that pressure altitude. **Humidity is
/// not an input.**
///
/// ```text
/// DA_ft = PA_ft + (120 * 5/9) * (T_F - ISA_T_F(PA))      ISA_T_F(PA) = 59 - 3.57 * PA_ft/1000
/// ```
///
/// This is the density altitude the engine has always reported —
/// `pdf_dope_card::calculate_density_altitude` delegates to it — and the one
/// [`resolve_atmosphere_for_density_altitude`], the density-altitude ENTRY mode, inverts
/// exactly.
///
/// It stays humidity-free for compatibility, not out of necessity. Every card printed to date
/// shows this figure, and the entry mode that inverts it has shipped since MBA-1366 taking a
/// temperature but no humidity (the run's humidity is applied separately, afterwards). A
/// humidity-aware forward could be inverted just as well given the run's humidity, the way
/// temperature already is; it would change every printed card and what an entered density
/// altitude resolves to.
///
/// It is not the National Weather Service's density altitude either, despite the NWS pressure
/// altitude inside it: the NWS calculator uses virtual temperature, so it includes humidity
/// and tracks [`density_matched_altitude_m`] to within a few meters.
///
/// And it is an approximation. The rule is a straight line where the real dependence curves,
/// and 120 ft/°C is slightly steeper than the exact slope at sea level (118.6 ft/°C for an
/// ideal gas, 119.0 for the CIPM density [`density_matched_altitude_m`] uses), so
/// it drifts from the exact figure as the temperature moves away from standard FOR THE
/// PRESSURE ALTITUDE: about +22 m at 30 °C and 1013.25 hPa in completely dry air. See
/// [`density_matched_altitude_m`] for the exact, humidity-inclusive figure, and why the bridge
/// reports both.
pub fn faa_rule_density_altitude_ft(station_pressure_hpa: f64, temperature_f: f64) -> f64 {
    let pressure_alt = nws_pressure_altitude_ft(station_pressure_hpa);

    // ISA temperature at pressure altitude (lapse rate: 3.57°F per 1000 ft)
    let isa_temp_f = DA_NWS_SEA_LEVEL_TEMP_F - (pressure_alt / 1000.0) * DA_NWS_LAPSE_F_PER_1000FT;

    pressure_alt + DA_FAA_RULE_FT_PER_DEG_F * (temperature_f - isa_temp_f)
}

/// ISA troposphere constants for [`density_matched_altitude_m`].
const ISA_SEA_LEVEL_TEMP_K: f64 = 288.15;
const ISA_LAPSE_K_PER_M: f64 = 0.0065;

/// The reference air that [`density_matched_altitude_m`] calls sea level: ICAO standard sea
/// level, perfectly dry, evaluated with the SAME CIPM-2007 model as the air being described.
const DENSITY_MATCH_REFERENCE_TEMP_C: f64 = 15.0;
const DENSITY_MATCH_REFERENCE_PRESSURE_HPA: f64 = 1013.25;

/// The altitude in the ICAO Standard Atmosphere whose standard density equals the ACTUAL
/// density of this air, humidity included — density altitude by its textbook definition.
///
/// ```text
/// ratio = rho_CIPM(T, P, RH) / rho_CIPM(15 C, 1013.25 hPa, 0%)
/// h_m   = (T0 / L) * (1 - ratio^(1/n))      T0 = 288.15 K, L = 0.0065 K/m,
///                                            n  = g0 / (R L) - 1  (= 4.25588)
/// ```
///
/// # Why this exists alongside [`faa_rule_density_altitude_ft`]
///
/// Two different quantities are both called "density altitude" in the field, and reporting
/// them under one name is how an external user came to ask why the engine read 0 m where
/// another tool read 36 m at ICAO standard with 50% humidity (2026-09). The FAA-rule figure
/// is humidity-free for compatibility (see its docs). This one answers the question a shooter
/// usually means — "which ISA altitude has the air density I actually have" — and so it must
/// include humidity: water vapor (18 g/mol) is lighter than dry air (~29 g/mol), so humid air
/// is thinner.
///
/// At 15 °C, 1013.25 hPa and 50% RH this reads 33.1 m where the FAA-rule figure reads 0. But
/// the difference is NOT purely humidity, and callers must not present it that way: the FAA
/// rule also drifts from this exact figure whenever the temperature is not standard for the
/// pressure altitude, so the two disagree by about 22 m at 30 °C in perfectly dry air.
///
/// # Anchoring
///
/// Density comes from [`calculate_atmosphere`] — the one CIPM-2007 path every solver uses —
/// and the ratio is taken against that SAME model's density for ICAO standard dry sea level,
/// not against the textbook 1.225 kg/m³. The model is accurate, but it is not built to
/// reproduce 1.225 exactly (its own anchor test allows ±0.002), and at sea level 0.002 kg/m³
/// is roughly 17 m of altitude. Anchoring to the model makes ICAO standard dry air read
/// exactly 0 and makes the humidity effect the model's own.
///
/// # Scope
///
/// The closed-form troposphere, with no geopotential correction — the same plain-height
/// convention as the FAA-rule figure and the DA entry mode. That covers every shooting
/// condition on Earth, but the arithmetic does not know where it stops: above 11 km it keeps
/// answering, wrongly. Callers bound it; the bridge command refuses conditions outside
/// -5 km..11 km (`crate::atmosphere_service`).
pub fn density_matched_altitude_m(
    temperature_c: f64,
    station_pressure_hpa: f64,
    humidity_percent: f64,
) -> f64 {
    let (density, _) = calculate_atmosphere(
        0.0,
        Some(temperature_c),
        Some(station_pressure_hpa),
        humidity_percent,
    );
    let (reference, _) = calculate_atmosphere(
        0.0,
        Some(DENSITY_MATCH_REFERENCE_TEMP_C),
        Some(DENSITY_MATCH_REFERENCE_PRESSURE_HPA),
        0.0,
    );
    let exponent = G_ACCEL_MPS2 / (R_AIR * ISA_LAPSE_K_PER_M) - 1.0;
    (ISA_SEA_LEVEL_TEMP_K / ISA_LAPSE_K_PER_M) * (1.0 - (density / reference).powf(1.0 / exponent))
}

/// 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 FAA-rule density-altitude model this engine
/// already uses to REPORT density altitude ([`faa_rule_density_altitude_ft`], which
/// `pdf_dope_card::calculate_density_altitude` delegates to):
///
/// ```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. The FAA-rule density altitude is
///   honored either way, but the DENSITY is not: the rule is a straight-line approximation, not
///   a curve of equal density, so the ISA-default case rebuilds different air — at 30 °C and
///   1013.25 hPa, a rebuilt 11.4 °C / 949 hPa, 0.3% denser, with a speed of sound 3% lower. Only
///   the real temperature reproduces the real air.
///
/// # Height convention
/// The FAA-rule formula being inverted ([`faa_rule_density_altitude_ft`]) performs no
/// geopotential correction — it treats its altitude as a plain height (the DOPE card, which
/// delegates to it, leaves its own `_altitude_ft` parameter unused 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 = DA_FAA_RULE_FT_PER_DEG_F;

    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);

    // 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. The cap is silent:
    // air whose vapor pressure would exceed its total pressure does not exist, and callers
    // that must refuse it rather than compute pure steam test
    // `water_vapor_mole_fraction_uncapped` first (the bridge's `atmosphere.density_altitude`).
    let x_v = vapor_mole_fraction(p_sv, f, humidity_percent, 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))
}

/// CIPM-2007 mole fraction of water vapor, `x_v = RH * f * p_sv / p`, UNCAPPED: the one
/// formula [`calculate_air_density_cimp`] uses (it caps the result at 1). A value at or above 1
/// means the stated humidity would need more vapor pressure than the air has total pressure.
pub(crate) fn water_vapor_mole_fraction_uncapped(
    temp_c: f64,
    pressure_hpa: f64,
    humidity_percent: f64,
) -> f64 {
    let t_k = temp_c + 273.15;
    let p_sv = enhanced_saturation_vapor_pressure(t_k);
    let pressure_pa = pressure_hpa * 100.0;
    let f = enhanced_enhancement_factor(pressure_pa, temp_c);
    vapor_mole_fraction(
        p_sv,
        f,
        humidity_percent,
        pressure_pa.max(f64::MIN_POSITIVE),
    )
}

#[inline(always)]
fn vapor_mole_fraction(p_sv_hpa: f64, f: f64, humidity_percent: f64, p_pa: f64) -> f64 {
    // 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_hpa;

    // 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.2216 — moist air is
    // LIGHTER than dry air).
    let p_v_pa = p_v * 100.0;
    p_v_pa / p_pa
}

/// 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) {
    get_local_atmosphere_humid_with_gravity(
        altitude_m,
        base_alt,
        base_temp_c,
        base_press_hpa,
        base_ratio,
        humidity_percent,
        Gravity::EARTH,
    )
}

#[allow(clippy::too_many_arguments)]
pub(crate) fn get_local_atmosphere_humid_with_gravity(
    altitude_m: f64,
    base_alt: f64,
    base_temp_c: f64,
    base_press_hpa: f64,
    base_ratio: f64,
    humidity_percent: f64,
    gravity: Gravity,
) -> (f64, f64) {
    let (temp_k, pressure_pa, density) = local_temp_pressure_density_with_gravity(
        altitude_m,
        base_alt,
        base_temp_c,
        base_press_hpa,
        base_ratio,
        gravity,
    );
    (
        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) {
    local_temp_pressure_density_with_gravity(
        altitude_m,
        base_alt,
        base_temp_c,
        base_press_hpa,
        base_ratio,
        Gravity::EARTH,
    )
}

fn local_temp_pressure_density_with_gravity(
    altitude_m: f64,
    base_alt: f64,
    base_temp_c: f64,
    base_press_hpa: f64,
    base_ratio: f64,
    gravity: Gravity,
) -> (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,
        gravity,
    );

    // 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,
    gravity: Gravity,
) -> (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,
                gravity,
            );
            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,
                gravity,
            );
            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,
    gravity: Gravity,
) -> (f64, f64) {
    let temp_k = base_temp_k + lapse_rate * height_diff;
    let pressure_hpa = if lapse_rate.abs() < 1e-10 {
        base_pressure_hpa * (-gravity.magnitude() * height_diff / (R_AIR * base_temp_k)).exp()
    } else {
        let temp_ratio = temp_k / base_temp_k;
        base_pressure_hpa * temp_ratio.powf(-gravity.magnitude() / (lapse_rate * R_AIR))
    };

    (temp_k, pressure_hpa)
}

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

    #[test]
    fn uniform_gravity_controls_hydrostatic_pressure_without_moving_the_station_anchor() {
        for magnitude in [0.0, 1.62, 9.80665, 20.0] {
            let gravity = Gravity::new(magnitude).unwrap();
            let (temperature, pressure) =
                integrate_local_atmosphere_segment(280.0, 900.0, 100.0, 0.0, gravity);
            assert_eq!(temperature, 280.0);
            assert!(
                (pressure / 900.0 - (-magnitude * 100.0 / (R_AIR * 280.0)).exp()).abs() < 1e-14
            );
            let station =
                local_temp_pressure_density_with_gravity(300.0, 300.0, 20.0, 900.0, 0.9, gravity);
            assert_eq!(
                station,
                local_temp_pressure_density(300.0, 300.0, 20.0, 900.0, 0.9)
            );
        }
        // Exercise both the lapse and isothermal branches while crossing layer boundaries.
        let (_, pressure, _) = local_temp_pressure_density_with_gravity(
            25_000.0,
            0.0,
            20.0,
            900.0,
            0.9,
            Gravity::new(0.0).unwrap(),
        );
        assert_eq!(pressure, 90_000.0);
    }
}

/// 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);
    }

    /// The band's edges are inside it, in both directions. At sea level the ICAO standard is
    /// exactly 101325 Pa, and 1.08 and 0.85 times it are exact in binary, so the division
    /// lands on the constants themselves rather than next to them.
    #[test]
    fn station_pressure_plausibility_band_edges_are_inside_the_band() {
        let at_sea_level = |pa| station_pressure_plausibility(pa, 0.0);
        assert_eq!(at_sea_level(101_325.0).icao_standard_pa, 101_325.0);

        let high_edge = at_sea_level(109_431.0);
        assert_eq!(high_edge.ratio, STATION_PRESSURE_MAX_PLAUSIBLE_RATIO);
        assert!(!high_edge.implausible, "{high_edge:?}");
        assert!(at_sea_level(109_432.0).implausible);

        let low_edge = at_sea_level(86_126.25);
        assert_eq!(low_edge.ratio, STATION_PRESSURE_MIN_PLAUSIBLE_RATIO);
        assert!(!low_edge.implausible, "{low_edge:?}");
        assert!(at_sea_level(86_126.0).implausible);

        // Away from sea level the edges hold to a part in a billion either side.
        for altitude_m in [-400.0, 1500.0, 3310.128] {
            let (_, icao_pa) = calculate_icao_standard_atmosphere(altitude_m);
            for (ratio, implausible) in [
                (STATION_PRESSURE_MAX_PLAUSIBLE_RATIO * (1.0 - 1e-9), false),
                (STATION_PRESSURE_MAX_PLAUSIBLE_RATIO * (1.0 + 1e-9), true),
                (STATION_PRESSURE_MIN_PLAUSIBLE_RATIO * (1.0 + 1e-9), false),
                (STATION_PRESSURE_MIN_PLAUSIBLE_RATIO * (1.0 - 1e-9), true),
            ] {
                let p = station_pressure_plausibility(ratio * icao_pa, altitude_m);
                assert_eq!(p.icao_standard_pa, icao_pa);
                assert_eq!(
                    p.implausible, implausible,
                    "{altitude_m} m, ratio {ratio}: {p:?}"
                );
            }
        }

        assert!(station_pressure_plausibility(f64::NAN, 0.0).implausible);
    }

    /// The record sea-level pressures the band is justified by stay inside it, and the reported
    /// placeholder at Cusco does not.
    #[test]
    fn station_pressure_plausibility_passes_record_weather_and_flags_a_placeholder() {
        for record_hpa in [870.0, 1083.8, 1089.1] {
            let p = station_pressure_plausibility(record_hpa * 100.0, 0.0);
            assert!(!p.implausible, "{record_hpa} hPa: {p:?}");
        }

        let placeholder = station_pressure_plausibility(101_325.0, 3310.128);
        assert!(placeholder.implausible, "{placeholder:?}");
        assert!((placeholder.ratio - 1.503).abs() < 5e-4, "{placeholder:?}");
        // A real Cusco station reading is close to standard.
        assert!(!station_pressure_plausibility(68_830.0, 3310.128).implausible);
        // The same reading at an omitted (0 m) altitude is far below standard.
        assert!(station_pressure_plausibility(68_830.0, 0.0).implausible);
        // 101325 Pa leaves the band between 640 m and 650 m.
        assert!(!station_pressure_plausibility(101_325.0, 640.0).implausible);
        assert!(station_pressure_plausibility(101_325.0, 650.0).implausible);
    }

    #[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 altimeter-setting relation.
        // NOTE the geopotential step — 1500 m geometric is 1499.646 m geopotential, and it is
        // the geopotential height that enters the formula:
        //   h_geo   = 6356766*1500/(6356766+1500)                    = 1499.6461... m
        //   n       = 1/5.255875601466713                            = 0.19026325503612...
        //   station = 1030.0 * (1 - (0.0065*h_geo/288.15)*(1013.25/1030.0)^n)^(1/n)
        //           = 860.0681114332783 hPa
        // The ratio form this replaced, 1030.0 * (1 - 0.0065*h_geo/288.15)^5.2558756, gave
        // 859.5753123926447 hPa: 0.49 hPa low, because it scales the standard pressure ratio
        // by QNH instead of adding the station height to the QNH's pressure altitude.
        let reduced = reduce_qnh_to_station_pressure(1030.0, 1500.0);
        assert!(
            (reduced - 860.068_111_433_278_3).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
            );
        }
    }

    /// The doc comment's comparison with the NWS/ASOS altimeter-setting algorithm, checked:
    /// with that algorithm's constants (n = 0.190284, k = 8.4228807e-5, the field elevation as
    /// entered) and without its 0.3 hPa sensor offset, its station pressure is about 0.1 hPa
    /// below this function's at 1500 m and 0.3 to 0.45 hPa below at 4000 to 5000 m, over QNH
    /// 950 to 1050 hPa.
    #[test]
    fn qnh_reduction_differs_from_nws_asos_by_the_documented_amounts() {
        let nws_station_hpa = |qnh_hpa: f64, elevation_m: f64| {
            (qnh_hpa.powf(0.190284) - 8.4228807e-5 * elevation_m).powf(1.0 / 0.190284)
        };
        for (altitude_m, low, high) in [
            (1500.0, 0.09, 0.11),
            (4000.0, 0.3, 0.45),
            (5000.0, 0.3, 0.45),
        ] {
            for qnh_hpa in [950.0, 980.0, 1013.25, 1030.0, 1050.0] {
                let gap = reduce_qnh_to_station_pressure(qnh_hpa, altitude_m)
                    - nws_station_hpa(qnh_hpa, altitude_m);
                assert!(
                    (low..=high).contains(&gap),
                    "QNH {qnh_hpa} at {altitude_m} m: this function is {gap} hPa above NWS/ASOS"
                );
            }
        }
    }

    #[test]
    fn qnh_reduction_adds_the_station_height_to_the_qnhs_pressure_altitude() {
        // What an altimeter setting means: an altimeter set to QNH reads field elevation on
        // the ground, so the station pressure's ISA pressure altitude is the QNH's own pressure
        // altitude plus the station's (geopotential) height. The troposphere's pressure
        // altitude is H(p) = (T0/-L) * (1 - (p/P0)^n), n = -L*R/g.
        let layer = &ICAO_LAYERS[0];
        let n = -layer.lapse_rate * R_AIR / G_ACCEL_MPS2;
        let p0_hpa = layer.base_pressure / 100.0;
        let pressure_altitude_m = |pressure_hpa: f64| {
            layer.base_temperature / -layer.lapse_rate * (1.0 - (pressure_hpa / p0_hpa).powf(n))
        };
        // The unit-test pin, Cusco (SPZO), La Paz (SLLP), low and high settings, a standard
        // setting, and a station below sea level.
        for (qnh_hpa, altitude_m) in [
            (1030.0, 1500.0),
            (1033.0, 3310.0),
            (1038.0, 4061.0),
            (1000.0, 2000.0),
            (980.0, 1500.0),
            (1040.0, 4500.0),
            (1013.25, 3000.0),
            (1030.0, -400.0),
        ] {
            let station_hpa = reduce_qnh_to_station_pressure(qnh_hpa, altitude_m);
            let added_m = pressure_altitude_m(station_hpa) - pressure_altitude_m(qnh_hpa);
            let height_m = geometric_to_geopotential_height_m(altitude_m);
            assert!(
                (added_m - height_m).abs() < 1e-6,
                "QNH {qnh_hpa} hPa at {altitude_m} m: station {station_hpa} hPa adds \
                 {added_m} m of pressure altitude, expected {height_m} m"
            );
        }
    }

    #[test]
    fn qnh_reduction_is_nan_where_no_station_pressure_exists() {
        // A station pressure exists only while the QNH's pressure altitude plus the station
        // height stays below T0/-L (44,331 m, where the troposphere's lapse rate reaches
        // absolute zero). Past that the result is NaN, which every caller's finite-pressure
        // check refuses, rather than a near-vacuum number that would be flown as real air.
        assert!(reduce_qnh_to_station_pressure(500.0, 40_000.0).is_nan());
        assert!(reduce_qnh_to_station_pressure(-1030.0, 1500.0).is_nan());
        // Everywhere real shooting happens the reduction is finite and positive, lowers the
        // pressure above sea level and raises it below.
        for qnh_step in 0..=24 {
            let qnh_hpa = 850.0 + 10.0 * f64::from(qnh_step);
            for altitude_step in 0..=20 {
                let altitude_m = -1000.0 + 500.0 * f64::from(altitude_step);
                let station_hpa = reduce_qnh_to_station_pressure(qnh_hpa, altitude_m);
                assert!(
                    station_hpa.is_finite() && station_hpa > 0.0,
                    "QNH {qnh_hpa} hPa at {altitude_m} m gave {station_hpa} hPa"
                );
                assert_eq!(
                    station_hpa.partial_cmp(&qnh_hpa),
                    0.0.partial_cmp(&altitude_m),
                    "QNH {qnh_hpa} hPa at {altitude_m} m gave {station_hpa} hPa"
                );
            }
        }
    }

    #[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 - 860.068_111_433_278_3).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}"
            );
        }
    }
}

/// Station-presence resolution (0.49.0): a supplied station reading is never read as omitted.
#[cfg(test)]
mod station_presence_tests {
    use super::*;

    /// The body of `resolve_station_conditions_with_pressure_mode` exactly as it stood at
    /// ab4eb0c, before presence existed. It calls only the primitives this change leaves
    /// alone, so comparing against it cannot pass by comparing the new code with itself.
    fn legacy_resolve_ab4eb0c(
        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)
    }

    const TEMPERATURES_C: [f64; 10] = [
        -40.0, 14.89, 14.9, 14.95, 15.0, 15.05, 15.0999, 15.1, 15.2, 30.0,
    ];
    const PRESSURES_HPA: [f64; 11] = [
        1012.7,
        1012.75,
        1012.8562053925472,
        1013.0,
        1013.2079,
        1013.25,
        1013.5,
        1013.75,
        1014.0,
        850.0,
        1031.0,
    ];
    const ALTITUDES_M: [f64; 13] = [
        -2.0, -1.0, -0.5, 0.0, 0.5, 1.0, 1.0000001, 3.0, 10.0, 150.0, 1500.0, 1524.0, 5000.0,
    ];
    const MODES: [PressureReferenceMode; 2] =
        [PressureReferenceMode::Absolute, PressureReferenceMode::Qnh];

    fn assert_bits(actual: (f64, f64), expected: (f64, f64), context: &str) {
        assert_eq!(
            actual.0.to_bits(),
            expected.0.to_bits(),
            "temperature {context}: {actual:?} vs {expected:?}"
        );
        assert_eq!(
            actual.1.to_bits(),
            expected.1.to_bits(),
            "pressure {context}: {actual:?} vs {expected:?}"
        );
    }

    #[test]
    fn neither_is_bit_identical_to_ab4eb0c() {
        for t in TEMPERATURES_C {
            for p in PRESSURES_HPA {
                for alt in ALTITUDES_M {
                    for mode in MODES {
                        let expected = legacy_resolve_ab4eb0c(t, p, alt, mode);
                        let context = format!("({t}, {p}, {alt}, {mode:?})");
                        assert_bits(
                            resolve_station_conditions_with_presence(
                                t,
                                p,
                                alt,
                                mode,
                                StationAtmospherePresence::NEITHER,
                            ),
                            expected,
                            &context,
                        );
                        assert_bits(
                            resolve_station_conditions_with_pressure_mode(t, p, alt, mode),
                            expected,
                            &context,
                        );
                    }
                }
            }
        }

        // Six points captured at ab4eb0c, so a change to a shared primitive cannot pass
        // unseen through the verbatim copy above.
        let pinned = [
            (
                (15.05, 1013.0, 150.0, PressureReferenceMode::Absolute),
                // captured at ab4eb0c
                (
                    f64::from_bits(0x402c0ccfd0c49b20),
                    f64::from_bits(0x408f1ae13f434f1f),
                ),
            ),
            (
                (15.0, 1013.25, 1500.0, PressureReferenceMode::Absolute),
                // captured at ab4eb0c
                (
                    f64::from_bits(0x4015025af8d136c0),
                    f64::from_bits(0x408a6cc6358cefbd),
                ),
            ),
            (
                (15.0, 1013.25, 10.0, PressureReferenceMode::Qnh),
                // captured at ab4eb0c
                (
                    f64::from_bits(0x402ddeb85559de00),
                    f64::from_bits(0x408fa064e4e5f101),
                ),
            ),
            (
                (15.0, 1013.2079, 1524.0, PressureReferenceMode::Qnh),
                // captured at ab4eb0c
                (
                    f64::from_bits(0x401462aff47e8bc0),
                    f64::from_bits(0x408a58982cf0204f),
                ),
            ),
            (
                (
                    14.95,
                    1012.8562053925472,
                    1.0000001,
                    PressureReferenceMode::Absolute,
                ),
                // captured at ab4eb0c
                (
                    f64::from_bits(0x402dfcac08345a20),
                    f64::from_bits(0x408fa909fb979713),
                ),
            ),
            (
                (15.0, 1013.25, 3.0, PressureReferenceMode::Absolute),
                // captured at ab4eb0c
                (
                    f64::from_bits(0x402df60418e281e0),
                    f64::from_bits(0x408fa71e04fe2409),
                ),
            ),
        ];
        for ((t, p, alt, mode), expected) in pinned {
            let context = format!("pinned ({t}, {p}, {alt}, {mode:?})");
            assert_bits(
                resolve_station_conditions_with_presence(
                    t,
                    p,
                    alt,
                    mode,
                    StationAtmospherePresence::NEITHER,
                ),
                expected,
                &context,
            );
            assert_bits(
                resolve_station_conditions_with_pressure_mode(t, p, alt, mode),
                expected,
                &context,
            );
        }
    }

    #[test]
    fn both_returns_the_supplied_pair() {
        for t in TEMPERATURES_C {
            for p in PRESSURES_HPA {
                for alt in ALTITUDES_M {
                    let context = format!("({t}, {p}, {alt})");
                    assert_bits(
                        resolve_station_conditions_with_presence(
                            t,
                            p,
                            alt,
                            PressureReferenceMode::Absolute,
                            StationAtmospherePresence::BOTH,
                        ),
                        (t, p),
                        &context,
                    );
                    assert_bits(
                        resolve_station_conditions_with_presence(
                            t,
                            p,
                            alt,
                            PressureReferenceMode::Qnh,
                            StationAtmospherePresence::BOTH,
                        ),
                        (t, reduce_qnh_to_station_pressure(p, alt)),
                        &context,
                    );
                }
            }
        }
        // The headline repro: QNH 1031 hPa at 150 m reduces into the band and still flies.
        let (_, p) = resolve_station_conditions_with_presence(
            20.0,
            1031.0,
            150.0,
            PressureReferenceMode::Qnh,
            StationAtmospherePresence::BOTH,
        );
        assert_eq!(p.to_bits(), 1012.8562053925472_f64.to_bits(), "{p:?}");
    }

    #[test]
    fn presence_is_per_field() {
        let (icao_t_k, icao_p_pa) = calculate_icao_standard_atmosphere(1500.0);
        let temperature_only = resolve_station_conditions_with_presence(
            15.0,
            1013.25,
            1500.0,
            PressureReferenceMode::Absolute,
            StationAtmospherePresence::new(true, false),
        );
        assert_bits(temperature_only, (15.0, icao_p_pa / 100.0), "(true, false)");
        let pressure_only = resolve_station_conditions_with_presence(
            15.0,
            1013.0,
            1500.0,
            PressureReferenceMode::Absolute,
            StationAtmospherePresence::new(false, true),
        );
        assert_bits(pressure_only, (icao_t_k - 273.15, 1013.0), "(false, true)");
    }

    #[test]
    fn presence_constructors_agree() {
        assert_eq!(
            StationAtmospherePresence::default(),
            StationAtmospherePresence::NEITHER
        );
        assert_eq!(
            StationAtmospherePresence::NEITHER,
            StationAtmospherePresence::new(false, false)
        );
        assert_eq!(
            StationAtmospherePresence::BOTH,
            StationAtmospherePresence::new(true, true)
        );
        assert_eq!(
            StationAtmospherePresence::from_readings(
                StationReading::supplied(15.0),
                StationReading::omitted(1013.25)
            ),
            StationAtmospherePresence::new(true, false)
        );
        assert_eq!(
            StationAtmospherePresence::from_readings(
                StationReading::omitted(59.0),
                StationReading::supplied(29.92)
            ),
            StationAtmospherePresence::new(false, true)
        );
        assert!(StationReading::omitted(29.92).is_omitted());
        assert!(!StationReading::supplied(29.92).is_omitted());
        assert_eq!(
            StationReading::omitted(29.92).value.to_bits(),
            29.92_f64.to_bits()
        );
    }

    #[derive(Debug, serde::Serialize, serde::Deserialize)]
    #[serde(deny_unknown_fields)]
    struct Probe {
        #[serde(
            default = "probe_default",
            skip_serializing_if = "StationReading::is_omitted"
        )]
        pressure_inhg: StationReading,
    }

    fn probe_default() -> StationReading {
        StationReading::omitted(29.92)
    }

    #[test]
    fn station_reading_serde() {
        let absent: Probe = serde_json::from_str("{}").unwrap();
        assert!(!absent.pressure_inhg.supplied);
        assert_eq!(absent.pressure_inhg.value.to_bits(), 29.92_f64.to_bits());

        let present: Probe = serde_json::from_str(r#"{"pressure_inhg": 29.92}"#).unwrap();
        assert!(present.pressure_inhg.supplied);
        assert_eq!(present.pressure_inhg.value.to_bits(), 29.92_f64.to_bits());

        assert_eq!(serde_json::to_string(&absent).unwrap(), "{}");
        assert_eq!(
            serde_json::to_string(&present).unwrap(),
            r#"{"pressure_inhg":29.92}"#
        );

        assert!(serde_json::from_str::<Probe>(r#"{"pressure_inhg": null}"#).is_err());
        assert!(serde_json::from_str::<Probe>(r#"{"pressure_inhg": "29.92"}"#).is_err());
    }
}