nyx_space/od/msr/
two_way.rs

1/*
2    Nyx, blazing fast astrodynamics
3    Copyright (C) 2018-onwards Christopher Rabotin <christopher.rabotin@gmail.com>
4
5    This program is free software: you can redistribute it and/or modify
6    it under the terms of the GNU Affero General Public License as published
7    by the Free Software Foundation, either version 3 of the License, or
8    (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU Affero General Public License for more details.
14
15    You should have received a copy of the GNU Affero General Public License
16    along with this program.  If not, see <https://www.gnu.org/licenses/>.
17*/
18
19use anise::constants::SPEED_OF_LIGHT_KM_S;
20use anise::constants::celestial_objects::SUN;
21use anise::constants::frames::SUN_J2000;
22use anise::constants::orientations::ICRS;
23use anise::errors::AlmanacPhysicsSnafu;
24use anise::frames::Frame;
25use anise::prelude::Almanac;
26use hifitime::{Duration, Epoch, TimeUnits};
27use nalgebra::Vector3;
28
29use crate::Spacecraft;
30use crate::md::prelude::Traj;
31use crate::od::ground_station::GroundStation;
32use crate::od::{ODAlmanacSnafu, ODError, ODPlanetaryDataSnafu, ODTrajSnafu};
33use snafu::ResultExt;
34
35const SUN_ICRS: Frame = Frame::new(SUN, ICRS);
36
37#[derive(Copy, Clone, Debug, PartialEq)]
38pub struct LightTimeLeg {
39    pub start_epoch: Epoch,
40    pub end_epoch: Epoch,
41    pub r_start_icrf_km: Vector3<f64>,
42    pub r_end_icrf_km: Vector3<f64>,
43    pub light_time: Duration,
44}
45
46impl LightTimeLeg {
47    pub fn range_km(&self) -> f64 {
48        (self.r_start_icrf_km - self.r_end_icrf_km).norm()
49    }
50}
51
52/// Solves the two-way tracking problem using the Picard fixed-point iteration like in JPL Moyer (2000).
53///
54/// This solves the problem strictly in ICRF because it is the only true inertial frame.
55///
56/// ```text
57/// t1 (Transmit)                 t2 (Bounce)                 t3 (Receive)
58/// Ground Station -------------> Spacecraft -------------> Ground Station
59///               [Uplink Leg]                [Downlink Leg]
60///               (Solved Second)              (Solved First)
61/// ```
62#[derive(Copy, Clone, Debug, PartialEq)]
63pub struct TwoWaySolution {
64    pub t1_transmit: Epoch,
65    pub t2_bounce: Epoch,
66    pub t3_receive: Epoch,
67    pub uplink: LightTimeLeg,
68    pub downlink: LightTimeLeg,
69}
70
71impl TwoWaySolution {
72    pub fn total_delay(&self) -> Duration {
73        self.uplink.light_time + self.downlink.light_time
74    }
75
76    pub fn range_km(&self) -> f64 {
77        0.5 * (self.t3_receive - self.t1_transmit).to_seconds() * SPEED_OF_LIGHT_KM_S
78        // (self.uplink.range_km() + self.downlink.range_km()) / 2.0
79    }
80}
81
82/// Computes the exact two-way light-time solution using Picard fixed-point iteration.
83/// IMPORTANT This does not check for elevation constraints. You must manually check that
84/// at the t3_bounce epoch.
85pub(crate) fn solve_two_way_picard(
86    t3: Epoch,
87    station: &GroundStation,
88    traj: &Traj<Spacecraft>,
89    almanac: &Almanac,
90) -> Result<TwoWaySolution, ODError> {
91    let sun_mu_km3_s2 = if station.relativistic_corrections {
92        Some(
93            almanac
94                .frame_info(SUN_J2000)
95                .context(ODPlanetaryDataSnafu {
96                    action: "fetching Sun grav param for Shapiro delay",
97                })?
98                .mu_km3_s2()
99                .context(AlmanacPhysicsSnafu {
100                    action: "Sun mu not defined",
101                })
102                .context(ODAlmanacSnafu {
103                    action: "fetching Sun grav param for Shapiro delay",
104                })?,
105        )
106    } else {
107        None
108    };
109    // Step 0: Anchor the reception state of the ground station in ICRF
110    let gs_rx_orbit = station.to_orbit(t3, almanac).context(ODAlmanacSnafu {
111        action: "building ground station orbit at t3",
112    })?;
113    let gs_rx_icrf = almanac
114        .transform_to(gs_rx_orbit, SUN_ICRS, None)
115        .context(ODAlmanacSnafu {
116            action: "transforming station at t3 to ICRF",
117        })?;
118    let r3_icrf_km = gs_rx_icrf.radius_km;
119
120    // Solve Downlink Leg (t3 -> t2)
121    // Find t2 such that c * (t3 - t2) = || r_sc(t2) - r_gs(t3) ||
122    let mut tau_down = Duration::ZERO;
123    let mut r2_icrf_km = Vector3::zeros();
124    let mut t2 = t3;
125
126    // Exactly 3 iterations converge to sub-millimeter precision in ICRF
127    for _ in 0..3 {
128        t2 = t3 - tau_down;
129        // NOTE Using with_context to lazy eval the format on the error
130        let sc_state = traj.at(t2).with_context(|_| ODTrajSnafu {
131            details: format!("interpolating spacecraft state at bounce epoch {t2}"),
132        })?;
133
134        // Transform spacecraft orbit to ICRF
135        let sc_icrf = almanac
136            .transform_to(sc_state.orbit, SUN_ICRS, None)
137            .context(ODAlmanacSnafu {
138                action: "transforming spacecraft at t2 to ICRF",
139            })?;
140
141        r2_icrf_km = sc_icrf.radius_km;
142        let dist_down_km = (r2_icrf_km - r3_icrf_km).norm();
143        // Geometric transit time + Solar Shapiro time dilation
144        let tau_geometric = dist_down_km / SPEED_OF_LIGHT_KM_S;
145        let tau_shapiro = if let Some(sun_mu) = sun_mu_km3_s2 {
146            shapiro_delay_s(&r2_icrf_km, &r3_icrf_km, dist_down_km, sun_mu)
147        } else {
148            0.0
149        };
150        tau_down = (tau_geometric + tau_shapiro).seconds();
151    }
152
153    let downlink = LightTimeLeg {
154        start_epoch: t2,
155        end_epoch: t3,
156        r_start_icrf_km: r2_icrf_km,
157        r_end_icrf_km: r3_icrf_km,
158        light_time: tau_down,
159    };
160
161    // Solve Uplink Leg (t2 -> t1)
162    // Spacecraft state r_sc(t2) is now fixed.
163    // Find t1 such that c * (t2 - t1) = || r_sc(t2) - r_gs(t1) ||
164    let mut tau_up = tau_down; // Good initial guess
165    let mut r1_icrf_km = Vector3::zeros();
166    let mut t1 = t2 - tau_up;
167
168    for _ in 0..3 {
169        t1 = t2 - tau_up;
170        let gs_tx_orbit = station.to_orbit(t1, almanac).context(ODAlmanacSnafu {
171            action: "building ground station orbit at t1",
172        })?;
173
174        let gs_tx_icrf =
175            almanac
176                .transform_to(gs_tx_orbit, SUN_ICRS, None)
177                .context(ODAlmanacSnafu {
178                    action: "transforming station at t1 to ICRF",
179                })?;
180
181        r1_icrf_km = gs_tx_icrf.radius_km;
182        let dist_up_km = (r2_icrf_km - r1_icrf_km).norm();
183        let tau_geometric = dist_up_km / SPEED_OF_LIGHT_KM_S;
184        let tau_shapiro = if let Some(sun_mu) = sun_mu_km3_s2 {
185            shapiro_delay_s(&r1_icrf_km, &r2_icrf_km, dist_up_km, sun_mu)
186        } else {
187            0.0
188        };
189        tau_up = (tau_geometric + tau_shapiro).seconds();
190    }
191
192    let uplink = LightTimeLeg {
193        start_epoch: t1,
194        end_epoch: t2,
195        r_start_icrf_km: r1_icrf_km,
196        r_end_icrf_km: r2_icrf_km,
197        light_time: tau_up,
198    };
199
200    // Package Two-Way Observables
201    Ok(TwoWaySolution {
202        t1_transmit: t1,
203        t2_bounce: t2,
204        t3_receive: t3,
205        uplink,
206        downlink,
207    })
208}
209
210/// Shapiro delay is computed when relativistic corrections are enabled.
211/// Signals passing near a massive object take slightly longer to travel to a target and longer to return than they
212/// would if the mass of the object were not present. The time delay is caused by time dilation, which increases
213/// the time it takes light to travel a given distance from the perspective of an outside observer.
214fn shapiro_delay_s(
215    r_start_icrf_km: &Vector3<f64>,
216    r_end_icrf_km: &Vector3<f64>,
217    rho_km: f64,
218    sun_mu_km3_s2: f64,
219) -> f64 {
220    let num = r_start_icrf_km.norm() + r_end_icrf_km.norm() + rho_km;
221    let denom = r_start_icrf_km.norm() + r_end_icrf_km.norm() - rho_km;
222
223    if denom <= 0.0 || num <= 0.0 {
224        return 0.0;
225    }
226
227    // PPN formulation (gamma = 1): 2 * mu / c^3 * ln(...)
228    let factor = (2.0 * sun_mu_km3_s2) / SPEED_OF_LIGHT_KM_S.powi(3);
229    factor * (num / denom).ln()
230}
nyx_space/od/msr/two_way.rs

nyx_space/od/msr/
two_way.rs
