lox_orbits/
ground.rs

1// SPDX-FileCopyrightText: 2024 Helge Eichhorn <git@helgeeichhorn.de>
2//
3// SPDX-License-Identifier: MPL-2.0
4
5use crate::orbits::{CartesianOrbit, TrajectorError, Trajectory};
6use crate::propagators::Propagator;
7use glam::{DMat3, DVec3};
8use lox_bodies::{DynOrigin, RotationalElements, Spheroid, TrySpheroid};
9use lox_core::coords::{Cartesian, LonLatAlt};
10use lox_core::types::units::Radians;
11use lox_core::units::Distance;
12use lox_frames::{DynFrame, Iau, ReferenceFrame};
13use lox_time::Time;
14use lox_time::deltas::TimeDelta;
15use lox_time::intervals::TimeInterval;
16use lox_time::time_scales::TimeScale;
17use thiserror::Error;
18
19/// Topocentric observation of a satellite from a ground location.
20#[derive(Clone, Debug)]
21#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
22pub struct Observables {
23    azimuth: Radians,
24    elevation: Radians,
25    range: f64,
26    range_rate: f64,
27}
28
29impl Observables {
30    /// Creates a new set of observables.
31    pub fn new(azimuth: Radians, elevation: Radians, range: f64, range_rate: f64) -> Self {
32        Observables {
33            azimuth,
34            elevation,
35            range,
36            range_rate,
37        }
38    }
39    /// Returns the azimuth angle in radians.
40    pub fn azimuth(&self) -> Radians {
41        self.azimuth
42    }
43
44    /// Returns the elevation angle in radians.
45    pub fn elevation(&self) -> Radians {
46        self.elevation
47    }
48
49    /// Returns the slant range in meters.
50    pub fn range(&self) -> f64 {
51        self.range
52    }
53
54    /// Returns the range rate in meters per second.
55    pub fn range_rate(&self) -> f64 {
56        self.range_rate
57    }
58}
59
60/// A location on the surface of a celestial body.
61#[derive(Clone, Debug)]
62#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
63pub struct GroundLocation<B: TrySpheroid> {
64    coordinates: LonLatAlt,
65    body: B,
66}
67
68/// Type alias for a ground location with a dynamic origin.
69pub type DynGroundLocation = GroundLocation<DynOrigin>;
70
71/// Infallible constructor — requires compile-time `Spheroid` guarantee.
72impl<B: Spheroid> GroundLocation<B> {
73    /// Creates a new ground location on a body that is guaranteed to be a spheroid.
74    pub fn new(coordinates: LonLatAlt, body: B) -> Self {
75        GroundLocation { coordinates, body }
76    }
77}
78
79/// Fallible constructor — for `DynOrigin` and other `TrySpheroid` types.
80impl<B: TrySpheroid> GroundLocation<B> {
81    /// Creates a new ground location, returning an error if the body has no spheroid.
82    pub fn try_new(coordinates: LonLatAlt, body: B) -> Result<Self, &'static str> {
83        if body.try_equatorial_radius().is_err() {
84            return Err("no spheroid");
85        }
86        Ok(GroundLocation { coordinates, body })
87    }
88}
89
90impl<B: TrySpheroid + Into<DynOrigin>> GroundLocation<B> {
91    /// Converts the ground location into a dynamic representation.
92    pub fn into_dyn(self) -> DynGroundLocation {
93        GroundLocation {
94            coordinates: self.coordinates,
95            body: self.body.into(),
96        }
97    }
98}
99
100impl<B: TrySpheroid> GroundLocation<B> {
101    /// Returns the central body.
102    pub fn origin(&self) -> B
103    where
104        B: Clone,
105    {
106        self.body.clone()
107    }
108
109    /// Returns the geodetic coordinates.
110    pub fn coordinates(&self) -> LonLatAlt {
111        self.coordinates
112    }
113
114    /// Returns the longitude in radians.
115    pub fn longitude(&self) -> f64 {
116        self.coordinates.lon().to_radians()
117    }
118
119    /// Returns the latitude in radians.
120    pub fn latitude(&self) -> f64 {
121        self.coordinates.lat().to_radians()
122    }
123
124    /// Returns altitude in km (for backward compatibility).
125    pub fn altitude(&self) -> f64 {
126        self.coordinates.alt().to_kilometers()
127    }
128
129    fn equatorial_radius(&self) -> Distance {
130        self.body
131            .try_equatorial_radius()
132            .expect("equatorial radius should be available")
133    }
134
135    fn flattening(&self) -> f64 {
136        self.body
137            .try_flattening()
138            .expect("flattening should be available")
139    }
140
141    /// Returns the body-fixed Cartesian position in meters.
142    pub fn body_fixed_position(&self) -> DVec3 {
143        self.coordinates
144            .to_body_fixed(self.equatorial_radius(), self.flattening())
145    }
146
147    /// Returns the rotation matrix from body-fixed to topocentric (SEZ) frame.
148    pub fn rotation_to_topocentric(&self) -> DMat3 {
149        self.coordinates.rotation_to_topocentric()
150    }
151
152    /// Computes topocentric observables from raw body-fixed position and velocity vectors.
153    pub fn compute_observables(&self, state_position: DVec3, state_velocity: DVec3) -> Observables {
154        let rot = self.rotation_to_topocentric();
155        let position = rot * (state_position - self.body_fixed_position());
156        let velocity = rot * state_velocity;
157        let range = position.length();
158        let range_rate = position.dot(velocity) / range;
159        let elevation = (position.z / range).asin();
160        let azimuth = position.y.atan2(-position.x);
161        Observables {
162            azimuth,
163            elevation,
164            range,
165            range_rate,
166        }
167    }
168
169    /// Computes topocentric observables from a Cartesian orbit in the body-fixed frame.
170    pub fn observables<T: TimeScale + Copy>(
171        &self,
172        state: CartesianOrbit<T, B, Iau<B>>,
173    ) -> Observables
174    where
175        B: RotationalElements + Copy,
176    {
177        self.compute_observables(state.position(), state.velocity())
178    }
179
180    /// Computes topocentric observables from a dynamic Cartesian orbit.
181    pub fn observables_dyn(&self, state: crate::orbits::DynCartesianOrbit) -> Observables {
182        self.compute_observables(state.position(), state.velocity())
183    }
184}
185
186/// Errors that can occur during ground propagation.
187#[derive(Debug, Error)]
188pub enum GroundPropagatorError {
189    /// A frame transformation failed.
190    #[error("frame transformation error: {0}")]
191    FrameTransformation(String),
192    /// A trajectory construction error occurred.
193    #[error(transparent)]
194    Trajectory(#[from] TrajectorError),
195}
196
197/// Propagator that produces a stationary body-fixed trajectory for a ground location.
198pub struct GroundPropagator<B: TrySpheroid, R: ReferenceFrame> {
199    location: GroundLocation<B>,
200    frame: R,
201    step: Option<TimeDelta>,
202}
203
204/// Type alias for a ground propagator with dynamic origin and frame.
205pub type DynGroundPropagator = GroundPropagator<DynOrigin, DynFrame>;
206
207/// Typed constructor -- for static bodies with `Spheroid + RotationalElements`.
208impl<B: Spheroid + RotationalElements> GroundPropagator<B, Iau<B>> {
209    /// Creates a new ground propagator in the body's IAU frame.
210    pub fn new(location: GroundLocation<B>) -> Self
211    where
212        B: Copy,
213    {
214        let frame = Iau::new(location.body);
215        GroundPropagator {
216            location,
217            frame,
218            step: None,
219        }
220    }
221}
222
223/// Fallible constructor for `DynOrigin`.
224impl GroundPropagator<DynOrigin, DynFrame> {
225    /// Creates a new ground propagator, returning an error if the body has no spheroid.
226    pub fn try_new(location: GroundLocation<DynOrigin>) -> Result<Self, &'static str> {
227        if location.body.try_equatorial_radius().is_err() {
228            return Err("no spheroid");
229        }
230        let frame = DynFrame::Iau(location.body);
231        Ok(GroundPropagator {
232            location,
233            frame,
234            step: None,
235        })
236    }
237}
238
239impl<B: TrySpheroid, R: ReferenceFrame> GroundPropagator<B, R> {
240    /// Sets the propagation time step.
241    pub fn with_step(mut self, step: TimeDelta) -> Self {
242        self.step = Some(step);
243        self
244    }
245
246    /// Returns a reference to the underlying ground location.
247    pub fn location(&self) -> &GroundLocation<B> {
248        &self.location
249    }
250
251    /// Compute the body-fixed state at a single time.
252    pub fn state_at<T: TimeScale + Copy>(&self, time: Time<T>) -> CartesianOrbit<T, B, R>
253    where
254        B: Copy,
255        R: Copy,
256    {
257        let pos = self.location.body_fixed_position();
258        CartesianOrbit::new(
259            Cartesian::from_vecs(pos, DVec3::ZERO),
260            time,
261            self.location.body,
262            self.frame,
263        )
264    }
265}
266
267/// Single `Propagator` impl covers both typed and Dyn paths.
268impl<T, B, R> Propagator<T, B> for GroundPropagator<B, R>
269where
270    T: TimeScale + Copy + PartialOrd,
271    B: TrySpheroid + lox_bodies::Origin + Copy,
272    R: ReferenceFrame + Copy,
273{
274    type Frame = R;
275    type Error = GroundPropagatorError;
276
277    fn state_at(&self, time: Time<T>) -> Result<CartesianOrbit<T, B, R>, GroundPropagatorError> {
278        Ok(self.state_at(time))
279    }
280
281    fn propagate(&self, interval: TimeInterval<T>) -> Result<Trajectory<T, B, R>, Self::Error> {
282        let pos = self.location.body_fixed_position();
283        let step = self.step.unwrap_or(TimeDelta::from_seconds(60));
284        let states: Vec<_> = interval
285            .step_by(step)
286            .map(|t| {
287                CartesianOrbit::new(
288                    Cartesian::from_vecs(pos, DVec3::ZERO),
289                    t,
290                    self.location.body,
291                    self.frame,
292                )
293            })
294            .collect();
295        Trajectory::try_new(states).map_err(Into::into)
296    }
297}
298
299#[cfg(test)]
300mod tests {
301    use lox_bodies::Earth;
302    use lox_core::coords::Cartesian;
303    use lox_frames::Icrf;
304    use lox_frames::providers::DefaultRotationProvider;
305    use lox_test_utils::assert_approx_eq;
306    use lox_time::intervals::Interval;
307    use lox_time::time_scales::Tdb;
308    use lox_time::utc::Utc;
309    use lox_time::{Time, time, utc};
310
311    use super::*;
312
313    #[test]
314    fn test_ground_location_to_body_fixed() {
315        let coords = LonLatAlt::from_degrees(-4.3676, 40.4527, 0.0).unwrap();
316        let location = GroundLocation::new(coords, Earth);
317        let expected = DVec3::new(4846130.017870638, -370132.8551351891, 4116364.272747229);
318        assert_approx_eq!(location.body_fixed_position(), expected);
319    }
320
321    #[test]
322    fn test_ground_location_rotation_to_topocentric() {
323        let coords = LonLatAlt::from_degrees(-4.3676, 40.4527, 0.0).unwrap();
324        let location = GroundLocation::new(coords, Earth);
325        let act = location.rotation_to_topocentric();
326        let exp = DMat3::from_cols(
327            DVec3::new(0.6469358921661584, 0.07615519584215287, 0.7587320591443464),
328            DVec3::new(
329                -0.049411020334552434,
330                0.9970959763965771,
331                -0.05794967578213965,
332            ),
333            DVec3::new(-0.7609418522440956, 0.0, 0.6488200809957448),
334        );
335        assert_approx_eq!(exp, act);
336    }
337
338    #[test]
339    fn test_ground_location_observables() {
340        let coords = LonLatAlt::from_degrees(-4.0, 41.0, 0.0).unwrap();
341        let location = GroundLocation::new(coords, Earth);
342        let position = DVec3::new(3359927.0, -2398072.0, 5153000.0);
343        let velocity = DVec3::new(5065.7, 5485.0, -744.0);
344        let time = time!(Tdb, 2012, 7, 1).unwrap();
345        let state = CartesianOrbit::new(
346            Cartesian::from_vecs(position, velocity),
347            time,
348            Earth,
349            Iau::new(Earth),
350        );
351        let observables = location.observables(state);
352        let expected_range = 2707700.0;
353        let expected_range_rate = -7160.0;
354        let expected_azimuth = -53.418f64.to_radians();
355        let expected_elevation = -7.077f64.to_radians();
356        assert_approx_eq!(observables.range, expected_range, rtol <= 1e-2);
357        assert_approx_eq!(observables.range_rate, expected_range_rate, rtol <= 1e-2);
358        assert_approx_eq!(observables.azimuth, expected_azimuth, rtol <= 1e-2);
359        assert_approx_eq!(observables.elevation, expected_elevation, rtol <= 1e-2);
360    }
361
362    #[test]
363    fn test_ground_propagator_body_fixed() {
364        let coords = LonLatAlt::from_degrees(-4.3676, 40.4527, 0.0).unwrap();
365        let location = GroundLocation::new(coords, Earth);
366        let propagator = GroundPropagator::new(location.clone());
367        let time = utc!(2022, 1, 31, 23).unwrap().to_time();
368        let t1 = time + TimeDelta::from_minutes(5);
369        let interval = Interval::new(time, t1);
370        let traj = propagator.propagate(interval).unwrap();
371        // All states should have the same body-fixed position
372        let expected = location.body_fixed_position();
373        for state in traj.states() {
374            assert_approx_eq!(state.position(), expected);
375            assert_approx_eq!(state.velocity(), DVec3::ZERO);
376        }
377    }
378
379    #[test]
380    fn test_ground_propagator_in_icrf() {
381        let coords = LonLatAlt::from_degrees(-4.3676, 40.4527, 0.0).unwrap();
382        let location = GroundLocation::new(coords, Earth);
383        let propagator = GroundPropagator::new(location);
384        let time = utc!(2022, 1, 31, 23).unwrap().to_time();
385        let t1 = time + TimeDelta::from_minutes(5);
386        let interval = Interval::new(time, t1);
387        let traj = propagator
388            .propagate(interval)
389            .unwrap()
390            .into_frame(Icrf, &DefaultRotationProvider)
391            .unwrap();
392        let state = traj.states()[0];
393        let expected = DVec3::new(-1765953.5510583583, 4524585.984442561, 4120189.198495323);
394        assert_approx_eq!(state.position(), expected);
395    }
396
397    #[test]
398    fn test_try_new_with_static_body() {
399        let coords = LonLatAlt::from_degrees(-4.3676, 40.4527, 0.0).unwrap();
400        let location = GroundLocation::try_new(coords, Earth).unwrap();
401        assert_approx_eq!(location.longitude(), -4.3676f64.to_radians());
402        assert_approx_eq!(location.latitude(), 40.4527f64.to_radians());
403        assert_approx_eq!(location.altitude(), 0.0);
404    }
405
406    #[test]
407    fn test_try_new_with_dyn_origin() {
408        let coords = LonLatAlt::from_degrees(-4.3676, 40.4527, 0.0).unwrap();
409        let location = GroundLocation::try_new(coords, DynOrigin::Earth).unwrap();
410        assert_eq!(location.origin(), DynOrigin::Earth);
411    }
412
413    #[test]
414    fn test_try_new_rejects_non_spheroid() {
415        let coords = LonLatAlt::from_degrees(0.0, 0.0, 0.0).unwrap();
416        let result = GroundLocation::try_new(coords, DynOrigin::SolarSystemBarycenter);
417        assert!(result.is_err());
418    }
419
420    #[test]
421    fn test_into_dyn_ground_location() {
422        let coords = LonLatAlt::from_degrees(-4.3676, 40.4527, 0.0).unwrap();
423        let location = GroundLocation::new(coords, Earth);
424        let dyn_location = location.into_dyn();
425        assert_eq!(dyn_location.origin(), DynOrigin::Earth);
426        assert_approx_eq!(dyn_location.longitude(), -4.3676f64.to_radians());
427        assert_approx_eq!(dyn_location.latitude(), 40.4527f64.to_radians());
428    }
429
430    #[test]
431    fn test_ground_propagator_try_new_with_dyn_origin() {
432        let coords = LonLatAlt::from_degrees(-4.3676, 40.4527, 0.0).unwrap();
433        let location = GroundLocation::try_new(coords, DynOrigin::Earth).unwrap();
434        let propagator = GroundPropagator::try_new(location).unwrap();
435        let time = utc!(2022, 1, 31, 23).unwrap().to_time();
436        let t1 = time + TimeDelta::from_minutes(5);
437        let interval = Interval::new(time, t1);
438        let traj = propagator
439            .propagate(interval)
440            .unwrap()
441            .into_frame(DynFrame::Icrf, &DefaultRotationProvider)
442            .unwrap();
443        let state = traj.states()[0];
444        // Same result as the static version
445        let expected = DVec3::new(-1765953.5510583583, 4524585.984442561, 4120189.198495323);
446        assert_approx_eq!(state.position(), expected);
447    }
448}
lox_orbits/ground.rs

lox_orbits/
ground.rs
