sidereon-core 0.9.0

//! Parity tests for SP3 interpolation.
//!
//! Two channels, two references:
//!
//! - **Position**: the sliding-window high-degree Lagrange (Neville) recipe is
//!   matched against **RTKLIB** `preceph.c` `interppol`/`pephpos` (the IGS
//!   standard), NOT scipy. The golden values were generated by building RTKLIB's
//!   `peph2pos` (opt = 0: no antenna offset, no transmit->receive Sagnac) against
//!   the same public IGS SP3 fixture this crate ships, dumping satellite ECEF
//!   positions at interior, near-day-boundary, and near-gap query epochs. The
//!   crate must match to sub-centimetre. RTKLIB's `peph2pos` applies the same
//!   OMGE per-node Earth-rotation correction this module does, so the comparison
//!   is frame-consistent (both positions are in the ECEF frame at the query
//!   instant).
//! - **Clock**: still the not-a-knot `scipy.interpolate.CubicSpline`; the small-
//!   node-count and J2000-second-exactness tests below cover the shared spline
//!   evaluator and the query conversion.

use super::{
    eval_cubic_spline_for_test as eval_spline, instant_to_j2000_seconds,
    interpolate_position_neville,
};
use crate::astro::constants::time::SECONDS_PER_DAY_I64;
use crate::astro::time::civil::{J2000_JULIAN_DAY_NUMBER, J2000_NOON_OFFSET_S};
use crate::astro::time::j2000_seconds_from_split;
use crate::astro::time::model::{Instant, TimeScale};
use crate::astro::time::scales::julian_day_number;
use crate::constants::J2000_JD;
use crate::Error;

/// Sub-centimetre agreement bound for the RTKLIB position parity (3D component
/// magnitude). RTKLIB carries SP3 positions in km internally and prints/holds
/// them in double precision; the dominant residual vs this crate is the last
/// few ULP of the Neville recurrence and the OMGE sin/cos, far under 1 cm.
const RTKLIB_POS_TOL_M: f64 = 0.01;

const KM_TO_M: f64 = 1_000.0;
const US_TO_S: f64 = 1.0e-6;

// --- shared cubic-spline evaluator coverage (clock channel) ---

#[test]
fn small_node_count_special_cases() {
    // n == 2: not-a-knot degenerates to the straight Hermite (both endpoint
    // slopes = secant slope), which is BLAS-free (no solve) and 0-ULP vs
    // scipy CubicSpline([0,900],[10,17.5])(437.0).
    let v2 = eval_spline(&[0.0, 900.0], &[10.0, 17.5], 437.0);
    assert_eq!(
        v2.to_bits(),
        0x402b488888888888,
        "n==2 not 0-ULP: {:#018x}",
        v2.to_bits()
    );

    // n == 3 not-a-knot: scipy solves a DENSE 3x3 parabola system via
    // `scipy.linalg.solve` -> LAPACK `gesv`. Held to a tight agreement bound
    // (the dense-solve BLAS caveat), not 0 ULP.
    let x3 = [0.0, 900.0, 1800.0];
    let y3 = [10.0, 17.5, 9.0];
    let v3a = eval_spline(&x3, &y3, 437.0);
    let v3b = eval_spline(&x3, &y3, 1500.0);
    let want_a = f64::from_bits(0x402f47adc1a12a3a);
    let want_b = f64::from_bits(0x402b38e38e38e38d);
    assert!(
        (v3a - want_a).abs() <= 4.0 * f64::EPSILON * want_a.abs(),
        "n==3 @437 off: got {} want {}",
        v3a,
        want_a
    );
    assert!(
        (v3b - want_b).abs() <= 4.0 * f64::EPSILON * want_b.abs(),
        "n==3 @1500 off: got {} want {}",
        v3b,
        want_b
    );
}

// --- RTKLIB position parity (interior + day-boundary) ---

/// One RTKLIB-reference position case.
///
/// `query` is the J2000 second (matching this crate's floored node axis), and
/// `pos_m` is RTKLIB `peph2pos(opt=0)` satellite ECEF position in meters.
struct Case {
    prn: &'static str,
    query: f64,
    pos_m: [f64; 3],
}

/// Interior + near-day-boundary cases, RTKLIB `peph2pos` over the full
/// `GRG0MGXFIN_20201760000_01D_15M_ORB.SP3` fixture for G01/G15/G32. The first
/// two queries per sat (00:07:30, 00:22:30) are the day-boundary intervals where
/// the global cubic spline erred by hundreds of metres; the rest are interior.
const INTERIOR_CASES: &[Case] = &[
    Case {
        prn: "G01",
        query: 646229250.0,
        pos_m: [-11125530.998355811, 19912984.19788936, -13542820.592182232],
    },
    Case {
        prn: "G01",
        query: 646230150.0,
        pos_m: [-12314137.386707624, 20651008.347774327, -11126395.854435727],
    },
    Case {
        prn: "G01",
        query: 646233075.0,
        pos_m: [-14503232.889903756, 21943065.43745991, -2176901.884074253],
    },
    Case {
        prn: "G01",
        query: 646272420.0,
        pos_m: [11256448.174644291, -19993852.434904076, -13306698.18210848],
    },
    Case {
        prn: "G01",
        query: 646274250.0,
        pos_m: [13359844.15738438, -21330281.774696253, -8151899.722851374],
    },
    Case {
        prn: "G01",
        query: 646293825.0,
        pos_m: [19940761.656224072, 11879948.063645441, 12987243.364720155],
    },
    Case {
        prn: "G01",
        query: 646312050.0,
        pos_m: [-4625770.090611601, 16746354.107113319, -20215259.682363085],
    },
    Case {
        prn: "G01",
        query: 646314750.0,
        pos_m: [-10099001.26139606, 19319559.36596251, -15160098.379524706],
    },
    Case {
        prn: "G15",
        query: 646229250.0,
        pos_m: [5706609.530972648, -21293396.701827034, 14224098.286349315],
    },
    Case {
        prn: "G15",
        query: 646230150.0,
        pos_m: [6460420.873947138, -19623115.117335513, 16165588.658165246],
    },
    Case {
        prn: "G15",
        query: 646233075.0,
        pos_m: [9909606.248824, -13334016.583217062, 20338698.406210665],
    },
    Case {
        prn: "G15",
        query: 646272420.0,
        pos_m: [-5776477.390518184, 21130639.11271612, 14435847.801910691],
    },
    Case {
        prn: "G15",
        query: 646274250.0,
        pos_m: [-7488450.286292078, 17536524.71907534, 18019350.08859665],
    },
    Case {
        prn: "G15",
        query: 646293825.0,
        pos_m: [-21790707.649836697, -5926051.980660642, -14590940.781290831],
    },
    Case {
        prn: "G15",
        query: 646312050.0,
        pos_m: [3853078.970594666, -25512367.920297455, 5075952.642399781],
    },
    Case {
        prn: "G15",
        query: 646314750.0,
        pos_m: [5237281.443057492, -22387368.612346206, 12648466.68057392],
    },
    Case {
        prn: "G32",
        query: 646229250.0,
        pos_m: [-14982710.33984932, -11416137.04370103, -18698169.967418224],
    },
    Case {
        prn: "G32",
        query: 646230150.0,
        pos_m: [
            -15202082.969047092,
            -13469926.879499504,
            -17106603.031938538,
        ],
    },
    Case {
        prn: "G32",
        query: 646233075.0,
        pos_m: [
            -16153963.273672031,
            -18576222.673382334,
            -10077096.844378414,
        ],
    },
    Case {
        prn: "G32",
        query: 646272420.0,
        pos_m: [15002030.462289223, 11632546.749006014, -18550835.263271578],
    },
    Case {
        prn: "G32",
        query: 646274250.0,
        pos_m: [15522666.139369572, 15547540.154455326, -14945081.460181221],
    },
    Case {
        prn: "G32",
        query: 646293825.0,
        pos_m: [-10949919.484357078, 15065869.848260796, 18966010.002153255],
    },
    Case {
        prn: "G32",
        query: 646312050.0,
        pos_m: [-15075567.757707063, -2590019.536678051, -21631778.62209533],
    },
    Case {
        prn: "G32",
        query: 646314750.0,
        pos_m: [-14877007.647514388, -9804548.32109721, -19658313.960636113],
    },
];

/// Near-gap cases. The crate is fed the GAPPED fixture (G01 missing the
/// 07:30..10:00 records); the golden is RTKLIB `peph2pos` run over the
/// pre-gap-only contiguous run (nodes 00:00..07:15). The crate must restrict its
/// window to that same contiguous run and reproduce the RTKLIB values, proving it
/// never interpolates across the gap.
const GAP_CASES: &[Case] = &[
    Case {
        prn: "G01",
        query: 646253550.0,
        pos_m: [-21780616.443619847, -14648637.054289348, 4885781.790914722],
    },
    Case {
        prn: "G01",
        query: 646254450.0,
        pos_m: [-21995073.802280508, -15046334.842752764, 2046724.44502548],
    },
    Case {
        prn: "G01",
        query: 646255350.0,
        pos_m: [-21971156.624103755, -15241796.37191432, -826930.205430558],
    },
];

fn id_for(prn: &str) -> crate::id::GnssSatelliteId {
    use crate::id::{GnssSatelliteId, GnssSystem};
    let n: u8 = prn[1..].parse().expect("prn number");
    GnssSatelliteId::new(GnssSystem::Gps, n).expect("valid satellite id")
}

fn assert_pos_matches(got: [f64; 3], want: [f64; 3], prn: &str, q: f64) {
    for axis in 0..3 {
        let d = (got[axis] - want[axis]).abs();
        assert!(
            d <= RTKLIB_POS_TOL_M,
            "{} axis {} at q={}: crate {} vs RTKLIB {} -> |diff| {:.6} m > {} m",
            prn,
            axis,
            q,
            got[axis],
            want[axis],
            d,
            RTKLIB_POS_TOL_M
        );
    }
}

/// Interior + near-day-boundary: the public `position()` API on the shipped IGS
/// fixture matches RTKLIB `pephpos` to sub-centimetre.
#[test]
fn position_matches_rtklib_interior_and_boundary() {
    use super::super::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 IGS SP3 fixture {path}: {e}"));
    let sp3 = Sp3::parse(&bytes).expect("parse IGS SP3");

    for case in INTERIOR_CASES {
        let id = id_for(case.prn);
        let state = sp3
            .position_at_j2000_seconds(id, case.query)
            .unwrap_or_else(|e| panic!("{} @ {}: position err {:?}", case.prn, case.query, e));
        assert_pos_matches(state.position.as_array(), case.pos_m, case.prn, case.query);
    }
}

#[test]
fn position_rejects_epoch_scale_mismatch() {
    use super::super::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 IGS SP3 fixture {path}: {e}"));
    let sp3 = Sp3::parse(&bytes).expect("parse IGS SP3");
    assert_eq!(sp3.header.time_scale, TimeScale::Gpst);

    let gps_epoch = sp3.epochs[0];
    let utc_epoch =
        Instant::from_julian_date(TimeScale::Utc, gps_epoch.julian_date().expect("split JD"));

    let err = sp3
        .position(id_for("G01"), utc_epoch)
        .expect_err("UTC-tagged epoch must not be evaluated on a GPST product axis");
    match err {
        Error::InvalidInput(msg) => {
            assert!(msg.contains("UTC"), "{msg}");
            assert!(msg.contains("GPST"), "{msg}");
        }
        other => panic!("expected InvalidInput for scale mismatch, got {other:?}"),
    }
}

#[test]
fn position_same_scale_matches_direct_j2000_query() {
    use super::super::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 IGS SP3 fixture {path}: {e}"));
    let sp3 = Sp3::parse(&bytes).expect("parse IGS SP3");
    let sat = id_for("G01");
    let epoch = sp3.epochs[4];

    let via_instant = sp3
        .position(sat, epoch)
        .expect("same-scale instant query succeeds");
    let query = instant_to_j2000_seconds(&epoch).expect("split JD converts to J2000 seconds");
    let via_seconds = sp3
        .position_at_j2000_seconds(sat, query)
        .expect("direct J2000 query succeeds");

    assert_eq!(via_instant, via_seconds);
}

/// Near-gap: the crate restricts the Lagrange window to the contiguous run and
/// matches RTKLIB run over that same run, never crossing the gap.
#[test]
fn position_matches_rtklib_near_gap() {
    use super::super::Sp3;

    let path = concat!(
        env!("CARGO_MANIFEST_DIR"),
        "/tests/fixtures/sp3/GAP_G01_20201760000_15M.sp3"
    );
    let bytes =
        std::fs::read(path).unwrap_or_else(|e| panic!("read gapped SP3 fixture {path}: {e}"));
    let sp3 = Sp3::parse(&bytes).expect("parse gapped SP3");

    for case in GAP_CASES {
        let id = id_for(case.prn);
        let state = sp3
            .position_at_j2000_seconds(id, case.query)
            .unwrap_or_else(|e| {
                panic!("{} @ {} (gap): position err {:?}", case.prn, case.query, e)
            });
        assert_pos_matches(state.position.as_array(), case.pos_m, case.prn, case.query);
    }
}

/// A query inside the coverage gap itself (no contiguous run brackets it within
/// the margin) is rejected, never silently interpolated across the gap.
#[test]
fn position_inside_gap_is_rejected() {
    use super::super::Sp3;

    let path = concat!(
        env!("CARGO_MANIFEST_DIR"),
        "/tests/fixtures/sp3/GAP_G01_20201760000_15M.sp3"
    );
    let bytes =
        std::fs::read(path).unwrap_or_else(|e| panic!("read gapped SP3 fixture {path}: {e}"));
    let sp3 = Sp3::parse(&bytes).expect("parse gapped SP3");
    let id = id_for("G01");
    // 08:45:00 GPST = 646260300 J2000 sec, deep inside the 07:30..10:00 G01 gap.
    let q = 646260300.0;
    assert!(
        sp3.position_at_j2000_seconds(id, q).is_err(),
        "query inside coverage gap must be rejected, not interpolated across"
    );
}

#[test]
fn position_gapped_product_edges_use_nominal_spacing() {
    use super::super::Sp3;

    let path = concat!(
        env!("CARGO_MANIFEST_DIR"),
        "/tests/fixtures/sp3/GAP_G01_20201760000_15M.sp3"
    );
    let bytes =
        std::fs::read(path).unwrap_or_else(|e| panic!("read gapped SP3 fixture {path}: {e}"));
    let sp3 = Sp3::parse(&bytes).expect("parse gapped SP3");
    let id = id_for("G01");
    let axis = sp3.epochs_j2000_seconds();
    let first = axis[0];
    let last = axis[axis.len() - 1];

    for query in [first - 950.0, last + 950.0] {
        assert!(
            matches!(
                sp3.position_at_j2000_seconds(id, query),
                Err(Error::EpochOutOfRange)
            ),
            "gapped product edge query {query} must be rejected beyond one nominal cadence"
        );
    }
}

#[test]
fn position_gap_window_uses_right_arc_before_second_edge() {
    let x = [0.0, 10.0, 20.0, 100.0, 110.0, 120.0];
    let kx = [0.0; 6];
    let ky = [0.0; 6];
    let kz = [-3.0, -3.0, -3.0, 7.0, 7.0, 7.0];

    let (_, _, left_z_m) = interpolate_position_neville(&x, &kx, &ky, &kz, 25.0);
    assert_eq!(left_z_m.to_bits(), (-3_000.0f64).to_bits());

    let (_, _, right_z_m) = interpolate_position_neville(&x, &kx, &ky, &kz, 95.0);
    assert_eq!(right_z_m.to_bits(), 7_000.0f64.to_bits());
}

#[test]
fn non_finite_direct_query_is_rejected() {
    use super::super::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 IGS SP3 fixture {path}: {e}"));
    let sp3 = Sp3::parse(&bytes).expect("parse IGS SP3");
    let id = id_for("G01");

    for query in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
        match sp3.position_at_j2000_seconds(id, query) {
            Err(crate::Error::InvalidInput(message)) => {
                assert!(message.contains("query_j2000_s"));
                assert!(message.contains("not finite"));
            }
            other => panic!("expected invalid non-finite query, got {other:?}"),
        }
    }
}

#[test]
fn duplicate_epochs_are_rejected_without_panic() {
    use super::super::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 IGS SP3 fixture {path}: {e}"));
    let mut sp3 = Sp3::parse(&bytes).expect("parse IGS SP3");
    sp3.epochs[1] = sp3.epochs[0];

    let result = std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
        sp3.position_at_j2000_seconds(id_for("G01"), 646229250.0)
    }));

    assert!(result.is_ok(), "duplicate SP3 epochs must not panic");
    match result.expect("duplicate SP3 epochs should not unwind") {
        Err(crate::Error::InvalidInput(message)) => {
            assert!(message.contains("strictly increasing"));
        }
        other => panic!("expected invalid duplicate epochs, got {other:?}"),
    }
}

/// Clock channel still interpolates (the cubic-spline path is unchanged): the
/// public API returns a clock at an interior epoch.
#[test]
fn clock_channel_still_interpolates() {
    use super::super::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 IGS SP3 fixture {path}: {e}"));
    let sp3 = Sp3::parse(&bytes).expect("parse IGS SP3");
    let state = sp3
        .position_at_j2000_seconds(id_for("G01"), 646272420.0)
        .expect("position+clock");
    assert!(state.clock_s.is_some(), "clock should interpolate");
}

/// `epochs_j2000_seconds` returns the exact node axis: one value per epoch, in
/// the file's time scale, ascending, equal to the per-epoch Julian-date
/// conversion. This is the public query grid the binding hands back to callers.
#[test]
fn epochs_j2000_seconds_is_the_exact_node_axis() {
    use super::super::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 IGS SP3 fixture {path}: {e}"));
    let sp3 = Sp3::parse(&bytes).expect("parse IGS SP3");

    let axis = sp3.epochs_j2000_seconds();
    assert_eq!(axis.len(), sp3.epoch_count(), "one value per epoch");
    assert!(
        axis.windows(2).all(|w| w[1] > w[0]),
        "epoch axis must be strictly ascending"
    );

    for (ep, &got) in sp3.epochs.iter().zip(axis.iter()) {
        let split = ep.julian_date().expect("SP3 epochs are Julian-date");
        let want = j2000_seconds_from_split(split.jd_whole, split.fraction);
        assert_eq!(got.to_bits(), want.to_bits(), "epoch axis not bit-exact");
    }
}

/// Direct exercise of the shared cubic-spline + unit-scaling path used by the
/// clock channel: a 1-D microsecond spline over a 9-node arc, * 1e-6 to seconds.
#[test]
fn clock_spline_unit_scaling() {
    const NODES: [f64; 9] = [
        646_228_800.0,
        646_229_700.0,
        646_230_600.0,
        646_231_500.0,
        646_232_400.0,
        646_233_300.0,
        646_234_200.0,
        646_235_100.0,
        646_236_000.0,
    ];
    const CLK_US: [f64; 9] = [
        123.456789, 123.966789, 124.496789, 125.046789, 125.616789, 126.206789, 126.816789,
        127.446789, 128.096789,
    ];
    let x = NODES.to_vec();
    let y = CLK_US.to_vec();
    // Interior query; the spline must lie between the bracketing node values.
    let c_s = eval_spline(&x, &y, 646_231_501.0) * US_TO_S;
    assert!(
        c_s > 125.0 * US_TO_S && c_s < 125.1 * US_TO_S,
        "clock spline out of expected interior range: {c_s}"
    );
    // km/us scaling constants are applied as single multiplies.
    let _ = KM_TO_M;
}

/// Pins the proven exactness bound of `instant_to_j2000_seconds` for the SP3
/// query path. The midnight term is an exact integer across the SP3 epoch
/// domain; the sub-second term round-trips to under 1 ULP at magnitude 86400.
#[test]
fn instant_to_j2000_seconds_query_exactness_bound() {
    const D: f64 = 86_400.0;

    let lo = julian_day_number(1994, 1, 1);
    let hi = julian_day_number(2100, 1, 1);
    assert!(
        (hi as f64) < 2f64.powi(53),
        "jdn must be exactly representable"
    );
    for jd in lo..=hi {
        let jd_whole = jd as f64 - 0.5;
        let days_whole = jd_whole - J2000_JD;
        let midnight = days_whole * D;
        assert_eq!(
            midnight.fract(),
            0.0,
            "midnight term not integer for jdn {jd}"
        );
        let exact = (jd - J2000_JULIAN_DAY_NUMBER) * SECONDS_PER_DAY_I64 - J2000_NOON_OFFSET_S;
        assert_eq!(
            midnight, exact as f64,
            "midnight term != exact i64 for jdn {jd}"
        );
    }

    const BOUND_S: f64 = 7.3e-12;
    let mut worst = 0.0f64;
    for s in 0..86_400u32 {
        let ds = s as f64;
        let e = ((ds / D) * D - ds).abs();
        if e > worst {
            worst = e;
        }
    }
    for grid in [2u32, 10, 100, 1000] {
        for i in 0..(86_400 * grid) {
            let ds = i as f64 / grid as f64;
            let e = ((ds / D) * D - ds).abs();
            if e > worst {
                worst = e;
            }
        }
    }
    assert!(
        worst <= BOUND_S,
        "sub-second round-trip error {worst:e} s exceeds proven bound {BOUND_S:e} s"
    );
}