astrodynamics-gnss 0.9.4

//! 0-ULP parity tests for the Saastamoinen + Niell tropospheric delay recipe.
//!
//! These assert the Rust port reproduces the canonical reference recipe
//! `parity/generator/troposphere.py` bit-for-bit, using the committed golden
//! fixture `parity/fixtures/troposphere_golden.json`. Values are serialised as
//! hex-float (Python `float.hex()`) so there is no decimal-parse ambiguity, and
//! parity is measured as ULP distance via the integer reinterpretation of the
//! IEEE-754 bit pattern, per the `skyfield_parity_test.exs` discipline.
//!
//! Every intermediate quantity (the water-vapour partial pressure, both zenith
//! delays, both mapping factors, the height correction, the seasonal phase) and
//! the final slant delay is checked per component, so a divergence is localised
//! to a single algorithm step rather than only seen at the output. The cases
//! span the full branch matrix: elevations from zenith down past the Niell low
//! range and the at/below-horizon gate; equator / mid / high / polar / southern
//! latitudes exercising the latitude-interpolation nodes, clamps, and the
//! southern seasonal-phase offset; day-of-year at the cosine-phase extremes and
//! a fractional mid-year value; sea level, high altitude, below sea level, the
//! Met-path height-gate edges; and humidity, temperature, and pressure extremes.
//! A separate standard-atmosphere family covers the pressure/temperature
//! synthesis and its height clamp.

use std::path::PathBuf;

use serde_json::Value;

use super::saastamoinen::{slant_components, standard_atmosphere};

/// Parse a C99 / Python `float.hex()` hex-float string into the exact `f64`.
///
/// Every hex frac digit is 4 mantissa bits and f64 has 52, so the reconstruction
/// of a 13-digit fraction is exact (no rounding).
fn parse_hex_float(s: &str) -> f64 {
    let s = s.trim();
    let (neg, rest) = match s.strip_prefix('-') {
        Some(r) => (true, r),
        None => (false, s),
    };
    let rest = rest
        .strip_prefix("0x")
        .or_else(|| rest.strip_prefix("0X"))
        .unwrap_or_else(|| panic!("not a hex float (missing 0x): {s:?}"));

    let (mantissa, exp_str) = rest
        .split_once(['p', 'P'])
        .unwrap_or_else(|| panic!("not a hex float (missing p exponent): {s:?}"));
    let exp2: i32 = exp_str
        .parse()
        .unwrap_or_else(|_| panic!("bad binary exponent in {s:?}"));

    let (int_part, frac_part) = match mantissa.split_once('.') {
        Some((i, f)) => (i, f),
        None => (mantissa, ""),
    };

    let int_val: f64 = i64::from_str_radix(int_part, 16)
        .unwrap_or_else(|_| panic!("bad integer hex digits in {s:?}"))
        as f64;

    let mut frac_val = 0.0f64;
    let mut scale = 1.0f64 / 16.0;
    for c in frac_part.chars() {
        let d = c
            .to_digit(16)
            .unwrap_or_else(|| panic!("bad hex frac digit {c:?} in {s:?}"));
        frac_val += (d as f64) * scale;
        scale /= 16.0;
    }

    let significand = int_val + frac_val;
    let val = significand * 2.0f64.powi(exp2);
    if neg {
        -val
    } else {
        val
    }
}

/// ULP distance between two `f64`, using the monotone signed-integer mapping of
/// the IEEE-754 bit pattern. Returns `u64::MAX` for any NaN so a NaN never
/// silently reads as 0 ULP.
fn ulp_distance(a: f64, b: f64) -> u64 {
    if a.is_nan() || b.is_nan() {
        return u64::MAX;
    }
    ordered_i64(a).abs_diff(ordered_i64(b))
}

/// Map an `f64` to a sign-magnitude-ordered `i64` so adjacent floats differ by 1.
fn ordered_i64(x: f64) -> i64 {
    let bits = x.to_bits() as i64;
    if bits < 0 {
        i64::MIN - bits
    } else {
        bits
    }
}

/// Render an `f64` as a Python-`float.hex()`-style string for diagnostics.
fn float_hex(x: f64) -> String {
    if x == 0.0 {
        return if x.is_sign_negative() {
            "-0x0.0p+0".into()
        } else {
            "0x0.0p+0".into()
        };
    }
    let bits = x.to_bits();
    let sign = if (bits >> 63) & 1 == 1 { "-" } else { "" };
    let exp = ((bits >> 52) & 0x7ff) as i64;
    let mantissa = bits & 0x000f_ffff_ffff_ffff;
    let unbiased = exp - 1023;
    if unbiased >= 0 {
        format!("{sign}0x1.{mantissa:013x}p+{unbiased}")
    } else {
        format!("{sign}0x1.{mantissa:013x}p{unbiased}")
    }
}

fn fixture_path() -> PathBuf {
    let crate_dir = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
    crate_dir
        .join("tests/fixtures/troposphere_golden.json")
        .canonicalize()
        .unwrap_or_else(|e| {
            panic!(
                "cannot locate tests/fixtures/troposphere_golden.json from {}: {e}",
                crate_dir.display()
            )
        })
}

fn hexf(v: &Value, key: &str) -> f64 {
    parse_hex_float(
        v[key]
            .as_str()
            .unwrap_or_else(|| panic!("missing/non-string {key}")),
    )
}

#[test]
fn troposphere_zero_ulp_full_branch_matrix() {
    let raw = std::fs::read_to_string(fixture_path()).expect("read troposphere_golden.json");
    let doc: Value = serde_json::from_str(&raw).expect("parse troposphere_golden.json");

    // Self-check the hex-float parser/serialiser round-trips a known bit pattern,
    // so a parser bug can never masquerade as parity.
    let probe = "0x1.921fb54442d18p+1"; // math.pi
    assert_eq!(
        float_hex(parse_hex_float(probe)),
        probe,
        "hex-float parser/serialiser round-trip is broken"
    );

    let mut failures: Vec<String> = Vec::new();
    let mut checks = 0usize;

    // --- Slant cases (zenith + mapping + composed delay). ---
    let slant_cases = doc["slant_cases"].as_array().expect("slant_cases array");
    assert!(
        slant_cases.len() >= 22,
        "expected the full slant branch matrix (>= 22 cases), found {}",
        slant_cases.len()
    );

    let slant_components_keys: &[&str] = &[
        "e", "zhd_m", "zwd_m", "mh", "mw", "dm", "cosy", "y", "slant_m",
    ];

    for case in slant_cases {
        let name = case["name"].as_str().unwrap_or("<unnamed>");
        let inp = &case["inputs"];
        let exp = &case["expect"];

        // Radian angles and the fractional day-of-year are taken straight from
        // the fixture, exactly as the reference recipe consumes them.
        let got = slant_components(
            hexf(inp, "el_rad"),
            hexf(inp, "lat_rad"),
            hexf(inp, "lon_rad"),
            hexf(inp, "height_m"),
            hexf(inp, "pressure_hpa"),
            hexf(inp, "temperature_k"),
            hexf(inp, "relative_humidity"),
            hexf(inp, "doy"),
        );

        let actual = |c: &str| -> f64 {
            match c {
                "e" => got.e,
                "zhd_m" => got.zhd_m,
                "zwd_m" => got.zwd_m,
                "mh" => got.mh,
                "mw" => got.mw,
                "dm" => got.dm,
                "cosy" => got.cosy,
                "y" => got.y,
                "slant_m" => got.slant_m,
                other => panic!("unknown component {other}"),
            }
        };

        for &c in slant_components_keys {
            let want = parse_hex_float(
                exp[c]
                    .as_str()
                    .unwrap_or_else(|| panic!("case {name}: missing expected component {c}")),
            );
            let a = actual(c);
            let ulp = ulp_distance(a, want);
            checks += 1;
            if ulp != 0 {
                failures.push(format!(
                    "slant {name}.{c}: {ulp} ULP (rust={} ref={})",
                    float_hex(a),
                    exp[c].as_str().unwrap()
                ));
            }
        }
    }

    // --- Standard-atmosphere helper family (P/T synthesis + height clamp). ---
    let std_cases = doc["standard_atmosphere_cases"]
        .as_array()
        .expect("standard_atmosphere_cases array");
    assert!(
        std_cases.len() >= 4,
        "expected the standard-atmosphere family (>= 4 cases), found {}",
        std_cases.len()
    );

    let std_keys: &[&str] = &["pressure_hpa", "temperature_k", "relative_humidity"];

    for case in std_cases {
        let name = case["name"].as_str().unwrap_or("<unnamed>");
        let inp = &case["inputs"];
        let exp = &case["expect"];

        let got = standard_atmosphere(hexf(inp, "height_m"), hexf(inp, "relative_humidity"));

        let actual = |c: &str| -> f64 {
            match c {
                "pressure_hpa" => got.pressure_hpa,
                "temperature_k" => got.temperature_k,
                "relative_humidity" => got.relative_humidity,
                other => panic!("unknown component {other}"),
            }
        };

        for &c in std_keys {
            let want = parse_hex_float(
                exp[c]
                    .as_str()
                    .unwrap_or_else(|| panic!("case {name}: missing expected component {c}")),
            );
            let a = actual(c);
            let ulp = ulp_distance(a, want);
            checks += 1;
            if ulp != 0 {
                failures.push(format!(
                    "std {name}.{c}: {ulp} ULP (rust={} ref={})",
                    float_hex(a),
                    exp[c].as_str().unwrap()
                ));
            }
        }
    }

    assert!(checks > 0, "no components were checked - fixture empty?");
    assert!(
        failures.is_empty(),
        "Troposphere Rust port diverged from the reference recipe on {} of {checks} components:\n  {}",
        failures.len(),
        failures.join("\n  ")
    );
}

/// Coverage for the PUBLIC `tropo_slant(epoch, ...)` wrapper, which the kernel
/// test above does not exercise (it feeds `doy` straight into
/// `slant_components`). The wrapper derives the day-of-year from `epoch`. Unlike
/// the ionosphere path there is no angle conversion (latitude stays in radians),
/// so when the epoch reconstructs the day-of-year exactly the public path is
/// bit-for-bit identical to the kernel/golden. This drives the `zenith_midlat`
/// case (day-of-year 28 -> 2000-01-28 00:00, reconstructed exactly) through the
/// public API and asserts 0 ULP against the golden, closing the coverage gap.
#[test]
fn tropo_public_wrapper_is_0_ulp_when_epoch_reconstructs_doy() {
    use astrodynamics::time::model::{Instant, JulianDateSplit, TimeScale};

    let raw = std::fs::read_to_string(fixture_path()).expect("read troposphere_golden.json");
    let doc: Value = serde_json::from_str(&raw).expect("parse troposphere_golden.json");
    let cases = doc["slant_cases"].as_array().expect("slant_cases array");

    let case = cases
        .iter()
        .find(|c| c["name"].as_str() == Some("zenith_midlat"))
        .expect("zenith_midlat case present");
    let inp = &case["inputs"];

    // This reconstruction is only exact for an integer day-of-year at midnight;
    // assert that precondition so a fixture change can't silently invalidate it.
    assert_eq!(
        hexf(inp, "doy"),
        28.0,
        "test assumes zenith_midlat is day-of-year 28"
    );

    // 2000-01-28 00:00:00: JDN(noon) = 2451572, midnight boundary = 2451571.5.
    // fractional_day_of_year() maps this to exactly 28.0.
    let epoch = Instant::from_julian_date(TimeScale::Gpst, JulianDateSplit::new(2_451_571.5, 0.0));

    let receiver = crate::frame::Wgs84Geodetic::new(
        hexf(inp, "lat_rad"),
        hexf(inp, "lon_rad"),
        hexf(inp, "height_m"),
    );
    let met = super::Met::new(
        hexf(inp, "pressure_hpa"),
        hexf(inp, "temperature_k"),
        hexf(inp, "relative_humidity"),
    );

    let got = super::tropo_slant(hexf(inp, "el_rad"), receiver, met, epoch);
    let want = parse_hex_float(case["expect"]["slant_m"].as_str().unwrap());

    assert_eq!(
        got.to_bits(),
        want.to_bits(),
        "public tropo_slant not 0-ULP: got={} want={}",
        float_hex(got),
        case["expect"]["slant_m"].as_str().unwrap()
    );
}