use sidereon_core::astro::time::model::JulianDateSplit;
use sidereon_core::astro::time::split_julian_date;
use sidereon_core::constants::{
    C_M_S, F_L1_HZ, SECONDS_PER_DAY, SECONDS_PER_HOUR, SECONDS_PER_MINUTE,
};
use sidereon_core::ephemeris::BroadcastEphemeris;
use sidereon_core::observables::{
    j2000_seconds_from_split, predict, ObservableEphemerisSource, ObservableState,
    ObservablesError, PredictOptions,
};
use sidereon_core::positioning::{
    solve, solve_doppler_velocity, solve_with_doppler_velocity, Corrections, DopplerObservation,
    DopplerVelocityInputs, EphemerisSource, KlobucharCoeffs, Observation, SolveInputs, SurfaceMet,
};
use sidereon_core::rinex::nav::parse_nav;
use sidereon_core::rinex::observations::{
    observation_values, ObsEpoch, ObsEpochTime, ObservationFilter, RinexObs,
};
use sidereon_core::rinex::qc::{repair_obs_text, RepairOptions};
use sidereon_core::velocity::range_rate_to_doppler;
use sidereon_core::{Error, GnssSatelliteId, GnssSystem};

#[derive(Debug, Clone, Copy)]
struct SatState {
    id: GnssSatelliteId,
    position0_m: [f64; 3],
    velocity_m_s: [f64; 3],
}

#[derive(Debug, Clone)]
struct LinearSource {
    t0_s: f64,
    sats: Vec<SatState>,
}

impl LinearSource {
    fn state(&self, sat: GnssSatelliteId, t_s: f64) -> Option<[f64; 3]> {
        let state = self.sats.iter().find(|state| state.id == sat)?;
        let dt = t_s - self.t0_s;
        Some([
            state.position0_m[0] + state.velocity_m_s[0] * dt,
            state.position0_m[1] + state.velocity_m_s[1] * dt,
            state.position0_m[2] + state.velocity_m_s[2] * dt,
        ])
    }
}

impl EphemerisSource for LinearSource {
    fn position_clock_at_j2000_s(
        &self,
        sat: GnssSatelliteId,
        t_j2000_s: f64,
    ) -> Option<([f64; 3], f64)> {
        Some((self.state(sat, t_j2000_s)?, 0.0))
    }
}

impl ObservableEphemerisSource for LinearSource {
    fn observable_state_at_j2000_s(
        &self,
        sat: GnssSatelliteId,
        t_j2000_s: f64,
    ) -> Result<ObservableState, ObservablesError> {
        let Some(position_ecef_m) = self.state(sat, t_j2000_s) else {
            return Err(ObservablesError::NoEphemeris);
        };
        Ok(ObservableState {
            position_ecef_m,
            clock_s: Some(0.0),
        })
    }
}

fn gps(prn: u8) -> GnssSatelliteId {
    GnssSatelliteId::new(GnssSystem::Gps, prn).expect("valid GPS satellite")
}

fn norm(v: [f64; 3]) -> f64 {
    (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt()
}

fn unit(v: [f64; 3]) -> [f64; 3] {
    let n = norm(v);
    [v[0] / n, v[1] / n, v[2] / n]
}

fn add_scaled(a: [f64; 3], b: [f64; 3], scale: f64) -> [f64; 3] {
    [
        a[0] + b[0] * scale,
        a[1] + b[1] * scale,
        a[2] + b[2] * scale,
    ]
}

fn dot(a: [f64; 3], b: [f64; 3]) -> f64 {
    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
}

fn synthetic_case() -> (LinearSource, [f64; 3], [f64; 3], f64) {
    let t0_s = 800_000.0;
    let receiver0_m = [6_378_137.0, 0.0, 0.0];
    let receiver_velocity_m_s = [8.0, -3.0, 2.0];
    let directions = [
        [0.75, 0.50, 0.43],
        [0.65, -0.60, 0.46],
        [0.85, 0.25, -0.46],
        [0.55, -0.10, 0.83],
        [0.70, 0.70, -0.14],
        [0.90, -0.35, 0.26],
    ];
    let sat_velocities = [
        [80.0, -320.0, 120.0],
        [-40.0, 260.0, -160.0],
        [140.0, 40.0, 210.0],
        [-180.0, -130.0, 70.0],
        [30.0, 180.0, -220.0],
        [110.0, -90.0, -140.0],
    ];
    let sats = directions
        .into_iter()
        .zip(sat_velocities)
        .enumerate()
        .map(|(index, (direction, velocity_m_s))| SatState {
            id: gps((index + 1) as u8),
            position0_m: add_scaled(receiver0_m, unit(direction), 23_000_000.0),
            velocity_m_s,
        })
        .collect();

    (
        LinearSource { t0_s, sats },
        receiver0_m,
        receiver_velocity_m_s,
        t0_s,
    )
}

fn receiver_position_at(
    receiver0_m: [f64; 3],
    velocity_m_s: [f64; 3],
    t0_s: f64,
    t_s: f64,
) -> [f64; 3] {
    add_scaled(receiver0_m, velocity_m_s, t_s - t0_s)
}

fn pseudorange_inputs(
    source: &LinearSource,
    receiver0_m: [f64; 3],
    receiver_velocity_m_s: [f64; 3],
    t0_s: f64,
    t_s: f64,
) -> SolveInputs {
    let receiver_m = receiver_position_at(receiver0_m, receiver_velocity_m_s, t0_s, t_s);
    let observations = source
        .sats
        .iter()
        .map(|sat| {
            let predicted = predict(source, sat.id, receiver_m, t_s, PredictOptions::default())
                .expect("synthetic range prediction");
            Observation {
                satellite_id: sat.id,
                pseudorange_m: predicted.geometric_range_m,
            }
        })
        .collect();
    SolveInputs {
        observations,
        t_rx_j2000_s: t_s,
        t_rx_second_of_day_s: t_s.rem_euclid(86_400.0),
        day_of_year: 100.0,
        initial_guess: [
            receiver_m[0] + 40.0,
            receiver_m[1] - 30.0,
            receiver_m[2] + 20.0,
            0.0,
        ],
        corrections: Corrections::NONE,
        klobuchar: KlobucharCoeffs {
            alpha: [0.0; 4],
            beta: [0.0; 4],
        },
        beidou_klobuchar: None,
        galileo_nequick: None,
        sbas_iono: None,
        glonass_channels: Default::default(),
        met: SurfaceMet::default(),
        robust: None,
    }
}

fn doppler_observations(
    source: &LinearSource,
    receiver0_m: [f64; 3],
    receiver_velocity_m_s: [f64; 3],
    t0_s: f64,
    clock_drift_s_s: f64,
) -> Vec<DopplerObservation> {
    source
        .sats
        .iter()
        .map(|sat| {
            let predicted = predict(source, sat.id, receiver0_m, t0_s, PredictOptions::default())
                .expect("synthetic Doppler prediction");
            let range_rate_m_s = predicted.range_rate_m_s
                - dot(predicted.los_unit, receiver_velocity_m_s)
                + C_M_S * clock_drift_s_s;
            DopplerObservation {
                satellite_id: sat.id,
                doppler_hz: range_rate_to_doppler(range_rate_m_s, F_L1_HZ)
                    .expect("range-rate to Doppler"),
                carrier_hz: F_L1_HZ,
                sat_clock_drift_s_s: 0.0,
            }
        })
        .collect()
}

fn synthetic_rinex_obs_epoch(
    source: &LinearSource,
    receiver0_m: [f64; 3],
    receiver_velocity_m_s: [f64; 3],
    t0_s: f64,
    clock_drift_s_s: f64,
) -> String {
    let mut text = String::new();
    text.push_str("3.04 OBSERVATION DATA M RINEX VERSION / TYPE\n");
    text.push_str(&obs_header_line("TEST", "MARKER NAME"));
    text.push_str(&obs_header_line(
        "PGM                 RUN                 DATE",
        "PGM / RUN BY / DATE",
    ));
    text.push_str(&obs_header_line(
        "OBS                 AGENCY",
        "OBSERVER / AGENCY",
    ));
    text.push_str(&obs_header_line(
        "REC                 TYPE                VERS",
        "REC # / TYPE / VERS",
    ));
    text.push_str(&obs_header_line("ANT                 TYPE", "ANT # / TYPE"));
    text.push_str(&obs_header_line(
        "  6378137.0000        0.0000        0.0000",
        "APPROX POSITION XYZ",
    ));
    text.push_str(&obs_header_line(
        "        0.0000        0.0000        0.0000",
        "ANTENNA: DELTA H/E/N",
    ));
    text.push_str("G 2 C1C D1C SYS / # / OBS TYPES\n");
    text.push_str(&obs_header_line(
        "  2020     1     1     0     0    0.0000000     GPS",
        "TIME OF FIRST OBS",
    ));
    text.push_str(&obs_header_line("", "END OF HEADER"));
    text.push_str(&format!(
        "> 2020 01 01 00 00  0.0000000  0{:3}\n",
        source.sats.len()
    ));

    for sat in &source.sats {
        let predicted = predict(source, sat.id, receiver0_m, t0_s, PredictOptions::default())
            .expect("synthetic RINEX prediction");
        let range_rate_m_s = predicted.range_rate_m_s
            - dot(predicted.los_unit, receiver_velocity_m_s)
            + C_M_S * clock_drift_s_s;
        let doppler_hz =
            range_rate_to_doppler(range_rate_m_s, F_L1_HZ).expect("range-rate to Doppler");
        text.push_str(&format!(
            "{:<3}{:14.3}  {:14.3}  \n",
            sat.id, predicted.geometric_range_m, doppler_hz
        ));
    }

    text
}

fn parsed_l1_rows(obs: &RinexObs) -> (Vec<Observation>, Vec<DopplerObservation>) {
    let codes = obs.obs_codes(GnssSystem::Gps).expect("GPS code table");
    let c1c = codes.iter().position(|code| code == "C1C").expect("C1C");
    let d1c = codes.iter().position(|code| code == "D1C").expect("D1C");
    let epoch = &obs.epochs()[0];
    let mut ranges = Vec::new();
    let mut doppler = Vec::new();
    for (sat, values) in &epoch.sats {
        ranges.push(Observation {
            satellite_id: *sat,
            pseudorange_m: values[c1c].value.expect("C1C value"),
        });
        doppler.push(DopplerObservation {
            satellite_id: *sat,
            doppler_hz: values[d1c].value.expect("D1C value"),
            carrier_hz: F_L1_HZ,
            sat_clock_drift_s_s: 0.0,
        });
    }
    (ranges, doppler)
}

fn fixture_text(parts: &[&str]) -> String {
    let mut path = std::path::PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("tests/fixtures");
    for part in parts {
        path.push(part);
    }
    std::fs::read_to_string(&path).unwrap_or_else(|err| panic!("read fixture {path:?}: {err}"))
}

fn civil_to_julian_split(epoch: ObsEpochTime) -> JulianDateSplit {
    let (jd_whole, fraction) = split_julian_date(
        epoch.year,
        i32::from(epoch.month),
        i32::from(epoch.day),
        i32::from(epoch.hour),
        i32::from(epoch.minute),
        epoch.second,
    );
    JulianDateSplit::new(jd_whole, fraction).expect("valid split Julian date")
}

fn j2000_seconds(epoch: ObsEpochTime) -> f64 {
    let split = civil_to_julian_split(epoch);
    j2000_seconds_from_split(split.jd_whole, split.fraction).expect("valid split Julian date")
}

fn second_of_day(epoch: ObsEpochTime) -> f64 {
    f64::from(epoch.hour) * SECONDS_PER_HOUR
        + f64::from(epoch.minute) * SECONDS_PER_MINUTE
        + epoch.second
}

fn esbc_epoch_rows(
    obs: &RinexObs,
    epoch: &ObsEpoch,
) -> (Vec<Observation>, Vec<DopplerObservation>) {
    let filter = ObservationFilter::from_entries([(
        GnssSystem::Gps,
        vec!["C1C".to_string(), "D1C".to_string()],
    )]);
    let values = observation_values(obs, epoch, &filter).expect("observation values");
    let mut ranges = Vec::new();
    let mut doppler = Vec::new();
    for (sat, rows) in values {
        if sat.system != GnssSystem::Gps {
            continue;
        }
        if let Some(pseudorange_m) = rows
            .iter()
            .find(|row| row.code == "C1C")
            .and_then(|row| row.value)
        {
            ranges.push(Observation {
                satellite_id: sat,
                pseudorange_m,
            });
        }
        if let Some(doppler_hz) = rows
            .iter()
            .find(|row| row.code == "D1C")
            .and_then(|row| row.value)
        {
            doppler.push(DopplerObservation {
                satellite_id: sat,
                doppler_hz,
                carrier_hz: F_L1_HZ,
                sat_clock_drift_s_s: 0.0,
            });
        }
    }
    (ranges, doppler)
}

fn esbc_solve_inputs(obs: &RinexObs, epoch: &ObsEpoch) -> SolveInputs {
    let (observations, _) = esbc_epoch_rows(obs, epoch);
    assert!(
        observations.len() >= 5,
        "need a redundant GPS set, got {}",
        observations.len()
    );
    let approx = obs.header().approx_position_m.expect("APPROX POSITION XYZ");
    let sod = second_of_day(epoch.epoch);
    SolveInputs {
        observations,
        t_rx_j2000_s: j2000_seconds(epoch.epoch),
        t_rx_second_of_day_s: sod,
        day_of_year: 177.0 + sod / SECONDS_PER_DAY,
        initial_guess: [approx[0], approx[1], approx[2], 0.0],
        corrections: Corrections {
            ionosphere: false,
            troposphere: true,
        },
        klobuchar: KlobucharCoeffs {
            alpha: [0.0; 4],
            beta: [0.0; 4],
        },
        beidou_klobuchar: None,
        galileo_nequick: None,
        sbas_iono: None,
        glonass_channels: Default::default(),
        met: SurfaceMet::default(),
        robust: None,
    }
}

#[test]
fn spp_doppler_velocity_matches_finite_difference_position_oracle() {
    let (source, receiver0_m, receiver_velocity_m_s, t0_s) = synthetic_case();
    let half_dt_s = 20.0;
    let before = pseudorange_inputs(
        &source,
        receiver0_m,
        receiver_velocity_m_s,
        t0_s,
        t0_s - half_dt_s,
    );
    let after = pseudorange_inputs(
        &source,
        receiver0_m,
        receiver_velocity_m_s,
        t0_s,
        t0_s + half_dt_s,
    );
    let before_solution = solve(&source, &before, false).expect("solve before epoch");
    let after_solution = solve(&source, &after, false).expect("solve after epoch");
    let before_pos = before_solution.position.as_array();
    let after_pos = after_solution.position.as_array();
    let finite_difference_m_s = [
        (after_pos[0] - before_pos[0]) / (2.0 * half_dt_s),
        (after_pos[1] - before_pos[1]) / (2.0 * half_dt_s),
        (after_pos[2] - before_pos[2]) / (2.0 * half_dt_s),
    ];

    let clock_drift_s_s = 2.0e-9;
    let doppler = doppler_observations(
        &source,
        receiver0_m,
        receiver_velocity_m_s,
        t0_s,
        clock_drift_s_s,
    );
    let velocity = solve_doppler_velocity(
        &source,
        &DopplerVelocityInputs {
            observations: doppler.clone(),
            receiver_ecef_m: receiver0_m,
            t_rx_j2000_s: t0_s,
            light_time: true,
            sagnac: true,
        },
    )
    .expect("solve Doppler velocity");

    let error_m_s = norm([
        velocity.velocity_m_s[0] - finite_difference_m_s[0],
        velocity.velocity_m_s[1] - finite_difference_m_s[1],
        velocity.velocity_m_s[2] - finite_difference_m_s[2],
    ]);
    assert!(error_m_s < 0.10, "velocity error {error_m_s}");
    assert!((velocity.clock_drift_s_s - clock_drift_s_s).abs() < 1.0e-12);
    assert!(velocity
        .state_covariance
        .iter()
        .flatten()
        .all(|value| value.is_finite()));
    for idx in 0..4 {
        assert!(velocity.state_covariance[idx][idx] > 0.0);
    }
}

#[test]
fn spp_wrapper_uses_doppler_parsed_from_rinex_obs_epoch() {
    let (source, receiver0_m, receiver_velocity_m_s, t0_s) = synthetic_case();
    let clock_drift_s_s = -1.5e-9;
    let rinex_text = synthetic_rinex_obs_epoch(
        &source,
        receiver0_m,
        receiver_velocity_m_s,
        t0_s,
        clock_drift_s_s,
    );
    let obs = RinexObs::parse(&rinex_text).expect("strict parse synthetic RINEX OBS");
    assert_eq!(
        obs.obs_codes(GnssSystem::Gps).expect("GPS code table"),
        ["C1C".to_string(), "D1C".to_string()]
    );
    let (ranges, doppler) = parsed_l1_rows(&obs);
    let mut inputs = pseudorange_inputs(&source, receiver0_m, receiver_velocity_m_s, t0_s, t0_s);
    inputs.observations = ranges;

    let fused = solve_with_doppler_velocity(&source, &inputs, &doppler, false)
        .expect("solve RINEX parsed Doppler epoch");
    let velocity = fused.velocity.expect("Doppler velocity solution");
    assert_eq!(fused.velocity_error, None);
    assert_eq!(
        fused.receiver.rx_clock_drift_s_s,
        Some(velocity.clock_drift_s_s)
    );

    let error_m_s = norm([
        velocity.velocity_m_s[0] - receiver_velocity_m_s[0],
        velocity.velocity_m_s[1] - receiver_velocity_m_s[1],
        velocity.velocity_m_s[2] - receiver_velocity_m_s[2],
    ]);
    assert!(error_m_s < 0.02, "RINEX parsed velocity error {error_m_s}");
    assert!((velocity.clock_drift_s_s - clock_drift_s_s).abs() < 2.0e-12);
}

#[test]
fn spp_wrapper_solves_real_rinex_obs_nav_doppler_against_position_difference() {
    let nav_text = fixture_text(&["nav", "ESBC00DNK_R_20201770000_01D_MN.rnx"]);
    let store = BroadcastEphemeris::from_nav(&nav_text).expect("parse ESBC broadcast NAV");
    let obs_text = fixture_text(&["obs", "ESBC00DNK_R_20201770000_01D_30S_MO_trim.rnx"]);
    let obs = RinexObs::parse(&obs_text).expect("parse ESBC OBS");
    assert!(
        obs.epochs().len() >= 2,
        "ESBC trim must carry two epochs for the finite difference"
    );
    let first = &obs.epochs()[0];
    let second = &obs.epochs()[1];
    let first_inputs = esbc_solve_inputs(&obs, first);
    let second_inputs = esbc_solve_inputs(&obs, second);
    let (_, doppler) = esbc_epoch_rows(&obs, first);
    assert!(
        doppler.len() >= 4,
        "need at least four real D1C rows, got {}",
        doppler.len()
    );

    let first_position = solve(&store, &first_inputs, false).expect("solve first SPP epoch");
    let second_position = solve(&store, &second_inputs, false).expect("solve second SPP epoch");
    let dt_s = second_inputs.t_rx_j2000_s - first_inputs.t_rx_j2000_s;
    assert!(dt_s > 0.0, "epochs must be increasing");
    let p0 = first_position.position.as_array();
    let p1 = second_position.position.as_array();
    let finite_difference_m_s = [
        (p1[0] - p0[0]) / dt_s,
        (p1[1] - p0[1]) / dt_s,
        (p1[2] - p0[2]) / dt_s,
    ];

    let fused = solve_with_doppler_velocity(&store, &first_inputs, &doppler, false)
        .expect("solve public SPP wrapper from real OBS plus NAV");
    let velocity = fused.velocity.expect("real Doppler velocity solution");
    assert_eq!(fused.velocity_error, None);
    assert_eq!(
        fused.receiver.rx_clock_drift_s_s,
        Some(velocity.clock_drift_s_s)
    );
    let error_m_s = norm([
        velocity.velocity_m_s[0] - finite_difference_m_s[0],
        velocity.velocity_m_s[1] - finite_difference_m_s[1],
        velocity.velocity_m_s[2] - finite_difference_m_s[2],
    ]);
    assert!(
        error_m_s < 0.75,
        "real OBS/NAV wrapper velocity differs from SPP position difference by {error_m_s} m/s"
    );
}

#[test]
fn receiver_solution_carries_optional_clock_drift_from_doppler_rows() {
    let (source, receiver0_m, receiver_velocity_m_s, t0_s) = synthetic_case();
    let inputs = pseudorange_inputs(&source, receiver0_m, receiver_velocity_m_s, t0_s, t0_s);
    let pseudorange_only = solve(&source, &inputs, false).expect("solve SPP");
    assert_eq!(pseudorange_only.rx_clock_drift_s_s, None);

    let drift_s_s = 2.0e-9;
    let doppler =
        doppler_observations(&source, receiver0_m, receiver_velocity_m_s, t0_s, drift_s_s);
    let fused =
        solve_with_doppler_velocity(&source, &inputs, &doppler, false).expect("solve fused epoch");
    let velocity = fused.velocity.expect("Doppler velocity solution");
    assert_eq!(fused.velocity_error, None);
    assert_eq!(
        fused.receiver.rx_clock_drift_s_s,
        Some(velocity.clock_drift_s_s)
    );
}

#[test]
fn receiver_solution_position_covariance_exposes_ecef_and_consistent_enu() {
    let (source, receiver0_m, receiver_velocity_m_s, t0_s) = synthetic_case();
    let inputs = pseudorange_inputs(&source, receiver0_m, receiver_velocity_m_s, t0_s, t0_s);
    let solution = solve(&source, &inputs, true).expect("solve SPP");
    let geodetic = solution.geodetic.expect("geodetic solution");
    let rotated = sidereon_core::dop::rotate_covariance_ecef_to_enu_m2(
        solution.position_covariance.ecef_m2,
        geodetic,
    )
    .expect("rotate covariance");
    for (row, (got_row, expected_row)) in solution
        .position_covariance
        .enu_m2
        .iter()
        .zip(rotated.iter())
        .enumerate()
    {
        for (col, (&got, &expected)) in got_row.iter().zip(expected_row.iter()).enumerate() {
            assert!(
                (got - expected).abs() <= 1.0e-9,
                "covariance[{row}][{col}] got {got} expected {expected}"
            );
        }
    }
}

fn obs_header_line(body: &str, label: &str) -> String {
    format!("{body:<60}{label}\n")
}

fn minimal_obs_text(version_line: &str, obs_types_line: &str, sat: &str) -> String {
    let mut text = String::new();
    text.push_str(version_line);
    text.push('\n');
    text.push_str(&obs_header_line("TEST", "MARKER NAME"));
    text.push_str(&obs_header_line(
        "PGM                 RUN                 DATE",
        "PGM / RUN BY / DATE",
    ));
    text.push_str(&obs_header_line(
        "OBS                 AGENCY",
        "OBSERVER / AGENCY",
    ));
    text.push_str(&obs_header_line(
        "REC                 TYPE                VERS",
        "REC # / TYPE / VERS",
    ));
    text.push_str(&obs_header_line("ANT                 TYPE", "ANT # / TYPE"));
    text.push_str(&obs_header_line(
        "  6378137.0000        0.0000        0.0000",
        "APPROX POSITION XYZ",
    ));
    text.push_str(&obs_header_line(
        "        0.0000        0.0000        0.0000",
        "ANTENNA: DELTA H/E/N",
    ));
    text.push_str(obs_types_line);
    text.push('\n');
    text.push_str(&obs_header_line(
        "  2020     1     1     0     0    0.0000000     GPS",
        "TIME OF FIRST OBS",
    ));
    text.push_str(&obs_header_line("", "END OF HEADER"));
    text.push_str("> 2020 01 01 00 00  0.0000000  0  1\n");
    text.push_str(&format!("{sat:<3}{:14.3}\n", 20_200_000.0));
    text
}

fn extended_count_obs_text() -> String {
    let codes = [
        "C1C", "L1C", "D1C", "S1C", "C2W", "L2W", "D2W", "S2W", "C5Q", "L5Q", "D5Q", "S5Q", "C1W",
    ];
    let mut text = String::new();
    text.push_str(
        "     3.04           OBSERVATION DATA    M                   RINEX VERSION / TYPE\n",
    );
    text.push_str(&obs_header_line("TEST", "MARKER NAME"));
    text.push_str(&obs_header_line(
        "PGM                 RUN                 DATE",
        "PGM / RUN BY / DATE",
    ));
    text.push_str(&obs_header_line(
        "OBS                 AGENCY",
        "OBSERVER / AGENCY",
    ));
    text.push_str(&obs_header_line(
        "REC                 TYPE                VERS",
        "REC # / TYPE / VERS",
    ));
    text.push_str(&obs_header_line("ANT                 TYPE", "ANT # / TYPE"));
    text.push_str(&obs_header_line(
        "  6378137.0000        0.0000        0.0000",
        "APPROX POSITION XYZ",
    ));
    text.push_str(&obs_header_line(
        "        0.0000        0.0000        0.0000",
        "ANTENNA: DELTA H/E/N",
    ));
    let mut obs_types = format!("G  {:>3}", codes.len());
    for code in &codes {
        obs_types.push_str(&format!(" {code:>3}"));
    }
    text.push_str(&obs_header_line(&obs_types, "SYS / # / OBS TYPES"));

    let mut scale_first = format!("G {:>4}  {:>2}", 10, codes.len());
    for code in &codes[..12] {
        scale_first.push_str(&format!(" {code:>3}"));
    }
    let mut scale_second = " ".repeat(10);
    for code in &codes[12..] {
        scale_second.push_str(&format!(" {code:>3}"));
    }
    text.push_str(&obs_header_line(&scale_first, "SYS / SCALE FACTOR"));
    text.push_str(&obs_header_line(&scale_second, "SYS / SCALE FACTOR"));

    let mut prn_first = "G01".to_string();
    for value in 1..=9 {
        prn_first.push_str(&format!("{value:6}"));
    }
    let mut prn_second = " ".repeat(3);
    for value in 10..=13 {
        prn_second.push_str(&format!("{value:6}"));
    }
    text.push_str(&obs_header_line(&prn_first, "PRN / # OF OBS"));
    text.push_str(&obs_header_line(&prn_second, "PRN / # OF OBS"));
    text.push_str(&obs_header_line(
        "  2020     1     1     0     0    0.0000000     GPS",
        "TIME OF FIRST OBS",
    ));
    text.push_str(&obs_header_line("", "END OF HEADER"));
    text.push_str("> 2020 01 01 00 00  0.0000000  0  1\n");
    text.push_str(&format!("G01{:14.3}\n", 20_200_000.0));
    text
}

fn assert_strict_header_columns(text: &str) {
    let labels = [
        "RINEX VERSION / TYPE",
        "PGM / RUN BY / DATE",
        "MARKER NAME",
        "SYS / # / OBS TYPES",
        "END OF HEADER",
    ];
    for line in text.lines() {
        if line.contains("END OF HEADER") {
            assert_eq!(line.get(60..).unwrap_or("").trim(), "END OF HEADER");
            break;
        }
        if labels.iter().any(|label| line.contains(label)) {
            assert!(line.len() >= 60, "{line:?}");
            let label = line.get(60..).unwrap_or("").trim();
            assert!(labels.contains(&label), "{line:?}");
        }
    }
}

#[test]
fn repair_obs_text_outputs_strict_parseable_headers_for_trial_whitespace_cases() {
    let cases = [
        minimal_obs_text(
            "     3.04           OBSERVATION DATA    M                   RINEX VERSION / TYPE",
            "G    1 C1C                                                  SYS / # / OBS TYPES",
            "G01",
        ),
        minimal_obs_text(
            "3.04 OBSERVATION DATA M RINEX VERSION / TYPE",
            "G 1 C1C SYS / # / OBS TYPES",
            "G 1",
        ),
        minimal_obs_text(
            "     3.04           OBSERVATION DATA    M                      RINEX VERSION / TYPE",
            "G    1 C1C                                                     SYS / # / OBS TYPES",
            "G01",
        ),
    ];

    for case in cases {
        let repair = repair_obs_text(&case, &RepairOptions::default()).expect("repair OBS text");
        let output = repair.repaired.to_rinex_string();
        assert_strict_header_columns(&output);
        RinexObs::parse(&output).expect("repaired output strict parses");
    }
}

#[test]
fn obs_writer_wraps_scale_factors_and_prn_counts_without_truncation() {
    let obs = RinexObs::parse(&extended_count_obs_text()).expect("parse extended count headers");
    assert_eq!(obs.header().scale_factors[0].codes.len(), 13);
    assert_eq!(
        obs.header()
            .prn_obs_counts
            .get(&gps(1))
            .expect("G01 counts")
            .len(),
        13
    );

    let output = obs.to_rinex_string();
    assert_eq!(
        output
            .lines()
            .filter(|line| line.get(60..).unwrap_or("").trim() == "SYS / SCALE FACTOR")
            .count(),
        2
    );
    assert_eq!(
        output
            .lines()
            .filter(|line| line.get(60..).unwrap_or("").trim() == "PRN / # OF OBS")
            .count(),
        2
    );
    let reparsed = RinexObs::parse(&output).expect("reparse wrapped headers");
    assert_eq!(reparsed.header().scale_factors[0].codes.len(), 13);
    assert_eq!(
        reparsed
            .header()
            .prn_obs_counts
            .get(&gps(1))
            .expect("G01 counts")
            .len(),
        13
    );
}

#[test]
fn obs_type_whitespace_fallback_handles_continuation_and_rejects_surplus() {
    let compact_continuation = minimal_obs_text(
        "3.04 OBSERVATION DATA M RINEX VERSION / TYPE",
        "G 3 C1C SYS / # / OBS TYPES\nL1C D1C SYS / # / OBS TYPES",
        "G01",
    );
    let obs = RinexObs::parse(&compact_continuation).expect("parse compact continuation");
    assert_eq!(
        obs.obs_codes(GnssSystem::Gps).expect("GPS codes"),
        ["C1C".to_string(), "L1C".to_string(), "D1C".to_string()]
    );

    let strict_surplus_line = obs_header_line("G    1 C1C D1C", "SYS / # / OBS TYPES");
    let strict_surplus = minimal_obs_text(
        "     3.04           OBSERVATION DATA    M                   RINEX VERSION / TYPE",
        strict_surplus_line.trim_end(),
        "G01",
    );
    let strict_error = RinexObs::parse(&strict_surplus).expect_err("strict surplus must fail");
    let compact_surplus = minimal_obs_text(
        "3.04 OBSERVATION DATA M RINEX VERSION / TYPE",
        "G 1 C1C D1C SYS / # / OBS TYPES",
        "G01",
    );
    let compact_error = RinexObs::parse(&compact_surplus).expect_err("compact surplus must fail");
    match (strict_error, compact_error) {
        (Error::Parse(strict), Error::Parse(compact)) => {
            assert!(strict.contains("lists more codes than declared"));
            assert!(compact.contains("lists more codes than declared"));
        }
        (strict, compact) => panic!("expected parse errors, got {strict:?} and {compact:?}"),
    }
}

#[test]
fn non_padded_satellite_ids_parse_equivalently_in_obs_and_nav() {
    let obs_padded = minimal_obs_text(
        "     3.04           OBSERVATION DATA    M                   RINEX VERSION / TYPE",
        "G    1 C1C                                                  SYS / # / OBS TYPES",
        "G01",
    );
    let obs_non_padded = minimal_obs_text(
        "     3.04           OBSERVATION DATA    M                   RINEX VERSION / TYPE",
        "G    1 C1C                                                  SYS / # / OBS TYPES",
        "G 1",
    );
    let padded = RinexObs::parse(&obs_padded).expect("parse padded OBS");
    let non_padded = RinexObs::parse(&obs_non_padded).expect("parse non-padded OBS");
    let padded_sat = padded.epochs()[0].sats.keys().next().copied();
    let non_padded_sat = non_padded.epochs()[0].sats.keys().next().copied();
    assert_eq!(padded_sat, Some(gps(1)));
    assert_eq!(non_padded_sat, padded_sat);

    let nav_padded = include_str!("fixtures/nav/BRD400DLR_S_20261800000_01H_MN_trim.rnx");
    let nav_non_padded = nav_padded.replacen("G01 ", "G 1 ", 1);
    let padded_nav = parse_nav(nav_padded).expect("parse padded NAV");
    let non_padded_nav = parse_nav(&nav_non_padded).expect("parse non-padded NAV");
    assert_eq!(padded_nav[0].satellite_id, gps(1));
    assert_eq!(non_padded_nav[0].satellite_id, padded_nav[0].satellite_id);
}
