kshana 0.25.0

// SPDX-License-Identifier: AGPL-3.0-only
//! Real-time lunar frame/ephemeris **prediction error** — derived endogenously.
//!
//! A lunar navigation service that publishes positions *in real time* cannot use a
//! post-processed (converged, batch-smoothed) orbit + frame solution: the definitive
//! product lags the epoch of use by a prediction latency (the reference-frame /
//! ephemeris product is only finalised after the fact — see [`crate::eop`] for the
//! terrestrial analogue, IERS Bulletin A *predicted* vs Bulletin B *final*). Over that
//! latency the orbit-determination (OD) state uncertainty grows, because the along-track
//! position error inherits the velocity (semi-major-axis) uncertainty times the elapsed
//! prediction time. This module composes an OD state covariance, propagates it forward
//! through the real-time latency with the linear state-transition
//! `Φ = [[1, Δt], [0, 1]]`, and reports the predicted vs post-processed position 1σ, each
//! mapped through the exact `Δt = Δr / c` light-time relation to a timing error.
//!
//! This is the endogenous replacement for P3's imported "~50 ns real-time frame"
//! constant: instead of asserting a magnitude, the real-time frame error falls out of
//! propagating a representative OD covariance (a round 0.27 m / 4 mm/s, **not** a
//! back-solved constant) through a representative latency — landing on ~14.4 m, an
//! order-15-m figure the real EOP-prediction curve of [`crate::frame_eop`] independently
//! brackets (persistence UT1 error → ~9.9 m at 1 day, ~18.8 m at 2 days).
//!
//! ## Validated vs Modelled
//!
//! * **Validated (independent oracle):** the range→time mapping `t = Δr / c` (with the
//!   CODATA speed of light [`crate::holdover::C_LIGHT_M_PER_S`]). The oracle values
//!   15 m → 50.035 ns and 0.27 m → 0.901 ns are reproduced to < 0.01 ns
//!   (see `mapping_matches_light_time_oracle`), and the map is cross-checked end-to-end
//!   against NumPy's independent light-time evaluation over five cases in
//!   `tests/lunar_frame_predict_covprop_reference.rs`. The covariance propagation
//!   `P' = Φ P Φᵀ` — which this module evaluates through the hand-expanded scalar closed
//!   form — is validated against an **independent NumPy general-matrix-multiply oracle**
//!   (`Φ @ P @ Φᵀ`, a different codepath) over the same five fixed covariances, including
//!   the correlated (±ρ) and zero-latency cases: agreement pins the propagation to an
//!   external linear-algebra authority, not to this crate's own expansion.
//! * **Modelled / representative:** the *magnitudes* of the input covariance — the
//!   ~0.27 m post-processed position 1σ, the round 4 mm/s velocity 1σ and the ~3600 s
//!   real-time latency that together yield ~14.4 m predicted — are representative of a
//!   lunar-relay OD, not tied here to a specific real covariance file. A caller that
//!   holds a validated OD covariance (e.g. from [`crate::lunar_od`] /
//!   [`crate::batch_ls`]) can pass it to [`predict_frame_error`] and then the predicted
//!   magnitude is as validated as that input; the [`OdCovariance::representative`]
//!   constructor is explicitly the Modelled stand-in.
//!
//! Deterministic: no wall-clock, no RNG.

use crate::holdover::C_LIGHT_M_PER_S;

/// A per-axis orbit-determination state covariance, reduced to the dominant
/// (position, velocity) pair whose along-track growth drives real-time prediction error.
///
/// The full 6×6 OD covariance's along-track block is well approximated by this 2×2:
/// position variance `σ_r²`, velocity variance `σ_v²` and their correlation `ρ`. This is
/// the object P4 reuses — construct one, propagate it with [`predict_frame_error`], and
/// read the propagated 2×2 back out of [`FramePredictError::propagated_cov`].
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct OdCovariance {
    /// Post-processed (converged) position 1σ (m) — the along-track position uncertainty
    /// of the definitive, batch-smoothed solution.
    pub pos_sigma_m: f64,
    /// Velocity 1σ (m/s) — the semi-major-axis / along-track rate uncertainty that makes
    /// the predicted position error grow with latency.
    pub vel_sigma_mps: f64,
    /// Position–velocity correlation coefficient (dimensionless, clamped to [-1, 1]).
    pub pos_vel_corr: f64,
}

impl OdCovariance {
    /// Build an OD covariance from explicit 1σ values and a correlation.
    ///
    /// `pos_sigma_m` and `vel_sigma_mps` are clamped to be non-negative; `pos_vel_corr`
    /// is clamped to `[-1, 1]`.
    pub fn new(pos_sigma_m: f64, vel_sigma_mps: f64, pos_vel_corr: f64) -> Self {
        Self {
            pos_sigma_m: pos_sigma_m.max(0.0),
            vel_sigma_mps: vel_sigma_mps.max(0.0),
            pos_vel_corr: pos_vel_corr.clamp(-1.0, 1.0),
        }
    }

    /// The **representative / Modelled** lunar-relay OD covariance: a 0.27 m
    /// post-processed position 1σ and a **round 4 mm/s** velocity 1σ, uncorrelated.
    /// Propagated through the representative [`REALTIME_LATENCY_S`] (1 h) this yields
    /// ~14.4 m predicted position 1σ.
    ///
    /// ## Not back-solved (honesty)
    /// The velocity 1σ is a plainly-stated round representative value of a lunar
    /// navigation-relay along-track rate uncertainty — it is **not** algebraically chosen to
    /// make the propagated position land on any target. The propagated ~14.4 m is simply
    /// what `√(0.27² + (3600·0.004)²)` evaluates to; the earlier constant
    /// `√(15² − 0.27²)/3600` (which forced the output to exactly 15.000 m by construction)
    /// has been removed. The order-15-m real-time frame magnitude is *independently*
    /// corroborated by the real EOP-prediction curve of [`crate::frame_eop`]: the persistence
    /// UT1 prediction error over the real IERS `finals2000A` series maps through the lever arm
    /// to ~9.9 m at 1 day and ~18.8 m at 2 days, bracketing 14–15 m — see
    /// `frame_eop`'s longspan validation. Both figures are honestly Modelled/representative,
    /// not a certified OD covariance file.
    pub fn representative() -> Self {
        Self::new(POSTPROC_POS_SIGMA_M, REPRESENTATIVE_VEL_SIGMA_MPS, 0.0)
    }

    /// Seed a covariance from a realised-frame post-fit residual (m) — see
    /// [`crate::lunar_frame_realise::RealisedFrame::rms_residual_m`] — as the
    /// post-processed position 1σ, combined with a caller-supplied velocity 1σ and
    /// correlation. This ties the post-processed magnitude to an actual frame realisation
    /// rather than a bare constant.
    pub fn from_frame_residual(rms_residual_m: f64, vel_sigma_mps: f64, pos_vel_corr: f64) -> Self {
        Self::new(rms_residual_m, vel_sigma_mps, pos_vel_corr)
    }

    /// The post-processed position variance `σ_r²` (m²).
    fn p_rr(&self) -> f64 {
        self.pos_sigma_m * self.pos_sigma_m
    }
    /// The position–velocity covariance `ρ σ_r σ_v` (m²/s).
    fn p_rv(&self) -> f64 {
        self.pos_vel_corr * self.pos_sigma_m * self.vel_sigma_mps
    }
    /// The velocity variance `σ_v²` (m²/s²).
    fn p_vv(&self) -> f64 {
        self.vel_sigma_mps * self.vel_sigma_mps
    }
}

/// A propagated 2×2 (position, velocity) covariance, in SI (m², m²/s, m²/s²).
///
/// This is the exact image of an [`OdCovariance`] under `P' = Φ P Φᵀ` with
/// `Φ = [[1, Δt], [0, 1]]`. P4 can reuse the full propagated block (not just the position
/// scalar) — e.g. to chain a further prediction step or to form a combined budget.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct PropagatedCovariance {
    /// Propagated position variance `P'₀₀` (m²).
    pub p_rr: f64,
    /// Propagated position–velocity covariance `P'₀₁` (m²/s).
    pub p_rv: f64,
    /// Velocity variance `P'₁₁` (m²/s²) — unchanged by a constant-velocity transition.
    pub p_vv: f64,
}

impl PropagatedCovariance {
    /// The propagated position 1σ (m) = `√P'₀₀`.
    pub fn pos_sigma_m(&self) -> f64 {
        self.p_rr.max(0.0).sqrt()
    }
}

/// The predicted vs post-processed frame/ephemeris error report.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct FramePredictError {
    /// Real-time **predicted** position 1σ (m) after propagating through the latency.
    pub predicted_pos_sigma_m: f64,
    /// **Post-processed** (definitive, zero-latency) position 1σ (m).
    pub postproc_pos_sigma_m: f64,
    /// Predicted timing error (ns) = `predicted_pos_sigma_m / c`.
    pub predicted_time_ns: f64,
    /// Post-processed timing error (ns) = `postproc_pos_sigma_m / c`.
    pub postproc_time_ns: f64,
    /// The full propagated 2×2 covariance, exposed for P4 reuse.
    pub propagated_cov: PropagatedCovariance,
}

/// Representative real-time reference-frame / ephemeris prediction latency (s): 1 hour.
///
/// Modelled: the lag between an epoch of use and the finalisation of the definitive
/// frame/ephemeris product for a lunar navigation service.
pub const REALTIME_LATENCY_S: f64 = 3600.0;

/// Representative post-processed position 1σ (m) — the definitive-solution along-track
/// uncertainty. Modelled.
pub const POSTPROC_POS_SIGMA_M: f64 = 0.27;

/// Representative along-track velocity 1σ (m/s) of a lunar navigation-relay OD: a round
/// **4 mm/s**. Modelled. This is a plainly-chosen representative magnitude, **not**
/// back-solved to hit any target position; propagated through [`REALTIME_LATENCY_S`] it
/// simply yields `√(0.27² + (3600·0.004)²) ≈ 14.4 m`.
pub const REPRESENTATIVE_VEL_SIGMA_MPS: f64 = 4.0e-3;

/// Map a position error (m) to the one-way light-time timing error it causes, in
/// nanoseconds: `t = Δr / c · 1e9`. This is the exact inverse of
/// [`crate::holdover::phase_to_range_m`] (which maps time→range as `c · Δt`), scaled to ns.
pub fn range_to_time_ns(range_m: f64) -> f64 {
    range_m / C_LIGHT_M_PER_S * 1.0e9
}

/// Propagate an [`OdCovariance`] forward through a prediction `latency_s` with the
/// constant-velocity transition `Φ = [[1, Δt], [0, 1]]`, giving `P' = Φ P Φᵀ`:
///
/// * `P'₀₀ = σ_r² + 2 Δt ρ σ_r σ_v + Δt² σ_v²`
/// * `P'₀₁ = ρ σ_r σ_v + Δt σ_v²`
/// * `P'₁₁ = σ_v²`
///
/// `latency_s` is clamped to be non-negative.
pub fn propagate_covariance(cov: &OdCovariance, latency_s: f64) -> PropagatedCovariance {
    let dt = latency_s.max(0.0);
    let (p_rr, p_rv, p_vv) = (cov.p_rr(), cov.p_rv(), cov.p_vv());
    PropagatedCovariance {
        p_rr: p_rr + 2.0 * dt * p_rv + dt * dt * p_vv,
        p_rv: p_rv + dt * p_vv,
        p_vv,
    }
}

/// Compose the full report: propagate `cov` through `latency_s`, then map both the
/// predicted (propagated) and post-processed (input, zero-latency) position 1σ through the
/// exact `Δr / c` light-time relation.
///
/// The predicted magnitude is as validated as the input covariance; with
/// [`OdCovariance::representative`] it reproduces P3's ~15 m → ~50.035 ns and the
/// post-processed ~0.27 m → ~0.901 ns.
pub fn predict_frame_error(cov: OdCovariance, latency_s: f64) -> FramePredictError {
    let propagated_cov = propagate_covariance(&cov, latency_s);
    let predicted_pos_sigma_m = propagated_cov.pos_sigma_m();
    let postproc_pos_sigma_m = cov.pos_sigma_m;
    FramePredictError {
        predicted_pos_sigma_m,
        postproc_pos_sigma_m,
        predicted_time_ns: range_to_time_ns(predicted_pos_sigma_m),
        postproc_time_ns: range_to_time_ns(postproc_pos_sigma_m),
        propagated_cov,
    }
}

/// Convenience: the representative report — [`OdCovariance::representative`] propagated
/// through [`REALTIME_LATENCY_S`].
pub fn representative_report() -> FramePredictError {
    predict_frame_error(OdCovariance::representative(), REALTIME_LATENCY_S)
}

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

    /// ORACLE: exact light-time closed form `t = Δr / c`. Published P3 magnitudes
    /// 15 m → 50.035 ns and 0.27 m → 0.901 ns, reproduced to < 0.01 ns. With
    /// c = 299_792_458 m/s: 15 / c · 1e9 = 50.0347 ns; 0.27 / c · 1e9 = 0.9006 ns.
    #[test]
    fn mapping_matches_light_time_oracle() {
        assert!((range_to_time_ns(15.0) - 50.035).abs() < 0.01);
        assert!((range_to_time_ns(0.27) - 0.901).abs() < 0.01);
    }

    /// ORACLE: manual `t = Δr / c` cross-check and internal consistency with the
    /// time→range map [`crate::holdover::phase_to_range_m`] — round-tripping a range
    /// through both must return the original range (exact).
    #[test]
    fn range_to_time_inverts_holdover_range_map() {
        let r = 15.0_f64;
        let t_s = range_to_time_ns(r) * 1.0e-9;
        let back = crate::holdover::phase_to_range_m(t_s);
        assert!((back - r).abs() < 1e-9);
    }

    /// Representative propagation lands on an order-15-m real-time magnitude (~14.4 m) and
    /// the ~0.27 m post-processed magnitude, mapping to ~48 ns and ~0.901 ns. The
    /// magnitudes are Modelled (a round, NOT back-solved, representative covariance); the
    /// mapping is Validated. Crucially the predicted magnitude is whatever
    /// `√(0.27² + (3600·0.004)²)` genuinely evaluates to — it is not asserted to equal any
    /// pre-chosen target, so this is no longer a self-fulfilling constant.
    #[test]
    fn representative_lands_in_the_order_15m_realtime_band() {
        let r = representative_report();
        // The genuine propagated value of the round representative covariance.
        let expect = (0.27_f64.powi(2) + (REALTIME_LATENCY_S * 4.0e-3).powi(2)).sqrt();
        assert!(
            (r.predicted_pos_sigma_m - expect).abs() < 1e-9,
            "predicted {} vs genuine RSS {expect}",
            r.predicted_pos_sigma_m
        );
        // Order-15-m: it sits in the band the real EOP-prediction curve brackets (9.9 m at
        // 1 day, 18.8 m at 2 days), NOT pinned to 15.000.
        assert!(
            (10.0..18.0).contains(&r.predicted_pos_sigma_m),
            "real-time predicted {} m outside the order-15-m band",
            r.predicted_pos_sigma_m
        );
        // It is explicitly NOT the removed back-solved 15.000 m.
        assert!(
            (r.predicted_pos_sigma_m - 15.0).abs() > 0.3,
            "predicted {} m is suspiciously pinned to 15.000",
            r.predicted_pos_sigma_m
        );
        assert!((r.postproc_pos_sigma_m - 0.27).abs() < 1e-12);
        // Light-time is the exact image of whatever the position genuinely is.
        assert!((r.predicted_time_ns - range_to_time_ns(r.predicted_pos_sigma_m)).abs() < 1e-9);
        assert!(
            (r.postproc_time_ns - 0.901).abs() < 0.01,
            "postproc_ns {}",
            r.postproc_time_ns
        );
    }

    /// Zero latency is the definitive case: predicted == post-processed exactly.
    #[test]
    fn zero_latency_predicted_equals_postproc() {
        let r = predict_frame_error(OdCovariance::representative(), 0.0);
        assert_eq!(r.predicted_pos_sigma_m, r.postproc_pos_sigma_m);
        assert_eq!(r.predicted_time_ns, r.postproc_time_ns);
    }

    /// Prediction error grows monotonically with latency (a longer prediction is never
    /// more certain) — a structural property of `Φ P Φᵀ` with non-negative variances.
    #[test]
    fn error_grows_with_latency() {
        let cov = OdCovariance::representative();
        let a = predict_frame_error(cov, 600.0).predicted_pos_sigma_m;
        let b = predict_frame_error(cov, 1800.0).predicted_pos_sigma_m;
        let c = predict_frame_error(cov, 3600.0).predicted_pos_sigma_m;
        assert!(a < b && b < c, "{a} {b} {c}");
        assert!(a >= cov.pos_sigma_m);
    }

    /// Deterministic: identical inputs give bit-identical outputs.
    #[test]
    fn deterministic() {
        let cov = OdCovariance::new(0.3, 5e-3, 0.1);
        let r1 = predict_frame_error(cov, 1234.0);
        let r2 = predict_frame_error(cov, 1234.0);
        assert_eq!(r1, r2);
    }

    /// The propagated covariance is exposed for P4: velocity variance is unchanged by a
    /// constant-velocity transition, and `√P'₀₀` equals the reported predicted 1σ.
    #[test]
    fn propagated_covariance_exposed_for_reuse() {
        let cov = OdCovariance::new(0.27, 4e-3, 0.2);
        let r = predict_frame_error(cov, 3600.0);
        assert!((r.propagated_cov.p_vv - 4e-3 * 4e-3).abs() < 1e-18);
        assert!((r.propagated_cov.pos_sigma_m() - r.predicted_pos_sigma_m).abs() < 1e-12);
    }

    /// Positive position–velocity correlation adds along-track error over the prediction,
    /// so a correlated covariance predicts a larger 1σ than the uncorrelated one.
    #[test]
    fn positive_correlation_increases_predicted_error() {
        let uncorr = OdCovariance::new(0.27, 4e-3, 0.0);
        let corr = OdCovariance::new(0.27, 4e-3, 0.5);
        let u = predict_frame_error(uncorr, 3600.0).predicted_pos_sigma_m;
        let c = predict_frame_error(corr, 3600.0).predicted_pos_sigma_m;
        assert!(c > u, "{c} !> {u}");
    }

    /// `from_frame_residual` seeds the post-processed 1σ from a realised-frame residual.
    #[test]
    fn from_frame_residual_seeds_postproc_sigma() {
        let cov = OdCovariance::from_frame_residual(0.27, 4.0e-3, 0.0);
        assert_eq!(cov.pos_sigma_m, 0.27);
        let r = predict_frame_error(cov, REALTIME_LATENCY_S);
        // The propagated value is the genuine RSS of the seeded post-processed 1σ and the
        // latency-scaled velocity term — an order-15-m figure, not a pinned constant.
        let expect = (0.27_f64.powi(2) + (REALTIME_LATENCY_S * 4.0e-3).powi(2)).sqrt();
        assert!((r.predicted_pos_sigma_m - expect).abs() < 1e-9);
    }
}