astrodynamics-gnss 0.9.7

//! 0-ULP parity tests for the SP3 CubicSpline interpolation recipe.
//!
//! These assert the Rust port reproduces `scipy.interpolate.CubicSpline(x,y)`
//! (defaults: not-a-knot, extrapolate=True) bit-for-bit, using the committed
//! golden fixture `parity/fixtures/sp3_interp.json` (generated by the canonical
//! recipe `parity/generator/sp3_interp.py`). Comparison is hex-float / ULP
//! distance per component, the `skyfield_parity_test.exs` discipline (spec
//! line 125).
//!
//! The fixture's raw node substrate (J2000 integer seconds, positions in km,
//! clocks in us) is embedded below so the test rebuilds the identical splines.
//! Position references are in METERS (km spline * 1000), clocks in SECONDS
//! (us spline * 1e-6) — the orbis API boundary, matching `position_m`/`clock_s`
//! in the recipe.

use super::eval_cubic_spline_for_test as eval_spline;

/// J2000 integer-second nodes (from the fixture).
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,
];

/// Position nodes in km, columns X/Y/Z (from the fixture `positions_km`).
const POS_KM: [[f64; 3]; 9] = [
    [26560.0, 0.0, 12000.0],
    [26527.106916090303, 1327.2467358292165, 12004.9],
    [26427.910629784365, 2651.1755461397556, 12009.6],
    [26262.659829981283, 3968.476798498796, 12014.1],
    [26031.76830746338, 5275.857425916825, 12018.4],
    [25735.813920634722, 6570.049157800129, 12022.5],
    [25375.537151176093, 7847.81668892518, 12026.4],
    [24951.839253226382, 9105.965766016789, 12030.1],
    [24465.78000071663, 10341.351171717757, 12033.6],
];

/// Clock nodes in microseconds (from the fixture `clocks_us`).
const CLK_US: [f64; 9] = [
    123.456789, 123.966789, 124.496789, 125.046789, 125.616789, 126.206789, 126.816789, 127.446789,
    128.096789,
];

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

/// One golden case: query epoch (J2000 s), expected position bits (m), clock bits (s).
struct Case {
    query: f64,
    pos_m_bits: [u64; 3],
    clk_s_bits: u64,
}

const CASES: [Case; 4] = [
    Case {
        query: 646_229_237.0,
        pos_m_bits: [0x41795280309f6e2c, 0x4123ac631b0b7677, 0x4166e48c86718662],
        clk_s_bits: 0x3f2036bf708e3a06,
    },
    Case {
        query: 646_231_501.0,
        pos_m_bits: [0x41790bba79e30898, 0x414e49c7c743418c, 0x4166ea431c70c434],
        clk_s_bits: 0x3f2063e51559f4a5,
    },
    Case {
        query: 646_235_222.5,
        pos_m_bits: [0x4177bc7eb396b50d, 0x4161b110722753ff, 0x4166f24f84b74f04],
        clk_s_bits: 0x3f20b755739c9dfc,
    },
    Case {
        query: 646_236_300.0,
        pos_m_bits: [0x41772a32b1eb52a0, 0x41647fd33e587e12, 0x4166f454471c71c6],
        clk_s_bits: 0x3f20d1a25e72abaa,
    },
];

fn col(axis: usize) -> [f64; 9] {
    let mut out = [0.0; 9];
    for (i, row) in POS_KM.iter().enumerate() {
        out[i] = row[axis];
    }
    out
}

#[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`. Unlike the n>=4 tridiagonal path
    // (`dgtsv`, which is 0-ULP — see `position_clock_interpolation_is_0_ulp_vs_scipy`),
    // a general dense `gesv` is the BLAS/LAPACK-bound carve-out the spec marks
    // as NOT independently 0-ULP-reproducible: Apple
    // Accelerate's blocked/vectorized `dgesv` last bits are owned by that build.
    // The n==3 case also never occurs in real SP3 (hundreds of epochs per arc),
    // so it is held to a tight agreement bound, not 0 ULP, and labeled as such.
    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);
    // Agreement: <= a few ULP (here 2), sub-femto relative — the dense-solve
    // BLAS caveat, not an algorithmic discrepancy.
    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
    );
}

/// Build a minimal SP3-c file carrying the fixture's nodes for one GPS sat, so
/// the full public `Sp3::position` path (parse -> J2000-second node axis ->
/// per-axis spline -> unit scaling) is exercised end-to-end and matched against
/// the same scipy golden values.
fn fixture_sp3_file() -> String {
    use std::fmt::Write;
    // Epochs: 2020-06-24 00:00:00 + i*900s, the fixture's datetime nodes.
    let mut s = String::new();
    s.push_str("#cP2020  6 24  0  0  0.00000000       9 ORBIT IGS14 FIT  TST\n");
    s.push_str("## 2111 432000.00000000   900.00000000 59024 0.0000000000000\n");
    s.push_str("+    1   G01  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0\n");
    s.push_str("++         0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0\n");
    s.push_str("%c G  cc GPS ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc\n");
    s.push_str("%c cc cc ccc ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc\n");
    s.push_str("%f  1.2500000  1.025000000  0.00000000000  0.000000000000000\n");
    s.push_str("%f  0.0000000  0.000000000  0.00000000000  0.000000000000000\n");
    s.push_str("%i    0    0    0    0      0      0      0      0         0\n");
    s.push_str("%i    0    0    0    0      0      0      0      0         0\n");
    // 9 epochs, 15-minute spacing within the hour, rolling minutes.
    let minutes = [0, 15, 30, 45, 60, 75, 90, 105, 120];
    for (i, &min) in minutes.iter().enumerate() {
        let hh = min / 60;
        let mm = min % 60;
        writeln!(s, "*  2020  6 24  {:>2} {:>2}  0.00000000", hh, mm).unwrap();
        // Position in km (full precision via {:.12}), clock in us.
        writeln!(
            s,
            "PG01{:>14.6}{:>14.6}{:>14.6}{:>14.6}",
            POS_KM[i][0], POS_KM[i][1], POS_KM[i][2], CLK_US[i]
        )
        .unwrap();
    }
    s.push_str("EOF\n");
    s
}

#[test]
fn public_position_api_end_to_end_matches_scipy() {
    use crate::id::{GnssSatelliteId, GnssSystem};
    use crate::sp3::Sp3;
    use astrodynamics::time::model::{Instant, InstantRepr, JulianDateSplit, TimeScale};

    let file = fixture_sp3_file();
    let sp3 = Sp3::parse(file.as_bytes()).expect("parse");
    let g01 = GnssSatelliteId::new(GnssSystem::Gps, 1);

    // The parser will re-derive km->m and us->s, then position() reverses them;
    // because the fixture writes only 6 decimals, the *position values stored*
    // are the rounded-to-6-decimal km — which differ from the full-precision
    // POS_KM constants. So this test verifies the node axis + spline wiring is
    // correct end-to-end (J2000 seconds, per-axis, unit scaling), comparing the
    // public API against an independent in-test spline over the SAME rounded
    // node values, rather than re-asserting the 6-decimal-lossy fixture bits.
    // (The exact-bit 0-ULP gate is the direct-spline test above, which uses the
    // full-precision nodes.)

    // First node J2000 seconds must equal the fixture's first node.
    let ep0 = sp3.epochs[0];
    let secs0 = match ep0.repr {
        InstantRepr::JulianDate(s) => {
            ((s.jd_whole - 2_451_545.0) * 86_400.0 + s.fraction * 86_400.0).floor()
        }
        _ => panic!("expected JD epoch"),
    };
    assert_eq!(secs0, NODES[0], "J2000-second node axis mismatch");

    // Query at an interior epoch and confirm the public API agrees with an
    // in-test spline over the parsed (rounded) node values.
    let query_secs = NODES[0] + 437.0;
    let jd_whole = 2_451_545.0 + (query_secs / 86_400.0).floor();
    let frac = (query_secs - (query_secs / 86_400.0).floor() * 86_400.0) / 86_400.0;
    let q = Instant::from_julian_date(TimeScale::Gpst, JulianDateSplit::new(jd_whole, frac));

    let state = sp3.position(g01, q).expect("position");

    // Rebuild the parsed (rounded-to-6-decimal) node values and spline them.
    let parsed_km: Vec<[f64; 3]> = (0..9)
        .map(|i| {
            let st = sp3.state(g01, i).unwrap();
            [
                st.position.x_m / 1000.0,
                st.position.y_m / 1000.0,
                st.position.z_m / 1000.0,
            ]
        })
        .collect();
    let x = NODES.to_vec();
    let cx: Vec<f64> = parsed_km.iter().map(|r| r[0]).collect();
    let cy: Vec<f64> = parsed_km.iter().map(|r| r[1]).collect();
    let cz: Vec<f64> = parsed_km.iter().map(|r| r[2]).collect();
    let want_x = eval_spline(&x, &cx, query_secs) * 1000.0;
    let want_y = eval_spline(&x, &cy, query_secs) * 1000.0;
    let want_z = eval_spline(&x, &cz, query_secs) * 1000.0;

    assert_eq!(state.position.x_m.to_bits(), want_x.to_bits());
    assert_eq!(state.position.y_m.to_bits(), want_y.to_bits());
    assert_eq!(state.position.z_m.to_bits(), want_z.to_bits());
    assert!(state.clock_s.is_some(), "clock should interpolate");
}

/// 0-ULP parity against goldens generated from a **named public IGS SP3 file**
/// — the CNES/CLS/GRGS MGEX final product `GRG0MGXFIN_20201760000_01D_15M_ORB.SP3`
/// (GPS week 2111, 2020-06-24, SP3-c, 15-minute) — read by gnssanalysis
/// (`gn_io.sp3.read_sp3`) and interpolated with the committed scipy
/// `CubicSpline(x,y)` recipe (not-a-knot, extrapolate=True), per Cartesian axis
/// for position and 1-D for the clock, over the full 96-epoch GPS arc.
///
/// This is the real-data counterpart to
/// `position_clock_interpolation_is_0_ulp_vs_scipy` (synthetic 9-node arc):
/// same recipe, same 0-ULP bar, but the node substrate is genuine IGS orbit/clock
/// values for satellites G01/G15/G32 and the nodes are real 96-point arcs. Bits
/// embedded in `super::sp3_golden_data` (f64::to_bits, exact). Per-component
/// hex-float/ULP-distance assertion (`skyfield_parity_test.exs` discipline).
#[test]
fn real_igs_sp3_interpolation_is_0_ulp_vs_gnssanalysis_scipy() {
    use super::sp3_golden_data::GOLDEN_SATS;

    let bits_to_f64 = f64::from_bits;
    let mut checked = 0usize;

    for sat in GOLDEN_SATS {
        let n = sat.nodes_seconds_j2000.len();
        assert_eq!(
            sat.positions_km_bits.len(),
            n,
            "{} node/pos length",
            sat.prn
        );
        assert_eq!(sat.clocks_us_bits.len(), n, "{} node/clk length", sat.prn);

        let x: Vec<f64> = sat
            .nodes_seconds_j2000
            .iter()
            .map(|&b| bits_to_f64(b))
            .collect();
        let cx: Vec<f64> = sat
            .positions_km_bits
            .iter()
            .map(|r| bits_to_f64(r[0]))
            .collect();
        let cy: Vec<f64> = sat
            .positions_km_bits
            .iter()
            .map(|r| bits_to_f64(r[1]))
            .collect();
        let cz: Vec<f64> = sat
            .positions_km_bits
            .iter()
            .map(|r| bits_to_f64(r[2]))
            .collect();
        let clk: Vec<f64> = sat.clocks_us_bits.iter().map(|&b| bits_to_f64(b)).collect();

        for case in sat.cases {
            let q = bits_to_f64(case.query_seconds_j2000);

            let x_m = eval_spline(&x, &cx, q) * KM_TO_M;
            let y_m = eval_spline(&x, &cy, q) * KM_TO_M;
            let z_m = eval_spline(&x, &cz, q) * KM_TO_M;
            let c_s = eval_spline(&x, &clk, q) * US_TO_S;

            assert_eq!(
                x_m.to_bits(),
                case.position_m_bits[0],
                "{} X mismatch at q={}: got {:#018x} ({}), want {:#018x}",
                sat.prn,
                q,
                x_m.to_bits(),
                x_m,
                case.position_m_bits[0]
            );
            assert_eq!(
                y_m.to_bits(),
                case.position_m_bits[1],
                "{} Y mismatch at q={}: got {:#018x} ({}), want {:#018x}",
                sat.prn,
                q,
                y_m.to_bits(),
                y_m,
                case.position_m_bits[1]
            );
            assert_eq!(
                z_m.to_bits(),
                case.position_m_bits[2],
                "{} Z mismatch at q={}: got {:#018x} ({}), want {:#018x}",
                sat.prn,
                q,
                z_m.to_bits(),
                z_m,
                case.position_m_bits[2]
            );
            assert_eq!(
                c_s.to_bits(),
                case.clock_s_bits,
                "{} clock mismatch at q={}: got {:#018x} ({}), want {:#018x}",
                sat.prn,
                q,
                c_s.to_bits(),
                c_s,
                case.clock_s_bits
            );
            checked += 1;
        }
    }

    // Guard against an empty/short golden silently passing.
    assert!(
        checked >= 15,
        "expected >= 15 real-file cases, checked {}",
        checked
    );
}

/// End-to-end 0-ULP gate: read the named public IGS SP3 file, run it through
/// the production `Sp3::parse`, and drive the public `Sp3::position` API — then
/// assert bit-equality against the gnssanalysis/scipy golden.
///
/// The sibling test above feeds gnssanalysis-extracted node bits straight into
/// the spline, so it proves the interpolation recipe but not the parser. This
/// test exercises the whole path: it first checks the parser's node substrate
/// (the J2000-second axis and the native km/us node values the interpolator
/// fits) is bit-identical to the golden, then checks the full parse->position
/// result is 0-ULP.
#[test]
fn real_igs_file_through_parse_is_0_ulp_end_to_end() {
    use super::super::Sp3;
    use super::instant_to_j2000_seconds_for_test as inst_secs;
    use super::sp3_golden_data::GOLDEN_SATS;
    use crate::id::{GnssSatelliteId, GnssSystem};
    use astrodynamics::time::model::{Instant, JulianDateSplit, TimeScale};

    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 real IGS SP3 fixture {path}: {e}"));
    let sp3 = Sp3::parse(&bytes).expect("parse real IGS SP3");

    let bits = f64::from_bits;
    let mut checked = 0usize;

    for sat in GOLDEN_SATS {
        let prn: u8 = sat.prn[1..].parse().expect("golden prn");
        let id = GnssSatelliteId::new(GnssSystem::Gps, prn);
        let n = sat.nodes_seconds_j2000.len();

        // (1) Parser fidelity: gather this sat's NATIVE-unit nodes (the exact
        // km the interpolator fits, straight from the ASCII) and prove they
        // equal the golden node substrate bit-for-bit. We assert the native km
        // here, NOT `position.x_m / 1000.0`: the public meters are `km * 1000`
        // and reversing them drifts up to 1 ULP by design — that round trip is
        // precisely why interpolation must use the native values.
        let mut secs = Vec::with_capacity(n);
        let mut kx = Vec::with_capacity(n);
        let mut ky = Vec::with_capacity(n);
        let mut kz = Vec::with_capacity(n);
        for (idx, ep) in sp3.epochs.iter().enumerate() {
            let (km, _clk_us) = match sp3.interp_node_for_test(id, idx) {
                Some(node) => node,
                None => continue,
            };
            secs.push(inst_secs(ep).expect("epoch -> J2000 sec"));
            kx.push(km[0]);
            ky.push(km[1]);
            kz.push(km[2]);
        }

        assert_eq!(
            secs.len(),
            n,
            "{}: parser yielded {} position nodes, golden has {}",
            sat.prn,
            secs.len(),
            n
        );
        for i in 0..n {
            assert_eq!(
                secs[i].to_bits(),
                sat.nodes_seconds_j2000[i],
                "{} J2000-second node[{}]: parser {} vs golden {}",
                sat.prn,
                i,
                secs[i],
                bits(sat.nodes_seconds_j2000[i])
            );
            assert_eq!(
                kx[i].to_bits(),
                sat.positions_km_bits[i][0],
                "{} km-X node[{}] after parse round-trip: got {:#018x} ({}), want {:#018x}",
                sat.prn,
                i,
                kx[i].to_bits(),
                kx[i],
                sat.positions_km_bits[i][0]
            );
            assert_eq!(
                ky[i].to_bits(),
                sat.positions_km_bits[i][1],
                "{} km-Y node[{}] after parse round-trip",
                sat.prn,
                i
            );
            assert_eq!(
                kz[i].to_bits(),
                sat.positions_km_bits[i][2],
                "{} km-Z node[{}] after parse round-trip",
                sat.prn,
                i
            );
        }

        // (2) End-to-end: drive the public position() API at each golden epoch.
        for case in sat.cases {
            let q = bits(case.query_seconds_j2000); // J2000 second (may be sub-second)
                                                    // Construct a split-JD Instant that converts to EXACTLY q (the query
                                                    // is evaluated un-floored; a sub-second case like 646272312.5 must
                                                    // survive the Instant boundary bit-for-bit, which it does for these
                                                    // golden epochs).
            let day = (q / 86_400.0).floor();
            let rem = q - day * 86_400.0;
            let epoch = Instant::from_julian_date(
                TimeScale::Gpst,
                JulianDateSplit::new(2_451_545.0 + day, rem / 86_400.0),
            );
            assert_eq!(
                inst_secs(&epoch).unwrap(),
                q,
                "{}: constructed query Instant does not reproduce golden second {} exactly",
                sat.prn,
                q
            );

            let state = sp3.position(id, epoch).expect("position via public API");
            assert_eq!(
                state.position.x_m.to_bits(),
                case.position_m_bits[0],
                "{} end-to-end X at q={}: got {:#018x} ({}), want {:#018x}",
                sat.prn,
                q,
                state.position.x_m.to_bits(),
                state.position.x_m,
                case.position_m_bits[0]
            );
            assert_eq!(
                state.position.y_m.to_bits(),
                case.position_m_bits[1],
                "{} end-to-end Y at q={}",
                sat.prn,
                q
            );
            assert_eq!(
                state.position.z_m.to_bits(),
                case.position_m_bits[2],
                "{} end-to-end Z at q={}",
                sat.prn,
                q
            );
            assert_eq!(
                state.clock_s.map(|c| c.to_bits()),
                Some(case.clock_s_bits),
                "{} end-to-end clock at q={}: got {:?}, want {:#018x}",
                sat.prn,
                q,
                state.clock_s.map(|c| c.to_bits()),
                case.clock_s_bits
            );
            checked += 1;
        }
    }

    assert!(
        checked >= 15,
        "expected >= 15 end-to-end cases, checked {}",
        checked
    );
}

#[test]
fn position_clock_interpolation_is_0_ulp_vs_scipy() {
    let x = NODES.to_vec();
    let cx = col(0).to_vec();
    let cy = col(1).to_vec();
    let cz = col(2).to_vec();
    let clk = CLK_US.to_vec();

    for case in &CASES {
        // Position: per-axis km spline, then single * 1000.0 to meters.
        let x_m = eval_spline(&x, &cx, case.query) * KM_TO_M;
        let y_m = eval_spline(&x, &cy, case.query) * KM_TO_M;
        let z_m = eval_spline(&x, &cz, case.query) * KM_TO_M;

        assert_eq!(
            x_m.to_bits(),
            case.pos_m_bits[0],
            "X mismatch at q={}: got {:#018x} ({}), want {:#018x}",
            case.query,
            x_m.to_bits(),
            x_m,
            case.pos_m_bits[0]
        );
        assert_eq!(
            y_m.to_bits(),
            case.pos_m_bits[1],
            "Y mismatch at q={}: got {:#018x} ({}), want {:#018x}",
            case.query,
            y_m.to_bits(),
            y_m,
            case.pos_m_bits[1]
        );
        assert_eq!(
            z_m.to_bits(),
            case.pos_m_bits[2],
            "Z mismatch at q={}: got {:#018x} ({}), want {:#018x}",
            case.query,
            z_m.to_bits(),
            z_m,
            case.pos_m_bits[2]
        );

        // Clock: us spline, then single * 1e-6 to seconds.
        let c_s = eval_spline(&x, &clk, case.query) * US_TO_S;
        assert_eq!(
            c_s.to_bits(),
            case.clk_s_bits,
            "clock mismatch at q={}: got {:#018x} ({}), want {:#018x}",
            case.query,
            c_s.to_bits(),
            c_s,
            case.clk_s_bits
        );
    }
}

/// Pins the proven exactness bound of `instant_to_j2000_seconds` for the SP3
/// query path.
///
/// The conversion is `q = (jd_whole - J2000_JD) * 86400 + fraction * 86400`,
/// where the parser builds `jd_whole = jdn - 0.5` (an exact integer minus 0.5)
/// and `fraction = day_seconds / 86400`. See `docs/gnss/time-recipe.md` §9.
///
/// Two facts, proven here for the SP3 epoch domain:
///
/// 1. The midnight term `(jd_whole - J2000_JD) * 86400` is an EXACT integer
///    (no rounding) for every civil day in the SP3 domain (1994..2100): the day
///    index fits well under 2^53, and `(k - 0.5) * 86400 = k*86400 - 43200` is
///    integer-valued and exactly representable.
/// 2. The sub-second term `(day_seconds / 86400) * 86400` reproduces
///    `day_seconds` to within ~7.28e-12 s (0.0073 ns) over the entire
///    `[0, 86400)` second-of-day range — well under one ULP at magnitude 86400.
///
/// The honest LIMIT (documented, not over-claimed): the *final* float sum
/// `midnight_integer + sub_second` is then rounded to the f64 grid at the
/// J2000-seconds magnitude (~1e8..3e9), whose ULP is ~0.03..0.48 us. That
/// granularity is inherent to representing seconds-since-J2000 as a single f64
/// and is NOT a defect of this conversion; it bounds any consumer of the
/// scalar-seconds form. Integer-second nodes are recovered exactly by the
/// `.floor()` on the node axis (proven by the end-to-end 0-ULP test above).
#[test]
fn instant_to_j2000_seconds_query_exactness_bound() {
    const D: f64 = 86_400.0;
    const J2000_JD: f64 = 2_451_545.0;

    // (1) Midnight term is an exact integer across the SP3 epoch domain.
    // jdn for 1994..2100; step by a few representative days incl. boundaries.
    fn jdn(y: i64, m: i64, d: i64) -> i64 {
        let a = (14 - m) / 12;
        let yy = y + 4800 - a;
        let mm = m + 12 * a - 3;
        d + (153 * mm + 2) / 5 + 365 * yy + yy / 4 - yy / 100 + yy / 400 - 32_045
    }
    let lo = jdn(1994, 1, 1);
    let hi = jdn(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;
        // Exact integer: fract() is bit-zero, and it equals the i64 path.
        assert_eq!(
            midnight.fract(),
            0.0,
            "midnight term not integer for jdn {jd}"
        );
        let exact = (jd - 2_451_545) * 86_400 - 43_200;
        assert_eq!(
            midnight, exact as f64,
            "midnight term != exact i64 for jdn {jd}"
        );
    }

    // (2) Sub-second round-trip bound over [0, 86400): exhaustive integer-second
    // grid plus fine sub-second grids. Worst case is ~0.0073 ns (< 1 ULP @ 86400).
    const BOUND_S: f64 = 7.3e-12; // proven empirical max is 7.275957614183426e-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"
    );
}