sgp4_rs/
lib.rs

1use chrono::prelude::*;
2use chrono::DateTime;
3use thiserror::Error;
4use uom::si::{angle, angular_velocity::radian_per_second, f64::*, length::kilometer};
5
6mod sgp4_sys;
7#[cfg(feature = "tlegen")]
8mod tlegen;
9
10#[derive(Debug, Error, PartialEq)]
11pub enum Error {
12    #[error("TLE was malformed: {0}")]
13    MalformedTwoLineElement(String),
14    #[error("{0}")]
15    UnknownError(String),
16    #[error(transparent)]
17    PropagationError(#[from] sgp4_sys::Error),
18    #[error("Optimization error: {0}")]
19    OptimizationError(String),
20}
21
22type Result<T> = std::result::Result<T, Error>;
23
24/// A state vector, in the TEME-ECI coordinate frame, for an orbiting body.
25///
26/// To obtain the state of an object at a specific time, use the propagation functions provided by
27/// [TwoLineElement].
28#[derive(Debug, Clone, Copy)]
29pub struct StateVector {
30    pub epoch: DateTime<Utc>,
31
32    /// The satellite position in km.
33    pub position: [f64; 3],
34
35    /// The satellite velocity in km/s.
36    pub velocity: [f64; 3],
37
38    pub coe: ClassicalOrbitalElements,
39}
40
41impl StateVector {
42    pub fn new(epoch: DateTime<Utc>, position: [f64; 3], velocity: [f64; 3]) -> Self {
43        Self {
44            epoch,
45            position,
46            velocity,
47            coe: sgp4_sys::to_classical_elements(&position, &velocity).into(),
48        }
49    }
50
51    pub fn semilatus_rectum(&self) -> Length {
52        self.coe.semilatus_rectum
53    }
54
55    pub fn semimajor_axis(&self) -> Length {
56        self.coe.semimajor_axis
57    }
58
59    pub fn inclination(&self) -> Angle {
60        self.coe.inclination
61    }
62
63    pub fn raan(&self) -> Angle {
64        self.coe.raan
65    }
66
67    pub fn mean_anomaly(&self) -> Angle {
68        self.coe.mean_anomaly
69    }
70
71    pub fn true_anomaly(&self) -> Angle {
72        self.coe.true_anomaly
73    }
74
75    pub fn eccentricity(&self) -> f64 {
76        self.coe.eccentricity
77    }
78
79    pub fn longitude_of_periapsis(&self) -> Angle {
80        self.coe.longitude_of_periapsis
81    }
82
83    pub fn true_longitude(&self) -> Angle {
84        self.coe.true_longitude
85    }
86
87    pub fn argument_of_perigee(&self) -> Angle {
88        self.coe.argument_of_perigee
89    }
90
91    pub fn argument_of_latitude(&self) -> Angle {
92        self.coe.argument_of_latitude
93    }
94}
95
96/// A Keplerian orbital element set.
97///
98/// This structure contains all of the "classical" orbital elements as derivable from a TLE. We
99/// lean on the `uom` crate to provide safe dimensional types which help to avoid bugs related to
100/// mixing units of measure.
101#[derive(Debug, Clone, Copy, PartialEq)]
102pub struct ClassicalOrbitalElements {
103    pub semilatus_rectum: Length,
104    pub semimajor_axis: Length,
105    pub eccentricity: f64,
106    pub inclination: Angle,
107    pub raan: Angle,
108    pub argument_of_perigee: Angle,
109    pub true_anomaly: Angle,
110    pub mean_anomaly: Angle,
111    pub argument_of_latitude: Angle,
112    pub true_longitude: Angle,
113    pub longitude_of_periapsis: Angle,
114}
115
116impl From<sgp4_sys::ClassicalOrbitalElements> for ClassicalOrbitalElements {
117    fn from(coe: sgp4_sys::ClassicalOrbitalElements) -> Self {
118        let semilatus_rectum = Length::new::<kilometer>(coe.p);
119        let semimajor_axis = Length::new::<kilometer>(coe.a);
120        let inclination = Angle::new::<angle::radian>(coe.incl);
121        let raan = Angle::new::<angle::radian>(coe.omega);
122        let mean_anomaly = Angle::new::<angle::radian>(coe.m);
123        let true_anomaly = Angle::new::<angle::radian>(coe.nu);
124        let eccentricity = coe.ecc;
125        let longitude_of_periapsis = Angle::new::<angle::radian>(coe.lonper);
126        let true_longitude = Angle::new::<angle::radian>(coe.truelon);
127        let argument_of_perigee = Angle::new::<angle::radian>(coe.argp);
128        let argument_of_latitude = Angle::new::<angle::radian>(coe.arglat);
129
130        Self {
131            semilatus_rectum,
132            semimajor_axis,
133            eccentricity,
134            inclination,
135            raan,
136            argument_of_perigee,
137            true_anomaly,
138            mean_anomaly,
139            argument_of_latitude,
140            true_longitude,
141            longitude_of_periapsis,
142        }
143    }
144}
145
146impl From<StateVector> for ClassicalOrbitalElements {
147    fn from(sv: StateVector) -> Self {
148        sv.coe
149    }
150}
151
152const TLE_LINE_LENGTH: usize = 69;
153
154/// A parsed, valid Two Line Element data set which can be used for orbital propagation.
155///
156/// Internally this uses SGP4's own structure representation. Various fields which are useful for
157/// analysis or simulation are exposed via methods like `raan()`/`set_raan()`
158/// which allow access to the values and direct modification of the underlying orbital element set
159/// in a type-safe manner. The `uom` crate provides dimensional analysis to help avoid
160/// unit-of-measure errors which can otherwise be quite difficult to detect.
161#[derive(Clone)]
162pub struct TwoLineElement {
163    elements: sgp4_sys::OrbitalElementSet,
164}
165
166impl TwoLineElement {
167    /// Create a validated TwoLineElement from a string.
168    pub fn new(line1: &str, line2: &str) -> Result<TwoLineElement> {
169        let line1 = line1.trim();
170        let line2 = line2.trim();
171
172        if line1.len() != TLE_LINE_LENGTH {
173            return Err(Error::MalformedTwoLineElement(format!(
174                "Line 1 is the wrong length. Expected {}, but got {}\n{}",
175                TLE_LINE_LENGTH,
176                line1.len(),
177                line1
178            )));
179        }
180
181        if line2.len() != TLE_LINE_LENGTH {
182            return Err(Error::MalformedTwoLineElement(format!(
183                "Line 2 is the wrong length. Expected {}, but got {}\n{}",
184                TLE_LINE_LENGTH,
185                line2.len(),
186                line2
187            )));
188        }
189
190        let elements = sgp4_sys::to_orbital_elements(
191            line1,
192            line2,
193            sgp4_sys::RunType::Verification,
194            sgp4_sys::OperationMode::Improved,
195            sgp4_sys::GravitationalConstant::Wgs84,
196        )
197        .map_err(|e| Error::MalformedTwoLineElement(e.to_string()))?;
198
199        Ok(TwoLineElement { elements })
200    }
201
202    /// Create a TwoLineElement from a string containing both lines, and optionally a header line.
203    pub fn from_lines(combined_lines: &str) -> Result<TwoLineElement> {
204        let lines: Vec<_> = {
205            let mut ls: Vec<_> = combined_lines
206                .split('\n')
207                .filter(|s| !s.is_empty())
208                .collect();
209            if ls.len() == 3 {
210                ls.split_off(1)
211            } else if ls.len() == 2 {
212                ls
213            } else {
214                return Err(Error::MalformedTwoLineElement(format!(
215                    "Expected two lines, got {}",
216                    ls.len()
217                )));
218            }
219        };
220        TwoLineElement::new(lines[0], lines[1])
221    }
222
223    /// Get the epoch of a TwoLineElement.
224    pub fn epoch(&self) -> Result<DateTime<Utc>> {
225        Ok(self.elements.epoch())
226    }
227
228    pub fn mean_motion(&self) -> AngularVelocity {
229        AngularVelocity::new::<radian_per_second>(self.elements.mean_motion() / 60.)
230    }
231
232    /// Propagate a TwoLineElement to the given time to obtain a state vector for the object.
233    pub fn propagate_to(&self, t: DateTime<Utc>) -> Result<StateVector> {
234        let tle_epoch = self.elements.epoch();
235        // TODO: determine correct behaviour for negative prop
236        // assert!(t >= tle_epoch);
237
238        let min_since_epoch = (t - tle_epoch).num_milliseconds() as f64 / 60_000.;
239
240        let (r, v) = sgp4_sys::run_sgp4(
241            self.elements,
242            sgp4_sys::GravitationalConstant::Wgs84,
243            min_since_epoch,
244        )?;
245
246        Ok(StateVector::new(t, r.to_owned(), v.to_owned()))
247    }
248}
249
250/// Wrapper type representing a Julian day.
251///
252/// This is the number of days since the start of the Julian astronomical calendar in 4713 BC, used
253/// to provide a consistent time reference for astronomical calculations.
254pub struct JulianDay(f64);
255
256impl From<DateTime<Utc>> for JulianDay {
257    fn from(d: DateTime<Utc>) -> Self {
258        JulianDay(sgp4_sys::datetime_to_julian_day(d))
259    }
260}
261
262impl From<JulianDay> for DateTime<Utc> {
263    fn from(jd: JulianDay) -> Self {
264        sgp4_sys::julian_day_to_datetime(jd.0)
265    }
266}
267
268/// Wrapper type representing the angular form of Greenwich Mean Sidereal Time.
269///
270/// This is primarily used to account for the Earth's rotation during conversion between fixed and
271/// inertial coordinate frames. Note that this is an angle measured in radians, and not a "time" as
272/// such. The value may range from 0 to 2π.
273pub struct GreenwichMeanSiderealTime(f64);
274
275impl GreenwichMeanSiderealTime {
276    pub fn as_radians(&self) -> f64 {
277        self.0
278    }
279}
280
281impl From<DateTime<Utc>> for GreenwichMeanSiderealTime {
282    fn from(d: DateTime<Utc>) -> Self {
283        GreenwichMeanSiderealTime(sgp4_sys::datetime_to_gstime(d))
284    }
285}
286
287#[cfg(test)]
288mod tests {
289    use super::*;
290
291    use chrono::Duration;
292    use float_cmp::approx_eq;
293
294    fn vecs_eq(l: &[f64; 3], r: &[f64; 3]) -> bool {
295        approx_eq!(f64, l[0], r[0]) && approx_eq!(f64, l[1], r[1]) && approx_eq!(f64, l[2], r[2])
296    }
297
298    #[test]
299    fn test_simple_propagation() -> Result<()> {
300        let line1 = "1 25544U 98067A   20148.21301450  .00001715  00000-0  38778-4 0  9992";
301        let line2 = "2 25544  51.6435  92.2789 0002570 358.0648 144.9972 15.49396855228767";
302
303        let tle = TwoLineElement::new(line1, line2)?;
304        let epoch = tle.epoch()?;
305
306        let s1 = tle.propagate_to(epoch)?;
307        let s2 = tle.propagate_to(epoch + Duration::hours(1))?;
308
309        assert!(!vecs_eq(&s1.position, &s2.position));
310        assert!(!vecs_eq(&s1.velocity, &s2.velocity));
311
312        Ok(())
313    }
314
315    #[test]
316    fn test_decay_error() -> Result<()> {
317        let line1 = "1 43051U 17071Q   22046.92182028  .07161566  12340-4  74927-3 0  9993";
318        let line2 = "2 43051  51.6207 236.5853 0009084 284.2762  75.7254 16.36736354237455";
319        let tle = TwoLineElement::new(line1, line2)?;
320        let epoch = tle.epoch()?;
321
322        let s2 = tle.propagate_to(epoch + Duration::days(30));
323        assert!(s2.is_err());
324        assert_eq!(
325            s2.unwrap_err(),
326            Error::PropagationError(sgp4_sys::Error::SatelliteDecay)
327        );
328        Ok(())
329    }
330
331    #[test]
332    fn mean_motion_round_trip() -> Result<()> {
333        let line1 = "1 25544U 98067A   20148.21301450  .00001715  00000-0  38778-4 0  9992";
334        let line2 = "2 25544  51.6435  92.2789 0002570 358.0648 144.9972 15.49396855228767";
335        let tle = TwoLineElement::new(line1, line2)?;
336        let mean_motion = tle.mean_motion().get::<radian_per_second>() * 24. * 60.0 * 60.0
337            / (2. * std::f64::consts::PI);
338        assert!(approx_eq!(f64, mean_motion, 15.493968, epsilon = 0.01));
339        Ok(())
340    }
341
342    #[test]
343    fn test_negative_time_propagation() -> Result<()> {
344        let line1 = "1 25544U 98067A   20148.21301450  .00001715  00000-0  38778-4 0  9992";
345        let line2 = "2 25544  51.6435  92.2789 0002570 358.0648 144.9972 15.49396855228767";
346
347        let tle = TwoLineElement::new(line1, line2)?;
348        let epoch = tle.epoch()?;
349
350        let s1 = tle.propagate_to(epoch)?;
351        let s2 = tle.propagate_to(epoch - Duration::days(30))?;
352
353        assert!(!vecs_eq(&s1.position, &s2.position));
354        assert!(!vecs_eq(&s1.velocity, &s2.velocity));
355
356        Ok(())
357    }
358
359    #[test]
360    fn test_tle_from_lines() -> Result<()> {
361        let lines = "1 25544U 98067A   20148.21301450  .00001715  00000-0  38778-4 0  9992
362                     2 25544  51.6435  92.2789 0002570 358.0648 144.9972 15.49396855228767";
363
364        let _tle = TwoLineElement::from_lines(lines)?;
365        Ok(())
366    }
367
368    #[test]
369    fn test_tle_from_lines_with_header() -> Result<()> {
370        let lines = "ISS (ZARYA)
371                     1 25544U 98067A   20148.21301450  .00001715  00000-0  38778-4 0  9992
372                     2 25544  51.6435  92.2789 0002570 358.0648 144.9972 15.49396855228767";
373
374        let _tle = TwoLineElement::from_lines(lines)?;
375        Ok(())
376    }
377
378    #[test]
379    fn test_tle_from_lines_with_surrounding_whitespace() -> Result<()> {
380        let lines = "\nISS (ZARYA)
381                     1 25544U 98067A   20148.21301450  .00001715  00000-0  38778-4 0  9992
382                     2 25544  51.6435  92.2789 0002570 358.0648 144.9972 15.49396855228767\n";
383
384        let _tle = TwoLineElement::from_lines(lines)?;
385
386        Ok(())
387    }
388
389    #[test]
390    fn test_julian_day_identity() {
391        let t = Utc.with_ymd_and_hms(2020, 1, 1, 0, 0, 0).unwrap();
392        assert_eq!(DateTime::<Utc>::from(JulianDay::from(t)), t);
393    }
394
395    #[test]
396    fn test_gmst_conversion() {
397        let t = Utc.with_ymd_and_hms(2020, 1, 1, 0, 0, 0).unwrap();
398        let a: f64 = 100.1218209532; // GMST for 2020-01-01T00:00:00 in degrees
399        let a_rad = a.to_radians();
400        assert!(sgp4_sys::close(
401            GreenwichMeanSiderealTime::from(t).as_radians(),
402            a_rad
403        ));
404    }
405}
sgp4_rs/lib.rs

sgp4_rs/
lib.rs
