lox_orbits/orbits/
keplerian.rs

1// SPDX-FileCopyrightText: 2025 Helge Eichhorn <git@helgeeichhorn.de>
2//
3// SPDX-License-Identifier: MPL-2.0
4
5use std::{
6    f64::consts::{PI, TAU},
7    iter::zip,
8};
9
10use super::{CartesianOrbit, KeplerianOrbit, Orbit, Trajectory};
11use lox_bodies::{Origin, PointMass, TryMeanRadius, TryPointMass, UndefinedOriginPropertyError};
12use lox_core::units::{AngleUnits, Distance};
13use lox_core::{
14    anomalies::{EccentricAnomaly, TrueAnomaly},
15    coords::{CartesianTrajectory, TimeStampedCartesian},
16    elements::{
17        ArgumentOfPeriapsis, Eccentricity, Inclination, Keplerian, LongitudeOfAscendingNode,
18    },
19    utils::Linspace,
20};
21use lox_frames::{NonQuasiInertialFrameError, QuasiInertial, ReferenceFrame, TryQuasiInertial};
22use lox_time::{deltas::TimeDelta, time_scales::TimeScale};
23use thiserror::Error;
24
25/// Errors that can occur when constructing a Keplerian orbit.
26#[derive(Debug, Clone, Error)]
27pub enum KeplerianOrbitError {
28    /// The reference frame is not quasi-inertial.
29    #[error(transparent)]
30    NonQuasiInertial(#[from] NonQuasiInertialFrameError),
31
32    /// The perigee radius is below the central body's mean radius.
33    #[error("perigee radius is below the origin mean radius")]
34    PerigeeCrossesBodyRadius,
35}
36
37impl<T, O, R> KeplerianOrbit<T, O, R>
38where
39    T: TimeScale,
40    O: Origin,
41    R: ReferenceFrame,
42{
43    /// Constructs a new Keplerian orbit in a quasi-inertial frame.
44    pub const fn new(keplerian: Keplerian, time: lox_time::Time<T>, origin: O, frame: R) -> Self
45    where
46        R: QuasiInertial,
47    {
48        Orbit::from_state(keplerian, time, origin, frame)
49    }
50
51    fn validate_keplerian(k: &Keplerian, origin: &O) -> Result<(), KeplerianOrbitError>
52    where
53        O: TryMeanRadius,
54    {
55        let a_km = k.semi_major_axis().to_kilometers();
56        let e = k.eccentricity().as_f64();
57        let perigee_km = a_km * (1.0 - e);
58
59        if let Ok(body_mean_radius) = origin.try_mean_radius()
60            && perigee_km < body_mean_radius.to_kilometers()
61        {
62            return Err(KeplerianOrbitError::PerigeeCrossesBodyRadius);
63        }
64
65        Ok(())
66    }
67
68    /// Constructs a Keplerian orbit with validation of the frame and perigee radius.
69    pub fn try_from_keplerian(
70        keplerian: Keplerian,
71        time: lox_time::Time<T>,
72        origin: O,
73        frame: R,
74    ) -> Result<Self, KeplerianOrbitError>
75    where
76        R: TryQuasiInertial,
77        O: TryMeanRadius,
78    {
79        frame.try_quasi_inertial()?;
80        Self::validate_keplerian(&keplerian, &origin)?;
81        Ok(Orbit::from_state(keplerian, time, origin, frame))
82    }
83
84    /// Returns the semi-major axis.
85    pub fn semi_major_axis(&self) -> Distance {
86        self.state().semi_major_axis()
87    }
88
89    /// Returns the eccentricity.
90    pub fn eccentricity(&self) -> Eccentricity {
91        self.state().eccentricity()
92    }
93
94    /// Returns the inclination.
95    pub fn inclination(&self) -> Inclination {
96        self.state().inclination()
97    }
98
99    /// Returns the longitude of the ascending node.
100    pub fn longitude_of_ascending_node(&self) -> LongitudeOfAscendingNode {
101        self.state().longitude_of_ascending_node()
102    }
103
104    /// Returns the argument of periapsis.
105    pub fn argument_of_periapsis(&self) -> ArgumentOfPeriapsis {
106        self.state().argument_of_periapsis()
107    }
108
109    /// Returns the true anomaly.
110    pub fn true_anomaly(&self) -> TrueAnomaly {
111        self.state().true_anomaly()
112    }
113
114    /// Converts this Keplerian orbit to Cartesian position and velocity.
115    pub fn to_cartesian(&self) -> CartesianOrbit<T, O, R>
116    where
117        T: Copy,
118        O: Copy + PointMass,
119        R: Copy,
120    {
121        Orbit::from_state(
122            self.state().to_cartesian(self.gravitational_parameter()),
123            self.time(),
124            self.origin(),
125            self.reference_frame(),
126        )
127    }
128
129    /// Converts this Keplerian orbit to Cartesian, returning an error if the gravitational parameter is undefined.
130    pub fn try_to_cartesian(&self) -> Result<CartesianOrbit<T, O, R>, UndefinedOriginPropertyError>
131    where
132        T: Copy,
133        O: Copy + TryPointMass,
134        R: Copy,
135    {
136        Ok(Orbit::from_state(
137            self.state()
138                .to_cartesian(self.try_gravitational_parameter()?),
139            self.time(),
140            self.origin(),
141            self.reference_frame(),
142        ))
143    }
144
145    /// Returns the orbital period, or `None` for non-elliptical orbits or undefined gravitational parameters.
146    pub fn orbital_period(&self) -> Option<TimeDelta>
147    where
148        O: TryPointMass,
149    {
150        self.state()
151            .orbital_period(self.try_gravitational_parameter().ok()?)
152    }
153
154    /// Generates a trajectory by tracing the full orbit with `n` evenly-spaced samples.
155    pub fn trace(&self, n: usize) -> Option<Trajectory<T, O, R>>
156    where
157        T: Copy,
158        O: TryPointMass + Copy,
159        R: Copy,
160    {
161        let period = self.orbital_period()?;
162        let mean_motion = TAU / period.to_seconds().to_f64();
163        let mean_anomaly_at_epoch = self.true_anomaly().to_mean(self.eccentricity()).ok()?;
164
165        let state = self.state();
166        let state_iter = state.iter_trace(self.try_gravitational_parameter().ok()?, n);
167
168        let data: CartesianTrajectory = zip(Linspace::new(-PI, PI, n), state_iter)
169            .map(|(ecc, state)| {
170                let mean_anomaly = EccentricAnomaly::new(ecc.rad()).to_mean(self.eccentricity());
171                let time_of_flight = (mean_anomaly - mean_anomaly_at_epoch).as_f64() / mean_motion;
172                TimeStampedCartesian {
173                    time: TimeDelta::from_seconds_f64(time_of_flight),
174                    state,
175                }
176            })
177            .collect();
178
179        Some(Trajectory::from_parts(
180            self.time(),
181            self.origin(),
182            self.reference_frame(),
183            data,
184        ))
185    }
186}
187
188#[cfg(test)]
189mod tests {
190    use lox_bodies::{Earth, MeanRadius};
191    use lox_frames::Icrf;
192    use lox_time::{Time, time_scales::Tai, utc::Utc};
193    use lox_units::DistanceUnits;
194
195    use super::*;
196
197    const JD1: f64 = 2458849.5;
198    const JD2: f64 = 49.78099017 - 1.0;
199
200    #[test]
201    fn test_valid_keplerian() {
202        let elements = Keplerian::new(
203            MeanRadius::mean_radius(&Earth) + 500.0.km(),
204            Eccentricity::try_new(0.0).unwrap(),
205            Inclination::try_new(97.0.deg()).unwrap(),
206            LongitudeOfAscendingNode::try_new(0.0.deg()).unwrap(),
207            ArgumentOfPeriapsis::try_new(0.0.deg()).unwrap(),
208            TrueAnomaly::new(0.0.deg()),
209        );
210
211        let epoch: Time<Tai> = Utc::from_delta(TimeDelta::from_two_part_julian_date(JD1, JD2))
212            .unwrap()
213            .to_time();
214
215        let result = KeplerianOrbit::try_from_keplerian(elements, epoch, Earth, Icrf);
216        assert!(result.is_ok());
217    }
218
219    #[test]
220    #[should_panic]
221    fn test_invalid_sma() {
222        let elements = Keplerian::new(
223            // negative altitude
224            MeanRadius::mean_radius(&Earth) - 500.0.km(),
225            Eccentricity::try_new(0.0).unwrap(),
226            Inclination::try_new(97.0.deg()).unwrap(),
227            LongitudeOfAscendingNode::try_new(0.0.deg()).unwrap(),
228            ArgumentOfPeriapsis::try_new(0.0.deg()).unwrap(),
229            TrueAnomaly::new(0.0.deg()),
230        );
231
232        let epoch: Time<Tai> = Utc::from_delta(TimeDelta::from_two_part_julian_date(JD1, JD2))
233            .unwrap()
234            .to_time();
235
236        KeplerianOrbit::try_from_keplerian(elements, epoch, Earth, Icrf).unwrap();
237    }
238}
lox_orbits/orbits/keplerian.rs

lox_orbits/orbits/
keplerian.rs
