astrodynamics-gnss 0.14.0

//! Two-track parity for the GPS L1 single-point-positioning pipeline.
//!
//! The reference recipe is `parity/generator/spp.py`; the committed fixtures
//! (vendored at `tests/fixtures/spp_trace_*.json`) record the synthesized
//! observations, the frozen-branch choices, the effective scipy options, and
//! the full iteration trace. Float values are serialized as the raw IEEE-754
//! bit pattern (`f64::to_bits`, a 16-hex-digit `0x...` literal) so there is no
//! decimal-parse ambiguity; parity is measured as ULP distance via the integer
//! reinterpretation of the bit pattern, per the existing SP3/DOP/atmosphere
//! parity discipline.
//!
//! Track 1 (0 ULP, libm/arithmetic-bound) is the trace-replay: for each
//! recorded state `x`, the per-satellite intermediates, the weighted residual
//! vector, and its 2-point finite-difference Jacobian are recomputed by the
//! Rust SPP substrate AT THAT x and asserted bit-for-bit. The Rust solver is
//! not run for this track. The ladder is built up by correction level (L0
//! geometry+clock+Sagnac, L1 +ionosphere, L2 +troposphere, L3 relativistic
//! no-op) so a miss localizes to the term added.
//!
//! Track 2 (sub-micron, BLAS-bound) is the independent-solve agreement: the
//! crate trust-region solver is run from the same inputs and the converged
//! position/clock is asserted to agree with both the recorded scipy solution
//! and the synthesized truth to a documented sub-micron bound. This is a
//! solver-agreement check, explicitly NOT a 0-ULP physics claim, because the
//! trust-region linear-algebra step is not bit-reproducible across BLAS builds.

use std::path::PathBuf;

use serde_json::Value;

use super::test_support;
use super::{
    solve, Corrections, KlobucharCoeffs, Observation, RejectionReason, SolveInputs, SppError,
    SurfaceMet,
};
use crate::id::{GnssSatelliteId, GnssSystem};
use astrodynamics::math::least_squares::{jacobian_2point, FD_REL_STEP_2POINT};
use nalgebra::DVector;

// ---------------------------------------------------------------------------
// Hex (f64::to_bits) helpers, ULP distance, NaN guard, parser self-check.
// ---------------------------------------------------------------------------

/// Parse a `0x...` raw-bits literal (Python `hexf` / Rust `f64::to_bits`) to f64.
fn bits(s: &str) -> f64 {
    let s = s.trim();
    let hex = s
        .strip_prefix("0x")
        .or_else(|| s.strip_prefix("0X"))
        .unwrap_or_else(|| panic!("not a 0x bits literal: {s:?}"));
    let u = u64::from_str_radix(hex, 16).unwrap_or_else(|_| panic!("bad hex bits in {s:?}"));
    f64::from_bits(u)
}

/// ULP distance between two f64; NaN on either side reads as `u64::MAX` so it can
/// never masquerade as 0 ULP.
fn ulp_distance(a: f64, b: f64) -> u64 {
    if a.is_nan() || b.is_nan() {
        return u64::MAX;
    }
    ordered(a).abs_diff(ordered(b))
}

fn ordered(x: f64) -> i64 {
    let b = x.to_bits() as i64;
    if b < 0 {
        i64::MIN - b
    } else {
        b
    }
}

fn fixture_path(name: &str) -> PathBuf {
    PathBuf::from(env!("CARGO_MANIFEST_DIR"))
        .join("tests/fixtures")
        .join(name)
}

fn read_fixture(name: &str) -> Value {
    let raw =
        std::fs::read_to_string(fixture_path(name)).unwrap_or_else(|e| panic!("read {name}: {e}"));
    serde_json::from_str(&raw).unwrap_or_else(|e| panic!("parse {name}: {e}"))
}

fn sp3() -> crate::sp3::Sp3 {
    let path = concat!(
        env!("CARGO_MANIFEST_DIR"),
        "/tests/fixtures/sp3/GRG0MGXFIN_20201760000_01D_15M_ORB.SP3"
    );
    let bytes = std::fs::read(path).unwrap_or_else(|e| panic!("read SP3 fixture {path}: {e}"));
    crate::sp3::Sp3::parse(&bytes).expect("parse real IGS SP3")
}

fn parse_prn(token: &str) -> GnssSatelliteId {
    let sys = GnssSystem::from_letter(token.chars().next().unwrap()).expect("known system letter");
    let prn: u8 = token[1..].parse().expect("prn digits");
    GnssSatelliteId::new(sys, prn)
}

fn arr3(v: &Value) -> [f64; 3] {
    let a = v.as_array().expect("array");
    [
        bits(a[0].as_str().unwrap()),
        bits(a[1].as_str().unwrap()),
        bits(a[2].as_str().unwrap()),
    ]
}

/// Read the level-independent inputs shared by every solve from a fixture.
struct Inputs {
    observations: Vec<Observation>,
    t_rx_j2000_s: f64,
    sod_s: f64,
    doy: f64,
    x0: [f64; 4],
    klobuchar: KlobucharCoeffs,
    met: SurfaceMet,
    corrections: Corrections,
}

fn corrections_for(level: &str) -> Corrections {
    match level {
        "L0_minimal" => Corrections::NONE,
        "L1_iono" => Corrections::IONO,
        "L2_tropo" | "L3_relativistic" => Corrections::IONO_TROPO,
        other => panic!("unknown level {other}"),
    }
}

fn load_inputs(doc: &Value, level: &str) -> Inputs {
    let f = &doc["fixture"];
    let inp = &f["inputs"];

    let observations = inp["observations"]
        .as_array()
        .expect("observations array")
        .iter()
        .map(|o| Observation {
            satellite_id: parse_prn(o["sat_id"].as_str().unwrap()),
            pseudorange_m: bits(o["p_meas_m"].as_str().unwrap()),
        })
        .collect();

    let alpha_v = inp["klobuchar_alpha"].as_array().unwrap();
    let beta_v = inp["klobuchar_beta"].as_array().unwrap();
    let klobuchar = KlobucharCoeffs {
        alpha: [
            bits(alpha_v[0].as_str().unwrap()),
            bits(alpha_v[1].as_str().unwrap()),
            bits(alpha_v[2].as_str().unwrap()),
            bits(alpha_v[3].as_str().unwrap()),
        ],
        beta: [
            bits(beta_v[0].as_str().unwrap()),
            bits(beta_v[1].as_str().unwrap()),
            bits(beta_v[2].as_str().unwrap()),
            bits(beta_v[3].as_str().unwrap()),
        ],
    };

    let met = SurfaceMet {
        pressure_hpa: bits(inp["met"]["pressure_hpa"].as_str().unwrap()),
        temperature_k: bits(inp["met"]["temperature_k"].as_str().unwrap()),
        relative_humidity: bits(inp["met"]["relative_humidity"].as_str().unwrap()),
    };

    let x0v = f["frozen"]["initial_guess_x0"].as_array().unwrap();
    let x0 = [
        bits(x0v[0].as_str().unwrap()),
        bits(x0v[1].as_str().unwrap()),
        bits(x0v[2].as_str().unwrap()),
        bits(x0v[3].as_str().unwrap()),
    ];

    Inputs {
        observations,
        t_rx_j2000_s: bits(inp["t_rx_j2000_s"].as_str().unwrap()),
        sod_s: bits(inp["t_rx_sod_s"].as_str().unwrap()),
        doy: bits(inp["doy"].as_str().unwrap()),
        x0,
        klobuchar,
        met,
        corrections: corrections_for(level),
    }
}

fn solve_inputs(i: &Inputs) -> SolveInputs {
    SolveInputs {
        observations: i.observations.clone(),
        t_rx_j2000_s: i.t_rx_j2000_s,
        t_rx_second_of_day_s: i.sod_s,
        day_of_year: i.doy,
        initial_guess: i.x0,
        corrections: i.corrections,
        klobuchar: i.klobuchar,
        beidou_klobuchar: None,
        met: i.met,
    }
}

const LEVELS: &[&str] = &["L0_minimal", "L1_iono", "L2_tropo", "L3_relativistic"];

fn fixture_name(level: &str) -> String {
    format!("spp_trace_{level}.json")
}

// ---------------------------------------------------------------------------
// Parser self-check: a known bit pattern round-trips. A parser bug must not be
// able to masquerade as parity.
// ---------------------------------------------------------------------------
#[test]
fn hex_bits_parser_round_trips() {
    let pi = std::f64::consts::PI;
    let hexed = format!("0x{:016x}", pi.to_bits());
    assert_eq!(
        bits(&hexed).to_bits(),
        pi.to_bits(),
        "bits parser round-trip broken"
    );
    // ULP distance of a value to itself is zero; to its neighbour is one.
    assert_eq!(ulp_distance(pi, pi), 0);
    let nxt = f64::from_bits(pi.to_bits() + 1);
    assert_eq!(ulp_distance(pi, nxt), 1);
    assert_eq!(ulp_distance(f64::NAN, 1.0), u64::MAX);
}

// ---------------------------------------------------------------------------
// Boundary cross-check: the meters/radians-native geodetic helper vs the core
// km/deg itrs_to_geodetic_compute. The 0-ULP claim is on the meters-native
// radians vs the recipe's recorded meters_native radians (both replicate the
// same Skyfield AU-internal algorithm); the core public deg API is checked
// against the recipe's recorded core_api_deg.
// ---------------------------------------------------------------------------
#[test]
fn geodetic_meters_native_zero_ulp() {
    let doc = read_fixture("spp_trace_L0_minimal.json");
    let cases = doc["geodetic_crosscheck"]["cases"]
        .as_array()
        .expect("cases");
    assert!(cases.len() >= 3, "expected >= 3 cross-check cases");

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

    for case in cases {
        let km = arr3(&case["ecef_km"]);
        let g = test_support::geodetic_from_ecef_m_for_test(
            km[0] * 1000.0,
            km[1] * 1000.0,
            km[2] * 1000.0,
        );

        let mn = &case["meters_native"];
        let want_lat = bits(mn["lat_rad"].as_str().unwrap());
        let want_lon = bits(mn["lon_rad"].as_str().unwrap());
        let want_h = bits(mn["height_m"].as_str().unwrap());
        for (label, got, want) in [
            ("lat_rad", g.lat_rad, want_lat),
            ("lon_rad", g.lon_rad, want_lon),
            ("height_m", g.height_m, want_h),
        ] {
            checks += 1;
            let u = ulp_distance(got, want);
            if u != 0 {
                failures.push(format!("meters_native.{label}: {u} ULP"));
            }
        }

        // Core public deg API vs the recipe's recorded core_api_deg.
        let (lat_deg, lon_deg, alt_km) =
            test_support::itrs_to_geodetic_core_km(km[0], km[1], km[2]);
        let want_lat_deg = bits(case["core_api_deg"]["lat_deg"].as_str().unwrap());
        let want_lon_deg = bits(case["core_api_deg"]["lon_deg"].as_str().unwrap());
        let want_alt_km = bits(case["core_internal"]["alt_km"].as_str().unwrap());
        for (label, got, want) in [
            ("core.lat_deg", lat_deg, want_lat_deg),
            ("core.lon_deg", lon_deg, want_lon_deg),
            ("core.alt_km", alt_km, want_alt_km),
        ] {
            checks += 1;
            let u = ulp_distance(got, want);
            if u != 0 {
                failures.push(format!("{label}: {u} ULP"));
            }
        }
    }

    assert!(checks > 0, "no cross-check components asserted");
    assert!(
        failures.is_empty(),
        "geodetic boundary cross-check diverged on {} of {checks} components:\n  {}",
        failures.len(),
        failures.join("\n  ")
    );
}

// ---------------------------------------------------------------------------
// TRACK 1 (0 ULP): trace-replay of the per-satellite intermediates, the
// weighted residual vector, and the 2-point FD Jacobian at each recorded state.
// ---------------------------------------------------------------------------

fn used_sats(doc: &Value) -> Vec<GnssSatelliteId> {
    doc["fixture"]["used_sats"]
        .as_array()
        .expect("used_sats")
        .iter()
        .map(|s| parse_prn(s.as_str().unwrap()))
        .collect()
}

/// The weighted residual closure (`sqrt(w) * (P_meas - P_hat)`) used by the FD
/// Jacobian, mirroring what scipy differences and what `LeastSquaresProblem::
/// with_weights` scales.
fn weighted_residual_at(
    sp3: &crate::sp3::Sp3,
    used: &[GnssSatelliteId],
    obs_by_id: &[(GnssSatelliteId, f64)],
    sqrt_w: &[f64],
    inputs: &Inputs,
    x: &[f64; 4],
) -> DVector<f64> {
    let rx = [x[0], x[1], x[2]];
    let b = x[3];
    let r: Vec<f64> = used
        .iter()
        .enumerate()
        .map(|(i, &sat)| {
            let p_meas = obs_by_id
                .iter()
                .find(|(id, _)| *id == sat)
                .map(|(_, p)| *p)
                .unwrap();
            let m = test_support::sat_model_for_test(
                sp3,
                sat,
                rx,
                b,
                p_meas,
                inputs.t_rx_j2000_s,
                inputs.sod_s,
                inputs.doy,
                inputs.corrections,
                &inputs.klobuchar,
                &inputs.met,
            )
            .expect("ephemeris present at trace state");
            sqrt_w[i] * (p_meas - m.p_hat_m)
        })
        .collect();
    DVector::from_vec(r)
}

fn trace_replay_level(level: &str) {
    let name = fixture_name(level);
    let doc = read_fixture(&name);
    let f = &doc["fixture"];
    let inputs = load_inputs(&doc, level);
    let sp3 = sp3();

    let used = used_sats(&doc);
    let obs_by_id: Vec<(GnssSatelliteId, f64)> = inputs
        .observations
        .iter()
        .map(|o| (o.satellite_id, o.pseudorange_m))
        .collect();

    // sqrt(weight) per used satellite, from the recorded frozen geometry.
    let geom = &f["used_sat_geometry"];
    let sqrt_w: Vec<f64> = used
        .iter()
        .map(|id| bits(geom[id.to_string()]["sqrt_weight"].as_str().unwrap()))
        .collect();

    let mut failures = Vec::new();
    let mut checks = 0usize;
    let mut check = |label: String, got: f64, want: f64, failures: &mut Vec<String>| {
        checks += 1;
        let u = ulp_distance(got, want);
        if u != 0 {
            failures.push(format!(
                "{label}: {u} ULP (rust=0x{:016x} ref=0x{:016x})",
                got.to_bits(),
                want.to_bits()
            ));
        }
    };

    let states = f["trace_states"].as_array().expect("trace_states");
    assert!(!states.is_empty(), "{level}: no trace states");

    for st in states {
        let ti = st["trace_index"].as_i64().unwrap();
        let xv = st["x"].as_array().unwrap();
        let x = [
            bits(xv[0].as_str().unwrap()),
            bits(xv[1].as_str().unwrap()),
            bits(xv[2].as_str().unwrap()),
            bits(xv[3].as_str().unwrap()),
        ];
        let rx = [x[0], x[1], x[2]];
        let b = x[3];

        // (1) Per-satellite named intermediates, bit-for-bit.
        let per_sat = st["per_sat"].as_array().unwrap();
        for (i, &sat) in used.iter().enumerate() {
            let ps = &per_sat[i];
            assert_eq!(
                ps["prn"].as_str().unwrap(),
                sat.to_string(),
                "{level}.state{ti}: per_sat order"
            );
            let p_meas = obs_by_id
                .iter()
                .find(|(id, _)| *id == sat)
                .map(|(_, p)| *p)
                .unwrap();
            let m = test_support::sat_model_for_test(
                &sp3,
                sat,
                rx,
                b,
                p_meas,
                inputs.t_rx_j2000_s,
                inputs.sod_s,
                inputs.doy,
                inputs.corrections,
                &inputs.klobuchar,
                &inputs.met,
            )
            .expect("ephemeris present");

            let pfx = format!("{level}.state{ti}.{sat}");
            check(
                format!("{pfx}.tau_s"),
                m.tau_s,
                bits(ps["tau_s"].as_str().unwrap()),
                &mut failures,
            );
            check(
                format!("{pfx}.t_tx_j2000_s"),
                m.t_tx_j2000_s,
                bits(ps["t_tx_j2000_s"].as_str().unwrap()),
                &mut failures,
            );
            let se = arr3(&ps["sat_ecef_m"]);
            check(
                format!("{pfx}.sat_ecef_x"),
                m.sat_ecef_m[0],
                se[0],
                &mut failures,
            );
            check(
                format!("{pfx}.sat_ecef_y"),
                m.sat_ecef_m[1],
                se[1],
                &mut failures,
            );
            check(
                format!("{pfx}.sat_ecef_z"),
                m.sat_ecef_m[2],
                se[2],
                &mut failures,
            );
            check(
                format!("{pfx}.dt_sat_s"),
                m.dt_sat_s,
                bits(ps["dt_sat_s"].as_str().unwrap()),
                &mut failures,
            );
            check(
                format!("{pfx}.theta_rad"),
                m.theta_rad,
                bits(ps["theta_rad"].as_str().unwrap()),
                &mut failures,
            );
            let sr = arr3(&ps["sat_rot_ecef_m"]);
            check(
                format!("{pfx}.sat_rot_x"),
                m.sat_rot_ecef_m[0],
                sr[0],
                &mut failures,
            );
            check(
                format!("{pfx}.sat_rot_y"),
                m.sat_rot_ecef_m[1],
                sr[1],
                &mut failures,
            );
            check(
                format!("{pfx}.sat_rot_z"),
                m.sat_rot_ecef_m[2],
                sr[2],
                &mut failures,
            );
            check(
                format!("{pfx}.rho_m"),
                m.rho_m,
                bits(ps["rho_m"].as_str().unwrap()),
                &mut failures,
            );
            check(
                format!("{pfx}.az_rad"),
                m.az_rad,
                bits(ps["az_rad"].as_str().unwrap()),
                &mut failures,
            );
            check(
                format!("{pfx}.el_rad"),
                m.el_rad,
                bits(ps["el_rad"].as_str().unwrap()),
                &mut failures,
            );
            check(
                format!("{pfx}.iono_m"),
                m.iono_m,
                bits(ps["iono_m"].as_str().unwrap()),
                &mut failures,
            );
            check(
                format!("{pfx}.tropo_m"),
                m.tropo_m,
                bits(ps["tropo_m"].as_str().unwrap()),
                &mut failures,
            );
            check(
                format!("{pfx}.p_hat_m"),
                m.p_hat_m,
                bits(ps["p_hat_m"].as_str().unwrap()),
                &mut failures,
            );
            // The weighted residual the solver sees.
            let r_w = sqrt_w[i] * (p_meas - m.p_hat_m);
            check(
                format!("{pfx}.residual_m"),
                r_w,
                bits(ps["residual_m"].as_str().unwrap()),
                &mut failures,
            );
        }

        // (2) Weighted residual vector (used-sat order).
        let r = weighted_residual_at(&sp3, &used, &obs_by_id, &sqrt_w, &inputs, &x);
        let res_v = st["residual"].as_array().unwrap();
        for (i, want) in res_v.iter().enumerate() {
            check(
                format!("{level}.state{ti}.residual[{i}]"),
                r[i],
                bits(want.as_str().unwrap()),
                &mut failures,
            );
        }

        // (3) 2-point FD Jacobian of the weighted residual at x. The fixture
        // records the FD rel_step; assert it equals the crate constant so the
        // step itself is pinned, then compare the assembled Jacobian.
        let fd = &st["fd_2point"];
        check(
            format!("{level}.state{ti}.fd_rel_step"),
            FD_REL_STEP_2POINT,
            bits(fd["rel_step"].as_str().unwrap()),
            &mut failures,
        );

        let f0 = r.clone();
        let x_vec = DVector::from_row_slice(&x);
        let resid_closure = |p: &DVector<f64>| -> DVector<f64> {
            let pa = [p[0], p[1], p[2], p[3]];
            weighted_residual_at(&sp3, &used, &obs_by_id, &sqrt_w, &inputs, &pa)
        };
        let jac = jacobian_2point(resid_closure, &x_vec, &f0);

        let jac_v = fd["jac"].as_array().unwrap();
        for (row, want_row) in jac_v.iter().enumerate() {
            let cols = want_row.as_array().unwrap();
            for (col, want) in cols.iter().enumerate() {
                check(
                    format!("{level}.state{ti}.jac[{row}][{col}]"),
                    jac[(row, col)],
                    bits(want.as_str().unwrap()),
                    &mut failures,
                );
            }
        }
    }

    assert!(checks > 0, "{level}: no components checked");
    assert!(
        failures.is_empty(),
        "{level}: SPP substrate diverged from the reference recipe on {} of {checks} components:\n  {}",
        failures.len(),
        failures.join("\n  ")
    );
}

#[test]
fn trace_replay_l0_minimal_zero_ulp() {
    trace_replay_level("L0_minimal");
}

#[test]
fn trace_replay_l1_iono_zero_ulp() {
    trace_replay_level("L1_iono");
}

#[test]
fn trace_replay_l2_tropo_zero_ulp() {
    trace_replay_level("L2_tropo");
}

#[test]
fn trace_replay_l3_relativistic_zero_ulp() {
    trace_replay_level("L3_relativistic");
}

/// The relativistic level is a documented no-op: SP3 precise clocks are used
/// as-is with no separate periodic term, so L3 must reproduce L2 bit-for-bit at
/// every recorded per-satellite predicted range.
#[test]
fn relativistic_level_equals_tropo_level() {
    let l2 = read_fixture("spp_trace_L2_tropo.json");
    let l3 = read_fixture("spp_trace_L3_relativistic.json");
    let p2 = &l2["fixture"]["trace_states"][0]["per_sat"];
    let p3 = &l3["fixture"]["trace_states"][0]["per_sat"];
    let a2 = p2.as_array().unwrap();
    let a3 = p3.as_array().unwrap();
    assert_eq!(a2.len(), a3.len(), "L2/L3 used-sat count differs");
    for (s2, s3) in a2.iter().zip(a3) {
        assert_eq!(
            s2["p_hat_m"].as_str().unwrap(),
            s3["p_hat_m"].as_str().unwrap(),
            "relativistic no-op changed p_hat for {}",
            s2["prn"].as_str().unwrap()
        );
    }
}

// ---------------------------------------------------------------------------
// TRACK 2 (sub-micron, BLAS-bound): independent-solve agreement. The crate
// solver runs from the same inputs; the converged position/clock is asserted to
// agree with both the recorded scipy solution and the synthesized truth to a
// documented sub-micron bound. This is a SOLVER-AGREEMENT check, never a 0-ULP
// physics claim.
// ---------------------------------------------------------------------------

/// Documented agreement bound: the converged receiver coordinates and the clock
/// length agree to better than one micron with both the recorded scipy solution
/// and the synthesized truth. The trust-region linear-algebra step is owned by
/// the platform BLAS, so the last bits are not independently reproducible;
/// sub-micron is the non-obtrusive bar from the spec's Solver-reality section.
///
/// In practice this implementation agrees to a few nanometres (observed
/// ~5e-9 m, ~7e-16 relative on the ~6.3e6 m position magnitude); the bound is
/// held at the spec's sub-micron figure with comfortable margin, and is a
/// SOLVER-AGREEMENT check, explicitly not a 0-ULP physics-parity claim.
const AGREEMENT_BOUND_M: f64 = 1.0e-6;

fn independent_solve_level(level: &str) {
    let name = fixture_name(level);
    let doc = read_fixture(&name);
    let f = &doc["fixture"];
    let inputs = load_inputs(&doc, level);
    let sp3 = sp3();

    let sol = solve(&sp3, &solve_inputs(&inputs), true).expect("solve converges");

    // used_sats / rejected_sats match the fixture exactly (deterministic order).
    let want_used = used_sats(&doc);
    assert_eq!(sol.used_sats, want_used, "{level}: used_sats order/content");

    let want_rej = f["rejected_sats"].as_array().unwrap();
    assert_eq!(
        sol.rejected_sats.len(),
        want_rej.len(),
        "{level}: rejected count"
    );
    for (got, want) in sol.rejected_sats.iter().zip(want_rej) {
        assert_eq!(
            got.satellite_id,
            parse_prn(want["id"].as_str().unwrap()),
            "{level}: rejected id"
        );
        let want_reason = match want["reason"].as_str().unwrap() {
            "no_ephemeris" => RejectionReason::NoEphemeris,
            "low_elevation" => RejectionReason::LowElevation,
            other => panic!("unexpected rejection reason {other}"),
        };
        assert_eq!(
            got.reason, want_reason,
            "{level}: rejected reason for {}",
            got.satellite_id
        );
    }

    // Converged position/clock vs the recorded scipy solution.
    let fs = &f["final_solution"];
    let scipy_x = fs["x"].as_array().unwrap();
    let sx = [
        bits(scipy_x[0].as_str().unwrap()),
        bits(scipy_x[1].as_str().unwrap()),
        bits(scipy_x[2].as_str().unwrap()),
        bits(scipy_x[3].as_str().unwrap()),
    ];
    let got = [
        sol.position.x_m,
        sol.position.y_m,
        sol.position.z_m,
        sol.rx_clock_s * super::C_M_S,
    ];
    for (k, (g, s)) in got.iter().zip(sx.iter()).enumerate() {
        assert!(
            (g - s).abs() <= AGREEMENT_BOUND_M,
            "{level}: component {k} disagrees with scipy: |{g} - {s}| = {} > {AGREEMENT_BOUND_M} m",
            (g - s).abs()
        );
    }

    // Converged position/clock vs the synthesized truth.
    let tx = fs["truth_x"].as_array().unwrap();
    let truth = [
        bits(tx[0].as_str().unwrap()),
        bits(tx[1].as_str().unwrap()),
        bits(tx[2].as_str().unwrap()),
        bits(tx[3].as_str().unwrap()),
    ];
    for (k, (g, t)) in got.iter().zip(truth.iter()).enumerate() {
        assert!(
            (g - t).abs() <= AGREEMENT_BOUND_M,
            "{level}: component {k} disagrees with truth: |{g} - {t}| = {} > {AGREEMENT_BOUND_M} m",
            (g - t).abs()
        );
    }

    // The clock-second boundary: rx_clock_s == b_m / c.
    let want_clock_s = bits(fs["truth_rx_clock_s"].as_str().unwrap());
    assert!(
        (sol.rx_clock_s - want_clock_s).abs() <= AGREEMENT_BOUND_M / super::C_M_S,
        "{level}: rx_clock_s off by {}",
        (sol.rx_clock_s - want_clock_s).abs()
    );

    assert!(sol.metadata.converged, "{level}: solver did not converge");
    assert!(sol.dop.is_some(), "{level}: DOP missing");
}

#[test]
fn independent_solve_l0_agreement() {
    independent_solve_level("L0_minimal");
}

#[test]
fn independent_solve_l1_agreement() {
    independent_solve_level("L1_iono");
}

#[test]
fn independent_solve_l2_agreement() {
    independent_solve_level("L2_tropo");
}

#[test]
fn independent_solve_l3_agreement() {
    independent_solve_level("L3_relativistic");
}

// ---------------------------------------------------------------------------
// DOP from the converged geometry agrees with the recipe's recorded DOP.
// (BLAS-agreement track: the DOP recipe is 0-ULP at a fixed geometry, but here
// it is recomputed at the independently-converged position, so it is a
// sub-micron-geometry agreement, not a 0-ULP claim.)
// ---------------------------------------------------------------------------
#[test]
fn dop_from_converged_geometry_agrees() {
    for &level in LEVELS {
        let doc = read_fixture(&fixture_name(level));
        let inputs = load_inputs(&doc, level);
        let sp3 = sp3();
        let sol = solve(&sp3, &solve_inputs(&inputs), false).expect("solve");
        let dop = sol.dop.expect("dop present");
        let want = &doc["fixture"]["dop"];
        for (label, got) in [
            ("gdop", dop.gdop),
            ("pdop", dop.pdop),
            ("hdop", dop.hdop),
            ("vdop", dop.vdop),
            ("tdop", dop.tdop),
        ] {
            let w = bits(want[label].as_str().unwrap());
            let rel = (got - w).abs() / w.max(1.0);
            assert!(
                rel <= 1e-9,
                "{level}: {label} disagrees: rust={got} ref={w} (rel {rel})"
            );
        }
    }
}

// ---------------------------------------------------------------------------
// Failure/rejection behavior.
// ---------------------------------------------------------------------------

/// Fewer than four usable satellites is the documented underdetermined failure.
#[test]
fn too_few_satellites_rejected() {
    let doc = read_fixture("spp_trace_L0_minimal.json");
    let mut inputs = load_inputs(&doc, "L0_minimal");
    // Keep only the first three used satellites' observations.
    let used = used_sats(&doc);
    let keep: Vec<_> = used.iter().take(3).copied().collect();
    inputs
        .observations
        .retain(|o| keep.contains(&o.satellite_id));
    let sp3 = sp3();
    match solve(&sp3, &solve_inputs(&inputs), false) {
        Err(SppError::TooFewSatellites { used, required }) => {
            // GPS-only here, so the requirement is the classic four.
            assert!(used < 4, "expected <4 usable, got {used}");
            assert_eq!(required, 4, "single-system solve requires 4 satellites");
        }
        other => panic!("expected TooFewSatellites, got {other:?}"),
    }
}

/// A satellite present in the observations but absent from the SP3 product is
/// rejected with `no_ephemeris` (the residual path must error/skip it, never
/// panic), and the solve still succeeds on the remaining satellites.
#[test]
fn no_ephemeris_satellite_is_rejected() {
    let doc = read_fixture("spp_trace_L0_minimal.json");
    let mut inputs = load_inputs(&doc, "L0_minimal");
    // GPS PRN 99 is not in the product, so it has no ephemeris at any epoch.
    let ghost = GnssSatelliteId::new(GnssSystem::Gps, 99);
    inputs.observations.push(Observation {
        satellite_id: ghost,
        pseudorange_m: 2.2e7,
    });

    let sp3 = sp3();
    let sol =
        solve(&sp3, &solve_inputs(&inputs), false).expect("solve succeeds on the real satellites");
    assert!(
        sol.rejected_sats
            .iter()
            .any(|r| r.satellite_id == ghost && r.reason == RejectionReason::NoEphemeris),
        "ghost satellite should be rejected with no_ephemeris; rejected = {:?}",
        sol.rejected_sats
    );
}

/// Duplicate observations for the same satellite are rejected deterministically
/// (by the smallest repeated id), so the result can never depend on input order.
#[test]
fn duplicate_observation_is_rejected() {
    let doc = read_fixture("spp_trace_L0_minimal.json");
    let mut inputs = load_inputs(&doc, "L0_minimal");
    let dup = inputs.observations[0];
    // Push a second observation for the same satellite with a different range.
    inputs.observations.push(Observation {
        satellite_id: dup.satellite_id,
        pseudorange_m: dup.pseudorange_m + 1234.5,
    });

    let sp3 = sp3();
    match solve(&sp3, &solve_inputs(&inputs), false) {
        Err(SppError::DuplicateObservation { satellite }) => {
            assert_eq!(satellite, dup.satellite_id)
        }
        other => panic!("expected DuplicateObservation, got {other:?}"),
    }
}

/// The residual path returns `Err(satellite)` instead of panicking when a used
/// satellite cannot be modeled at the query state. This is the condition
/// `solve()` records in its closure and surfaces as `SppError::EphemerisLost`
/// if it occurs during a solver probe (the harder case where a satellite
/// survives selection and then drops). With real SP3 coverage a
/// selection-surviving satellite does not actually drop across the bounded
/// transmit-time probes, so the `Err` path itself is exercised directly here to
/// lock in the no-panic guarantee the solver relies on.
#[test]
fn residual_errs_instead_of_panicking_on_unmodelable_satellite() {
    let doc = read_fixture("spp_trace_L0_minimal.json");
    let si = solve_inputs(&load_inputs(&doc, "L0_minimal"));
    let sp3 = sp3();

    // A satellite with no ephemeris in the product, placed in the `used` set as
    // if it had survived selection; the residual must error on it, not panic.
    let ghost = GnssSatelliteId::new(GnssSystem::Gps, 99);
    let used = [ghost];
    let obs_by_id = [(ghost, 2.2e7)];
    let r = super::residual_unweighted(&sp3, &used, &obs_by_id, &si.initial_guess, &si);
    assert_eq!(
        r,
        Err(ghost),
        "residual must return Err for an unmodelable used satellite, never panic"
    );
}

/// The solve's rejected set, reasons, and order match the fixture exactly. This
/// L0 geometry exercises the `low_elevation` mask; the `no_ephemeris` branch is
/// covered separately by `no_ephemeris_satellite_is_rejected`.
#[test]
fn rejection_reasons_match_fixture() {
    let doc = read_fixture("spp_trace_L0_minimal.json");
    let inputs = load_inputs(&doc, "L0_minimal");
    let sp3 = sp3();
    let sol = solve(&sp3, &solve_inputs(&inputs), false).expect("solve");

    let want: Vec<(GnssSatelliteId, RejectionReason)> = doc["fixture"]["rejected_sats"]
        .as_array()
        .unwrap()
        .iter()
        .map(|r| {
            let id = parse_prn(r["id"].as_str().unwrap());
            let reason = match r["reason"].as_str().unwrap() {
                "no_ephemeris" => RejectionReason::NoEphemeris,
                "low_elevation" => RejectionReason::LowElevation,
                other => panic!("unexpected reason {other}"),
            };
            (id, reason)
        })
        .collect();

    let got: Vec<(GnssSatelliteId, RejectionReason)> = sol
        .rejected_sats
        .iter()
        .map(|r| (r.satellite_id, r.reason))
        .collect();
    assert_eq!(
        got, want,
        "rejected set/reasons/order diverged from fixture"
    );
    // At least one low-elevation rejection exists in this geometry.
    assert!(
        want.iter()
            .any(|(_, r)| *r == RejectionReason::LowElevation),
        "fixture should exercise the low-elevation mask"
    );
}

/// A degenerate geometry — several satellites at coincident positions, so every
/// line-of-sight row is identical and the design matrix is rank-deficient — must
/// be handled gracefully, never a panic. The Levenberg-damped solver regularizes
/// the rank deficiency and converges, but the DOP cofactor inverse detects the
/// singular geometry and reports no DOP. (`SppError::Singular` is the defensive
/// path for a subproblem the damping cannot step through, and is not reachable
/// from this robustly-regularized geometry.)
#[test]
fn degenerate_geometry_is_handled_gracefully() {
    let bytes = std::fs::read(fixture_path("sp3/degenerate_coincident_5sat.sp3"))
        .expect("read degenerate SP3");
    let sp3 = crate::sp3::Sp3::parse(&bytes).expect("parse degenerate SP3");

    // Identical pseudorange for every satellite, so the Sagnac rotation is the
    // same for all of them and the rows stay coincident (rank-deficient).
    let p = 20_181_863.0;
    let observations = (1..=5)
        .map(|prn| Observation {
            satellite_id: GnssSatelliteId::new(GnssSystem::Gps, prn),
            pseudorange_m: p,
        })
        .collect();
    // A receive epoch inside the product's [00:00, 00:15] window.
    let inputs = SolveInputs {
        observations,
        t_rx_j2000_s: 646_229_000.0,
        t_rx_second_of_day_s: 200.0,
        day_of_year: 176.0,
        initial_guess: [6_378_137.0, 0.0, 0.0, 0.0],
        corrections: Corrections::NONE,
        klobuchar: KlobucharCoeffs {
            alpha: [0.0; 4],
            beta: [0.0; 4],
        },
        beidou_klobuchar: None,
        met: SurfaceMet {
            pressure_hpa: 1013.25,
            temperature_k: 288.15,
            relative_humidity: 0.5,
        },
    };

    match solve(&sp3, &inputs, false) {
        // Damped solver converges, but the rank-deficient geometry yields no DOP.
        Ok(sol) => assert!(
            sol.dop.is_none(),
            "a rank-deficient geometry must not report a DOP, got {:?}",
            sol.dop
        ),
        // Also acceptable if the subproblem genuinely cannot be stepped.
        Err(SppError::Singular(_)) => {}
        other => {
            panic!("degenerate geometry mishandled (expected Ok/no-DOP or Singular): {other:?}")
        }
    }
}

/// The satellite ordering the solve pins (`GnssSatelliteId` `Ord`) matches a
/// zero-padded PRN string sort for GPS, so the Rust ordering agrees with the
/// fixture generator's string-keyed sort. This is the GPS-only v1 assumption;
/// multi-GNSS would need the same property to hold across the system letters.
#[test]
fn gnss_satellite_id_orders_like_zero_padded_prn_strings() {
    let mut ids: Vec<GnssSatelliteId> = (1..=12u8)
        .rev()
        .map(|prn| GnssSatelliteId::new(GnssSystem::Gps, prn))
        .collect();
    ids.sort();
    let by_ord: Vec<String> = ids.iter().map(|s| s.to_string()).collect();

    let mut by_string = by_ord.clone();
    by_string.sort();

    assert_eq!(
        by_ord, by_string,
        "GnssSatelliteId Ord must match the zero-padded PRN string order"
    );
    assert_eq!(by_ord.first().map(String::as_str), Some("G01"));
    assert_eq!(by_ord.last().map(String::as_str), Some("G12"));
}