siderust 0.9.1

High-precision astronomy and satellite mechanics in Rust.
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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164

//! # RTN orbit-to-orbit comparison
//!
//! ## Scientific scope
//!
//! Comparing an estimated trajectory against a reference orbit is most
//! interpretable in the local radial-transverse-normal frame, where along-
//! track, radial, and cross-track errors separate naturally. This module
//! computes that frame-aware difference for matched epochs.
//!
//! It assumes both trajectories are already sampled at identical epochs and
//! expressed in compatible inertial coordinates. Interpolation and frame
//! reconciliation are intentionally left to caller code.
//!
//! ## Technical scope
//!
//! The public items are `RtnDiff`, `RtnSummary`, `rtn_diff`, and
//! `rtn_summary`. Callers pass slices of `OrbitState` values and receive
//! per-epoch RTN differences plus a condensed summary for QC reporting.
//!
//! This module does not estimate the orbit or ingest external product files
//! directly.
//!
//! ## References
//!
//! - Tapley, B. D., Schutz, B. E., & Born, G. H. (2004). Statistical Orbit
//!   Determination. Elsevier Academic Press.
//! - Vallado, D. A. (2013). Fundamentals of Astrodynamics and Applications
//!   (4th ed.). Microcosm Press.
use crate::astro::dynamics::frames::RTN;
use crate::astro::dynamics::OrbitState;
use affn::cartesian::Displacement;
use affn::frames::GCRS;
use qtty::unit::Kilometer;

/// Per-epoch RTN difference between an estimated and reference state.
#[derive(Debug, Clone)]
pub struct RtnDiff {
    /// Julian Date (TT) of the comparison epoch.
    pub jd_tt: f64,
    /// Radial difference (m).
    pub r_m: f64,
    /// Along-track difference (m).
    pub t_m: f64,
    /// Cross-track difference (m).
    pub n_m: f64,
}

impl RtnDiff {
    /// Position difference as a typed displacement in the RTN frame (kilometres).
    pub fn position_rtn_km(&self) -> Displacement<RTN, Kilometer> {
        Displacement::new(self.r_m * 1e-3, self.t_m * 1e-3, self.n_m * 1e-3)
    }
}

/// Aggregate RTN statistics over a comparison run.
#[derive(Debug, Clone)]
pub struct RtnSummary {
    /// Number of compared epochs.
    pub n: usize,
    /// 3-D RMS, metres.
    pub rms_3d_m: f64,
    /// RMS of radial component, metres.
    pub rms_r_m: f64,
    /// RMS of along-track, metres.
    pub rms_t_m: f64,
    /// RMS of cross-track, metres.
    pub rms_n_m: f64,
}

/// Pairwise RTN difference between estimated and reference states; epochs are
/// matched by index (caller is responsible for alignment and time-tag check).
pub fn rtn_diff(estimated: &[OrbitState], reference: &[OrbitState]) -> Vec<RtnDiff> {
    let n = estimated.len().min(reference.len());
    let mut out = Vec::with_capacity(n);
    for i in 0..n {
        let e = estimated[i];
        let r = reference[i];
        // Position difference in GCRS (typed Displacement<GCRS, Km>).
        let d_gcrs: Displacement<GCRS, Kilometer> = e.position - r.position;
        let frame = crate::astro::dynamics::frames::LocalOrbitalFrame::<RTN>::try_from_state(&r)
            .expect("RTN frame from reference state");
        let d_rtn = frame.to_local(d_gcrs);
        out.push(RtnDiff {
            jd_tt: e.epoch.to::<tempoch::JD>().value(),
            r_m: d_rtn.x().value() * 1000.0,
            t_m: d_rtn.y().value() * 1000.0,
            n_m: d_rtn.z().value() * 1000.0,
        });
    }
    out
}

/// Aggregate RMS statistics over a list of RTN differences.
pub fn rtn_summary(diffs: &[RtnDiff]) -> RtnSummary {
    let n = diffs.len();
    if n == 0 {
        return RtnSummary {
            n: 0,
            rms_3d_m: 0.0,
            rms_r_m: 0.0,
            rms_t_m: 0.0,
            rms_n_m: 0.0,
        };
    }
    let mut sr = 0.0;
    let mut st = 0.0;
    let mut sn = 0.0;
    for d in diffs {
        sr += d.r_m * d.r_m;
        st += d.t_m * d.t_m;
        sn += d.n_m * d.n_m;
    }
    let rms_r = (sr / n as f64).sqrt();
    let rms_t = (st / n as f64).sqrt();
    let rms_n = (sn / n as f64).sqrt();
    let rms_3d = (rms_r * rms_r + rms_t * rms_t + rms_n * rms_n).sqrt();
    RtnSummary {
        n,
        rms_3d_m: rms_3d,
        rms_r_m: rms_r,
        rms_t_m: rms_t,
        rms_n_m: rms_n,
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::astro::dynamics::{Position, Velocity};
    use crate::time::JulianDate;

    #[test]
    fn rtn_zero_for_identical_orbits() {
        let s = OrbitState::new(
            JulianDate::new(2_451_545.0).to_j2000s(),
            Position::new(7000.0, 0.0, 0.0),
            Velocity::new(0.0, 7.5, 0.0),
        );
        let d = rtn_diff(&[s], &[s]);
        assert_eq!(d.len(), 1);
        assert!(d[0].r_m.abs() < 1e-9);
        let summary = rtn_summary(&d);
        assert!(summary.rms_3d_m < 1e-9);
        assert_eq!(summary.n, 1);
    }

    #[test]
    fn radial_offset_only_appears_in_r_component() {
        let r = OrbitState::new(
            JulianDate::new(2_451_545.0).to_j2000s(),
            Position::new(7000.0, 0.0, 0.0),
            Velocity::new(0.0, 7.5, 0.0),
        );
        let e = OrbitState::new(
            JulianDate::new(2_451_545.0).to_j2000s(),
            Position::new(7000.001, 0.0, 0.0),
            Velocity::new(0.0, 7.5, 0.0),
        );
        let d = rtn_diff(&[e], &[r]);
        assert!((d[0].r_m - 1.0).abs() < 1e-9);
        assert!(d[0].t_m.abs() < 1e-9);
        assert!(d[0].n_m.abs() < 1e-9);
    }
}