kshana 0.37.0

// SPDX-License-Identifier: AGPL-3.0-only
//! Two-way satellite time and frequency transfer (TWSTFT): the hardware-delay equation and
//! its uncertainty budget, a two-way link with flicker phase noise and a diurnal
//! non-reciprocity, and the stability of real `UTC − UTC(k)` and IGS station-clock series.
//!
//! ## The two-way equation
//!
//! Each station's time-interval counter is started by its own clock tick and stopped by the
//! other station's signal. With `T_A`, `T_B` the two clocks, the counter readings are
//!
//! ```text
//!   TI(B) = (T_A − T_B) + TX_A + ρ_AS/c + up_A + SAT_AB + ρ_SB/c + down_B + RX_B + S(A→S→B)
//!   TI(A) = (T_B − T_A) + TX_B + ρ_BS/c + up_B + SAT_BA + ρ_SA/c + down_A + RX_A + S(B→S→A)
//! ```
//!
//! The terms are:
//! - `TX` and `RX`: transmit and receive equipment delays;
//! - `up` and `down`: frequency-dependent atmospheric delays of each station's uplink and
//!   downlink;
//! - `SAT`: transponder delays;
//! - `S(·)`: the Sagnac light-time term of [`crate::timegeo::sagnac_correction`] along each
//!   path.
//!
//! The geometric ranges are reciprocal and cancel in the difference, which gives
//!
//! ```text
//!   T_A − T_B = ½[TI(B) − TI(A)]
//!             − ½[(TX_A − RX_A) − (TX_B − RX_B)]
//!             − ½[(up_A − down_A) − (up_B − down_B)]
//!             − ½[SAT_AB − SAT_BA]
//!             − ½[S(A→S→B) − S(B→S→A)]
//! ```
//!
//! This is the structure of Recommendation ITU-R TF.1153, written in this module's own sign
//! convention. Two further pieces:
//!
//! - **Ionosphere.** Its term is `40.3·TEC·(1/f_up² − 1/f_down²)/c` per station, because the
//!   uplink and downlink use different frequencies. The troposphere is frequency-independent
//!   at Ku band and cancels.
//! - **Calibration.** The station-delay difference `(TX_A − RX_A) − (TX_B − RX_B)` is what a
//!   travelling calibration station measures. It enters as one calibrated quantity with its
//!   uncertainty.
//!
//! ## Uncertainty budget
//!
//! Each term carries a standard uncertainty. The budget combines them in quadrature, with the
//! statistical (type A) term taken from the link's time deviation at the averaging time.
//!
//! ## Link noise
//!
//! [`two_way_series`] adds three things to white phase noise:
//!
//! - a **flicker phase-noise** floor, generated by the Kasdin–Walter fractional-difference
//!   filter (α = 1 in phase);
//! - a **diurnal non-reciprocity**, a sinusoid of stated amplitude and one-day period, the
//!   signature of temperature-driven delay variation in the ground equipment;
//! - a constant **reciprocity residual**.
//!
//! The constant residual shifts the offset. It is a bias, and the TDEV cannot see it.
//!
//! HONEST SCOPE (MODELLED): the equation and the budget are bookkeeping with stated inputs.
//! No real TWSTFT data or calibration report has been compared. The noise model is a
//! parametric stand-in, not a fit to a real link.

use crate::allan::{overlapping_adev, time_deviation};
use crate::portable_distr::{Distribution, StandardNormal};
use crate::portable_math::PortableFloat;
use crate::timegeo::{sagnac_correction, C_M_PER_S};
use rand::SeedableRng;
use rand_chacha::ChaCha8Rng;

type Vec3 = [f64; 3];

/// Ionospheric group delay (s) for `tec_tecu` total electron content (TECU) at `f_hz`.
pub fn iono_delay_s(tec_tecu: f64, f_hz: f64) -> f64 {
    40.3 * tec_tecu * 1e16 / (f_hz * f_hz) / C_M_PER_S
}

/// One ground station of a TWSTFT link.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Station {
    /// Earth-fixed position (m).
    pub r_m: Vec3,
    /// Transmit equipment delay (s).
    pub tx_delay_s: f64,
    /// Receive equipment delay (s).
    pub rx_delay_s: f64,
    /// Uplink carrier (Hz).
    pub uplink_hz: f64,
    /// Downlink carrier (Hz).
    pub downlink_hz: f64,
    /// Slant TEC of this station's path to the satellite (TECU).
    pub slant_tec_tecu: f64,
    /// Slant tropospheric delay (s); frequency-independent, so it cancels.
    pub tropo_s: f64,
}

/// The satellite of a TWSTFT link.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Relay {
    /// Earth-fixed position (m).
    pub r_m: Vec3,
    /// Transponder delay for the A→B path (s).
    pub delay_ab_s: f64,
    /// Transponder delay for the B→A path (s).
    pub delay_ba_s: f64,
}

/// `up − down` for station `s`: the uplink is at `uplink_hz`, the downlink at `downlink_hz`,
/// through the same slant TEC; the troposphere cancels.
fn path_asymmetry(s: &Station) -> f64 {
    iono_delay_s(s.slant_tec_tecu, s.uplink_hz) - iono_delay_s(s.slant_tec_tecu, s.downlink_hz)
}

fn norm(a: Vec3, b: Vec3) -> f64 {
    ((a[0] - b[0]).ppowi(2) + (a[1] - b[1]).ppowi(2) + (a[2] - b[2]).ppowi(2)).sqrt()
}

/// The two counter readings `(TI(A), TI(B))` for a true offset `T_A − T_B` (forward model).
pub fn counter_readings(a: &Station, b: &Station, sat: &Relay, offset_s: f64) -> (f64, f64) {
    let rho_a = norm(a.r_m, sat.r_m) / C_M_PER_S;
    let rho_b = norm(b.r_m, sat.r_m) / C_M_PER_S;
    let up = |s: &Station| iono_delay_s(s.slant_tec_tecu, s.uplink_hz) + s.tropo_s;
    let down = |s: &Station| iono_delay_s(s.slant_tec_tecu, s.downlink_hz) + s.tropo_s;
    let s_ab = sagnac_correction(a.r_m, sat.r_m) + sagnac_correction(sat.r_m, b.r_m);
    let s_ba = sagnac_correction(b.r_m, sat.r_m) + sagnac_correction(sat.r_m, a.r_m);
    let ti_b = offset_s
        + a.tx_delay_s
        + rho_a
        + up(a)
        + sat.delay_ab_s
        + rho_b
        + down(b)
        + b.rx_delay_s
        + s_ab;
    let ti_a = -offset_s
        + b.tx_delay_s
        + rho_b
        + up(b)
        + sat.delay_ba_s
        + rho_a
        + down(a)
        + a.rx_delay_s
        + s_ba;
    (ti_a, ti_b)
}

/// Each term of the two-way equation (s), for reporting.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct TwoWayTerms {
    /// `½[TI(B) − TI(A)]`.
    pub half_counter_difference: f64,
    /// `−½[(TX_A − RX_A) − (TX_B − RX_B)]`.
    pub station_delays: f64,
    /// `−½[(up_A − down_A) − (up_B − down_B)]` (ionosphere).
    pub ionosphere: f64,
    /// `−½[SAT_AB − SAT_BA]`.
    pub transponder: f64,
    /// `−½[S(A→S→B) − S(B→S→A)]`.
    pub sagnac: f64,
}

impl TwoWayTerms {
    /// `T_A − T_B`, the sum of the terms.
    pub fn offset(&self) -> f64 {
        self.half_counter_difference
            + self.station_delays
            + self.ionosphere
            + self.transponder
            + self.sagnac
    }
}

/// Solve the two-way equation for `T_A − T_B` from the counter readings and the calibrated
/// delays.
pub fn solve_offset(a: &Station, b: &Station, sat: &Relay, ti_a: f64, ti_b: f64) -> TwoWayTerms {
    let s_ab = sagnac_correction(a.r_m, sat.r_m) + sagnac_correction(sat.r_m, b.r_m);
    let s_ba = sagnac_correction(b.r_m, sat.r_m) + sagnac_correction(sat.r_m, a.r_m);
    TwoWayTerms {
        half_counter_difference: 0.5 * (ti_b - ti_a),
        station_delays: -0.5 * ((a.tx_delay_s - a.rx_delay_s) - (b.tx_delay_s - b.rx_delay_s)),
        ionosphere: -0.5 * (path_asymmetry(a) - path_asymmetry(b)),
        transponder: -0.5 * (sat.delay_ab_s - sat.delay_ba_s),
        sagnac: -0.5 * (s_ab - s_ba),
    }
}

/// One line of the uncertainty budget.
#[derive(Clone, Debug, PartialEq)]
pub struct BudgetLine {
    /// What the term is.
    pub name: String,
    /// Standard uncertainty (s).
    pub u_s: f64,
    /// `"A"` (statistical) or `"B"` (other).
    pub kind: &'static str,
}

/// Combine a budget in quadrature (s).
pub fn combined_uncertainty(lines: &[BudgetLine]) -> f64 {
    lines.iter().map(|l| l.u_s * l.u_s).sum::<f64>().sqrt()
}

/// Two-way link noise.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct LinkNoise {
    /// White phase noise per measurement, 1-σ (s).
    pub white_s: f64,
    /// Flicker phase noise: standard deviation of the unit-variance-driven filter output scale (s).
    pub flicker_s: f64,
    /// Diurnal non-reciprocity amplitude (s).
    pub diurnal_amp_s: f64,
    /// Constant reciprocity residual (s).
    pub reciprocity_bias_s: f64,
}

/// Flicker (α = 1) phase noise of `n` samples by the Kasdin–Walter fractional-difference
/// filter: `h_0 = 1`, `h_k = h_{k−1}(k − 1 + α/2)/k`, convolved with unit white noise and scaled
/// by `scale`.
pub fn flicker_phase(n: usize, scale: f64, seed: u64) -> Vec<f64> {
    let mut rng = ChaCha8Rng::seed_from_u64(seed);
    let w: Vec<f64> = (0..n).map(|_| StandardNormal.sample(&mut rng)).collect();
    let mut h = vec![1.0; n];
    for k in 1..n {
        h[k] = h[k - 1] * ((k as f64 - 1.0) + 0.5) / k as f64;
    }
    (0..n)
        .map(|i| scale * (0..=i).map(|k| h[k] * w[i - k]).sum::<f64>())
        .collect()
}

/// A simulated `T_A − T_B` measurement series around `offset_s`, sampled every `step_s`: white
/// and flicker phase noise, a diurnal sinusoid and a constant reciprocity residual.
pub fn two_way_series(
    offset_s: f64,
    noise: &LinkNoise,
    step_s: f64,
    n: usize,
    seed: u64,
) -> Vec<f64> {
    let mut rng = ChaCha8Rng::seed_from_u64(seed);
    let fl = if noise.flicker_s > 0.0 {
        flicker_phase(n, noise.flicker_s, seed ^ 0xf11c)
    } else {
        vec![0.0; n]
    };
    (0..n)
        .map(|i| {
            let t = i as f64 * step_s;
            let w: f64 = StandardNormal.sample(&mut rng);
            offset_s
                + noise.reciprocity_bias_s
                + noise.white_s * w
                + fl[i]
                + noise.diurnal_amp_s * (2.0 * std::f64::consts::PI * t / 86_400.0).psin()
        })
        .collect()
}

/// `(τ, TDEV)` at octave averaging factors of a phase series.
pub fn tdev_curve(x: &[f64], step_s: f64) -> Vec<(f64, f64)> {
    let mut out = Vec::new();
    let mut m = 1usize;
    while x.len() >= 3 * m {
        let td = time_deviation(x, step_s, m);
        if td.is_finite() {
            out.push((m as f64 * step_s, td));
        }
        m *= 2;
    }
    out
}

/// Stability summary of a time-offset series `(t_s, offset_s)` on a regular grid: the least-
/// squares rate (s/s), the residual RMS (s) after removing mean and rate, and the TDEV and ADEV
/// curves of the residual at octave averaging factors.
#[derive(Clone, Debug, PartialEq)]
pub struct SeriesStability {
    pub n: usize,
    pub step_s: f64,
    pub rate: f64,
    pub residual_rms_s: f64,
    pub tdev: Vec<(f64, f64)>,
    pub adev: Vec<(f64, f64)>,
}

/// The stability of a regularly sampled offset series. Errors when the samples are fewer than
/// eight or not on a regular grid (within 1e-6 of a step).
pub fn series_stability(series: &[(f64, f64)]) -> Result<SeriesStability, String> {
    if series.len() < 8 {
        return Err("at least eight samples are needed".into());
    }
    let step = series[1].0 - series[0].0;
    if !(step.is_finite() && step > 0.0) {
        return Err("samples must increase in time".into());
    }
    for w in series.windows(2) {
        if ((w[1].0 - w[0].0) - step).abs() > 1e-6 * step {
            return Err("samples must be on a regular grid".into());
        }
    }
    let n = series.len() as f64;
    let tm = series.iter().map(|p| p.0).sum::<f64>() / n;
    let xm = series.iter().map(|p| p.1).sum::<f64>() / n;
    let sxx = series.iter().map(|p| (p.0 - tm).ppowi(2)).sum::<f64>();
    let sxy = series.iter().map(|p| (p.0 - tm) * (p.1 - xm)).sum::<f64>();
    let rate = sxy / sxx;
    let resid: Vec<f64> = series
        .iter()
        .map(|p| p.1 - xm - rate * (p.0 - tm))
        .collect();
    let rms = (resid.iter().map(|r| r * r).sum::<f64>() / n).sqrt();
    let mut adev = Vec::new();
    let mut m = 1usize;
    while series.len() > 2 * m {
        let a = overlapping_adev(&resid, step, m);
        if a.is_finite() {
            adev.push((m as f64 * step, a));
        }
        m *= 2;
    }
    Ok(SeriesStability {
        n: series.len(),
        step_s: step,
        rate,
        residual_rms_s: rms,
        tdev: tdev_curve(&resid, step),
        adev,
    })
}

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

    fn link() -> (Station, Station, Relay) {
        let a = Station {
            r_m: [4_075_580.0, 931_853.0, 4_801_568.0], // central Europe
            tx_delay_s: 412.3e-9,
            rx_delay_s: 398.7e-9,
            uplink_hz: 14.0e9,
            downlink_hz: 11.0e9,
            slant_tec_tecu: 35.0,
            tropo_s: 8.1e-9,
        };
        let b = Station {
            r_m: [1_130_000.0, -4_831_000.0, 3_994_000.0], // eastern North America
            tx_delay_s: 501.9e-9,
            rx_delay_s: 455.2e-9,
            uplink_hz: 14.0e9,
            downlink_hz: 11.0e9,
            slant_tec_tecu: 22.0,
            tropo_s: 9.4e-9,
        };
        let s = Relay {
            r_m: [33_000_000.0, -24_000_000.0, 0.0],
            delay_ab_s: 812.0e-9,
            delay_ba_s: 809.5e-9,
        };
        (a, b, s)
    }

    #[test]
    fn the_two_way_equation_recovers_the_offset_exactly() {
        let (a, b, s) = link();
        for off in [-3.7e-6, 0.0, 12.5e-9, 1.0e-3] {
            let (ta, tb) = counter_readings(&a, &b, &s, off);
            let terms = solve_offset(&a, &b, &s, ta, tb);
            assert!(
                (terms.offset() - off).abs() < 1e-15,
                "{} vs {off}",
                terms.offset()
            );
        }
    }

    #[test]
    fn the_sagnac_term_is_the_bipm_closed_form() {
        let (a, b, s) = link();
        let terms = solve_offset(&a, &b, &s, 0.0, 0.0);
        // −½[S(A→S→B) − S(B→S→A)] = −S(A→S→B) by antisymmetry, which is minus the loop
        // A→S→B→A less its B→A leg; compare with the module's loop closed form.
        let loop_ab = crate::timetransfer_adv::twstft_sagnac(a.r_m, s.r_m, b.r_m)
            - sagnac_correction(b.r_m, a.r_m);
        assert!((terms.sagnac + loop_ab).abs() < 1e-18);
    }

    #[test]
    fn the_ionospheric_term_follows_the_frequency_split() {
        let (mut a, mut b, s) = link();
        a.slant_tec_tecu = 50.0;
        b.slant_tec_tecu = 0.0;
        let t = solve_offset(&a, &b, &s, 0.0, 0.0);
        let want =
            -0.5 * 40.3 * 50e16 * (1.0 / 14e9f64.ppowi(2) - 1.0 / 11e9f64.ppowi(2)) / C_M_PER_S;
        assert!((t.ionosphere - want).abs() < 1e-20);
        // About 0.1 ns for 50 TECU at 14/11 GHz.
        assert!(
            t.ionosphere > 0.5e-10 && t.ionosphere < 2e-10,
            "{}",
            t.ionosphere
        );
        // The troposphere cancels.
        b.tropo_s += 5e-9;
        a.tropo_s -= 3e-9;
        let t2 = solve_offset(&a, &b, &s, 0.0, 0.0);
        assert_eq!(t.ionosphere, t2.ionosphere);
    }

    #[test]
    fn flicker_phase_has_a_flat_tdev_and_white_falls() {
        let n = 2048;
        let fl = flicker_phase(n, 1e-10, 7);
        let td = tdev_curve(&fl, 1.0);
        // Flicker PM: TDEV is nearly flat; white PM falls as τ^-1/2.
        let first = td[0].1;
        let late = td[5].1;
        assert!(late / first > 0.5 && late / first < 2.0, "{first} {late}");
        let mut rng = ChaCha8Rng::seed_from_u64(9);
        let wh: Vec<f64> = (0..n)
            .map(|_| {
                1e-10 * {
                    let z: f64 = StandardNormal.sample(&mut rng);
                    z
                }
            })
            .collect();
        let tw = tdev_curve(&wh, 1.0);
        assert!(tw[5].1 / tw[0].1 < 0.3);
    }

    #[test]
    fn diurnal_non_reciprocity_lifts_tdev_near_half_a_day() {
        let quiet = LinkNoise {
            white_s: 100e-12,
            flicker_s: 0.0,
            diurnal_amp_s: 0.0,
            reciprocity_bias_s: 0.0,
        };
        let warm = LinkNoise {
            diurnal_amp_s: 500e-12,
            ..quiet
        };
        let step = 3600.0;
        let n = 24 * 20;
        let a = tdev_curve(&two_way_series(0.0, &quiet, step, n, 1), step);
        let b = tdev_curve(&two_way_series(0.0, &warm, step, n, 1), step);
        let at = |c: &[(f64, f64)], tau: f64| c.iter().find(|p| p.0 == tau).unwrap().1;
        assert!(at(&b, 8.0 * 3600.0) > 3.0 * at(&a, 8.0 * 3600.0));
    }

    #[test]
    fn a_reciprocity_bias_shifts_the_offset_but_not_the_tdev() {
        let n = LinkNoise {
            white_s: 100e-12,
            flicker_s: 0.0,
            diurnal_amp_s: 0.0,
            reciprocity_bias_s: 0.0,
        };
        let biased = LinkNoise {
            reciprocity_bias_s: 2e-9,
            ..n
        };
        let x = two_way_series(0.0, &n, 300.0, 512, 3);
        let y = two_way_series(0.0, &biased, 300.0, 512, 3);
        for (a, b) in x.iter().zip(&y) {
            assert!((b - a - 2e-9).abs() < 1e-18);
        }
        for (p, q) in tdev_curve(&x, 300.0).iter().zip(&tdev_curve(&y, 300.0)) {
            assert_eq!(p.0, q.0);
            assert!((p.1 - q.1).abs() <= 1e-12 * p.1, "{p:?} {q:?}");
        }
    }

    #[test]
    fn series_stability_removes_a_rate_and_reports_white_phase_noise() {
        let mut rng = ChaCha8Rng::seed_from_u64(4);
        let s: Vec<(f64, f64)> = (0..512)
            .map(|i| {
                let t = i as f64 * 432_000.0;
                (
                    t,
                    3e-9 + 1e-15 * t
                        + 1e-9 * {
                            let z: f64 = StandardNormal.sample(&mut rng);
                            z
                        },
                )
            })
            .collect();
        let st = series_stability(&s).unwrap();
        assert!((st.rate - 1e-15).abs() < 2e-17, "{}", st.rate);
        assert!((st.residual_rms_s - 1e-9).abs() < 1.5e-10);
        // White PM: TDEV at m = 1 is σ/√3... a looser check: it is of the order of σ.
        assert!(st.tdev[0].1 > 3e-10 && st.tdev[0].1 < 1.2e-9);
        assert!(series_stability(&s[..5]).is_err());
    }

    #[test]
    fn budget_combines_in_quadrature() {
        let l = vec![
            BudgetLine {
                name: "a".into(),
                u_s: 3e-9,
                kind: "B",
            },
            BudgetLine {
                name: "b".into(),
                u_s: 4e-9,
                kind: "A",
            },
        ];
        assert!((combined_uncertainty(&l) - 5e-9).abs() < 1e-21);
    }
}