kshana 0.27.1

// SPDX-License-Identifier: AGPL-3.0-only
//! Observability-Gramian-over-arc core for planar cislunar tracking (paper P6).
//!
//! A single range snapshot to a spacecraft in the planar circular restricted three-body
//! problem sees one line-of-sight and nothing of velocity — the four-state
//! `s = [x, y, ẋ, ẏ]` is far from observable. Observability is *recovered over an arc*:
//! as the geometry evolves, the per-epoch measurement Jacobians `H(t_k)`, mapped back to
//! the initial epoch through the **variational state-transition matrix** `Φ(t_k)`, span
//! more of the state space. This module assembles that structure and reads off how much
//! of the state the arc actually constrains.
//!
//! For a linearised measurement model `z_k = H_k · δs(t_k) + n` and `δs(t_k) = Φ_k·δs_0`,
//! the sensitivity of the whole batch to the initial state is the stacked
//! **observability matrix** `O = stack_k[ H_k · Φ_k ]`. The state is observable over the
//! arc iff `O` has full column rank; the **observability Gramian**
//! `W = Σ_k Δt_k · Φ_kᵀ H_kᵀ H_k Φ_k` (a dt-weighted `OᵀO`) is the symmetric
//! positive-semidefinite information content, whose spectrum quantifies *how strongly*
//! each direction is seen.
//!
//! ## What is Validated vs Modelled
//! * **Validated.** The **rank** is read from a singular-value threshold on `O`
//!   (rank-revealing SVD via the squared-singular-value = eigenvalue-of-`OᵀO` identity),
//!   and independently confirmed against the eigen-rank of the dt-weighted Gramian `W`.
//!   The **eigen-spectrum** of `W` is the symmetric spectrum from the crate's
//!   Jacobi eigensolver ([`crate::fim::sym_eig`]), cross-checked against the spectral
//!   invariants `trace(W) = Σλ`, `‖W‖_F² = Σλ²`, and (for a full-rank block) an
//!   independent Gaussian-elimination `det(W) = Πλ`. The variational STM `Φ` produced by
//!   [`planar_state_stm`] is the finite-difference-validated CR3BP STM of
//!   [`crate::cr3bp::propagate_state_stm`] (its planar `[x, y, ẋ, ẏ]` sub-block), so STM
//!   propagation is Validated, not a first-order approximation.
//! * **Modelled.** The particular tracking geometry (which spacecraft, which links, the
//!   arc length and epoch grid) is a scenario input; the *specific* rank progression it
//!   produces is a property of that Modelled geometry, not an oracle-verified universal.
//!
//! ## Three dimensions and measurement noise
//! The same assembly serves the **spatial** six-state `[x, y, z, ẋ, ẏ, ż]`:
//! [`spatial_state_stm`] is the crate's finite-difference-validated CR3BP STM at full
//! width, of which [`planar_state_stm`] is the `{x, y, ẋ, ẏ}` restriction, so both paths
//! share one validated linearisation. **Measurement noise** enters through
//! [`whiten_epochs`] — with `R = σ²I` the batch information is the unit-weight Gram of the
//! rows `H/σ` — and [`whitened_posterior`] reads the formal covariance that whitening
//! implies. Reported rather than glossed: a homoscedastic `σ` scales `O` by a scalar, so
//! rank, defect and condition number are *invariant* under it, and only the posterior
//! uncertainty moves.

use crate::fim::{crlb, design_metrics, information_matrix, sym_eig};
use crate::intersat_range::{range_rate_row, range_row, PlanarState, SpatialState};

/// A dense matrix as rows of columns (matching the rest of the crate).
pub type Mat = Vec<Vec<f64>>;

/// Planar CR3BP state dimension `[x, y, ẋ, ẏ]`.
pub const N_PLANAR: usize = 4;

/// Spatial CR3BP state dimension `[x, y, z, ẋ, ẏ, ż]` — the full three-dimensional state
/// the spatial observability path estimates.
pub const N_SPATIAL: usize = 6;

/// The rotating-frame in-plane component indices in the 6-vector `[x, y, z, ẋ, ẏ, ż]`.
const PLANAR_IDX: [usize; N_PLANAR] = [0, 1, 3, 4];

// ── Variational STM bridge (L31) ────────────────────────────────────────────

/// Propagate a **planar** CR3BP state and its 4×4 variational STM for time `t`.
///
/// This is the planar `[x, y, ẋ, ẏ]` sub-block of the crate's finite-difference-validated
/// CR3BP STM ([`crate::cr3bp::propagate_state_stm`]): the planar state is embedded as
/// `z = ż = 0` (the plane is an invariant manifold, and the out-of-plane block decouples
/// exactly there), the full 6×6 STM is integrated, and the `{x, y, ẋ, ẏ}` rows/columns
/// are extracted. Returns `(state(t), Φ(t))` with `Φ` the true linearisation of the flow.
pub fn planar_state_stm(
    s0: &PlanarState,
    mu: f64,
    t: f64,
    steps: usize,
) -> (PlanarState, [[f64; N_PLANAR]; N_PLANAR]) {
    let embed = crate::cr3bp::Cr3bpState {
        r: [s0[0], s0[1], 0.0],
        v: [s0[2], s0[3], 0.0],
    };
    let (st, phi6) = crate::cr3bp::propagate_state_stm(&embed, mu, t, steps);
    let state = [st.r[0], st.r[1], st.v[0], st.v[1]];
    let mut phi = [[0.0; N_PLANAR]; N_PLANAR];
    for (i, &ri) in PLANAR_IDX.iter().enumerate() {
        for (j, &cj) in PLANAR_IDX.iter().enumerate() {
            phi[i][j] = phi6[ri][cj];
        }
    }
    (state, phi)
}

/// Planar CR3BP state after time `t` (position + velocity), without the STM — the flow
/// used to place the reference spacecraft along the arc.
pub fn planar_propagate(s0: &PlanarState, mu: f64, t: f64, steps: usize) -> PlanarState {
    let embed = crate::cr3bp::Cr3bpState {
        r: [s0[0], s0[1], 0.0],
        v: [s0[2], s0[3], 0.0],
    };
    let st = crate::cr3bp::propagate_cr3bp(embed, mu, t, steps);
    [st.r[0], st.r[1], st.v[0], st.v[1]]
}

/// Propagate a **spatial** CR3BP state and its 6×6 variational STM for time `t`.
///
/// This is the crate's finite-difference-validated CR3BP STM
/// ([`crate::cr3bp::propagate_state_stm`]) used at full width — no sub-block extraction,
/// no re-derivation. Returns `(state(t), Φ(t))` with `Φ` the true linearisation of the
/// three-dimensional flow. [`planar_state_stm`] is the `{x, y, ẋ, ẏ}` restriction of this
/// same matrix, so the two paths share one validated linearisation.
pub fn spatial_state_stm(
    s0: &SpatialState,
    mu: f64,
    t: f64,
    steps: usize,
) -> (SpatialState, [[f64; N_SPATIAL]; N_SPATIAL]) {
    let embed = crate::cr3bp::Cr3bpState {
        r: [s0[0], s0[1], s0[2]],
        v: [s0[3], s0[4], s0[5]],
    };
    let (st, phi) = crate::cr3bp::propagate_state_stm(&embed, mu, t, steps);
    ([st.r[0], st.r[1], st.r[2], st.v[0], st.v[1], st.v[2]], phi)
}

/// Spatial CR3BP state after time `t` (position + velocity), without the STM — the flow
/// used to place the reference spacecraft along a three-dimensional arc.
pub fn spatial_propagate(s0: &SpatialState, mu: f64, t: f64, steps: usize) -> SpatialState {
    let embed = crate::cr3bp::Cr3bpState {
        r: [s0[0], s0[1], s0[2]],
        v: [s0[3], s0[4], s0[5]],
    };
    let st = crate::cr3bp::propagate_cr3bp(embed, mu, t, steps);
    [st.r[0], st.r[1], st.r[2], st.v[0], st.v[1], st.v[2]]
}

// ── Observability assembly (L27) ─────────────────────────────────────────────

/// One tracking epoch: the measurement Jacobian rows `H_k` (each row a length-`n`
/// partial), the variational STM `Φ_k` mapping the initial state to this epoch, and the
/// integration weight `Δt_k` folded into the Gramian.
#[derive(Clone, Debug)]
pub struct ObsEpoch {
    /// Measurement Jacobian rows at this epoch (`m_k × n`).
    pub h: Mat,
    /// Variational STM `Φ(t_k)` from the initial epoch (`n × n`).
    pub phi: Mat,
    /// Integration weight (the sub-arc length this epoch represents).
    pub dt: f64,
}

/// A single stacked observability row `h · Φ` (row vector times matrix).
#[allow(clippy::needless_range_loop)]
fn row_times_matrix(h: &[f64], phi: &Mat) -> Vec<f64> {
    let n = phi.len();
    let mut out = vec![0.0; n];
    for c in 0..h.len() {
        let hc = h[c];
        if hc == 0.0 {
            continue;
        }
        for j in 0..n {
            out[j] += hc * phi[c][j];
        }
    }
    out
}

/// Assemble the stacked **observability matrix** `O = stack_k[ H_k · Φ_k ]` and the
/// per-row integration weights (each row inherits its epoch's `Δt`). Returns
/// `(O, weights)`.
pub fn observability_matrix(epochs: &[ObsEpoch]) -> (Mat, Vec<f64>) {
    let mut o = Vec::new();
    let mut w = Vec::new();
    for ep in epochs {
        for h in &ep.h {
            o.push(row_times_matrix(h, &ep.phi));
            w.push(ep.dt);
        }
    }
    (o, w)
}

/// The dt-weighted observability **Gramian** `W = Σ_k Δt_k · Φ_kᵀ H_kᵀ H_k Φ_k`.
///
/// This is exactly the weighted Gram matrix `Σ_row w_row · o_rowᵀ o_row` of the stacked
/// observability rows, so it is assembled through the crate's Fisher-information kernel
/// [`crate::fim::information_matrix`].
pub fn gramian(epochs: &[ObsEpoch]) -> Mat {
    let (o, w) = observability_matrix(epochs);
    information_matrix(&o, &w)
}

/// Singular values of `O` in **descending** order, via the rank-revealing identity
/// `σ_i = √λ_i(OᵀO)` (the crate's symmetric Jacobi eigensolver on the Gram matrix). This
/// is the SVD spectrum the observability rank is thresholded from.
pub fn singular_values(o: &Mat) -> Vec<f64> {
    if o.is_empty() {
        return vec![];
    }
    // Unweighted Gram OᵀO; its eigenvalues are the squared singular values of O.
    let ones = vec![1.0; o.len()];
    let gram = information_matrix(o, &ones);
    let e = sym_eig(&gram);
    let mut sv: Vec<f64> = e.values.iter().map(|&l| l.max(0.0).sqrt()).collect();
    sv.sort_by(|a, b| b.total_cmp(a));
    sv
}

/// Count of singular values that clear the observability rank threshold, given the sorted
/// (descending) singular-value spectrum of `O`.
///
/// **The one, internally-consistent rank definition (a singular-value read on `O`).** A singular
/// value counts as an observed direction iff `σ_i > rel_tol · σ_max`. This is the SVD read of
/// observability that `numpy.linalg.matrix_rank(O)` implements, and it is applied **identically**
/// wherever P6 reports an observable rank — the rank-vs-arc table ([`rank_vs_arc`]), the full-arc
/// [`observable_rank`], the SRIF cross-validation ([`crate::cislunar_srif`]) **and** the
/// Gramian-spectrum rank ([`gramian_spectrum`]) — so the committed table cannot be
/// self-contradictory (one column reading a `σ`-threshold while another reads a mismatched
/// `λ`-threshold).
///
/// **Consistency with the eigenvalue side.** The crate forms these singular values as
/// `σ_i = √λ_i(OᵀO)` (eigenvalues of the Gram matrix, via [`crate::fim::sym_eig`]). Because the
/// Gramian `W = OᵀO` squares the data, the equivalent eigenvalue floor for `σ > rel_tol · σ_max`
/// is `λ = σ² > rel_tol² · λ_max`. [`gramian_spectrum`] therefore reads its rank through **this
/// same function** (on `√λ(W)`), rather than through the generic `rel_tol · λ_max` eigenvalue floor
/// of [`crate::fim::design_metrics`] (which is calibrated for un-squared information matrices):
/// the two would otherwise disagree on a weakly-observed direction whose `λ`-ratio sits between
/// `rel_tol²` and `rel_tol`. Reporting one `rel_tol` on the `σ` scale everywhere removes that
/// ambiguity.
pub fn rank_from_singular_values(sv_descending: &[f64], rel_tol: f64) -> usize {
    let smax = sv_descending.first().copied().unwrap_or(0.0);
    if smax <= 0.0 {
        return 0;
    }
    let thr = rel_tol * smax;
    sv_descending.iter().filter(|&&s| s > thr).count()
}

/// The relative singular-value tolerance below which a rank read stops being a rank.
///
/// The `σ > rel_tol·σ_max` count is taken on singular values reconstructed as `σ = √λ(OᵀO)`,
/// so on the eigenvalue side the floor is `rel_tol²·λ_max`. The Gram matrix itself is only
/// accurate to about `f64::EPSILON·λ_max`, so once `rel_tol² < f64::EPSILON` — that is, once
/// `rel_tol < √f64::EPSILON ≈ 1.49e-8` — the floor sits *below* the rounding noise of the
/// matrix being decomposed and the count starts including directions that are pure
/// arithmetic residue. [`bounded_rank_from_singular_values`] says so in its stated reason
/// rather than silently returning the inflated count.
///
/// The constant is `f64::EPSILON.sqrt()` written out (`sqrt` is not a `const fn`); the
/// module's unit tests assert the two are bit-identical.
pub const RANK_TOLERANCE_NOISE_FLOOR: f64 = 1.490_116_119_384_765_6e-8;

/// A numerical rank read together with the algebraic bound that limited it.
///
/// A matrix of `m` rows and `n` columns cannot have rank above `min(m, n)` — Sylvester's
/// bound, an identity, not a tolerance. A singular-value count at a tight `rel_tol` can
/// exceed it anyway (the reconstructed spectrum carries f64 residue below
/// [`RANK_TOLERANCE_NOISE_FLOOR`]), and then the number reported is not a rank at all. This
/// carries both the reported rank and what it would have been, so the clamp is visible.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct BoundedRank {
    /// The reported rank: the singular-value count, never above [`Self::shape_bound`].
    pub rank: usize,
    /// The raw `σ > rel_tol·σ_max` count, before the algebraic bound was applied.
    pub counted: usize,
    /// `min(rows, cols)` — the largest rank the matrix shape admits.
    pub shape_bound: usize,
    /// Why the reported rank is below the raw count. `Some` exactly when the count was
    /// clamped; never a silent clamp.
    pub reason: Option<String>,
}

/// The rank of an `n_rows × n_cols` matrix from its descending singular-value spectrum,
/// with the algebraic bound `rank ≤ min(rows, cols)` applied and the clamp stated.
///
/// The threshold convention is unchanged — `σ_i > rel_tol · σ_max`, the one P6 convention of
/// [`rank_from_singular_values`], which this calls. What is added is the shape bound: the
/// count is an *estimate* of the rank and can exceed what the matrix dimensions allow once
/// `rel_tol` drops below [`RANK_TOLERANCE_NOISE_FLOOR`], where the equivalent eigenvalue
/// floor `rel_tol²·λ_max` sinks under the f64 noise of the Gram matrix. Above that floor the
/// clamp never fires and the read is bit-identical to the unbounded one.
pub fn bounded_rank_from_singular_values(
    sv_descending: &[f64],
    rel_tol: f64,
    n_rows: usize,
    n_cols: usize,
) -> BoundedRank {
    let counted = rank_from_singular_values(sv_descending, rel_tol);
    let shape_bound = n_rows.min(n_cols);
    if counted <= shape_bound {
        return BoundedRank {
            rank: counted,
            counted,
            shape_bound,
            reason: None,
        };
    }
    let floor_note = if rel_tol < RANK_TOLERANCE_NOISE_FLOOR {
        format!(
            " The tolerance is below the f64 rank-read noise floor {RANK_TOLERANCE_NOISE_FLOOR:.3e} \
             (= sqrt(f64::EPSILON)): its equivalent eigenvalue floor rel_tol^2 = {:.3e} sits under \
             the rounding noise of the Gram matrix, so the extra directions are arithmetic \
             residue, not observability.",
            rel_tol * rel_tol
        )
    } else {
        String::new()
    };
    BoundedRank {
        rank: shape_bound,
        counted,
        shape_bound,
        reason: Some(format!(
            "rank clamped from the singular-value count {counted} to min(rows, cols) = \
             min({n_rows}, {n_cols}) = {shape_bound} at rel_tol {rel_tol:.3e}: a {n_rows}x{n_cols} \
             matrix cannot have rank {counted}.{floor_note}"
        )),
    }
}

/// A stated reason when `rel_tol` sits below [`RANK_TOLERANCE_NOISE_FLOOR`], or `None` when
/// it does not — the document-level companion to [`bounded_rank_from_singular_values`].
///
/// The tolerance is *not* altered: flooring it would move every rank read taken below the
/// floor, which is a change to a published convention rather than a repair. The read is left
/// exactly as asked for and the caller is told what it is reading.
pub fn rank_tolerance_note(rel_tol: f64) -> Option<String> {
    if rel_tol >= RANK_TOLERANCE_NOISE_FLOOR || rel_tol <= 0.0 {
        return None;
    }
    Some(format!(
        "rel_tol {rel_tol:.3e} is below the f64 rank-read noise floor \
         {RANK_TOLERANCE_NOISE_FLOOR:.3e} (= sqrt(f64::EPSILON)). Ranks are read as \
         sigma > rel_tol*sigma_max on sigma = sqrt(lambda(O^T O)), so the equivalent eigenvalue \
         floor is rel_tol^2 = {:.3e}; below sqrt(f64::EPSILON) that floor sits under the rounding \
         noise of the Gram matrix and the count can include directions that are arithmetic \
         residue. The tolerance asked for is used unchanged; every read is additionally held to \
         rank <= min(rows, cols), and any read that hits that bound says so.",
        rel_tol * rel_tol
    ))
}

/// Numerical **rank** of `O` from the singular-value threshold — the rank-revealing SVD read of
/// observability. See [`rank_from_singular_values`] for the one threshold convention (the same
/// `rel_tol · σ_max` floor used by every P6 rank read, matching `numpy.linalg.matrix_rank(O)`),
/// and [`bounded_rank_from_singular_values`] for the algebraic `rank ≤ min(rows, cols)` bound
/// this applies on top of it (a no-op at every tolerance above
/// [`RANK_TOLERANCE_NOISE_FLOOR`]).
pub fn observable_rank(o: &Mat, rel_tol: f64) -> usize {
    bounded_observable_rank(o, rel_tol).rank
}

/// [`observable_rank`] with the clamp made visible: the reported rank, the raw
/// singular-value count, and a stated reason whenever the two differ.
pub fn bounded_observable_rank(o: &Mat, rel_tol: f64) -> BoundedRank {
    let sv = singular_values(o);
    let n_cols = o.first().map(|r| r.len()).unwrap_or(0);
    bounded_rank_from_singular_values(&sv, rel_tol, o.len(), n_cols)
}

/// The symmetric spectrum and conditioning of an observability Gramian `W`.
#[derive(Clone, Debug)]
pub struct GramianSpectrum {
    /// Eigenvalues of `W` in ascending order (the symmetric spectrum).
    pub eigenvalues: Vec<f64>,
    /// Smallest eigenvalue `λ_min` — the worst-observed direction.
    pub min_eigenvalue: f64,
    /// Largest eigenvalue `λ_max`.
    pub max_eigenvalue: f64,
    /// `trace(W) = Σλ` — total information (an eigen-invariant cross-check anchor).
    pub trace: f64,
    /// Condition number `λ_max / λ_min` over the observable subspace (`+∞` if singular).
    pub condition: f64,
    /// Numerical rank (observable directions) of `W`.
    pub rank: usize,
    /// Datum-defect dimension `n − rank` (unobservable directions).
    pub defect: usize,
}

/// Eigen-spectrum, `λ_min`, condition number and rank of a Gramian `W = OᵀO`, read off the
/// crate's symmetric eigensolver and experiment-design metrics.
///
/// **Rank, defect AND condition all use the one observability rank convention**
/// ([`rank_from_singular_values`]): the eigenvalues are mapped back to singular values `σ = √λ` and
/// thresholded at `σ > rel_tol · σ_max`, so this Gramian rank agrees exactly with the SVD rank of
/// `O` ([`observable_rank`] / [`rank_vs_arc`]) and with `numpy.linalg.matrix_rank(O)`. Reading the
/// rank off the raw `λ > rel_tol · λ_max` floor of [`crate::fim::design_metrics`] instead would
/// apply a *different* effective floor (`rel_tol²` on the `σ` scale) and could report a different
/// rank on a weakly-observed direction — the self-contradiction this reconciliation removes. The
/// condition number is `λ_max / λ_min` over the **same** observable subspace this rank defines (the
/// top-`rank` eigenvalues), so a full-rank Gramian's condition is the true `λ_max / λ_min`, matching
/// `numpy.linalg.cond(W)`; a rank-deficient one's condition is taken over the observed modes only.
pub fn gramian_spectrum(w: &Mat, rel_tol: f64) -> GramianSpectrum {
    let e = sym_eig(w);
    let min_eigenvalue = e.values.first().copied().unwrap_or(0.0);
    let max_eigenvalue = e.values.last().copied().unwrap_or(0.0);
    let trace = e.values.iter().sum();
    // The one rank convention: σ = √λ thresholded at rel_tol·σ_max (= λ > rel_tol²·λ_max), so the
    // Gramian rank matches the SVD rank of O and numpy.linalg.matrix_rank(O).
    let mut sv: Vec<f64> = e.values.iter().map(|&l| l.max(0.0).sqrt()).collect();
    sv.sort_by(|a, b| b.total_cmp(a));
    let rank = rank_from_singular_values(&sv, rel_tol);
    let n = w.len();
    let defect = n.saturating_sub(rank);
    // Condition over the SAME observable subspace the rank defines: λ_max / (smallest of the
    // top-`rank` eigenvalues). `e.values` is ascending, so the observed modes are the last `rank`,
    // and the smallest observed eigenvalue is `e.values[n - rank]`. This keeps the condition
    // consistent with the reported rank (a full-rank W's condition is the true λ_max/λ_min, so it
    // matches numpy.linalg.cond(W)); an empty/zero subspace is `+∞`.
    let condition = if rank == 0 {
        f64::INFINITY
    } else {
        let lmin_obs = e.values[n - rank];
        if lmin_obs > 0.0 {
            max_eigenvalue / lmin_obs
        } else {
            f64::INFINITY
        }
    };
    GramianSpectrum {
        eigenvalues: e.values,
        min_eigenvalue,
        max_eigenvalue,
        trace,
        condition,
        rank,
        defect,
    }
}

/// One point of the rank-vs-arc-length table: the observable rank over the tracking arc
/// truncated at epoch `epoch_index` (arc time `arc_time`).
#[derive(Clone, Debug)]
pub struct RankArcPoint {
    /// Index of the last epoch in this prefix.
    pub epoch_index: usize,
    /// Elapsed arc time from the first epoch (normalised rotating-frame time units).
    pub arc_time: f64,
    /// Total stacked measurement rows accumulated up to and including this epoch.
    pub n_rows: usize,
    /// Numerical observable rank of the arc so far (SVD threshold).
    pub rank: usize,
    /// Largest singular value of the stacked `O` so far.
    pub sigma_max: f64,
    /// Smallest singular value of the stacked `O` so far.
    pub sigma_min: f64,
    /// Why [`Self::rank`] is below the raw singular-value count at this prefix — `Some`
    /// exactly when the algebraic `rank ≤ min(rows, cols)` bound clamped the read (see
    /// [`bounded_rank_from_singular_values`]), so the clamp is never silent.
    pub rank_limited_by: Option<String>,
}

/// The **rank-vs-arc-length** table: for each growing prefix of the epoch sequence, the
/// numerical observable rank of the accumulated observability matrix. As the arc extends,
/// the rank grows toward full observability (paper P6 Table 1).
pub fn rank_vs_arc(epochs: &[ObsEpoch], rel_tol: f64) -> Vec<RankArcPoint> {
    let mut o: Mat = Vec::new();
    let mut arc = 0.0;
    let mut out = Vec::with_capacity(epochs.len());
    for (k, ep) in epochs.iter().enumerate() {
        arc += ep.dt;
        for h in &ep.h {
            o.push(row_times_matrix(h, &ep.phi));
        }
        let sv = singular_values(&o);
        let sigma_max = sv.first().copied().unwrap_or(0.0);
        let sigma_min = sv.last().copied().unwrap_or(0.0);
        // The one, eigenvalue-consistent rank definition (see `rank_from_singular_values`): the
        // σ-floor is rel_tol·σ_max, i.e. the eigenvalue floor rel_tol²·λ_max the Gramian read
        // uses. (This comment previously named √rel_tol·σ_max and rel_tol·λ_max, contradicting
        // the convention every other doc-comment in the module states; the code always did the
        // above.)
        // …and the read is additionally held to the algebraic bound rank <= min(rows, cols):
        // a prefix of two measurement rows cannot observe three directions however tight the
        // tolerance. The clamp states itself in `rank_limited_by` rather than happening quietly.
        let n_cols = o.first().map(|r| r.len()).unwrap_or(0);
        let read = bounded_rank_from_singular_values(&sv, rel_tol, o.len(), n_cols);
        out.push(RankArcPoint {
            epoch_index: k,
            arc_time: arc,
            n_rows: o.len(),
            rank: read.rank,
            sigma_max,
            sigma_min,
            rank_limited_by: read.reason,
        });
    }
    out
}

// ── Measurement noise: whitening and the formal posterior (R3) ───────────────
//
// A rank test knows nothing about measurement noise. The standard way noise enters a
// least-squares observability analysis is **whitening**: with a measurement covariance
// `R = σ²I` the information is `Σ Hᵀ R⁻¹ H = Σ (H/σ)ᵀ(H/σ)`, i.e. every Jacobian row is
// divided by its measurement sigma before the Gram matrix is formed. Two consequences,
// both reported rather than glossed:
//
// 1. For a **homoscedastic** σ the whitened observability matrix is `Õ = O/σ`, a *scalar
//    multiple* of `O`. Its singular-value spectrum is uniformly scaled, so its RELATIVE
//    spectrum — and therefore its numerical rank, datum defect and condition number — is
//    bit-for-bit invariant. A rank-based arc-length threshold cannot move with noise; any
//    reported movement would be an artefact, not physics.
// 2. What noise does move is the **formal posterior uncertainty** `P = (ÕᵀÕ)⁻¹ = σ²(OᵀO)⁻¹`,
//    which scales as `σ²` in variance (`σ` in standard deviation). That is the quantity an
//    estimator designer actually has to clear, so it is the criterion under which an
//    arc-length threshold is noise-dependent at all.

/// Divide every measurement Jacobian row of an epoch sequence by the measurement sigma —
/// the **whitening** step that puts measurement noise into the Gramian.
///
/// With a measurement covariance `R = σ²I` the Fisher information of the batch is
/// `Σ Hᵀ R⁻¹ H`, which is exactly the unit-weight Gram matrix of the rows `H/σ`. A
/// non-finite or non-positive `sigma` is treated as the **noise-free** case (unit weight),
/// so `sigma = 0` reproduces the un-whitened epochs exactly.
///
/// Because a homoscedastic `sigma` scales every row identically, this changes the SCALE of
/// the observability spectrum but not its shape: rank, defect and condition number are
/// invariant (asserted in this module's unit tests). The posterior covariance it implies
/// is not.
pub fn whiten_epochs(epochs: &[ObsEpoch], sigma: f64) -> Vec<ObsEpoch> {
    let scale = if sigma.is_finite() && sigma > 0.0 {
        1.0 / sigma
    } else {
        1.0
    };
    epochs
        .iter()
        .map(|ep| ObsEpoch {
            h: ep
                .h
                .iter()
                .map(|row| row.iter().map(|v| v * scale).collect())
                .collect(),
            phi: ep.phi.clone(),
            dt: ep.dt,
        })
        .collect()
}

/// The formal posterior uncertainty of a batch least-squares estimate of the initial
/// state, read from a **noise-whitened** observability matrix `Õ` (see [`whiten_epochs`]).
///
/// `P = (ÕᵀÕ)⁻¹` is the covariance of the batch estimator under the whitened measurement
/// model; the fields below are its standard-deviation summaries in the state's own
/// (normalised, nondimensional) units. A rank-deficient geometry has no finite `P`, so the
/// position/velocity summaries are `None` rather than a fabricated number — the same
/// honesty rule [`cislunar_gdop`] applies to a singular geometry.
#[derive(Clone, Debug)]
pub struct WhitenedPosterior {
    /// Numerical rank of `Õ` (the one P6 convention, `σ > rel_tol·σ_max`).
    pub rank: usize,
    /// Datum-defect dimension `n − rank` (unobservable directions).
    pub defect: usize,
    /// Condition number `λ_max/λ_min` of `ÕᵀÕ` over the observable subspace (`+∞` if empty).
    pub condition: f64,
    /// Per-state 1σ standard deviations from the diagonal of `P` (or of the Moore–Penrose
    /// pseudo-inverse when rank-deficient, where the null directions are meaningless).
    pub sigma_state: Vec<f64>,
    /// Root-sum-square of the position-block 1σ values (nondimensional length units), or
    /// `None` when the geometry is rank-deficient and no finite covariance exists.
    pub sigma_position: Option<f64>,
    /// Root-sum-square of the velocity-block 1σ values (nondimensional velocity units), or
    /// `None` when the geometry is rank-deficient.
    pub sigma_velocity: Option<f64>,
    /// Orthonormal basis of the unobservable directions as columns (`n × defect`).
    pub null_space: Mat,
    /// Why [`Self::rank`] is below the raw eigenvalue count — `Some` exactly when the
    /// algebraic `rank ≤ min(rows, cols)` bound clamped the read (see
    /// [`bounded_rank_from_singular_values`]). The posterior is then recomputed over the
    /// bounded subspace, so `defect`, `null_space` and the σ summaries all agree with the
    /// reported rank rather than with the inflated count.
    pub rank_limited_by: Option<String>,
}

/// Formal posterior uncertainty of the initial state from a noise-whitened observability
/// matrix `o`, with the first `n_pos` state components taken as the position block.
///
/// The rank threshold is the one P6 convention (`σ > rel_tol·σ_max`), transported to the
/// eigenvalue side of the Gram matrix as `λ > rel_tol²·λ_max`, so this rank agrees exactly
/// with [`observable_rank`], [`rank_vs_arc`] and [`gramian_spectrum`] on the same matrix.
pub fn whitened_posterior(o: &Mat, n_pos: usize, rel_tol: f64) -> WhitenedPosterior {
    if o.is_empty() || o[0].is_empty() {
        return WhitenedPosterior {
            rank: 0,
            defect: 0,
            condition: f64::INFINITY,
            sigma_state: vec![],
            sigma_position: None,
            sigma_velocity: None,
            null_space: vec![],
            rank_limited_by: None,
        };
    }
    let n = o[0].len();
    let ones = vec![1.0; o.len()];
    let gram = information_matrix(o, &ones);
    // λ-floor rel_tol² ⇔ the σ-floor rel_tol the rest of P6 reads its rank at.
    let c = crlb(&gram, rel_tol * rel_tol);
    // …then the algebraic bound: a batch of `o.len()` scalar measurement rows cannot determine
    // more than `min(rows, n)` directions of the state, however tight the tolerance. When the
    // eigenvalue count exceeds it the whole posterior — not just the rank — was read off
    // directions that are f64 residue, so the pseudo-inverse is rebuilt over the bounded
    // subspace and the clamp states itself.
    let shape_bound = o.len().min(n);
    let mut rank_limited_by: Option<String> = None;
    let c = if c.rank > shape_bound {
        // Reconstruct the σ-scale count so the stated reason speaks the one P6 convention.
        let mut sv: Vec<f64> = c.eigenvalues.iter().map(|&l| l.max(0.0).sqrt()).collect();
        sv.sort_by(|a, b| b.total_cmp(a));
        rank_limited_by = bounded_rank_from_singular_values(&sv, rel_tol, o.len(), n).reason;
        // A relative λ-floor that keeps exactly `shape_bound` directions: the geometric
        // midpoint between the smallest kept eigenvalue and the largest dropped one. Nothing
        // is tuned — the bound is the matrix shape, and this is only the threshold that
        // realises it.
        let lmax = c.eigenvalues.last().copied().unwrap_or(0.0);
        // `shape_bound < n` holds here (the branch needs `c.rank > shape_bound` and
        // `c.rank <= n`), so both indices are in range.
        let smallest_kept = c.eigenvalues[n - shape_bound];
        let largest_dropped = c.eigenvalues[n - shape_bound - 1];
        let rel = if lmax > 0.0 {
            if largest_dropped > 0.0 {
                (smallest_kept * largest_dropped).sqrt() / lmax
            } else {
                0.5 * smallest_kept / lmax
            }
        } else {
            rel_tol * rel_tol
        };
        crlb(&gram, rel)
    } else {
        c
    };
    let lmax = c.eigenvalues.last().copied().unwrap_or(0.0);
    let condition = if c.rank == 0 {
        f64::INFINITY
    } else {
        let lmin_obs = c.eigenvalues[n - c.rank];
        if lmin_obs > 0.0 {
            lmax / lmin_obs
        } else {
            f64::INFINITY
        }
    };
    let full = c.defect == 0;
    let rss = |lo: usize, hi: usize| -> Option<f64> {
        if !full {
            return None;
        }
        let s: f64 = c.crlb_diag[lo..hi].iter().map(|v| v.max(0.0)).sum();
        Some(s.sqrt())
    };
    let n_pos = n_pos.min(n);
    WhitenedPosterior {
        rank: c.rank,
        defect: c.defect,
        condition,
        sigma_state: c.crlb_std.clone(),
        sigma_position: rss(0, n_pos),
        sigma_velocity: rss(n_pos, n),
        null_space: c.null_space,
        rank_limited_by,
    }
}

// ── Range-rate design lever (L30) ────────────────────────────────────────────

/// The instantaneous (single-epoch) observability lever of adding Doppler.
#[derive(Clone, Debug)]
pub struct RankLever {
    /// Number of inter-satellite links in the snapshot.
    pub n_links: usize,
    /// Rank of a range-only measurement stack at one epoch (velocity columns are zero, so
    /// this can never exceed the position dimension).
    pub rank_range_only: usize,
    /// Rank of a range **and** range-rate stack at one epoch — Doppler's non-zero velocity
    /// columns lift the rank toward the full four-state.
    pub rank_range_rate: usize,
}

/// Compare the **instantaneous** rank of range-only vs range+range-rate measurements from
/// a chief spacecraft to a set of reference spacecraft at a single epoch. Range-only rows
/// have zero velocity columns (rank capped at the position dimension); the range-rate rows
/// add non-zero velocity columns, so the combined stack observes more of the state.
pub fn range_vs_range_rate_rank(
    chief: &PlanarState,
    refs: &[PlanarState],
    rel_tol: f64,
) -> RankLever {
    let mut h_range: Mat = Vec::new();
    let mut h_both: Mat = Vec::new();
    for r in refs {
        let (_rho, rr) = range_row(chief, r);
        h_range.push(rr.to_vec());
        h_both.push(rr.to_vec());
        let (_rd, rrr) = range_rate_row(chief, r);
        h_both.push(rrr.to_vec());
    }
    RankLever {
        n_links: refs.len(),
        rank_range_only: observable_rank(&h_range, rel_tol),
        rank_range_rate: observable_rank(&h_both, rel_tol),
    }
}

// ── GDOP-singular reporting (L33) ────────────────────────────────────────────

/// A geometric-dilution report for a cislunar snapshot: either a finite value or an
/// explicit *undefined* verdict for a rank-deficient (singular) geometry.
#[derive(Clone, Debug, PartialEq)]
pub enum CislunarGdop {
    /// A well-posed geometry: the geometric dilution of precision `√trace(M⁻¹)`.
    Defined {
        /// The geometric dilution of precision.
        gdop: f64,
        /// Numerical rank of the information matrix (full rank ⇒ observable).
        rank: usize,
    },
    /// A rank-deficient / singular geometry: GDOP is **undefined** (not a bogus finite
    /// number). Carries the numerical rank, datum defect and a human-readable reason.
    Undefined {
        /// Numerical rank of the information matrix.
        rank: usize,
        /// Datum-defect dimension `n − rank` (unobservable directions).
        defect: usize,
        /// Why the value is undefined.
        reason: String,
    },
}

/// Report GDOP for a cislunar geometry, or flag it **undefined** when the geometry is
/// rank-deficient — the honest analogue of [`crate::pvt::solve_spp`]'s singular guard
/// (its `invert4(GᵀG)` returns `None` for a singular geometry, yielding no fix). Here the
/// same singular geometry is reported explicitly through [`crate::fim::design_metrics`]:
/// a non-zero datum defect or an infinite condition number means no finite dilution
/// exists, so a value is never fabricated.
pub fn cislunar_gdop(rows: &Mat, rel_tol: f64) -> CislunarGdop {
    if rows.is_empty() {
        return CislunarGdop::Undefined {
            rank: 0,
            defect: 0,
            reason: "no measurement rows: geometry is empty".to_string(),
        };
    }
    let n = rows[0].len();
    let weights = vec![1.0; rows.len()];
    let m = information_matrix(rows, &weights);
    let dm = design_metrics(&m, rel_tol);
    if dm.defect > 0 || !dm.condition.is_finite() {
        return CislunarGdop::Undefined {
            rank: dm.rank,
            defect: n - dm.rank,
            reason: format!(
                "GDOP undefined (rank-deficient / singular geometry): rank {} of {} \
                 states, datum defect {}, condition {}",
                dm.rank,
                n,
                n - dm.rank,
                if dm.condition.is_finite() {
                    format!("{:.3e}", dm.condition)
                } else {
                    "inf".to_string()
                }
            ),
        };
    }
    // Full rank: GDOP = √trace(M⁻¹), the sum of the per-state variance lower bounds.
    let c = crate::fim::crlb(&m, rel_tol);
    let trace_inv: f64 = c.crlb_diag.iter().sum();
    CislunarGdop::Defined {
        gdop: trace_inv.max(0.0).sqrt(),
        rank: dm.rank,
    }
}

// ── Oracle helper: an independent determinant ────────────────────────────────

/// Determinant of a square matrix by Gaussian elimination with partial pivoting — an
/// **independent** route to `det(W)` (a different algorithm from the eigen-product
/// `Πλ`), used to cross-check the eigensolver's spectrum on a full-rank block.
// Dense Gaussian-elimination kernel: explicit (col, r, c) index arithmetic is the natural
// form here (and matches the crate's other matrix code) — iterator rewrites obscure it.
#[allow(clippy::needless_range_loop)]
pub fn determinant(m: &Mat) -> f64 {
    let n = m.len();
    if n == 0 {
        return 1.0;
    }
    let mut a: Vec<Vec<f64>> = m.to_vec();
    let mut det = 1.0;
    for col in 0..n {
        // Partial pivot: largest magnitude in this column at or below the diagonal.
        let mut piv = col;
        let mut best = a[col][col].abs();
        for r in (col + 1)..n {
            let v = a[r][col].abs();
            if v > best {
                best = v;
                piv = r;
            }
        }
        if best == 0.0 {
            return 0.0;
        }
        if piv != col {
            a.swap(piv, col);
            det = -det;
        }
        det *= a[col][col];
        let pivot = a[col][col];
        for r in (col + 1)..n {
            let factor = a[r][col] / pivot;
            if factor != 0.0 {
                for c in col..n {
                    a[r][c] -= factor * a[col][c];
                }
            }
        }
    }
    det
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::cr3bp::EARTH_MOON_MU;

    fn frob_sq(m: &Mat) -> f64 {
        m.iter().flat_map(|r| r.iter()).map(|v| v * v).sum()
    }

    // ── L31: the planar variational STM matches a central finite difference ──────
    /// ORACLE (Validated): the 4×4 planar STM equals a central finite-difference STM of
    /// the CR3BP flow (a different code path — plain state propagation) to tolerance.
    #[test]
    fn planar_stm_matches_finite_difference() {
        let s0: PlanarState = [1.08, 0.03, 0.10, -0.50];
        let (t, steps) = (0.20, 4000);
        let (_st, phi) = planar_state_stm(&s0, EARTH_MOON_MU, t, steps);
        let eps = 1e-6;
        for j in 0..N_PLANAR {
            let mut sp = s0;
            let mut sm = s0;
            sp[j] += eps;
            sm[j] -= eps;
            let ep = planar_propagate(&sp, EARTH_MOON_MU, t, steps);
            let em = planar_propagate(&sm, EARTH_MOON_MU, t, steps);
            for i in 0..N_PLANAR {
                let fd = (ep[i] - em[i]) / (2.0 * eps);
                assert!(
                    (phi[i][j] - fd).abs() < 1e-5,
                    "planar STM[{i}][{j}] = {} vs finite-diff {fd}",
                    phi[i][j]
                );
            }
        }
    }

    #[test]
    #[allow(clippy::needless_range_loop)]
    fn planar_stm_is_identity_at_zero_time() {
        let s0: PlanarState = [1.1, 0.0, 0.0, -0.5];
        let (_st, phi) = planar_state_stm(&s0, EARTH_MOON_MU, 0.0, 10);
        for i in 0..N_PLANAR {
            for j in 0..N_PLANAR {
                let want = if i == j { 1.0 } else { 0.0 };
                assert!((phi[i][j] - want).abs() < 1e-12);
            }
        }
    }

    // ── L27: eigen-spectrum cross-checks against spectral invariants ─────────────
    /// ORACLE (Validated): the eigenvalues of a symmetric-positive-definite Gramian obey
    /// `trace = Σλ`, `‖W‖_F² = Σλ²`, and `det = Πλ` (an independent Gaussian-elimination
    /// determinant) — three invariants pinning the returned spectrum to the true one.
    #[test]
    fn gramian_spectrum_satisfies_spectral_invariants() {
        // Build a genuine arc Gramian from two epochs so W is full rank.
        let epochs = sample_arc();
        let w = gramian(&epochs);
        let spec = gramian_spectrum(&w, 1e-9);
        let sum: f64 = spec.eigenvalues.iter().sum();
        let sum_sq: f64 = spec.eigenvalues.iter().map(|l| l * l).sum();
        let prod: f64 = spec.eigenvalues.iter().product();
        let tr: f64 = (0..w.len()).map(|i| w[i][i]).sum();
        assert!(
            (sum - tr).abs() <= 1e-9 * (1.0 + tr.abs()),
            "trace {tr} vs Σλ {sum}"
        );
        assert!(
            (sum_sq - frob_sq(&w)).abs() <= 1e-9 * (1.0 + frob_sq(&w)),
            "Frobenius² {} vs Σλ² {sum_sq}",
            frob_sq(&w)
        );
        let det = determinant(&w);
        assert!(
            (prod - det).abs() <= 1e-8 * (1.0 + det.abs()),
            "det {det} vs Πλ {prod}"
        );
        // A symmetric PSD Gramian: every eigenvalue is non-negative.
        assert!(spec.min_eigenvalue >= -1e-12);
    }

    /// The SVD rank of O agrees with the eigen-rank of the dt-weighted Gramian W — two
    /// independent rank-revealing routes on the same observable subspace.
    #[test]
    fn svd_rank_matches_gramian_eigen_rank() {
        let epochs = sample_arc();
        let (o, _w) = observability_matrix(&epochs);
        let svd_rank = observable_rank(&o, 1e-9);
        let w = gramian(&epochs);
        let spec = gramian_spectrum(&w, 1e-9);
        assert_eq!(svd_rank, spec.rank, "SVD rank vs Gramian eigen-rank");
    }

    // A short two-epoch single-link arc that is full-rank observable (used by the
    // invariant + rank cross-check tests).
    fn sample_arc() -> Vec<ObsEpoch> {
        let chief: PlanarState = [1.10, 0.02, 0.05, -0.50];
        let reference: PlanarState = [1.02, -0.03, -0.06, -0.55];
        let mu = EARTH_MOON_MU;
        let ts = [0.02_f64, 0.05_f64];
        let mut out = Vec::new();
        let mut prev = 0.0;
        for &t in &ts {
            let (cs, phi) = planar_state_stm(&chief, mu, t, 2000);
            let rs = planar_propagate(&reference, mu, t, 2000);
            let (_rho, r_row) = range_row(&cs, &rs);
            let (_rd, rr_row) = range_rate_row(&cs, &rs);
            out.push(ObsEpoch {
                h: vec![r_row.to_vec(), rr_row.to_vec()],
                phi: phi.iter().map(|r| r.to_vec()).collect(),
                dt: t - prev,
            });
            prev = t;
        }
        out
    }

    // ── L30: the range-rate lever raises instantaneous rank ─────────────────────
    #[test]
    fn range_rate_raises_instantaneous_rank() {
        let chief: PlanarState = [1.10, 0.02, 0.05, -0.50];
        let refs = [
            [1.02, -0.03, -0.06, -0.55],
            [1.15, 0.05, 0.10, -0.45],
            [1.05, 0.06, 0.18, -0.40],
        ];
        let lever = range_vs_range_rate_rank(&chief, &refs, 1e-9);
        assert!(
            lever.rank_range_rate > lever.rank_range_only,
            "range+rate rank {} must exceed range-only rank {}",
            lever.rank_range_rate,
            lever.rank_range_only
        );
        // Range-only can never exceed the planar position dimension (velocity blind).
        assert!(lever.rank_range_only <= 2);
    }

    #[test]
    fn range_only_single_link_is_rank_one() {
        let chief: PlanarState = [1.10, 0.02, 0.05, -0.50];
        let refs = [[1.02, -0.03, -0.06, -0.55]];
        let lever = range_vs_range_rate_rank(&chief, &refs, 1e-9);
        assert_eq!(lever.rank_range_only, 1, "one range snapshot is rank-1");
        assert!(lever.rank_range_rate >= 2, "range+rate sees velocity too");
    }

    // ── L33: rank-deficient geometry flags GDOP undefined ───────────────────────
    #[test]
    fn rank_deficient_geometry_flags_gdop_undefined() {
        // Range-only rows at a single epoch: velocity columns are all zero, so the
        // 4-state geometry is rank-deficient and GDOP must be undefined, not finite.
        let chief: PlanarState = [1.10, 0.02, 0.05, -0.50];
        let refs = [[1.02, -0.03, -0.06, -0.55], [1.15, 0.05, 0.10, -0.45]];
        let mut rows: Mat = Vec::new();
        for r in &refs {
            let (_rho, rr) = range_row(&chief, r);
            rows.push(rr.to_vec());
        }
        match cislunar_gdop(&rows, 1e-9) {
            CislunarGdop::Undefined { defect, .. } => assert!(defect >= 1),
            CislunarGdop::Defined { gdop, .. } => {
                panic!("rank-deficient geometry must not yield a finite GDOP {gdop}")
            }
        }
    }

    #[test]
    fn full_rank_geometry_yields_finite_gdop() {
        // Range + range-rate to three references spans the full four-state, so GDOP is a
        // finite, positive number.
        let chief: PlanarState = [1.10, 0.02, 0.05, -0.50];
        let refs = [
            [1.02, -0.03, -0.06, -0.55],
            [1.15, 0.05, 0.10, -0.45],
            [1.05, 0.06, 0.18, -0.40],
        ];
        let mut rows: Mat = Vec::new();
        for r in &refs {
            let (_rho, rr) = range_row(&chief, r);
            rows.push(rr.to_vec());
            let (_rd, rrr) = range_rate_row(&chief, r);
            rows.push(rrr.to_vec());
        }
        match cislunar_gdop(&rows, 1e-9) {
            CislunarGdop::Defined { gdop, rank } => {
                assert_eq!(rank, N_PLANAR);
                assert!(gdop.is_finite() && gdop > 0.0, "GDOP {gdop}");
            }
            CislunarGdop::Undefined { reason, .. } => panic!("expected finite GDOP: {reason}"),
        }
    }

    #[test]
    fn determinant_matches_known_values() {
        // Identity → 1; a 2×2 with a known determinant; a singular row → 0.
        let id: Mat = vec![vec![1.0, 0.0], vec![0.0, 1.0]];
        assert!((determinant(&id) - 1.0).abs() < 1e-12);
        let m: Mat = vec![vec![4.0, 3.0], vec![6.0, 3.0]];
        assert!((determinant(&m) - (4.0 * 3.0 - 3.0 * 6.0)).abs() < 1e-12);
        let sing: Mat = vec![vec![1.0, 2.0], vec![2.0, 4.0]];
        assert!(determinant(&sing).abs() < 1e-12);
    }

    // ── Spatial STM bridge ───────────────────────────────────────────────────

    /// ORACLE (Validated): the 6×6 spatial STM equals a central finite-difference STM of
    /// the CR3BP flow (a different code path — plain state propagation) to tolerance.
    #[test]
    fn spatial_stm_matches_finite_difference() {
        let s0: SpatialState = [1.05, 0.02, -0.10, 0.10, 0.20, -0.05];
        let (t, steps) = (0.20, 4000);
        let (_st, phi) = spatial_state_stm(&s0, EARTH_MOON_MU, t, steps);
        let eps = 1e-6;
        for j in 0..N_SPATIAL {
            let mut sp = s0;
            let mut sm = s0;
            sp[j] += eps;
            sm[j] -= eps;
            let ep = spatial_propagate(&sp, EARTH_MOON_MU, t, steps);
            let em = spatial_propagate(&sm, EARTH_MOON_MU, t, steps);
            for i in 0..N_SPATIAL {
                let fd = (ep[i] - em[i]) / (2.0 * eps);
                assert!(
                    (phi[i][j] - fd).abs() < 1e-5,
                    "spatial STM[{i}][{j}] = {} vs finite-diff {fd}",
                    phi[i][j]
                );
            }
        }
    }

    /// The planar STM is EXACTLY the `{x, y, ẋ, ẏ}` sub-block of the spatial STM in the
    /// `z = ż = 0` embedding — one validated linearisation, two views of it.
    #[test]
    fn planar_stm_is_the_spatial_stm_sub_block() {
        let s4: PlanarState = [1.08, 0.03, 0.10, -0.50];
        let s6: SpatialState = [s4[0], s4[1], 0.0, s4[2], s4[3], 0.0];
        let (t, steps) = (0.05, 2000);
        let (_a, phi4) = planar_state_stm(&s4, EARTH_MOON_MU, t, steps);
        let (_b, phi6) = spatial_state_stm(&s6, EARTH_MOON_MU, t, steps);
        let idx = [0usize, 1, 3, 4];
        for (i, &ri) in idx.iter().enumerate() {
            for (j, &cj) in idx.iter().enumerate() {
                assert_eq!(phi4[i][j], phi6[ri][cj], "sub-block [{i}][{j}]");
            }
        }
    }

    /// In the `z = 0` plane the out-of-plane block DECOUPLES exactly: the STM's
    /// in-plane↔out-of-plane cross terms are identically zero. This is the dynamical half
    /// of why a wholly planar constellation cannot observe a six-state (the measurement
    /// half is the zero `û_z` column of a coplanar range row).
    #[test]
    fn out_of_plane_block_decouples_at_z_zero() {
        let s6: SpatialState = [1.08, 0.03, 0.0, 0.10, -0.50, 0.0];
        let (_st, phi) = spatial_state_stm(&s6, EARTH_MOON_MU, 0.05, 2000);
        let inplane = [0usize, 1, 3, 4];
        let outplane = [2usize, 5];
        for &i in &inplane {
            for &j in &outplane {
                assert_eq!(
                    phi[i][j], 0.0,
                    "Φ[{i}][{j}] couples out-of-plane into plane"
                );
                assert_eq!(
                    phi[j][i], 0.0,
                    "Φ[{j}][{i}] couples plane into out-of-plane"
                );
            }
        }
    }

    // ── Measurement noise: whitening ─────────────────────────────────────────

    /// Homoscedastic whitening is a SCALAR multiple of the observability matrix, so the
    /// numerical rank, the datum defect and the condition number are invariant — the
    /// reason a rank-based arc-length threshold cannot move with measurement noise.
    #[test]
    fn homoscedastic_whitening_leaves_rank_and_condition_invariant() {
        let epochs = sample_arc();
        let (o0, _) = observability_matrix(&epochs);
        let spec0 = gramian_spectrum(&gramian(&epochs), 1e-9);
        for &sigma in &[1e-3_f64, 1e-6, 1e-9, 1e-12] {
            let w = whiten_epochs(&epochs, sigma);
            let (o1, _) = observability_matrix(&w);
            assert_eq!(
                observable_rank(&o0, 1e-9),
                observable_rank(&o1, 1e-9),
                "rank moved under whitening at σ = {sigma}"
            );
            let spec1 = gramian_spectrum(&gramian(&w), 1e-9);
            assert_eq!(spec0.rank, spec1.rank);
            assert_eq!(spec0.defect, spec1.defect);
            let ratio = spec1.condition / spec0.condition;
            // Invariant to eigensolver precision: the Jacobi rotations are not bit-exact
            // under a uniform rescale, but the condition number agrees to ~1e-8 relative.
            assert!(
                (ratio - 1.0).abs() < 1e-6,
                "condition moved under whitening at σ = {sigma}: ratio {ratio}"
            );
        }
        // σ = 0 is the noise-free case: the epochs come back untouched.
        let free = whiten_epochs(&epochs, 0.0);
        let (o_free, _) = observability_matrix(&free);
        assert_eq!(o_free, o0);
    }

    /// The posterior standard deviation scales EXACTLY linearly with the measurement
    /// sigma (`P = σ²(OᵀO)⁻¹`) — the quantity under which an arc-length threshold is
    /// genuinely noise-dependent.
    #[test]
    fn posterior_sigma_scales_linearly_with_measurement_sigma() {
        let epochs = sample_arc();
        let (o1, _) = observability_matrix(&whiten_epochs(&epochs, 1e-9));
        let (o2, _) = observability_matrix(&whiten_epochs(&epochs, 2e-9));
        let p1 = whitened_posterior(&o1, 2, 1e-9);
        let p2 = whitened_posterior(&o2, 2, 1e-9);
        assert_eq!(p1.rank, p2.rank);
        let (s1, s2) = (p1.sigma_position.unwrap(), p2.sigma_position.unwrap());
        assert!(
            ((s2 / s1) - 2.0).abs() < 1e-9,
            "posterior σ ratio {} is not 2 (σ doubled)",
            s2 / s1
        );
    }

    /// A rank-deficient geometry yields NO finite posterior summary — `None`, never a
    /// fabricated number — and reports its unobservable directions.
    #[test]
    fn rank_deficient_posterior_is_none_with_a_null_space() {
        // A single instantaneous range row: rank 1 of 4.
        let chief: PlanarState = [1.10, 0.02, 0.05, -0.50];
        let reference: PlanarState = [1.02, -0.03, -0.06, -0.55];
        let (_rho, row) = range_row(&chief, &reference);
        let o: Mat = vec![row.to_vec()];
        let p = whitened_posterior(&o, 2, 1e-9);
        assert_eq!(p.rank, 1);
        assert_eq!(p.defect, 3);
        assert!(p.sigma_position.is_none() && p.sigma_velocity.is_none());
        assert_eq!(p.null_space.len(), N_PLANAR);
        assert_eq!(p.null_space[0].len(), 3);
        // The condition number is taken over the SAME observable subspace the rank
        // defines (one mode here), exactly as `gramian_spectrum` does — so it is 1, and
        // the honest "no finite covariance" verdict is carried by the `None` summaries
        // and the non-zero defect, not by an infinite condition.
        assert!(
            (p.condition - 1.0).abs() < 1e-12,
            "condition {}",
            p.condition
        );
    }

    #[test]
    fn rank_vs_arc_grows_and_reaches_full_rank() {
        // A single-link range-only arc: instantaneously rank-1, growing to full rank 4
        // as the arc lengthens and the STM couples position into velocity.
        let chief: PlanarState = [1.10, 0.02, 0.05, -0.50];
        let reference: PlanarState = [1.02, -0.03, -0.06, -0.55];
        let mu = EARTH_MOON_MU;
        let n_epochs = 24;
        let arc = 0.06_f64; // ~6 rotating-frame hours of coupling
        let mut epochs = Vec::new();
        let mut prev = 0.0;
        for k in 0..n_epochs {
            let t = arc * (k as f64) / ((n_epochs - 1) as f64);
            let (cs, phi) = planar_state_stm(&chief, mu, t, 3000);
            let rs = planar_propagate(&reference, mu, t, 3000);
            let (_rho, r_row) = range_row(&cs, &rs);
            epochs.push(ObsEpoch {
                h: vec![r_row.to_vec()],
                phi: phi.iter().map(|r| r.to_vec()).collect(),
                dt: t - prev,
            });
            prev = t;
        }
        let table = rank_vs_arc(&epochs, 1e-6);
        // First epoch (a single instantaneous range) is rank 1.
        assert_eq!(table[0].rank, 1, "single instantaneous range is rank-1");
        // Rank is non-decreasing along the arc.
        for w in table.windows(2) {
            assert!(
                w[1].rank >= w[0].rank,
                "rank must not decrease along the arc"
            );
        }
        // By the end of the arc the four-state is fully observable.
        assert_eq!(
            table.last().unwrap().rank,
            N_PLANAR,
            "full observability over arc"
        );
    }

    // ── The rank read can no longer report more rank than the matrix has rows ──────

    /// The noise-floor constant is exactly `f64::EPSILON.sqrt()`, written out because `sqrt`
    /// is not a `const fn`. If the two ever drift apart, every message that quotes the floor
    /// is quoting a different number from the one the code compares against.
    #[test]
    fn the_rank_tolerance_noise_floor_is_the_square_root_of_machine_epsilon() {
        assert_eq!(RANK_TOLERANCE_NOISE_FLOOR, f64::EPSILON.sqrt());
    }

    /// A two-row, four-column observability matrix cannot have rank 3, whatever the
    /// tolerance. Built so the raw singular-value count DOES exceed the shape bound at a
    /// tight `rel_tol` (the reconstructed spectrum carries two f64-residue directions below
    /// the noise floor), so this exercises the clamp rather than asserting a tautology.
    #[test]
    fn a_two_row_matrix_cannot_report_rank_three() {
        let o: Mat = vec![vec![1.0, 2.0, 3.0, 4.0], vec![4.0, 3.0, 2.0, 1.0]];
        let sv = singular_values(&o);
        // The true rank is 2: the Gram has exactly two non-zero eigenvalues, and the other
        // two come back as rounding residue rather than exact zeros.
        let counted = rank_from_singular_values(&sv, 1e-14);
        assert!(
            counted > 2,
            "setup: the unbounded count must exceed the row count for this to test anything \
             (counted {counted}, spectrum {sv:?})"
        );
        let read = bounded_rank_from_singular_values(&sv, 1e-14, o.len(), 4);
        assert_eq!(read.rank, 2, "rank is held to min(rows, cols) = 2");
        assert_eq!(read.counted, counted);
        assert_eq!(read.shape_bound, 2);
        let reason = read.reason.expect("a clamp must state its reason");
        assert!(
            reason.contains("min(2, 4) = 2") && reason.contains("noise floor"),
            "the stated reason must name the bound and the floor: {reason}"
        );
        assert_eq!(observable_rank(&o, 1e-14), 2);
        // The posterior read is held to the same bound, and its defect, null space and the
        // sigma summaries are rebuilt over the bounded subspace rather than left disagreeing
        // with the reported rank.
        let post = whitened_posterior(&o, 2, 1e-14);
        assert_eq!(post.rank, 2);
        assert_eq!(post.defect, 2);
        assert_eq!(post.null_space[0].len(), 2, "null space matches the defect");
        assert!(
            post.sigma_position.is_none() && post.sigma_velocity.is_none(),
            "an underdetermined batch has no finite covariance"
        );
        assert!(post.rank_limited_by.is_some(), "the clamp states itself");
    }

    /// R1 — above the f64 noise floor the bound never binds, so every rank read is
    /// bit-identical to the unbounded one it replaced. Asserted over the tolerances the
    /// repository actually ships at, on a growing single-link arc.
    #[test]
    fn the_bound_changes_nothing_at_any_shipped_tolerance() {
        let chief: PlanarState = [1.10, 0.02, 0.05, -0.50];
        let reference: PlanarState = [1.02, -0.03, -0.06, -0.55];
        let mu = EARTH_MOON_MU;
        let n_epochs = 12;
        let arc = 0.06_f64;
        let mut epochs = Vec::new();
        let mut prev = 0.0;
        for k in 0..n_epochs {
            let t = arc * (k as f64) / ((n_epochs - 1) as f64);
            let (cs, phi) = planar_state_stm(&chief, mu, t, 1500);
            let rs = planar_propagate(&reference, mu, t, 1500);
            let (_rho, r_row) = range_row(&cs, &rs);
            epochs.push(ObsEpoch {
                h: vec![r_row.to_vec()],
                phi: phi.iter().map(|r| r.to_vec()).collect(),
                dt: t - prev,
            });
            prev = t;
        }
        for rel_tol in [1e-4_f64, 1e-5, 1e-6, 1e-7, 1e-8] {
            assert!(
                rel_tol >= RANK_TOLERANCE_NOISE_FLOOR || rel_tol == 1e-8,
                "the shipped tolerances sit at or above the noise floor"
            );
            let table = rank_vs_arc(&epochs, rel_tol);
            for (k, point) in table.iter().enumerate() {
                let (o, _w) = observability_matrix(&epochs[..=k]);
                let unbounded = rank_from_singular_values(&singular_values(&o), rel_tol);
                assert_eq!(
                    point.rank, unbounded,
                    "rel_tol {rel_tol:e}, prefix {k}: the bound moved a shipped value"
                );
                assert!(
                    point.rank_limited_by.is_none(),
                    "rel_tol {rel_tol:e}, prefix {k}: nothing should be clamped here"
                );
            }
        }
    }

    /// The tolerance itself is never re-floored — a read below the noise floor is still
    /// taken at the tolerance asked for, and is reported as such.
    #[test]
    fn a_tolerance_below_the_noise_floor_is_reported_not_rewritten() {
        assert!(rank_tolerance_note(1e-6).is_none());
        assert!(rank_tolerance_note(RANK_TOLERANCE_NOISE_FLOOR).is_none());
        let note = rank_tolerance_note(1e-12).expect("a sub-floor tolerance must be named");
        assert!(
            note.contains("1.000e-12") && note.contains("1.490e-8"),
            "the note must quote both the tolerance and the floor: {note}"
        );
        // Rank reads at a loose tolerance are unaffected by the note's existence.
        let o: Mat = vec![vec![1.0, 0.0], vec![0.0, 1.0]];
        assert_eq!(observable_rank(&o, 1e-6), 2);
    }
}