astrodynamics-gnss 0.9.2

//! Multi-source SP3 combination: clock-datum alignment across analysis centers.
//!
//! Precise clock products from different analysis centers are referenced to
//! different station/ensemble clocks, so their raw clock values differ by a
//! per-epoch common offset — the reference-clock difference — that drifts over
//! the day. Before clocks from two centers can be compared or combined, that
//! datum must be removed. [`clock_reference_offset`] estimates it robustly (the
//! median, over the satellites both products report at each epoch, of
//! `other - reference`); subtract it from `other`'s clocks to put both products
//! on `reference`'s datum.
//!
//! Orbit positions need no such treatment: every center reports ITRF
//! center-of-mass coordinates, so cross-center position differences are already
//! directly comparable.

use std::collections::{BTreeMap, BTreeSet};

use astrodynamics::time::model::Instant;

use super::interp::instant_to_j2000_seconds;
use super::{RawNode, Sp3, Sp3DataType, Sp3Flags, Sp3Header, Sp3State};
use crate::frame::ItrfPositionM;
use crate::id::GnssSatelliteId;
use crate::{Error, Result};

/// One epoch's reference-clock offset of `other` relative to `reference`.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct ClockReferenceOffset {
    /// The matched epoch.
    pub epoch: Instant,
    /// `other - reference` clock datum at this epoch, in seconds. Positive means
    /// `other`'s clock datum runs ahead of `reference`'s; subtract it from
    /// `other`'s clocks to align them to `reference`.
    pub offset_s: f64,
    /// Number of satellites that contributed to the (median) estimate.
    pub satellites: usize,
}

/// Estimate the per-epoch reference-clock offset of `other` relative to
/// `reference`.
///
/// For each epoch present in both products, the offset is the median over the
/// satellites both report (each with a finite clock) of
/// `other_clock - reference_clock`. The median makes the estimate robust to a
/// single satellite whose clock one center has wrong — but only with enough
/// satellites, so `min_common` is the minimum number of common clocked
/// satellites required to emit an offset for an epoch (a sound robust median
/// wants at least three, so one outlier can be outvoted). Epochs with fewer
/// common clocks are omitted rather than reported as a fragile one- or
/// two-satellite estimate.
///
/// Epochs are matched by their J2000 second floored to a whole second (the same
/// node-axis convention the interpolator uses). Non-finite clock differences are
/// skipped. Epochs present in only one product, or below `min_common`, are
/// omitted from the result.
pub fn clock_reference_offset(
    reference: &Sp3,
    other: &Sp3,
    min_common: usize,
) -> Vec<ClockReferenceOffset> {
    let mut other_index: std::collections::HashMap<i64, usize> = std::collections::HashMap::new();
    for (idx, epoch) in other.epochs.iter().enumerate() {
        if let Some(seconds) = instant_to_j2000_seconds(epoch) {
            other_index.insert(seconds.floor() as i64, idx);
        }
    }

    let mut offsets = Vec::new();

    for (ref_idx, epoch) in reference.epochs.iter().enumerate() {
        let Some(ref_seconds) = instant_to_j2000_seconds(epoch) else {
            continue;
        };
        let Some(&other_idx) = other_index.get(&(ref_seconds.floor() as i64)) else {
            continue;
        };

        let (Ok(ref_states), Ok(other_states)) =
            (reference.states_at(ref_idx), other.states_at(other_idx))
        else {
            continue;
        };

        let mut diffs: Vec<f64> = Vec::new();
        for (sat, ref_state) in ref_states.iter() {
            let Some(ref_clock) = ref_state.clock_s else {
                continue;
            };
            if let Some(other_state) = other_states.get(sat) {
                if let Some(other_clock) = other_state.clock_s {
                    let diff = other_clock - ref_clock;
                    // SP3 should not carry NaN/inf clocks, but the parser can
                    // accept them; merge infrastructure must not panic on data.
                    if diff.is_finite() {
                        diffs.push(diff);
                    }
                }
            }
        }

        if diffs.len() >= min_common.max(1) {
            if let Some(offset_s) = median(&mut diffs) {
                offsets.push(ClockReferenceOffset {
                    epoch: *epoch,
                    offset_s,
                    satellites: diffs.len(),
                });
            }
        }
    }

    offsets
}

fn median(values: &mut [f64]) -> Option<f64> {
    if values.is_empty() {
        return None;
    }

    // Inputs are pre-filtered to finite values; total_cmp never panics regardless.
    values.sort_by(f64::total_cmp);

    let n = values.len();
    if n % 2 == 1 {
        Some(values[n / 2])
    } else {
        Some((values[n / 2 - 1] + values[n / 2]) / 2.0)
    }
}

// ===========================================================================
// Multi-source merge
// ===========================================================================

/// How the agreeing (consensus) sources for a cell are combined into the merged
/// value.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum MergeCombine {
    /// Arithmetic mean of the consensus sources. The clustering step has already
    /// removed outliers, so the mean uses every agreeing measurement. Default.
    Mean,
    /// Component-wise median of the consensus sources.
    Median,
    /// The value from the highest-precedence (earliest-listed) consensus source.
    Precedence,
}

/// Options for [`merge`].
#[derive(Debug, Clone, PartialEq)]
pub struct MergeOptions {
    /// Maximum 3D position difference (meters) for two sources to be in
    /// agreement.
    pub position_tolerance_m: f64,
    /// Maximum clock difference (seconds, after datum alignment) for two sources
    /// to be in agreement.
    pub clock_tolerance_s: f64,
    /// Minimum number of mutually-agreeing sources required to accept a cell that
    /// has two or more sources. A cell with a single source is always carried
    /// through (gap fill, recorded as `single_source`); a cell with several
    /// sources but no agreeing subset this large is quarantined rather than
    /// averaged across disagreeing centers.
    pub min_agree: usize,
    /// Minimum common clocked satellites for the per-epoch clock-datum estimate
    /// between two sources (see [`clock_reference_offset`]).
    pub clock_min_common: usize,
    /// How to combine the agreeing sources.
    pub combine: MergeCombine,
}

impl Default for MergeOptions {
    /// Defaults tuned for the common case of ~3 analysis centers: agreement is a
    /// 2-of-3 majority (`min_agree = 2`); combine the agreeing subset by mean.
    fn default() -> Self {
        Self {
            position_tolerance_m: 0.5,
            clock_tolerance_s: 5.0e-9,
            min_agree: 2,
            clock_min_common: 5,
            combine: MergeCombine::Mean,
        }
    }
}

/// One (epoch, satellite) cell the merge handled with a caveat. Nothing is
/// dropped or averaged silently — every such cell is recorded here.
#[derive(Debug, Clone, PartialEq)]
pub struct MergeFlag {
    /// The epoch.
    pub epoch: Instant,
    /// The satellite.
    pub satellite: GnssSatelliteId,
    /// The source indices (into the input slice) this flag refers to: for
    /// `single_source`, the lone contributor; for `quarantined`, all sources
    /// that disagreed; for `position_outliers`, the sources rejected from an
    /// otherwise-accepted consensus.
    pub sources: Vec<usize>,
}

/// Audit trail for a [`merge`].
#[derive(Debug, Clone, Default, PartialEq)]
pub struct MergeReport {
    /// Cells where two or more sources disagreed beyond tolerance with no
    /// agreeing subset of `min_agree` — omitted from the merged product.
    pub quarantined: Vec<MergeFlag>,
    /// Cells carried from a single source (no cross-check was possible).
    pub single_source: Vec<MergeFlag>,
    /// Cells accepted by consensus where one or more sources were rejected as
    /// position outliers.
    pub position_outliers: Vec<MergeFlag>,
}

/// Merge several SP3 products from different analysis centers into one
/// consistent precise-ephemeris dataset.
///
/// Orthogonal to time-stitching: this combines providers at the **same** epochs.
/// For every (epoch, satellite) cell in the union of the inputs:
///
/// - **Union coverage.** A satellite present in any input is present in the
///   output for that epoch, filling a single center's dropouts.
/// - **Position consensus.** With one source the value is carried through
///   (`single_source`). With several, the largest subset of sources mutually
///   within `position_tolerance_m` is found; if it has at least `min_agree`
///   members it is combined per `combine` and any sources outside it are recorded
///   as `position_outliers`. If no such subset exists the cell is `quarantined`
///   (omitted) — never averaged across disagreeing centers.
/// - **Clock consensus.** Clocks are first put on a common datum (each source
///   aligned to the first via [`clock_reference_offset`]), then combined by the
///   same agreement rule; a cell with no clock consensus carries no clock. A
///   non-reference source whose datum cannot be estimated at an epoch (below
///   `clock_min_common` common clocks) contributes **no** clock there rather than
///   an unaligned one — its position is still merged.
///
/// All inputs must share a time scale and coordinate system (epochs are matched
/// across products in that scale); mismatches are rejected. The merged record
/// flags are the union (OR) of the contributing sources' flags — in particular a
/// `clock_event` on any clock-consensus member is preserved, so the interpolator
/// still splits the clock arc. The merged header is **synthetic**: its first-epoch
/// fields describe the union's first epoch and its data type is position-only.
///
/// Pure and deterministic: order the inputs by center precedence and ties (equal
/// cluster sizes, `Precedence` combine) resolve to the earliest-listed source.
/// The merged product's interpolation nodes are the consensus values, so it
/// samples and interpolates like any other [`Sp3`] (it is a derived combination,
/// not a byte-faithful copy of any one center).
pub fn merge(sources: &[Sp3], opts: &MergeOptions) -> Result<(Sp3, MergeReport)> {
    if sources.is_empty() {
        return Err(Error::InvalidInput(
            "merge requires at least one SP3 product".into(),
        ));
    }

    // Inputs must be combinable: epochs are matched in a common time scale, and
    // positions are only comparable in a common coordinate system / frame.
    let base = &sources[0].header;
    for s in &sources[1..] {
        if s.header.time_scale != base.time_scale {
            return Err(Error::InvalidInput(format!(
                "merge inputs have mismatched time scales ({:?} vs {:?})",
                base.time_scale, s.header.time_scale
            )));
        }
        if s.header.coordinate_system != base.coordinate_system {
            return Err(Error::InvalidInput(format!(
                "merge inputs have mismatched coordinate systems ({:?} vs {:?})",
                base.coordinate_system, s.header.coordinate_system
            )));
        }
    }

    // floored-J2000-second -> epoch index, per source.
    let epoch_index: Vec<BTreeMap<i64, usize>> = sources
        .iter()
        .map(|s| {
            s.epochs
                .iter()
                .enumerate()
                .filter_map(|(i, ep)| {
                    instant_to_j2000_seconds(ep).map(|sec| (sec.floor() as i64, i))
                })
                .collect()
        })
        .collect();

    // Per-source per-epoch clock-datum offset relative to source 0. Source 0 is
    // the datum, so its offset is identically zero.
    let clock_offset: Vec<BTreeMap<i64, f64>> = sources
        .iter()
        .enumerate()
        .map(|(idx, s)| {
            if idx == 0 {
                BTreeMap::new()
            } else {
                clock_reference_offset(&sources[0], s, opts.clock_min_common)
                    .into_iter()
                    .filter_map(|o| {
                        instant_to_j2000_seconds(&o.epoch)
                            .map(|sec| (sec.floor() as i64, o.offset_s))
                    })
                    .collect()
            }
        })
        .collect();

    // Union of epochs (by floored second), keeping a representative Instant.
    let mut epoch_keys: BTreeMap<i64, Instant> = BTreeMap::new();
    for s in sources {
        for ep in &s.epochs {
            if let Some(sec) = instant_to_j2000_seconds(ep) {
                epoch_keys.entry(sec.floor() as i64).or_insert(*ep);
            }
        }
    }

    let mut out_epochs: Vec<Instant> = Vec::with_capacity(epoch_keys.len());
    let mut out_states: Vec<BTreeMap<GnssSatelliteId, Sp3State>> =
        Vec::with_capacity(epoch_keys.len());
    let mut out_raw: Vec<BTreeMap<GnssSatelliteId, RawNode>> = Vec::with_capacity(epoch_keys.len());
    let mut report = MergeReport::default();
    let mut all_sats: BTreeSet<GnssSatelliteId> = BTreeSet::new();

    for (&key, &epoch) in &epoch_keys {
        out_epochs.push(epoch);
        let mut states: BTreeMap<GnssSatelliteId, Sp3State> = BTreeMap::new();
        let mut raws: BTreeMap<GnssSatelliteId, RawNode> = BTreeMap::new();

        // Satellites present at this epoch in any source.
        let mut sats: BTreeSet<GnssSatelliteId> = BTreeSet::new();
        for (idx, s) in sources.iter().enumerate() {
            if let Some(&ei) = epoch_index[idx].get(&key) {
                if let Ok(map) = s.states_at(ei) {
                    sats.extend(map.keys().copied());
                }
            }
        }

        for sat in sats {
            // (source_idx, position_m, flags) and (source_idx, datum-aligned
            // clock_s, flags). A non-reference source contributes a clock only
            // when its datum offset could be estimated at this epoch; otherwise
            // its clock would be unaligned, so it is omitted (the position is
            // still gathered).
            let mut pos: Vec<(usize, [f64; 3], Sp3Flags)> = Vec::new();
            let mut clk: Vec<(usize, f64, Sp3Flags)> = Vec::new();
            for (idx, s) in sources.iter().enumerate() {
                let Some(&ei) = epoch_index[idx].get(&key) else {
                    continue;
                };
                let Ok(map) = s.states_at(ei) else { continue };
                let Some(state) = map.get(&sat) else { continue };
                pos.push((idx, state.position.as_array(), state.flags));
                if let Some(c) = state.clock_s {
                    let offset = if idx == 0 {
                        Some(0.0)
                    } else {
                        clock_offset[idx].get(&key).copied()
                    };
                    if let Some(off) = offset {
                        let aligned = c - off;
                        if aligned.is_finite() {
                            clk.push((idx, aligned, state.flags));
                        }
                    }
                }
            }

            let flag = |srcs: Vec<usize>| MergeFlag {
                epoch,
                satellite: sat,
                sources: srcs,
            };

            // Position consensus -> the merged position and the indices (into
            // `pos`) of the sources that contributed it.
            let (position_m, pos_members) = if pos.len() == 1 {
                report.single_source.push(flag(vec![pos[0].0]));
                (pos[0].1, vec![0usize])
            } else {
                let pts: Vec<[f64; 3]> = pos.iter().map(|(_, p, _)| *p).collect();
                let cluster = largest_within(&pts, |a, b| dist3(a, b) <= opts.position_tolerance_m);
                if cluster.len() >= opts.min_agree {
                    let rejected: Vec<usize> = (0..pos.len())
                        .filter(|i| !cluster.contains(i))
                        .map(|i| pos[i].0)
                        .collect();
                    if !rejected.is_empty() {
                        report.position_outliers.push(flag(rejected));
                    }
                    let members: Vec<(usize, [f64; 3])> =
                        cluster.iter().map(|&i| (pos[i].0, pos[i].1)).collect();
                    (combine3(&members, opts.combine), cluster)
                } else {
                    report
                        .quarantined
                        .push(flag(pos.iter().map(|(i, _, _)| *i).collect()));
                    continue;
                }
            };

            // Clock consensus, independent of position -> the merged clock and the
            // indices (into `clk`) of the sources that contributed it.
            let (clock_s, clk_members): (Option<f64>, Vec<usize>) = if clk.is_empty() {
                (None, Vec::new())
            } else if clk.len() == 1 {
                (Some(clk[0].1), vec![0usize])
            } else {
                let vals: Vec<f64> = clk.iter().map(|(_, c, _)| *c).collect();
                let cluster = largest_within(&vals, |a, b| (a - b).abs() <= opts.clock_tolerance_s);
                if cluster.len() >= opts.min_agree {
                    let members: Vec<(usize, f64)> =
                        cluster.iter().map(|&i| (clk[i].0, clk[i].1)).collect();
                    (Some(combine_axis(&members, opts.combine)), cluster)
                } else {
                    (None, Vec::new())
                }
            };

            // Preserve record flags: OR the orbit flags across the position
            // members and the clock flags across the clock members, so a
            // `clock_event` (clock reset) or maneuver on any contributing source
            // survives into the merged product.
            let mut flags = Sp3Flags::default();
            for &i in &pos_members {
                flags.maneuver |= pos[i].2.maneuver;
                flags.orbit_predicted |= pos[i].2.orbit_predicted;
            }
            for &i in &clk_members {
                flags.clock_event |= clk[i].2.clock_event;
                flags.clock_predicted |= clk[i].2.clock_predicted;
            }

            all_sats.insert(sat);
            states.insert(
                sat,
                Sp3State {
                    position: ItrfPositionM::new(position_m[0], position_m[1], position_m[2]),
                    clock_s,
                    velocity: None,
                    clock_rate_s_s: None,
                    flags,
                },
            );
            raws.insert(
                sat,
                RawNode {
                    km: [
                        position_m[0] / 1000.0,
                        position_m[1] / 1000.0,
                        position_m[2] / 1000.0,
                    ],
                    clock_us: clock_s.map(|c| c * 1.0e6),
                    clock_event: flags.clock_event,
                },
            );
        }

        out_states.push(states);
        out_raw.push(raws);
    }

    // Base the header on the source that owns the union's first epoch, so its
    // first-epoch fields (week / seconds-of-week / MJD) describe `out_epochs[0]`
    // rather than being stale from source 0. The merged product carries no
    // velocity records, so it is position-only regardless of the inputs' type.
    let base_idx = sources
        .iter()
        .enumerate()
        .filter_map(|(i, s)| {
            s.epochs
                .first()
                .and_then(instant_to_j2000_seconds)
                .map(|sec| (sec, i))
        })
        .min_by(|a, b| a.0.total_cmp(&b.0).then(a.1.cmp(&b.1)))
        .map(|(_, i)| i)
        .unwrap_or(0);

    // Recompute the nominal epoch spacing from the merged epochs rather than
    // inheriting it: the union of mixed-cadence products has its own (finer)
    // spacing, so the base header's interval would be stale. Use the smallest
    // positive gap between consecutive epochs (the finest cadence present),
    // falling back to the base header for a degenerate single-epoch product.
    let epoch_interval_s = {
        let secs: Vec<f64> = out_epochs
            .iter()
            .filter_map(instant_to_j2000_seconds)
            .collect();
        secs.windows(2)
            .map(|w| w[1] - w[0])
            .filter(|g| *g > 0.0)
            .min_by(f64::total_cmp)
            .unwrap_or(sources[base_idx].header.epoch_interval_s)
    };

    let header = Sp3Header {
        num_epochs: out_epochs.len() as u64,
        satellites: all_sats.into_iter().collect(),
        data_type: Sp3DataType::Position,
        epoch_interval_s,
        ..sources[base_idx].header.clone()
    };

    let merged = Sp3 {
        header,
        epochs: out_epochs,
        states: out_states,
        interp_raw: out_raw,
        comments: vec![format!("MERGED from {} SP3 products", sources.len())],
    };

    Ok((merged, report))
}

fn dist3(a: &[f64; 3], b: &[f64; 3]) -> f64 {
    let dx = a[0] - b[0];
    let dy = a[1] - b[1];
    let dz = a[2] - b[2];
    (dx * dx + dy * dy + dz * dz).sqrt()
}

/// Indices of the largest subset of `items` whose members are *mutually* within
/// `within`. Brute-force max-clique over the agreement graph; the item count is
/// the number of input products (a handful), so the `2^n` enumeration is fine.
/// Ties resolve to the lowest-indexed subset (precedence).
fn largest_within<T>(items: &[T], within: impl Fn(&T, &T) -> bool) -> Vec<usize> {
    let n = items.len();
    if n <= 1 {
        return (0..n).collect();
    }
    let mut best: Vec<usize> = vec![0];
    for mask in 1u32..(1u32 << n) {
        let members: Vec<usize> = (0..n).filter(|i| mask & (1 << i) != 0).collect();
        if members.len() <= best.len() {
            continue;
        }
        let consistent = members.iter().enumerate().all(|(mi, &i)| {
            members[mi + 1..]
                .iter()
                .all(|&j| within(&items[i], &items[j]))
        });
        if consistent {
            best = members;
        }
    }
    best
}

fn combine3(members: &[(usize, [f64; 3])], how: MergeCombine) -> [f64; 3] {
    [0usize, 1, 2].map(|axis| {
        let axis_members: Vec<(usize, f64)> = members.iter().map(|(s, v)| (*s, v[axis])).collect();
        combine_axis(&axis_members, how)
    })
}

fn combine_axis(members: &[(usize, f64)], how: MergeCombine) -> f64 {
    match how {
        MergeCombine::Mean => members.iter().map(|(_, v)| *v).sum::<f64>() / members.len() as f64,
        MergeCombine::Median => {
            let mut vals: Vec<f64> = members.iter().map(|(_, v)| *v).collect();
            median(&mut vals).expect("consensus cluster is non-empty")
        }
        MergeCombine::Precedence => members
            .iter()
            .min_by_key(|(s, _)| *s)
            .map(|(_, v)| *v)
            .expect("consensus cluster is non-empty"),
    }
}

/// Return a copy of `other` with its clocks shifted onto `reference`'s clock
/// datum.
///
/// This applies the per-epoch reference-clock offset from
/// [`clock_reference_offset`]: at each epoch where the offset could be estimated
/// (at least `min_common` common clocked satellites), every clocked satellite's
/// offset has the datum subtracted, so the result's clocks are directly
/// comparable to `reference`'s. Positions are untouched (already comparable).
///
/// Epochs where the offset could not be estimated are left unchanged — they are
/// *not* on `reference`'s datum, so a caller mixing aligned and unaligned epochs
/// should consult [`clock_reference_offset`] to see which epochs were aligned.
/// The returned product interpolates like any other [`Sp3`].
pub fn align_clock_reference(reference: &Sp3, other: &Sp3, min_common: usize) -> Sp3 {
    let offsets: BTreeMap<i64, f64> = clock_reference_offset(reference, other, min_common)
        .into_iter()
        .filter_map(|o| {
            instant_to_j2000_seconds(&o.epoch).map(|sec| (sec.floor() as i64, o.offset_s))
        })
        .collect();

    let mut aligned = other.clone();
    for ei in 0..aligned.epochs.len() {
        let Some(sec) = instant_to_j2000_seconds(&aligned.epochs[ei]) else {
            continue;
        };
        let Some(&off) = offsets.get(&(sec.floor() as i64)) else {
            continue;
        };
        for state in aligned.states[ei].values_mut() {
            if let Some(c) = state.clock_s.as_mut() {
                *c -= off;
            }
        }
        for node in aligned.interp_raw[ei].values_mut() {
            if let Some(us) = node.clock_us.as_mut() {
                *us -= off * 1.0e6;
            }
        }
    }
    aligned
}

#[cfg(test)]
mod tests {
    use super::super::Sp3;
    use super::{align_clock_reference, clock_reference_offset, merge, MergeOptions};
    use crate::id::{GnssSatelliteId, GnssSystem};

    fn gps(prn: u8) -> GnssSatelliteId {
        GnssSatelliteId::new(GnssSystem::Gps, prn)
    }

    // Single-epoch SP3-c from explicit `(satellite, [x,y,z] km, clock us, flag
    // suffix)` records under coordinate system `cs` (5 chars, e.g. `"IGS14"`).
    // `flags` is appended verbatim after the 60-column record body, so a test can
    // place an SP3 flag (e.g. `"              E"` -> the `E` clock-event flag at
    // column 75). A `None` clock writes the SP3 bad-clock sentinel.
    fn sp3_build(records: &[(&str, [f64; 3], Option<f64>, &str)], cs: &str) -> Sp3 {
        let n = records.len();
        let mut sats = String::new();
        for (sat, _, _, _) in records {
            sats.push_str(sat);
        }
        for _ in n..17 {
            sats.push_str("  0");
        }
        let mut body = String::new();
        body.push_str(&format!(
            "#cP2020  6 25  0  0  0.00000000       1 ORBIT {cs} FIT  TST\n"
        ));
        body.push_str("## 2111 432000.00000000   900.00000000 59025 0.0000000000000\n");
        body.push_str(&format!("+   {n:2}   {sats}\n"));
        body.push_str("++         0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0\n");
        body.push_str("%c G  cc GPS ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc\n");
        body.push_str("%c cc cc ccc ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc\n");
        body.push_str("%f  1.2500000  1.025000000  0.00000000000  0.000000000000000\n");
        body.push_str("%f  0.0000000  0.000000000  0.00000000000  0.000000000000000\n");
        body.push_str("%i    0    0    0    0      0      0      0      0         0\n");
        body.push_str("%i    0    0    0    0      0      0      0      0         0\n");
        body.push_str("/* TEST SP3-c FIXTURE\n");
        body.push_str("*  2020  6 25  0  0  0.00000000\n");
        for (sat, p, clk, flags) in records {
            let c = clk.unwrap_or(999_999.999_999);
            body.push_str(&format!(
                "P{sat}{:14.6}{:14.6}{:14.6}{c:14.6}{flags}\n",
                p[0], p[1], p[2]
            ));
        }
        body.push_str("EOF\n");
        Sp3::parse(body.as_bytes()).expect("parse test sp3")
    }

    // The common case: IGS14, no flags.
    fn sp3_records(records: &[(&str, [f64; 3], Option<f64>)]) -> Sp3 {
        let full: Vec<(&str, [f64; 3], Option<f64>, &str)> =
            records.iter().map(|(s, p, c)| (*s, *p, *c, "")).collect();
        sp3_build(&full, "IGS14")
    }

    #[test]
    fn merge_unions_coverage_when_one_center_misses_a_satellite() {
        // Center A reports G01/G02/G03; center B is missing G03. The merged
        // product must still cover G03 at that epoch (filled from A).
        let a = sp3_records(&[
            ("G01", [15000.0, -20000.0, 5000.0], Some(100.0)),
            ("G02", [16000.0, -21000.0, 6000.0], Some(200.0)),
            ("G03", [17000.0, -22000.0, 7000.0], Some(300.0)),
        ]);
        let b = sp3_records(&[
            ("G01", [15000.0, -20000.0, 5000.0], Some(100.0)),
            ("G02", [16000.0, -21000.0, 6000.0], Some(200.0)),
        ]);

        let (merged, report) = merge(&[a, b], &MergeOptions::default()).expect("merge");

        let states = merged.states_at(0).expect("epoch 0");
        assert!(
            states.contains_key(&gps(3)),
            "merged output must cover G03 from the center that has it"
        );
        assert_eq!(states.len(), 3, "union is G01/G02/G03");
        // G01 agreed across both centers -> consensus clock is their value.
        let g01 = states[&gps(1)];
        assert!((g01.clock_s.unwrap() - 100.0e-6).abs() < 1.0e-15);
        // G03 had a single source -> carried through, recorded, not quarantined.
        assert!(report.quarantined.is_empty());
        assert_eq!(report.single_source.len(), 1);
        assert_eq!(report.single_source[0].satellite, gps(3));
    }

    #[test]
    fn merge_combines_two_of_three_agreeing_sources_and_rejects_the_outlier() {
        // A and B agree on G01; C is 10 m off in X (> the default 0.5 m tolerance).
        let a = sp3_records(&[("G01", [15000.0, -20000.0, 5000.0], Some(100.0))]);
        let b = sp3_records(&[("G01", [15000.0, -20000.0, 5000.0], Some(100.0))]);
        let c = sp3_records(&[("G01", [15000.010, -20000.0, 5000.0], Some(100.0))]);

        let (merged, report) = merge(&[a, b, c], &MergeOptions::default()).expect("merge");

        let states = merged.states_at(0).expect("epoch 0");
        let g01 = states[&gps(1)];
        // Consensus is A/B (15000 km == 1.5e7 m); not dragged toward C.
        assert!(
            (g01.position.as_array()[0] - 15_000_000.0).abs() < 1.0e-3,
            "got {}",
            g01.position.as_array()[0]
        );
        // C is source index 2 -> recorded as the rejected position outlier.
        assert_eq!(report.position_outliers.len(), 1);
        assert_eq!(report.position_outliers[0].sources, vec![2]);
        assert!(report.quarantined.is_empty());
    }

    #[test]
    fn merge_quarantines_a_satellite_all_centers_disagree_on() {
        // Three sources, mutually beyond tolerance on G01: no 2-of-3 consensus.
        let a = sp3_records(&[("G01", [15000.000, -20000.0, 5000.0], Some(100.0))]);
        let b = sp3_records(&[("G01", [15000.010, -20000.0, 5000.0], Some(100.0))]);
        let c = sp3_records(&[("G01", [15000.020, -20000.0, 5000.0], Some(100.0))]);

        let (merged, report) = merge(&[a, b, c], &MergeOptions::default()).expect("merge");

        assert!(
            merged.states_at(0).expect("epoch 0").is_empty(),
            "no consensus -> G01 omitted, not averaged across disagreeing centers"
        );
        assert_eq!(report.quarantined.len(), 1);
        assert_eq!(report.quarantined[0].satellite, gps(1));
    }

    #[test]
    fn merge_rejects_an_empty_input() {
        assert!(merge(&[], &MergeOptions::default()).is_err());
    }

    #[test]
    fn merge_omits_an_unalignable_secondary_clock() {
        // Only 3 common satellites, but the default clock datum needs 5, so
        // center B's clocks cannot be put on A's datum. They must be dropped
        // rather than emitted raw, and a B-only satellite gets a position but no
        // clock.
        let a = sp3_records(&[
            ("G01", [15000.0, -20000.0, 5000.0], Some(100.0)),
            ("G02", [16000.0, -21000.0, 6000.0], Some(200.0)),
            ("G03", [17000.0, -22000.0, 7000.0], Some(300.0)),
        ]);
        let b = sp3_records(&[
            ("G01", [15000.0, -20000.0, 5000.0], Some(150.0)),
            ("G02", [16000.0, -21000.0, 6000.0], Some(250.0)),
            ("G03", [17000.0, -22000.0, 7000.0], Some(350.0)),
            ("G04", [18000.0, -23000.0, 8000.0], Some(450.0)),
        ]);

        let (merged, _) = merge(&[a, b], &MergeOptions::default()).expect("merge");
        let states = merged.states_at(0).expect("epoch 0");

        // G04 is B-only (gap fill): position carried, clock unalignable -> dropped.
        assert!(states.contains_key(&gps(4)));
        assert!(
            states[&gps(4)].clock_s.is_none(),
            "an unalignable secondary clock must be dropped, not emitted raw"
        );
        // G01's clock comes from the reference (source 0), which is on its own datum.
        let g01_clock = states[&gps(1)]
            .clock_s
            .expect("G01 carries the reference clock");
        assert!((g01_clock - 100.0e-6).abs() < 1.0e-12, "got {g01_clock}");
    }

    #[test]
    fn merge_rejects_mismatched_coordinate_systems() {
        let a = sp3_build(
            &[("G01", [15000.0, -20000.0, 5000.0], Some(100.0), "")],
            "IGS14",
        );
        let b = sp3_build(
            &[("G01", [15000.0, -20000.0, 5000.0], Some(100.0), "")],
            "IGS20",
        );

        assert!(merge(&[a, b], &MergeOptions::default()).is_err());
    }

    #[test]
    fn merge_preserves_a_clock_event_flag() {
        // Source A carries an `E` clock-event flag on G01 (column 75); the merged
        // product must keep it so the interpolator still splits the clock arc.
        let a = sp3_build(
            &[(
                "G01",
                [15000.0, -20000.0, 5000.0],
                Some(100.0),
                "              E",
            )],
            "IGS14",
        );
        let b = sp3_build(
            &[("G01", [15000.0, -20000.0, 5000.0], Some(100.0), "")],
            "IGS14",
        );

        let (merged, _) = merge(&[a, b], &MergeOptions::default()).expect("merge");
        let g01 = merged.states_at(0).expect("epoch 0")[&gps(1)];

        assert!(
            g01.flags.clock_event,
            "merged cell must preserve a contributing source's clock-event flag"
        );
    }

    #[test]
    fn merge_recomputes_epoch_interval_from_the_union() {
        // The header DECLARES a 300 s interval, but the two epochs are 15 min
        // (900 s) apart. The synthetic merged header must report the spacing of
        // the actual merged epochs, not inherit the stale declared value.
        let body = "#cP2020  6 25  0  0  0.00000000       2 ORBIT IGS14 FIT  TST\n\
            ## 2111 432000.00000000   300.00000000 59025 0.0000000000000\n\
            +    1   G01  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0\n\
            ++         0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0\n\
            %c G  cc GPS ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc\n\
            %c cc cc ccc ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc\n\
            %f  1.2500000  1.025000000  0.00000000000  0.000000000000000\n\
            %f  0.0000000  0.000000000  0.00000000000  0.000000000000000\n\
            %i    0    0    0    0      0      0      0      0         0\n\
            %i    0    0    0    0      0      0      0      0         0\n\
            /* TEST SP3-c FIXTURE\n\
            *  2020  6 25  0  0  0.00000000\n\
            PG01  15000.000000 -20000.000000   5000.000000    100.000000\n\
            *  2020  6 25  0 15  0.00000000\n\
            PG01  15001.000000 -20001.000000   5001.000000    101.000000\n\
            EOF\n";
        let a = Sp3::parse(body.as_bytes()).expect("parse test sp3");

        let (merged, _) = merge(&[a], &MergeOptions::default()).expect("merge");

        assert!(
            (merged.header.epoch_interval_s - 900.0).abs() < 1.0e-6,
            "got {}",
            merged.header.epoch_interval_s
        );
    }

    #[test]
    fn align_clock_reference_puts_other_on_the_reference_datum() {
        // `other`'s clocks all run +50 us ahead; after alignment they should sit
        // on `reference`'s datum (G01: 150 us - 50 us = 100 us = 1e-4 s).
        let reference = sp3([100.0, 200.0, 300.0]);
        let other = sp3([150.0, 250.0, 350.0]);

        let aligned = align_clock_reference(&reference, &other, 3);

        let g01 = aligned.states_at(0).expect("epoch 0")[&gps(1)];
        assert!(
            (g01.clock_s.unwrap() - 100.0e-6).abs() < 1.0e-15,
            "got {}",
            g01.clock_s.unwrap()
        );
        // Positions are untouched by clock alignment.
        let original = other.states_at(0).expect("epoch 0")[&gps(1)];
        assert_eq!(g01.position.as_array(), original.position.as_array());
    }

    // Minimal single-epoch SP3-c with three satellites; each `clocks_us` entry is
    // that satellite's clock in microseconds (positions are arbitrary but non-zero
    // so they parse as valid records).
    fn sp3(clocks_us: [f64; 3]) -> Sp3 {
        let body = format!(
            "#cP2020  6 25  0  0  0.00000000       1 ORBIT IGS14 FIT  TST\n\
             ## 2111 432000.00000000   900.00000000 59025 0.0000000000000\n\
             +    3   G01G02G03  0  0  0  0  0  0  0  0  0  0  0  0  0  0\n\
             ++         0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0\n\
             %c G  cc GPS ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc\n\
             %c cc cc ccc ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc\n\
             %f  1.2500000  1.025000000  0.00000000000  0.000000000000000\n\
             %f  0.0000000  0.000000000  0.00000000000  0.000000000000000\n\
             %i    0    0    0    0      0      0      0      0         0\n\
             %i    0    0    0    0      0      0      0      0         0\n\
             /* TEST SP3-c FIXTURE\n\
             *  2020  6 25  0  0  0.00000000\n\
             PG01  15000.000000 -20000.000000   5000.000000 {:13.6}\n\
             PG02  -1234.567890   2345.678901  -3456.789012 {:13.6}\n\
             PG03   8000.000000  12000.000000 -19000.000000 {:13.6}\n\
             EOF\n",
            clocks_us[0], clocks_us[1], clocks_us[2]
        );
        Sp3::parse(body.as_bytes()).expect("parse test sp3")
    }

    #[test]
    fn recovers_a_uniform_datum_shift() {
        // every `other` clock is +50 us (= 5e-5 s) from `reference`.
        let reference = sp3([100.0, 200.0, 300.0]);
        let other = sp3([150.0, 250.0, 350.0]);

        let offsets = clock_reference_offset(&reference, &other, 3);

        assert_eq!(offsets.len(), 1);
        assert_eq!(offsets[0].satellites, 3);
        assert!(
            (offsets[0].offset_s - 5.0e-5).abs() < 1.0e-12,
            "got {}",
            offsets[0].offset_s
        );
    }

    #[test]
    fn median_rejects_a_single_outlier_clock() {
        // Two satellites agree (+50 us); one is a wild outlier (+9000 us). The
        // median over the three tracks the consensus instead of being dragged out.
        let reference = sp3([100.0, 200.0, 300.0]);
        let other = sp3([150.0, 250.0, 9_300.0]);

        let offsets = clock_reference_offset(&reference, &other, 3);

        assert_eq!(offsets.len(), 1);
        assert!(
            (offsets[0].offset_s - 5.0e-5).abs() < 1.0e-12,
            "got {}",
            offsets[0].offset_s
        );
    }

    #[test]
    fn omits_epochs_below_min_common() {
        // Three common clocked satellites, but require four: the fragile estimate
        // is omitted rather than reported.
        let reference = sp3([100.0, 200.0, 300.0]);
        let other = sp3([150.0, 250.0, 350.0]);

        assert!(clock_reference_offset(&reference, &other, 4).is_empty());
    }
}