astrodynamics-gnss 0.14.0

//! RINEX NAV parser coverage against the committed multi-GNSS fixture.
//!
//! Parsing is deterministic byte-to-record translation, so these are round-trip
//! and schema assertions (counts, field ranges, message classification, a
//! physical sanity check via the evaluator), NOT a 0-ULP parity claim.

use super::*;
use crate::broadcast::satellite_state;

fn fixture_text() -> String {
    let path = concat!(
        env!("CARGO_MANIFEST_DIR"),
        "/tests/fixtures/nav/ESBC00DNK_R_20201770000_01D_MN.rnx"
    );
    std::fs::read_to_string(path).unwrap_or_else(|e| panic!("read NAV fixture {path}: {e}"))
}

fn records() -> Vec<BroadcastRecord> {
    parse_nav(&fixture_text()).expect("parse NAV fixture")
}

fn v4_fixture_text() -> String {
    let path = concat!(
        env!("CARGO_MANIFEST_DIR"),
        "/tests/fixtures/nav/KMS300DNK_R_20221591000_01H_MN.rnx"
    );
    std::fs::read_to_string(path).unwrap_or_else(|e| panic!("read v4 NAV fixture {path}: {e}"))
}

fn glonass_fixture_text() -> String {
    let path = concat!(
        env!("CARGO_MANIFEST_DIR"),
        "/tests/fixtures/nav/ESBC00DNK_R_20201770000_01D_RN.rnx"
    );
    std::fs::read_to_string(path).unwrap_or_else(|e| panic!("read GLONASS fixture {path}: {e}"))
}

#[test]
fn parses_and_evaluates_glonass_records() {
    use crate::spp::EphemerisSource;

    let text = glonass_fixture_text();
    let recs = parse_glonass(&text).expect("parse GLONASS records");
    assert_eq!(recs.len(), 510, "GLONASS record count");
    assert_eq!(
        parse_leap_seconds(&text),
        Some(18.0),
        "GPS-UTC leap seconds"
    );

    // Every record's broadcast state sits on the GLONASS orbit (~25,510 km).
    for r in &recs {
        let radius_km =
            (r.pos_m[0].powi(2) + r.pos_m[1].powi(2) + r.pos_m[2].powi(2)).sqrt() / 1000.0;
        assert!(
            (25_000.0..26_000.0).contains(&radius_km),
            "{:?} GLONASS radius {radius_km} km out of band",
            r.satellite_id
        );
    }

    // The store evaluates a GLONASS satellite through the RK4 propagator. At the
    // record's own reference epoch (tk = 0) the position is the broadcast state,
    // so the radius is the GLONASS orbit radius.
    let store = BroadcastStore::from_nav(&text).expect("parse GLONASS NAV");
    assert_eq!(store.glonass_records().len(), 510);
    let r0 = store.glonass_records()[0];
    let t_toe_gpst = r0.toe_utc_j2000_s + 18.0; // leap seconds for 2020
    let (pos, _clk) = store
        .position_clock_at_j2000_s(r0.satellite_id, t_toe_gpst)
        .expect("GLONASS position at its toe");
    let radius_km = (pos[0].powi(2) + pos[1].powi(2) + pos[2].powi(2)).sqrt() / 1000.0;
    assert!(
        (25_000.0..26_000.0).contains(&radius_km),
        "evaluated GLONASS radius {radius_km} km out of band"
    );
    // tk = 0 means no integration, so the evaluated position equals the state.
    assert_eq!(
        [pos[0], pos[1], pos[2]],
        r0.pos_m,
        "tk=0 returns the broadcast state"
    );

    // A query far outside the product's coverage (a day before any record) has
    // no record within the +/-15 min validity window, so no ephemeris. (A query
    // an hour later would instead be served by the next half-hourly record.)
    assert!(
        store
            .position_clock_at_j2000_s(r0.satellite_id, t_toe_gpst - 86_400.0)
            .is_none(),
        "a query a day before any record is outside every validity window"
    );
}

#[test]
fn spp_solves_from_broadcast_glonass() {
    use crate::spp::{
        solve, test_support, Corrections, KlobucharCoeffs, Observation, SolveInputs, SurfaceMet,
        ELEVATION_MASK_RAD,
    };

    // GLONASS-only store from the RN fixture.
    let store = BroadcastStore::from_nav(&glonass_fixture_text()).expect("parse GLONASS NAV");

    // 2020-06-25 12:00 GPST, mid-day so GLONASS satellites have a near-epoch
    // record. The ionosphere correction is unsupported for GLONASS (no modeled
    // single-frequency carrier), so this geometry-only solve leaves it off.
    let t_rx = 646_358_400.0_f64;
    let sod = 12.0 * 3600.0;
    let doy = 177.0;
    let x_true = [3_512_900.0, 780_500.0, 5_248_700.0, 0.0];
    let corr = Corrections::NONE;
    let kl = KlobucharCoeffs {
        alpha: [0.0; 4],
        beta: [0.0; 4],
    };
    let met = SurfaceMet {
        pressure_hpa: 1013.25,
        temperature_k: 288.15,
        relative_humidity: 0.5,
    };

    let mut sats: Vec<_> = store
        .glonass_records()
        .iter()
        .map(|r| r.satellite_id)
        .collect();
    sats.sort_unstable();
    sats.dedup();
    let mut observations = Vec::new();
    for sat in sats {
        if let Some(m) = test_support::sat_model_for_test(
            &store,
            sat,
            [x_true[0], x_true[1], x_true[2]],
            x_true[3],
            20_000_000.0,
            t_rx,
            sod,
            doy,
            corr,
            &kl,
            &met,
        ) {
            if m.el_rad >= ELEVATION_MASK_RAD {
                observations.push(Observation {
                    satellite_id: sat,
                    pseudorange_m: m.p_hat_m,
                });
            }
        }
    }
    assert!(
        observations.len() >= 4,
        "need >=4 visible GLONASS sats, got {}",
        observations.len()
    );

    let inputs = SolveInputs {
        observations,
        t_rx_j2000_s: t_rx,
        t_rx_second_of_day_s: sod,
        day_of_year: doy,
        initial_guess: [
            x_true[0] + 1000.0,
            x_true[1] - 1000.0,
            x_true[2] + 1000.0,
            0.0,
        ],
        corrections: corr,
        klobuchar: kl,
        beidou_klobuchar: None,
        met,
    };

    let sol = solve(&store, &inputs, true).expect("GLONASS broadcast SPP solve");
    let p = sol.position;
    let err =
        ((p.x_m - x_true[0]).powi(2) + (p.y_m - x_true[1]).powi(2) + (p.z_m - x_true[2]).powi(2))
            .sqrt();
    assert!(err < 1.0e-3, "recovered position off by {err} m");
    // A single-system GLONASS solve carries one receiver clock.
    assert_eq!(sol.system_clocks_s.len(), 1, "one GLONASS clock");
    assert_eq!(sol.system_clocks_s[0].0, GnssSystem::Glonass);
}

#[test]
fn beidou_uses_its_own_klobuchar_coefficients() {
    use crate::spp::{
        solve, test_support, Corrections, KlobucharCoeffs, Observation, SolveInputs, SurfaceMet,
        ELEVATION_MASK_RAD,
    };

    // BeiDou-only store.
    let store = BroadcastStore::new(
        records()
            .into_iter()
            .filter(|r| r.satellite_id.system == GnssSystem::BeiDou)
            .collect(),
    );
    let t_rx = 646_358_400.0_f64;
    let sod = 12.0 * 3600.0;
    let doy = 177.0;
    let x_true = [3_512_900.0, 780_500.0, 5_248_700.0];
    // The broadcast BeiDou Klobuchar-8 set (BDSA/BDSB).
    let bds = KlobucharCoeffs {
        alpha: [1.1180e-08, 2.9800e-08, -4.1720e-07, 6.5570e-07],
        beta: [1.4130e05, -5.2430e05, 1.6380e06, -4.5880e05],
    };
    let met = SurfaceMet {
        pressure_hpa: 1013.25,
        temperature_k: 288.15,
        relative_humidity: 0.5,
    };

    // Synthesize BeiDou observations with the ionosphere applied using the BeiDou
    // coefficients (sat_model scales the L1 delay to B1I for BeiDou).
    let mut sats: Vec<_> = store.records().iter().map(|r| r.satellite_id).collect();
    sats.sort_unstable();
    sats.dedup();
    let mut observations = Vec::new();
    for sat in sats {
        if let Some(m) = test_support::sat_model_for_test(
            &store,
            sat,
            x_true,
            0.0,
            22_000_000.0,
            t_rx,
            sod,
            doy,
            Corrections::IONO,
            &bds,
            &met,
        ) {
            if m.el_rad >= ELEVATION_MASK_RAD {
                observations.push(Observation {
                    satellite_id: sat,
                    pseudorange_m: m.p_hat_m,
                });
            }
        }
    }
    assert!(
        observations.len() >= 4,
        "need >=4 BeiDou sats, got {}",
        observations.len()
    );

    let base = |beidou_klobuchar| SolveInputs {
        observations: observations.clone(),
        t_rx_j2000_s: t_rx,
        t_rx_second_of_day_s: sod,
        day_of_year: doy,
        initial_guess: [
            x_true[0] + 1000.0,
            x_true[1] - 1000.0,
            x_true[2] + 1000.0,
            0.0,
        ],
        corrections: Corrections::IONO,
        // Zero GPS-side coefficients: if BeiDou wrongly used these, no ionosphere
        // would be applied and the synthesized delay would bias the solution.
        klobuchar: KlobucharCoeffs {
            alpha: [0.0; 4],
            beta: [0.0; 4],
        },
        beidou_klobuchar,
        met,
    };

    // With the BeiDou coefficients supplied, BeiDou uses them and the truth is
    // recovered (the applied ionosphere matches the synthesized one).
    let sol = solve(&store, &base(Some(bds)), false).expect("BeiDou-native iono solve");
    let p = sol.position;
    let err =
        ((p.x_m - x_true[0]).powi(2) + (p.y_m - x_true[1]).powi(2) + (p.z_m - x_true[2]).powi(2))
            .sqrt();
    assert!(
        err < 1.0e-3,
        "with BDSA/BDSB the solve recovers; off by {err} m"
    );

    // Without them, BeiDou falls back to the (zero) shared set, so the modelled
    // ionosphere is missing and the solution is biased — proving the per-system
    // coefficients are actually used.
    let sol0 = solve(&store, &base(None), false).expect("fallback solve");
    let p0 = sol0.position;
    let err0 = ((p0.x_m - x_true[0]).powi(2)
        + (p0.y_m - x_true[1]).powi(2)
        + (p0.z_m - x_true[2]).powi(2))
    .sqrt();
    assert!(
        err0 > 0.1,
        "without BeiDou coeffs the unmodelled ionosphere biases the fix; off by {err0} m"
    );
}

fn brdc_gop_text() -> String {
    let path = concat!(
        env!("CARGO_MANIFEST_DIR"),
        "/tests/fixtures/nav/BRDC00GOP_R_20210010000_01D_MN.rnx"
    );
    std::fs::read_to_string(path).unwrap_or_else(|e| panic!("read BRDC00GOP fixture {path}: {e}"))
}

#[test]
fn parses_broadcast_ionosphere_coefficients() {
    // The main fixture carries GPS (and Galileo NeQuick) coefficients but no
    // BeiDou set.
    let esbc = parse_iono_corrections(&fixture_text());
    let gps = esbc.gps.expect("ESBC has GPSA/GPSB");
    assert!(
        (gps.alpha[0] - 4.6566e-09).abs() < 1e-19,
        "GPSA a0 {}",
        gps.alpha[0]
    );
    assert!(
        (gps.beta[0] - 8.1920e04).abs() < 1e-3,
        "GPSB b0 {}",
        gps.beta[0]
    );
    assert!(esbc.beidou.is_none(), "ESBC has no BDSA/BDSB");

    // The merged BRDC header carries the BeiDou Klobuchar-8 set.
    let brdc = parse_iono_corrections(&brdc_gop_text());
    let bds = brdc.beidou.expect("BRDC00GOP has BDSA/BDSB");
    assert!(
        (bds.alpha[0] - 1.1180e-08).abs() < 1e-18,
        "BDSA a0 {}",
        bds.alpha[0]
    );
    assert!(
        (bds.alpha[2] - -4.1720e-07).abs() < 1e-17,
        "BDSA a2 {}",
        bds.alpha[2]
    );
    assert!(
        (bds.beta[0] - 1.4130e05).abs() < 1e-3,
        "BDSB b0 {}",
        bds.beta[0]
    );
    assert!(
        (bds.beta[1] - -5.2430e05).abs() < 1e-3,
        "BDSB b1 {}",
        bds.beta[1]
    );
    assert!(brdc.gps.is_some(), "BRDC00GOP also has GPSA/GPSB");
}

#[test]
fn broadcast_store_exposes_header_ionosphere_coefficients() {
    // from_nav captures the header coefficients; new() leaves them empty.
    let store = BroadcastStore::from_nav(&brdc_gop_text()).expect("parse BRDC00GOP");
    assert!(
        store.iono_corrections().beidou.is_some(),
        "BeiDou coeffs from header"
    );

    let bare = BroadcastStore::new(vec![]);
    assert_eq!(
        bare.iono_corrections(),
        Default::default(),
        "new() has no coeffs"
    );
}

#[test]
fn parses_a_real_rinex_v4_file() {
    let recs = parse_nav(&v4_fixture_text()).expect("parse v4 NAV fixture");
    let count = |sys| recs.iter().filter(|r| r.satellite_id.system == sys).count();
    let msg = |m| recs.iter().filter(|r| r.message == m).count();

    // Supported Keplerian records only: GPS LNAV, Galileo I/NAV + F/NAV, BeiDou
    // D1 + D2. GLONASS (FDMA), QZSS, SBAS, STO and ION frames are skipped.
    assert_eq!(count(GnssSystem::Gps), 30, "GPS LNAV count");
    assert_eq!(count(GnssSystem::Galileo), 108, "Galileo count");
    assert_eq!(count(GnssSystem::BeiDou), 36, "BeiDou count");
    assert_eq!(recs.len(), 174, "only G/E/C are parsed");
    assert_eq!(
        count(GnssSystem::Glonass) + count(GnssSystem::Qzss) + count(GnssSystem::Sbas),
        0,
        "GLONASS/QZSS/SBAS must be skipped"
    );

    // Message type comes from the v4 marker token.
    assert_eq!(msg(NavMessage::GpsLnav), 30);
    assert_eq!(msg(NavMessage::GalileoInav), 55);
    assert_eq!(msg(NavMessage::GalileoFnav), 53);
    assert_eq!(msg(NavMessage::BeidouD1), 33);
    assert_eq!(msg(NavMessage::BeidouD2), 3);

    // Parsed records evaluate to physical orbit radii (parser-to-evaluator sanity
    // on real v4 bytes), MEO/IGSO/GEO bands across the constellations.
    for sys in [GnssSystem::Gps, GnssSystem::Galileo, GnssSystem::BeiDou] {
        let r = recs.iter().find(|r| r.satellite_id.system == sys).unwrap();
        let st = satellite_state(
            &r.elements,
            &r.clock,
            &r.constants(),
            r.elements.toe_sow,
            r.group_delay_s,
            crate::rinex_nav::is_beidou_geo(r.satellite_id),
        );
        let p = st.orbit.position();
        let radius_km = (p.x_m * p.x_m + p.y_m * p.y_m + p.z_m * p.z_m).sqrt() / 1000.0;
        assert!(
            (20_000.0..50_000.0).contains(&radius_km),
            "{sys:?} v4 radius {radius_km} km out of band"
        );
    }
}

#[test]
fn parses_gps_galileo_and_beidou_records() {
    let recs = records();
    let count = |sys| recs.iter().filter(|r| r.satellite_id.system == sys).count();
    let gps = count(GnssSystem::Gps);
    let gal = count(GnssSystem::Galileo);
    let bds = count(GnssSystem::BeiDou);
    // The committed fixture is filtered to GPS + Galileo + BeiDou.
    assert_eq!(gps, 257, "GPS record count");
    assert_eq!(gal, 1602, "Galileo record count");
    assert_eq!(bds, 357, "BeiDou record count");
    assert_eq!(
        recs.len(),
        gps + gal + bds,
        "only GPS+Galileo+BeiDou are returned"
    );
}

#[test]
fn gps_record_fields_are_in_range() {
    let recs = records();
    let g01 = recs
        .iter()
        .find(|r| r.satellite_id == GnssSatelliteId::new(GnssSystem::Gps, 1))
        .expect("a G01 record");

    assert_eq!(g01.message, NavMessage::GpsLnav);
    assert_eq!(g01.week, 2111, "GPS week 2111 for this product");
    // GPS semi-major axis ~26560 km => sqrt(a) ~ 5153.6 sqrt(m).
    assert!(
        (5100.0..5200.0).contains(&g01.elements.sqrt_a),
        "sqrt_a {}",
        g01.elements.sqrt_a
    );
    assert!(
        (0.0..0.05).contains(&g01.elements.e),
        "e {}",
        g01.elements.e
    );
    // For this record the clock and ephemeris reference epochs coincide.
    assert_eq!(g01.clock.toc_sow, g01.elements.toe_sow);
    assert_eq!(g01.sv_health, 0.0, "G01 is healthy");
    assert!(g01.group_delay_s.abs() < 1.0e-6, "TGD is a small delay");
}

#[test]
fn galileo_messages_are_classified() {
    let recs = records();
    let gal: Vec<_> = recs
        .iter()
        .filter(|r| r.satellite_id.system == GnssSystem::Galileo)
        .collect();
    let inav = gal
        .iter()
        .filter(|r| r.message == NavMessage::GalileoInav)
        .count();
    let fnav = gal
        .iter()
        .filter(|r| r.message == NavMessage::GalileoFnav)
        .count();
    assert_eq!(inav, 821, "Galileo I/NAV record count");
    assert_eq!(fnav, 781, "Galileo F/NAV record count");
    assert_eq!(inav + fnav, gal.len(), "every Galileo record is classified");
}

#[test]
fn parsed_records_evaluate_to_physical_orbit_radii() {
    let recs = records();
    // Evaluate each constellation's first record at its toe and check the ECEF
    // radius is in the expected MEO band (parser-to-evaluator sanity).
    for (system, lo_km, hi_km) in [
        (GnssSystem::Gps, 25_000.0, 27_500.0),
        (GnssSystem::Galileo, 29_000.0, 30_500.0),
    ] {
        let r = recs
            .iter()
            .find(|r| r.satellite_id.system == system)
            .expect("a record");
        let state = satellite_state(
            &r.elements,
            &r.clock,
            &r.constants(),
            r.elements.toe_sow,
            r.group_delay_s,
            false,
        );
        let p = state.orbit.position();
        let radius_km = (p.x_m * p.x_m + p.y_m * p.y_m + p.z_m * p.z_m).sqrt() / 1000.0;
        assert!(
            (lo_km..hi_km).contains(&radius_km),
            "{system:?} radius {radius_km} km out of band"
        );
    }
}

#[test]
fn spp_solves_from_broadcast_gps() {
    use crate::spp::{
        solve, test_support, Corrections, KlobucharCoeffs, Observation, SolveInputs, SurfaceMet,
        ELEVATION_MASK_RAD,
    };

    // GPS-only store (avoids any Galileo I/NAV vs F/NAV selection ambiguity).
    let store = BroadcastStore::new(
        records()
            .into_iter()
            .filter(|r| r.satellite_id.system == GnssSystem::Gps)
            .collect(),
    );

    // 2020-06-25 12:00 GPST (DOY 177 noon), as a J2000 second; mid-day so every
    // GPS satellite has a near-toe record.
    let t_rx = 646_358_400.0_f64;
    let sod = 12.0 * 3600.0;
    let doy = 177.0;
    // A true receiver on the ground near the ESBC station (Esbjerg, Denmark).
    let x_true = [3_512_900.0, 780_500.0, 5_248_700.0, 0.0];
    let corr = Corrections::NONE;
    let kl = KlobucharCoeffs {
        alpha: [0.0; 4],
        beta: [0.0; 4],
    };
    let met = SurfaceMet {
        pressure_hpa: 1013.25,
        temperature_k: 288.15,
        relative_humidity: 0.5,
    };

    // Synthesize one pseudorange per visible GPS satellite with the same forward
    // model the solver inverts, so the true state is the zero-residual solution.
    let mut sats: Vec<_> = store.records().iter().map(|r| r.satellite_id).collect();
    sats.sort_unstable();
    sats.dedup();
    let mut observations = Vec::new();
    for sat in sats {
        if let Some(m) = test_support::sat_model_for_test(
            &store,
            sat,
            [x_true[0], x_true[1], x_true[2]],
            x_true[3],
            22_000_000.0,
            t_rx,
            sod,
            doy,
            corr,
            &kl,
            &met,
        ) {
            if m.el_rad >= ELEVATION_MASK_RAD {
                observations.push(Observation {
                    satellite_id: sat,
                    pseudorange_m: m.p_hat_m,
                });
            }
        }
    }
    assert!(
        observations.len() >= 4,
        "need >=4 visible GPS sats, got {}",
        observations.len()
    );

    let inputs = SolveInputs {
        observations,
        t_rx_j2000_s: t_rx,
        t_rx_second_of_day_s: sod,
        day_of_year: doy,
        initial_guess: [
            x_true[0] + 1000.0,
            x_true[1] - 1000.0,
            x_true[2] + 1000.0,
            0.0,
        ],
        corrections: corr,
        klobuchar: kl,
        beidou_klobuchar: None,
        met,
    };

    let sol = solve(&store, &inputs, true).expect("broadcast SPP solve");
    let p = sol.position;
    let err =
        ((p.x_m - x_true[0]).powi(2) + (p.y_m - x_true[1]).powi(2) + (p.z_m - x_true[2]).powi(2))
            .sqrt();
    assert!(err < 1.0e-3, "recovered position off by {err} m");
}

#[test]
fn rejects_a_non_navigation_header() {
    let bogus = "     3.05           OBSERVATION DATA   M                   RINEX VERSION / TYPE\n\
                 END OF HEADER\n";
    assert!(matches!(
        parse_nav(bogus),
        Err(NavParseError::UnsupportedHeader(_))
    ));
}

#[test]
fn reports_missing_header_end() {
    let truncated =
        "     3.05           NAVIGATION DATA     M                   RINEX VERSION / TYPE\n";
    assert_eq!(parse_nav(truncated), Err(NavParseError::MissingHeaderEnd));
}

// An exact GPS LNAV record block copied from the committed fixture (the v3 parser
// is already proven on this data), reused to build inline v3 and v4 inputs so the
// v4 path can be cross-checked against the v3 result. Continuation lines keep
// their fixed-column leading spaces.
const G01_LINES: &[&str] = &[
    "G01 2020 06 25 04 00 00 1.604342833161e-05 7.048583938740e-12 0.000000000000e+00",
    "     5.800000000000e+01-3.968750000000e+01 4.304822170265e-09 6.342094507864e-01",
    "    -2.177432179451e-06 1.000394229777e-02 1.937150955200e-06 5.153707128525e+03",
    "     3.600000000000e+05-1.508742570877e-07 2.572838528869e+00 1.359730958939e-07",
    "     9.806518601091e-01 3.539687500000e+02 7.941703015008e-01-8.384634967987e-09",
    "    -5.714523747137e-11 1.000000000000e+00 2.111000000000e+03 0.000000000000e+00",
    "     2.000000000000e+00 0.000000000000e+00 5.122274160385e-09 5.800000000000e+01",
    "     3.561060000000e+05 4.000000000000e+00",
];

// An exact Galileo record block whose data-source word (orbit-5 field 2 = 258,
// bit 0x100 set) infers F/NAV under the v3 rule, used to show the v4 marker token
// is authoritative over that inference.
const E01_LINES: &[&str] = &[
    "E01 2020 06 24 23 30 00-8.846927667037e-04-7.972289495228e-12 0.000000000000e+00",
    "     6.100000000000e+01 1.865625000000e+01 2.656539226950e-09-1.832282909549e+00",
    "     8.568167686462e-07 9.650341235101e-05 1.049041748047e-05 5.440602037430e+03",
    "     3.438000000000e+05 1.862645149231e-09 2.123282284601e-01-1.452863216400e-07",
    "     9.828296477370e-01 1.298750000000e+02-2.778709093141e+00-5.216288707934e-09",
    "    -6.996720012901e-10 2.580000000000e+02 2.111000000000e+03",
    "     3.120000000000e+00 0.000000000000e+00-1.862645149231e-09 0.000000000000e+00",
    "     3.445400000000e+05",
];

const V4_NAV_HEADER: &str =
    "     4.00           NAVIGATION DATA     M                   RINEX VERSION / TYPE\n\
     XXX                                                         END OF HEADER\n";

const V3_NAV_HEADER: &str =
    "     3.05           NAVIGATION DATA     M                   RINEX VERSION / TYPE\n\
     XXX                                                         END OF HEADER\n";

fn join(lines: &[&str]) -> String {
    let mut s = lines.join("\n");
    s.push('\n');
    s
}

#[test]
fn parses_rinex_v4_eph_frames_and_skips_the_rest() {
    // One GPS LNAV EPH frame, plus a CNAV EPH frame and an STO frame that must be
    // skipped (CNAV reorders the orbit columns; STO carries no ephemeris).
    let mut text = String::from(V4_NAV_HEADER);
    text.push_str("> EPH G01 LNAV\n");
    text.push_str(&join(G01_LINES));
    text.push_str("> EPH G03 CNAV\n");
    text.push_str(&join(G01_LINES)); // body content is irrelevant; it must be skipped, not parsed
    text.push_str("> STO G01 LNAV\n");
    text.push_str("    2020 06 25 00 00 00 GPUT 0.0 0.0 0 0\n");

    let recs = parse_nav(&text).expect("parse v4 NAV");
    assert_eq!(recs.len(), 1, "only the LNAV EPH frame is parsed");
    assert_eq!(
        recs[0].satellite_id,
        GnssSatelliteId::new(GnssSystem::Gps, 1)
    );
    assert_eq!(recs[0].message, NavMessage::GpsLnav);

    // The v4 record must equal the same block parsed as v3, field for field.
    let v3_text = format!("{V3_NAV_HEADER}{}", join(G01_LINES));
    let v3 = parse_nav(&v3_text).expect("parse v3 NAV");
    assert_eq!(v3.len(), 1);
    assert_eq!(recs[0].elements, v3[0].elements, "elements differ v4 vs v3");
    assert_eq!(recs[0].clock, v3[0].clock, "clock differs v4 vs v3");
    assert_eq!(recs[0].week, v3[0].week);
    assert_eq!(recs[0].fit_interval_s, v3[0].fit_interval_s);
}

#[test]
fn rinex_v4_message_type_comes_from_the_marker() {
    // The v3 data-source-word rule infers F/NAV for this Galileo block...
    let v3_text = format!("{V3_NAV_HEADER}{}", join(E01_LINES));
    let v3 = parse_nav(&v3_text).expect("parse v3 NAV");
    assert_eq!(
        v3[0].message,
        NavMessage::GalileoFnav,
        "v3 infers F/NAV here"
    );

    // ...but a v4 marker that says INAV is authoritative.
    let v4_text = format!("{V4_NAV_HEADER}> EPH E01 INAV\n{}", join(E01_LINES));
    let v4 = parse_nav(&v4_text).expect("parse v4 NAV");
    assert_eq!(v4.len(), 1);
    assert_eq!(
        v4[0].message,
        NavMessage::GalileoInav,
        "v4 message must come from the marker token, not the data-source word"
    );
}

#[test]
fn accepts_v4_nav_header_rejects_v4_non_nav() {
    // A 4.00 NAV header with one frame parses.
    let ok = format!("{V4_NAV_HEADER}> EPH G01 LNAV\n{}", join(G01_LINES));
    assert_eq!(parse_nav(&ok).expect("v4 NAV header accepted").len(), 1);

    // A 4.00 header that is not a navigation file (column 20 != 'N') is rejected.
    let bogus = "     4.00           OBSERVATION DATA   M                   RINEX VERSION / TYPE\n\
                 END OF HEADER\n";
    assert!(matches!(
        parse_nav(bogus),
        Err(NavParseError::UnsupportedHeader(_))
    ));
}

#[test]
fn from_nav_keeps_only_healthy_supported_messages() {
    let store = BroadcastStore::from_nav(&fixture_text()).expect("parse NAV");
    let recs = store.records();
    assert!(!recs.is_empty());
    // Every kept record is healthy and a supported single-frequency message:
    // GPS LNAV, Galileo I/NAV, or BeiDou D1/D2. Galileo F/NAV and unhealthy
    // satellites are dropped.
    assert!(
        recs.iter().all(|r| r.sv_health == 0.0),
        "unhealthy record kept"
    );
    assert!(
        recs.iter().all(|r| matches!(
            r.message,
            NavMessage::GpsLnav
                | NavMessage::GalileoInav
                | NavMessage::BeidouD1
                | NavMessage::BeidouD2
        )),
        "an unsupported message type was kept"
    );
    for sys in [GnssSystem::Gps, GnssSystem::Galileo, GnssSystem::BeiDou] {
        assert!(
            recs.iter().any(|r| r.satellite_id.system == sys),
            "no {sys:?} records kept"
        );
    }
    assert!(
        recs.iter().all(|r| r.message != NavMessage::GalileoFnav),
        "Galileo F/NAV must be excluded"
    );
    // The fixture's BeiDou set includes the geostationary C05 (a D2 message).
    assert!(
        recs.iter().any(
            |r| r.satellite_id == GnssSatelliteId::new(GnssSystem::BeiDou, 5)
                && r.message == NavMessage::BeidouD2
        ),
        "expected the geostationary C05 (D2) record"
    );
}

#[test]
fn a_wrong_week_epoch_has_no_ephemeris() {
    use crate::spp::EphemerisSource;
    let store = BroadcastStore::from_nav(&fixture_text()).expect("parse NAV");
    let sat = store.records()[0].satellite_id;

    // 2020-06-25 12:00 GPST as a J2000 second: a usable epoch for this product.
    let t_ok = 646_358_400.0_f64;
    assert!(
        store.position_clock_at_j2000_s(sat, t_ok).is_some(),
        "expected ephemeris at a valid epoch"
    );

    // The same wall-clock one week earlier: the nearest record is a week stale,
    // so the store must report no ephemeris rather than extrapolating a wrong
    // week's elements.
    let t_wrong_week = t_ok - 604_800.0;
    assert!(
        store.position_clock_at_j2000_s(sat, t_wrong_week).is_none(),
        "a wrong-week epoch must not silently produce an ephemeris"
    );
}

/// The J2000 second at which a record's reference epoch (`toe`) occurs, given the
/// satellite's timescale. BeiDou runs on BDT (= GPST - 14 s) with its week epoch
/// 1356 weeks after the GPS epoch; GPS/Galileo are GPST-aligned.
fn toe_as_j2000_s(rec: &BroadcastRecord) -> f64 {
    let toe_continuous = f64::from(rec.week) * 604_800.0 + rec.elements.toe_sow;
    let gps_epoch_to_j2000 = 630_763_200.0;
    if rec.satellite_id.system == GnssSystem::BeiDou {
        toe_continuous + 14.0 + 1356.0 * 604_800.0 - gps_epoch_to_j2000
    } else {
        toe_continuous - gps_epoch_to_j2000
    }
}

#[test]
fn broadcast_store_evaluates_beidou_including_geo() {
    use crate::spp::EphemerisSource;

    let store = BroadcastStore::from_nav(&fixture_text()).expect("parse NAV");
    // The geostationary C05 and a MEO (C19+) BeiDou satellite, evaluated at each
    // one's own reference epoch through the store's BDT timescale mapping.
    let geo = GnssSatelliteId::new(GnssSystem::BeiDou, 5);
    let meo = store
        .records()
        .iter()
        .map(|r| r.satellite_id)
        .find(|s| s.system == GnssSystem::BeiDou && s.prn >= 19)
        .expect("a BeiDou MEO satellite");

    for (sat, lo_km, hi_km) in [(geo, 41_000.0, 43_000.0), (meo, 27_000.0, 29_000.0)] {
        let rec = store
            .records()
            .iter()
            .find(|r| r.satellite_id == sat)
            .unwrap();
        let t = toe_as_j2000_s(rec);
        let (pos, _clk) = store
            .position_clock_at_j2000_s(sat, t)
            .unwrap_or_else(|| panic!("{sat:?} should evaluate at its toe"));
        let radius_km = (pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]).sqrt() / 1000.0;
        assert!(
            (lo_km..hi_km).contains(&radius_km),
            "{sat:?} radius {radius_km} km out of band"
        );
    }
    // The geostationary satellite sits near the equatorial plane.
    let c05 = store
        .records()
        .iter()
        .find(|r| r.satellite_id == geo)
        .unwrap();
    let (geo_pos, _) = store
        .position_clock_at_j2000_s(geo, toe_as_j2000_s(c05))
        .unwrap();
    let radius = (geo_pos[0].powi(2) + geo_pos[1].powi(2) + geo_pos[2].powi(2)).sqrt();
    assert!(
        geo_pos[2].abs() / radius < 0.2,
        "GEO should be near-equatorial"
    );
}

#[test]
fn broadcast_store_rejects_unsupported_systems() {
    use crate::broadcast::{ClockPolynomial, KeplerianElements};
    use crate::spp::EphemerisSource;

    // `BroadcastStore::new` accepts arbitrary records; a GLONASS satellite (a
    // non-Keplerian state-vector model) must report no ephemeris rather than be
    // evaluated with the wrong model.
    let sat = GnssSatelliteId::new(GnssSystem::Glonass, 1);
    let rec = BroadcastRecord {
        satellite_id: sat,
        message: NavMessage::GpsLnav,
        week: 2111,
        elements: KeplerianElements {
            sqrt_a: 5153.0,
            e: 0.001,
            m0: 0.0,
            delta_n: 0.0,
            omega0: 0.0,
            i0: 0.9,
            omega: 0.0,
            omega_dot: 0.0,
            idot: 0.0,
            cuc: 0.0,
            cus: 0.0,
            crc: 0.0,
            crs: 0.0,
            cic: 0.0,
            cis: 0.0,
            toe_sow: 0.0,
        },
        clock: ClockPolynomial {
            af0: 0.0,
            af1: 0.0,
            af2: 0.0,
            toc_sow: 0.0,
        },
        group_delay_s: 0.0,
        sv_health: 0.0,
        sv_accuracy_m: 0.0,
        fit_interval_s: None,
    };
    let store = BroadcastStore::new(vec![rec]);
    assert!(
        store
            .position_clock_at_j2000_s(sat, 646_358_400.0)
            .is_none(),
        "an unsupported system must report no ephemeris"
    );
}

#[test]
fn gps_fit_interval_bounds_record_validity() {
    use crate::broadcast::{ClockPolynomial, KeplerianElements};
    use crate::spp::EphemerisSource;

    // A minimal but evaluable record; only the system and fit interval matter to
    // selection (the orbit values just have to produce a finite state).
    let make = |system, fit_interval_s| BroadcastRecord {
        satellite_id: GnssSatelliteId::new(system, 1),
        message: NavMessage::GpsLnav,
        week: 2111,
        elements: KeplerianElements {
            sqrt_a: 5153.0,
            e: 0.001,
            m0: 0.0,
            delta_n: 0.0,
            omega0: 0.0,
            i0: 0.9,
            omega: 0.0,
            omega_dot: 0.0,
            idot: 0.0,
            cuc: 0.0,
            cus: 0.0,
            crc: 0.0,
            crs: 0.0,
            cic: 0.0,
            cis: 0.0,
            toe_sow: 0.0,
        },
        clock: ClockPolynomial {
            af0: 0.0,
            af1: 0.0,
            af2: 0.0,
            toc_sow: 0.0,
        },
        group_delay_s: 0.0,
        sv_health: 0.0,
        sv_accuracy_m: 0.0,
        fit_interval_s,
    };

    // `toe` at GPS week 2111, second-of-week 0, expressed as a J2000 second.
    let toe_j2000 = 2111.0 * 604_800.0 - 630_763_200.0;

    // GPS with a four-hour fit interval is valid within +/-2 h of toe only.
    let g = GnssSatelliteId::new(GnssSystem::Gps, 1);
    let gps = BroadcastStore::new(vec![make(GnssSystem::Gps, Some(4.0 * 3600.0))]);
    assert!(
        gps.position_clock_at_j2000_s(g, toe_j2000 + 3600.0)
            .is_some(),
        "1 h after toe is inside the 4 h fit interval"
    );
    assert!(
        gps.position_clock_at_j2000_s(g, toe_j2000 + 3.0 * 3600.0)
            .is_none(),
        "3 h after toe is outside the 4 h fit interval"
    );

    // A record with no broadcast fit interval (Galileo/BeiDou) falls back to the
    // coarse age bound, so the same 3 h offset is still accepted.
    let e = GnssSatelliteId::new(GnssSystem::Galileo, 1);
    let gal = BroadcastStore::new(vec![make(GnssSystem::Galileo, None)]);
    assert!(
        gal.position_clock_at_j2000_s(e, toe_j2000 + 3.0 * 3600.0)
            .is_some(),
        "without a fit interval the coarse 4 h bound applies"
    );
}

#[test]
fn select_prefers_a_valid_farther_record_over_an_expired_nearer_one() {
    use crate::broadcast::{ClockPolynomial, KeplerianElements};
    use crate::spp::EphemerisSource;

    let elements = |toe_sow| KeplerianElements {
        sqrt_a: 5153.0,
        e: 0.001,
        m0: 0.0,
        delta_n: 0.0,
        omega0: 0.0,
        i0: 0.9,
        omega: 0.0,
        omega_dot: 0.0,
        idot: 0.0,
        cuc: 0.0,
        cus: 0.0,
        crc: 0.0,
        crs: 0.0,
        cic: 0.0,
        cis: 0.0,
        toe_sow,
    };
    let rec = |toe_sow, fit_interval_s| BroadcastRecord {
        satellite_id: GnssSatelliteId::new(GnssSystem::Gps, 1),
        message: NavMessage::GpsLnav,
        week: 2111,
        elements: elements(toe_sow),
        clock: ClockPolynomial {
            af0: 0.0,
            af1: 0.0,
            af2: 0.0,
            toc_sow: toe_sow,
        },
        group_delay_s: 0.0,
        sv_health: 0.0,
        sv_accuracy_m: 0.0,
        fit_interval_s,
    };

    // Two records for one satellite: a nearer one with the nominal 4 h fit
    // (valid +/-2 h) and a farther one (toe 3 h earlier) with an extended 26 h
    // fit (valid +/-13 h).
    let near = rec(10_800.0, Some(4.0 * 3600.0));
    let far = rec(0.0, Some(26.0 * 3600.0));
    let store = BroadcastStore::new(vec![near, far]);

    let g = GnssSatelliteId::new(GnssSystem::Gps, 1);
    // 3 h after the nearer record's toe: outside its +/-2 h window, but inside
    // the farther record's +/-13 h window. Selecting nearest-then-checking would
    // wrongly reject this; filtering by validity first serves it from the farther
    // record.
    let near_toe_j2000 = 2111.0 * 604_800.0 + 10_800.0 - 630_763_200.0;
    let q = near_toe_j2000 + 3.0 * 3600.0;
    assert!(
        store.position_clock_at_j2000_s(g, q).is_some(),
        "a query past the nearest record's fit interval must fall back to a \
         farther record whose own window still covers it"
    );
}

#[test]
fn gps_fit_interval_field_distinguishes_blank_zero_value_and_malformed() {
    // Place a value in ORBIT-7 field 2 (columns 23..42): 23 leading blanks then
    // the field, so `field(line, 23, 42)` reads exactly the value.
    let with_field2 = |val: &str| format!("{:23}{:<19}", "", val);

    // Blank/absent -> the nominal four hours.
    assert_eq!(gps_fit_interval_s(&with_field2("")), Ok(4.0 * 3600.0));
    // Explicit zero -> the nominal four hours.
    assert_eq!(
        gps_fit_interval_s(&with_field2("0.000000000000e+00")),
        Ok(4.0 * 3600.0)
    );
    // A nonzero interval is taken verbatim (hours -> seconds).
    assert_eq!(
        gps_fit_interval_s(&with_field2("6.000000000000e+00")),
        Ok(6.0 * 3600.0)
    );
    // Present but non-numeric -> an error, not a silent nominal substitution.
    assert!(gps_fit_interval_s(&with_field2("garbage")).is_err());
}

#[test]
fn mixed_constellation_solve_recovers_the_receiver() {
    use crate::spp::{
        solve, test_support, Corrections, KlobucharCoeffs, Observation, SolveInputs, SurfaceMet,
        ELEVATION_MASK_RAD,
    };

    // The default store carries both GPS LNAV and Galileo I/NAV (healthy).
    let store = BroadcastStore::from_nav(&fixture_text()).expect("parse NAV");
    let t_rx = 646_358_400.0_f64;
    let sod = 12.0 * 3600.0;
    let doy = 177.5;
    let x_true = [3_512_900.0, 780_500.0, 5_248_700.0];
    let corr = Corrections::NONE;
    let kl = KlobucharCoeffs {
        alpha: [0.0; 4],
        beta: [0.0; 4],
    };
    let met = SurfaceMet {
        pressure_hpa: 1013.25,
        temperature_k: 288.15,
        relative_humidity: 0.5,
    };

    let mut sats: Vec<_> = store.records().iter().map(|r| r.satellite_id).collect();
    sats.sort_unstable();
    sats.dedup();

    let mut observations = Vec::new();
    let (mut have_gps, mut have_gal) = (false, false);
    for sat in sats {
        if let Some(m) = test_support::sat_model_for_test(
            &store,
            sat,
            x_true,
            0.0,
            22_000_000.0,
            t_rx,
            sod,
            doy,
            corr,
            &kl,
            &met,
        ) {
            if m.el_rad >= ELEVATION_MASK_RAD {
                observations.push(Observation {
                    satellite_id: sat,
                    pseudorange_m: m.p_hat_m,
                });
                have_gps |= sat.system == GnssSystem::Gps;
                have_gal |= sat.system == GnssSystem::Galileo;
            }
        }
    }
    assert!(
        have_gps && have_gal,
        "fixture must yield both GPS and Galileo observations"
    );

    let inputs = SolveInputs {
        observations,
        t_rx_j2000_s: t_rx,
        t_rx_second_of_day_s: sod,
        day_of_year: doy,
        initial_guess: [
            x_true[0] + 1000.0,
            x_true[1] - 1000.0,
            x_true[2] + 1000.0,
            0.0,
        ],
        corrections: corr,
        klobuchar: kl,
        beidou_klobuchar: None,
        met,
    };

    // The combined GPS+Galileo solve carries a per-system clock (a reference
    // clock plus the GPS/Galileo inter-system bias), so it recovers the receiver
    // from the mixed set. The geometry also yields a multi-system DOP.
    let sol = solve(&store, &inputs, true).expect("mixed-constellation solve");
    let p = sol.position;
    let err =
        ((p.x_m - x_true[0]).powi(2) + (p.y_m - x_true[1]).powi(2) + (p.z_m - x_true[2]).powi(2))
            .sqrt();
    assert!(
        err < 1.0e-3,
        "mixed solve recovered position off by {err} m"
    );

    let used_gps = sol.used_sats.iter().any(|s| s.system == GnssSystem::Gps);
    let used_gal = sol
        .used_sats
        .iter()
        .any(|s| s.system == GnssSystem::Galileo);
    assert!(
        used_gps && used_gal,
        "the solve must use both constellations"
    );
    let dop = sol
        .dop
        .expect("multi-system DOP present for the mixed solve");
    for (v, name) in [
        (dop.gdop, "GDOP"),
        (dop.pdop, "PDOP"),
        (dop.hdop, "HDOP"),
        (dop.vdop, "VDOP"),
        (dop.tdop, "TDOP"),
    ] {
        assert!(
            v.is_finite() && v > 0.0,
            "multi-system {name} not finite/positive: {v}"
        );
    }
}

#[test]
fn mixed_constellation_solve_recovers_a_nonzero_inter_system_bias() {
    use crate::spp::{
        solve, test_support, Corrections, KlobucharCoeffs, Observation, SolveInputs, SurfaceMet,
        C_M_S, ELEVATION_MASK_RAD,
    };

    let store = BroadcastStore::from_nav(&fixture_text()).expect("parse NAV");
    let t_rx = 646_358_400.0_f64;
    let sod = 12.0 * 3600.0;
    let doy = 177.5;
    let x_true = [3_512_900.0, 780_500.0, 5_248_700.0];
    // The Galileo receiver clock leads the GPS one by a real inter-system bias.
    let gal_bias_m = 50.0_f64;
    let corr = Corrections::NONE;
    let kl = KlobucharCoeffs {
        alpha: [0.0; 4],
        beta: [0.0; 4],
    };
    let met = SurfaceMet {
        pressure_hpa: 1013.25,
        temperature_k: 288.15,
        relative_humidity: 0.5,
    };

    let mut sats: Vec<_> = store.records().iter().map(|r| r.satellite_id).collect();
    sats.sort_unstable();
    sats.dedup();

    // Synthesize each pseudorange at the true position with the receiver clock
    // its own system sees: 0 for GPS (the reference), gal_bias_m for Galileo.
    let mut observations = Vec::new();
    let (mut have_gps, mut have_gal) = (false, false);
    for sat in sats {
        let b = if sat.system == GnssSystem::Galileo {
            gal_bias_m
        } else {
            0.0
        };
        if let Some(m) = test_support::sat_model_for_test(
            &store,
            sat,
            x_true,
            b,
            22_000_000.0,
            t_rx,
            sod,
            doy,
            corr,
            &kl,
            &met,
        ) {
            if m.el_rad >= ELEVATION_MASK_RAD {
                observations.push(Observation {
                    satellite_id: sat,
                    pseudorange_m: m.p_hat_m,
                });
                have_gps |= sat.system == GnssSystem::Gps;
                have_gal |= sat.system == GnssSystem::Galileo;
            }
        }
    }
    assert!(
        have_gps && have_gal,
        "need both GPS and Galileo observations"
    );

    let inputs = SolveInputs {
        observations,
        t_rx_j2000_s: t_rx,
        t_rx_second_of_day_s: sod,
        day_of_year: doy,
        initial_guess: [
            x_true[0] + 1000.0,
            x_true[1] - 1000.0,
            x_true[2] + 1000.0,
            0.0,
        ],
        corrections: corr,
        klobuchar: kl,
        beidou_klobuchar: None,
        met,
    };

    let sol = solve(&store, &inputs, false).expect("mixed solve with inter-system bias");

    // Position is still recovered despite the inter-system bias.
    let p = sol.position;
    let err =
        ((p.x_m - x_true[0]).powi(2) + (p.y_m - x_true[1]).powi(2) + (p.z_m - x_true[2]).powi(2))
            .sqrt();
    assert!(err < 1.0e-3, "recovered position off by {err} m");

    // The per-system clocks are recovered: GPS ~ 0, Galileo ~ the injected bias.
    let clk = |sys| {
        sol.system_clocks_s
            .iter()
            .find(|(s, _)| *s == sys)
            .map(|(_, c)| *c * C_M_S)
            .unwrap_or_else(|| panic!("no {sys:?} clock"))
    };
    assert!(
        clk(GnssSystem::Gps).abs() < 1.0e-3,
        "GPS clock {} m",
        clk(GnssSystem::Gps)
    );
    assert!(
        (clk(GnssSystem::Galileo) - gal_bias_m).abs() < 1.0e-3,
        "Galileo clock {} m, expected ~{gal_bias_m}",
        clk(GnssSystem::Galileo)
    );
}

#[test]
fn mixed_solve_recovers_with_gps_galileo_and_beidou() {
    use crate::spp::{
        solve, test_support, Corrections, KlobucharCoeffs, Observation, SolveInputs, SurfaceMet,
        C_M_S, ELEVATION_MASK_RAD,
    };

    let store = BroadcastStore::from_nav(&fixture_text()).expect("parse NAV");
    let t_rx = 646_358_400.0_f64;
    let sod = 12.0 * 3600.0;
    let doy = 177.5;
    let x_true = [3_512_900.0, 780_500.0, 5_248_700.0];
    let corr = Corrections::NONE;
    let kl = KlobucharCoeffs {
        alpha: [0.0; 4],
        beta: [0.0; 4],
    };
    let met = SurfaceMet {
        pressure_hpa: 1013.25,
        temperature_k: 288.15,
        relative_humidity: 0.5,
    };

    // A distinct receiver-clock bias per system (GPS is the reference).
    let bias_m = |sys| match sys {
        GnssSystem::Galileo => 50.0,
        GnssSystem::BeiDou => 120.0,
        _ => 0.0,
    };

    let mut sats: Vec<_> = store.records().iter().map(|r| r.satellite_id).collect();
    sats.sort_unstable();
    sats.dedup();

    let mut observations = Vec::new();
    let (mut g, mut e, mut c) = (false, false, false);
    for sat in sats {
        if let Some(m) = test_support::sat_model_for_test(
            &store,
            sat,
            x_true,
            bias_m(sat.system),
            22_000_000.0,
            t_rx,
            sod,
            doy,
            corr,
            &kl,
            &met,
        ) {
            if m.el_rad >= ELEVATION_MASK_RAD {
                observations.push(Observation {
                    satellite_id: sat,
                    pseudorange_m: m.p_hat_m,
                });
                g |= sat.system == GnssSystem::Gps;
                e |= sat.system == GnssSystem::Galileo;
                c |= sat.system == GnssSystem::BeiDou;
            }
        }
    }
    assert!(
        g && e && c,
        "need GPS, Galileo, and BeiDou observations (got {g} {e} {c})"
    );

    let inputs = SolveInputs {
        observations,
        t_rx_j2000_s: t_rx,
        t_rx_second_of_day_s: sod,
        day_of_year: doy,
        initial_guess: [
            x_true[0] + 1000.0,
            x_true[1] - 1000.0,
            x_true[2] + 1000.0,
            0.0,
        ],
        corrections: corr,
        klobuchar: kl,
        beidou_klobuchar: None,
        met,
    };

    let sol = solve(&store, &inputs, false).expect("three-constellation solve");
    let p = sol.position;
    let err =
        ((p.x_m - x_true[0]).powi(2) + (p.y_m - x_true[1]).powi(2) + (p.z_m - x_true[2]).powi(2))
            .sqrt();
    assert!(err < 1.0e-3, "recovered position off by {err} m");

    let clk = |sys| {
        sol.system_clocks_s
            .iter()
            .find(|(s, _)| *s == sys)
            .map(|(_, v)| *v * C_M_S)
            .unwrap_or_else(|| panic!("no {sys:?} clock"))
    };
    assert!(
        clk(GnssSystem::Gps).abs() < 1.0e-3,
        "GPS clock {}",
        clk(GnssSystem::Gps)
    );
    assert!(
        (clk(GnssSystem::Galileo) - 50.0).abs() < 1.0e-3,
        "GAL clock {}",
        clk(GnssSystem::Galileo)
    );
    assert!(
        (clk(GnssSystem::BeiDou) - 120.0).abs() < 1.0e-3,
        "BDS clock {}",
        clk(GnssSystem::BeiDou)
    );
}

#[test]
fn ionosphere_correction_is_applied_to_beidou_b1i() {
    use crate::spp::{
        solve, test_support, Corrections, KlobucharCoeffs, Observation, SolveInputs, SurfaceMet,
        ELEVATION_MASK_RAD,
    };

    let store = BroadcastStore::from_nav(&fixture_text()).expect("parse NAV");
    let t_rx = 646_358_400.0_f64;
    let sod = 12.0 * 3600.0;
    let doy = 177.5;
    let x_true = [3_512_900.0, 780_500.0, 5_248_700.0];
    // Ionosphere on. The broadcast Klobuchar L1 delay is scaled to each carrier
    // by (f_L1/f)^2 — exactly 1 for GPS L1 / Galileo E1, and scaled for BeiDou
    // B1I — so a BeiDou-bearing iono-corrected solve is now supported (not
    // rejected) and recovers the truth from observations synthesized with the
    // same frequency-aware model.
    let corr = Corrections::IONO;
    let kl = KlobucharCoeffs {
        alpha: [1.0e-8, 0.0, 0.0, 0.0],
        beta: [9.0e4, 0.0, 0.0, 0.0],
    };
    let met = SurfaceMet {
        pressure_hpa: 1013.25,
        temperature_k: 288.15,
        relative_humidity: 0.5,
    };

    let mut sats: Vec<_> = store.records().iter().map(|r| r.satellite_id).collect();
    sats.sort_unstable();
    sats.dedup();

    let mut observations = Vec::new();
    let mut saw_beidou = false;
    for sat in sats {
        if let Some(m) = test_support::sat_model_for_test(
            &store,
            sat,
            x_true,
            0.0,
            22_000_000.0,
            t_rx,
            sod,
            doy,
            corr,
            &kl,
            &met,
        ) {
            if m.el_rad >= ELEVATION_MASK_RAD {
                saw_beidou |= sat.system == GnssSystem::BeiDou;
                observations.push(Observation {
                    satellite_id: sat,
                    pseudorange_m: m.p_hat_m,
                });
            }
        }
    }
    assert!(
        saw_beidou,
        "the iono-corrected set must include a BeiDou satellite"
    );

    let inputs = SolveInputs {
        observations,
        t_rx_j2000_s: t_rx,
        t_rx_second_of_day_s: sod,
        day_of_year: doy,
        initial_guess: [
            x_true[0] + 1000.0,
            x_true[1] - 1000.0,
            x_true[2] + 1000.0,
            0.0,
        ],
        corrections: corr,
        klobuchar: kl,
        beidou_klobuchar: None,
        met,
    };

    let sol = solve(&store, &inputs, false).expect("BeiDou-bearing iono-corrected solve");
    let p = sol.position;
    let err =
        ((p.x_m - x_true[0]).powi(2) + (p.y_m - x_true[1]).powi(2) + (p.z_m - x_true[2]).powi(2))
            .sqrt();
    assert!(err < 1.0e-3, "recovered position off by {err} m");
}