astrodynamics-gnss 0.9.4

//! 0-ULP parity tests for the IONEX slant ionospheric delay pipeline.
//!
//! These assert the Rust port reproduces the canonical reference recipe
//! `parity/generator/ionex.py` bit-for-bit, using the committed golden fixture
//! `parity/fixtures/ionex_golden.json` and the synthetic IONEX product the
//! fixture was generated from. 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.
//!
//! Two things are checked. First the Rust IONEX parser is run on the committed
//! synthetic product and its grid (latitude/longitude node axes, every per-map
//! TEC value, and the J2000-second epoch axis) is asserted bit-for-bit against
//! the parser substrate recorded in the golden, so the byte/record reader is
//! pinned. Second the float pipeline is run per case and every intermediate
//! (pierce point, per-map bilinear VTEC, time-blended VTEC, STEC, meters) is
//! checked per component, so a divergence is localised to a single algorithm
//! step. The cases span the branch matrix: longitude wrap at both seams,
//! descending-latitude bracket, EXPONENT scaling, pierce-point latitude clamp at
//! the grid edge, epoch coincident with a map versus interior, endpoint hold
//! before the first and after the last map, and the L1/L2/L5 frequency scaling.

use std::path::PathBuf;

use serde_json::Value;

use super::grid::Ionex;
use super::slant::slant_delay_components;

/// 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 fixtures_dir() -> PathBuf {
    let crate_dir = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
    crate_dir
        .join("tests/fixtures")
        .canonicalize()
        .unwrap_or_else(|e| {
            panic!(
                "cannot locate tests/fixtures 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}")),
    )
}

/// Compare a parsed `f64` against a golden hex-float, recording any nonzero ULP.
fn check(failures: &mut Vec<String>, label: String, got: f64, want_hex: &str) {
    let want = parse_hex_float(want_hex);
    let ulp = ulp_distance(got, want);
    if ulp != 0 {
        failures.push(format!(
            "{label}: {ulp} ULP (rust={} ref={want_hex})",
            float_hex(got)
        ));
    }
}

#[test]
fn ionex_slant_zero_ulp_full_branch_matrix() {
    let fx = fixtures_dir();
    let golden_path = fx.join("ionex_golden.json");
    let raw = std::fs::read_to_string(&golden_path)
        .unwrap_or_else(|e| panic!("read {}: {e}", golden_path.display()));
    let doc: Value = serde_json::from_str(&raw).expect("parse ionex_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"
    );

    // Parse the committed synthetic IONEX product with the Rust parser.
    let file_meta = &doc["ionex_file"];
    let ionex_name = file_meta["name"].as_str().expect("ionex file name");
    let ionex_path = fx.join("ionex").join(ionex_name);
    let ionex_bytes =
        std::fs::read(&ionex_path).unwrap_or_else(|e| panic!("read {}: {e}", ionex_path.display()));
    let ionex = Ionex::parse(&ionex_bytes).expect("parse synthetic IONEX product");

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

    // ---- Parser parity: node axes, epochs, and every TEC value ----
    let lat_ref = file_meta["lat_arr"].as_array().expect("lat_arr");
    let lon_ref = file_meta["lon_arr"].as_array().expect("lon_arr");
    assert_eq!(
        ionex.lat_nodes_deg().len(),
        lat_ref.len(),
        "parsed latitude node count"
    );
    assert_eq!(
        ionex.lon_nodes_deg().len(),
        lon_ref.len(),
        "parsed longitude node count"
    );
    for (i, want) in lat_ref.iter().enumerate() {
        check(
            &mut failures,
            format!("lat_arr[{i}]"),
            ionex.lat_nodes_deg()[i],
            want.as_str().unwrap(),
        );
    }
    for (j, want) in lon_ref.iter().enumerate() {
        check(
            &mut failures,
            format!("lon_arr[{j}]"),
            ionex.lon_nodes_deg()[j],
            want.as_str().unwrap(),
        );
    }

    let exp_ref = file_meta["exponent"].as_f64().expect("exponent") as i32;
    assert_eq!(ionex.exponent(), exp_ref, "parsed EXPONENT");

    let epochs_ref = file_meta["map_epochs_s"].as_array().expect("map_epochs_s");
    assert_eq!(
        ionex.map_epochs_s().len(),
        epochs_ref.len(),
        "parsed map count"
    );
    for (m, want) in epochs_ref.iter().enumerate() {
        assert_eq!(
            ionex.map_epochs_s()[m],
            want.as_i64().expect("epoch int"),
            "parsed map epoch[{m}] (J2000 seconds)"
        );
    }

    let maps_ref = file_meta["maps_vtec"].as_array().expect("maps_vtec");
    assert_eq!(
        ionex.tec_maps().len(),
        maps_ref.len(),
        "parsed TEC map count"
    );
    for (m, map_ref) in maps_ref.iter().enumerate() {
        let rows = map_ref.as_array().unwrap();
        for (i, row_ref) in rows.iter().enumerate() {
            let cols = row_ref.as_array().unwrap();
            for (j, want) in cols.iter().enumerate() {
                check(
                    &mut failures,
                    format!("maps_vtec[{m}][{i}][{j}]"),
                    ionex.tec_maps()[m][i][j],
                    want.as_str().unwrap(),
                );
            }
        }
    }

    // ---- Pipeline parity: every intermediate, per case ----
    let cases = doc["cases"].as_array().expect("cases array");
    assert!(
        cases.len() >= 12,
        "expected the full branch matrix (>= 12 cases), found {}",
        cases.len()
    );

    let lat_arr = ionex.lat_nodes_deg();
    let lon_arr = ionex.lon_nodes_deg();
    let dlat = ionex.dlat_deg();
    let dlon = ionex.dlon_deg();
    let re = ionex.base_radius_km();
    let h = ionex.shell_height_km();
    let epochs = ionex.map_epochs_s();
    let maps = ionex.tec_maps();

    let mut checks = 0usize;

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

        let epoch_s = inp["epoch_s"].as_i64().expect("epoch_s int");

        let got = slant_delay_components(
            hexf(inp, "lat_rad"),
            hexf(inp, "lon_rad"),
            hexf(inp, "az_rad"),
            hexf(inp, "el_rad"),
            hexf(inp, "frequency_hz"),
            re,
            h,
            epoch_s,
            epochs,
            maps,
            lat_arr,
            lon_arr,
            dlat,
            dlon,
        );

        // The temporal-bracket index is a discrete branch outcome, not a float.
        assert_eq!(
            got.map_index as i64,
            case["map_index"].as_i64().expect("map_index"),
            "case {name}: temporal bracket index"
        );

        let components: &[(&str, f64)] = &[
            ("s", got.s),
            ("psi", got.psi),
            ("phi_ipp_deg", got.phi_ipp_deg),
            ("lambda_ipp_deg_raw", got.lambda_ipp_deg_raw),
            ("lambda_ipp_deg", got.lambda_ipp_deg),
            ("w", got.w),
            ("vtec0", got.vtec0),
            ("vtec1", got.vtec1),
            ("p0", got.p0),
            ("q0", got.q0),
            ("vtec", got.vtec),
            ("m", got.m),
            ("stec", got.stec),
            ("delay_m", got.delay_m),
        ];

        for &(c, value) in components {
            let want_hex = exp[c]
                .as_str()
                .unwrap_or_else(|| panic!("case {name}: missing expected component {c}"));
            check(&mut failures, format!("{name}.{c}"), value, want_hex);
            checks += 1;
        }
    }

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

/// Regression: a single-map IONEX product must not panic in the temporal
/// bracket. The multi-map path computes `nmaps - 2`, which underflows (usize)
/// for one map and then indexes a second, non-existent map. A one-map product
/// has no interval to interpolate, so it holds that single map; querying it must
/// return the same value the equivalent two-map product returns at its first
/// epoch (where the temporal weight is 0).
#[test]
fn ionex_single_map_does_not_panic_and_holds_the_map() {
    let fx = fixtures_dir();
    let two_map_path = fx.join("ionex/synthetic_2map_7x7.20i");
    let full = std::fs::read_to_string(&two_map_path)
        .unwrap_or_else(|e| panic!("read {}: {e}", two_map_path.display()));
    let lines: Vec<&str> = full.lines().collect();
    let hdr_end = lines
        .iter()
        .position(|l| l.contains("END OF HEADER"))
        .expect("END OF HEADER");
    let first_map_end = lines
        .iter()
        .position(|l| l.contains("END OF TEC MAP"))
        .expect("END OF TEC MAP");

    // Build a one-map product by reusing the real file's exact formatting:
    // header (with the maps-count digit forced to 1) + the first map block only.
    let mut single = String::new();
    for l in &lines[..=hdr_end] {
        if l.contains("# OF MAPS IN FILE") {
            single.push_str(&l.replacen('2', "1", 1));
        } else {
            single.push_str(l);
        }
        single.push('\n');
    }
    for l in &lines[(hdr_end + 1)..=first_map_end] {
        single.push_str(l);
        single.push('\n');
    }

    let one = super::Ionex::parse(single.as_bytes()).expect("parse single-map IONEX");
    assert_eq!(one.map_epochs_s().len(), 1, "expected exactly one map");
    let two = super::Ionex::parse(full.as_bytes()).expect("parse two-map IONEX");

    let receiver = crate::frame::Wgs84Geodetic::new(30.0_f64.to_radians(), 0.0, 0.0);
    let el = 45.0_f64.to_radians();
    let az = 90.0_f64.to_radians();
    let f_l1 = 1_575_420_000.0;
    let epoch0 = one.map_epochs_s()[0];

    // Must not panic, and must be a finite positive delay.
    let d_one = super::ionex_slant_delay(&one, receiver, el, az, epoch0, f_l1);
    assert!(
        d_one.is_finite() && d_one > 0.0,
        "single-map delay not finite/positive: {d_one}"
    );

    // At its first epoch the two-map product weights the first map only (w == 0),
    // so the single-map hold must reproduce it bit-for-bit.
    let d_two = super::ionex_slant_delay(&two, receiver, el, az, epoch0, f_l1);
    assert_eq!(
        d_one.to_bits(),
        d_two.to_bits(),
        "single-map delay {d_one} != two-map-at-first-epoch {d_two}"
    );
}

/// A degenerate grid with fewer than two nodes on an axis must be rejected at
/// parse time, not accepted and then panicked on at evaluation (bilinear
/// interpolation brackets a cell with `node[i+1]`). Build a one-latitude-node
/// product by forcing `LAT2 == LAT1` and keeping a single band, reusing the real
/// fixture's column layout otherwise.
#[test]
fn ionex_degenerate_single_node_axis_is_rejected_at_parse() {
    let fx = fixtures_dir();
    let path = fx.join("ionex/synthetic_2map_7x7.20i");
    let full =
        std::fs::read_to_string(&path).unwrap_or_else(|e| panic!("read {}: {e}", path.display()));
    let lines: Vec<&str> = full.lines().collect();
    let hdr_end = lines
        .iter()
        .position(|l| l.contains("END OF HEADER"))
        .unwrap();
    let map_start = hdr_end + 1; // START OF TEC MAP
    let map_end = lines
        .iter()
        .position(|l| l.contains("END OF TEC MAP"))
        .unwrap();

    let mut s = String::new();
    for l in &lines[..=hdr_end] {
        if l.contains("LAT1 / LAT2 / DLAT") {
            // 60.0 / -60.0 / -20.0 -> 60.0 / 60.0 / -20.0 == a single latitude node.
            s.push_str(&l.replacen("-60.0", " 60.0", 1));
        } else if l.contains("# OF MAPS IN FILE") {
            s.push_str(&l.replacen('2', "1", 1));
        } else {
            s.push_str(l);
        }
        s.push('\n');
    }
    // One map, first band only (latitude 60.0 with its seven longitude values).
    for l in &[
        lines[map_start],
        lines[map_start + 1],
        lines[map_start + 2],
        lines[map_start + 3],
        lines[map_end],
    ] {
        s.push_str(l);
        s.push('\n');
    }

    let parsed = super::Ionex::parse(s.as_bytes());
    assert!(
        parsed.is_err(),
        "degenerate single-node grid should be rejected"
    );
    let msg = format!("{:?}", parsed.err().unwrap()).to_lowercase();
    assert!(
        msg.contains("node"),
        "expected a node-count parse error, got: {msg}"
    );
}