kshana 0.11.0

Open, reproducible PNT-resilience simulator with quantum-sensor performance models
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89

// SPDX-License-Identifier: Apache-2.0
//! Geometric time-transfer corrections: the Sagnac effect and GNSS common-view.
//!
//! [`crate::timetransfer`] models the *stochastic* error of a two-way link (jitter,
//! flicker, random walk). This module adds the two *deterministic* geometric effects a
//! real comparison must account for:
//!
//! - **Sagnac correction** — because the signal propagates over a rotating Earth, a path
//!   with an east–west component accrues an extra delay `Δt = (ω_E/c²)·(x₁y₂ − x₂y₁)`
//!   (equivalently `2·A·ω_E/c²` for the equatorial-projected area `A` of the triangle
//!   Earth-centre–`r₁`–`r₂`). It is tens of nanoseconds for continental baselines and must
//!   be removed before two clocks can be compared.
//! - **GNSS common-view** — two ground stations observing the *same* satellite at the same
//!   instant difference their range-corrected pseudoranges; the satellite-clock error (and
//!   other common-mode errors) cancels exactly, leaving the inter-station clock offset.
//!
//! Both are exact closed forms. Scope (honest): a full TWSTFT transponder/hardware-delay
//! budget and a PPP ionosphere-free time-transfer solution are follow-ons (see `ROADMAP.md`).

/// Earth rotation rate (rad/s), IERS.
pub const OMEGA_EARTH: f64 = 7.292_115_9e-5;
/// Speed of light (m/s).
pub const C_M_PER_S: f64 = 299_792_458.0;

type Vec3 = [f64; 3];

/// Sagnac time correction (s) for a signal propagating from `r1` to `r2` (ECEF metres):
/// `Δt = (ω_E/c²)·(x₁·y₂ − x₂·y₁)`. Positive for an eastward-projected path; zero for a
/// purely radial or polar path. Antisymmetric in the path direction.
pub fn sagnac_correction(r1: Vec3, r2: Vec3) -> f64 {
    OMEGA_EARTH * (r1[0] * r2[1] - r2[0] * r1[1]) / (C_M_PER_S * C_M_PER_S)
}

/// Common-view single-difference estimate of the inter-station clock offset
/// `(clock_A − clock_B)` (s), from each station's pseudorange to the *same* satellite with
/// its computed geometric range removed: `[(ρ_A − R_A) − (ρ_B − R_B)] / c`. The satellite
/// clock error — common to both observations — cancels exactly.
pub fn common_view_offset(pr_a: f64, range_a: f64, pr_b: f64, range_b: f64) -> f64 {
    ((pr_a - range_a) - (pr_b - range_b)) / C_M_PER_S
}

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

    const RE: f64 = 6_378_137.0;

    #[test]
    fn sagnac_equatorial_quarter_turn_is_about_33_ns() {
        // r1 on the +x equator, r2 on the +y equator (90° apart): x1·y2 − x2·y1 = Re².
        // Δt = ω_E·Re²/c² ≈ 33.0 ns.
        let dt = sagnac_correction([RE, 0.0, 0.0], [0.0, RE, 0.0]);
        assert!((dt - 3.3006e-8).abs() < 1e-11, "Δt = {} s", dt);
    }

    #[test]
    fn sagnac_is_antisymmetric_and_zero_for_radial_or_polar_paths() {
        let a = [RE, 1.0e6, 2.0e6];
        let b = [-2.0e6, RE, 5.0e5];
        assert!((sagnac_correction(a, b) + sagnac_correction(b, a)).abs() < 1e-20);
        // Purely radial path (same x,y direction): no east–west component ⇒ 0.
        assert_eq!(sagnac_correction([RE, 0.0, 0.0], [2.0 * RE, 0.0, 0.0]), 0.0);
        // Polar path (x = y = 0): 0.
        assert_eq!(sagnac_correction([0.0, 0.0, RE], [0.0, 0.0, 2.0 * RE]), 0.0);
    }

    #[test]
    fn common_view_cancels_the_satellite_clock_and_recovers_the_offset() {
        // True inter-station offset to recover.
        let dt_a = 1.0e-7; // station A clock (s)
        let dt_b = -3.0e-8; // station B clock (s)
        let range_a = 2.10e7;
        let range_b = 2.13e7;
        // Pseudorange = geometric range + c·(rx_clock − sv_clock).
        let pr = |range: f64, rx: f64, sv: f64| range + C_M_PER_S * (rx - sv);
        // The estimate must equal (dt_a − dt_b) regardless of the satellite clock error.
        for &sv in &[0.0, 1.0e-6, -5.0e-7] {
            let est = common_view_offset(
                pr(range_a, dt_a, sv),
                range_a,
                pr(range_b, dt_b, sv),
                range_b,
            );
            // Femtosecond tolerance: the residual is pure f64 roundoff of c·dt at the
            // ~2×10⁷ m pseudorange magnitude, far below any real clock measurement.
            assert!((est - (dt_a - dt_b)).abs() < 1e-15, "sv {sv}: est {est}");
        }
    }
}