lox_orbits/propagators/
sgp4.rs

1// SPDX-FileCopyrightText: 2024 Helge Eichhorn <git@helgeeichhorn.de>
2//
3// SPDX-License-Identifier: MPL-2.0
4
5use glam::DVec3;
6pub use sgp4::Elements;
7use sgp4::{Constants, ElementsError, MinutesSinceEpoch};
8use thiserror::Error;
9
10use lox_bodies::Earth;
11use lox_core::f64::consts::SECONDS_PER_MINUTE;
12use lox_time::Time;
13use lox_time::deltas::TimeDelta;
14use lox_time::time_scales::Tai;
15use lox_time::utc::{Utc, UtcError};
16
17use crate::propagators::Propagator;
18use crate::states::State;
19use crate::trajectories::TrajectoryError;
20use lox_frames::Icrf;
21
22#[derive(Debug, Clone, Error)]
23pub enum Sgp4Error {
24    #[error(transparent)]
25    ElementsError(#[from] ElementsError),
26    #[error(transparent)]
27    TrajectoryError(#[from] TrajectoryError),
28    #[error(transparent)]
29    Sgp4(#[from] sgp4::Error),
30    #[error(transparent)]
31    Utc(#[from] UtcError),
32}
33
34pub struct Sgp4 {
35    constants: Constants,
36    time: Time<Tai>,
37}
38
39impl Sgp4 {
40    pub fn new(initial_state: Elements) -> Result<Self, Sgp4Error> {
41        let epoch = initial_state.epoch();
42        let time = Utc::from_delta(TimeDelta::from_julian_years(epoch))?.to_time();
43        let constants = Constants::from_elements(&initial_state)?;
44        Ok(Self { constants, time })
45    }
46
47    pub fn time(&self) -> Time<Tai> {
48        self.time
49    }
50}
51
52impl Propagator<Tai, Earth, Icrf> for Sgp4 {
53    type Error = Sgp4Error;
54
55    fn propagate(&self, time: Time<Tai>) -> Result<State<Tai, Earth, Icrf>, Self::Error> {
56        let dt = (time - self.time).as_seconds_f64() / SECONDS_PER_MINUTE;
57        let prediction = self.constants.propagate(MinutesSinceEpoch(dt))?;
58        let position = DVec3::from_array(prediction.position);
59        let velocity = DVec3::from_array(prediction.velocity);
60        Ok(State::new(time, position, velocity, Earth, Icrf))
61    }
62}
63
64#[cfg(test)]
65mod tests {
66    use lox_test_utils::assert_approx_eq;
67
68    use super::*;
69
70    #[test]
71    fn test_sgp4() {
72        let tle = Elements::from_tle(
73            Some("ISS (ZARYA)".to_string()),
74            "1 25544U 98067A   24170.37528350  .00016566  00000+0  30244-3 0  9996".as_bytes(),
75            "2 25544  51.6410 309.3890 0010444 339.5369 107.8830 15.49495945458731".as_bytes(),
76        )
77        .unwrap();
78        let sgp4 = Sgp4::new(tle).unwrap();
79        let orbital_period = 92.821;
80        let t1 = sgp4.time() + TimeDelta::from_minutes(orbital_period);
81        let s1 = sgp4.propagate(t1).unwrap();
82        let k1 = s1.to_keplerian();
83        assert_approx_eq!(
84            k1.orbital_period().as_seconds_f64() / SECONDS_PER_MINUTE,
85            orbital_period,
86            rtol <= 1e-4
87        );
88    }
89}