Skip to main content

lox_frames/rotations/
to_icrf.rs

1// SPDX-FileCopyrightText: 2026 Helge Eichhorn <git@helgeeichhorn.de>
2//
3// SPDX-License-Identifier: MPL-2.0
4
5//! Rotations between reference frames, composed through ICRF.
6//!
7//! Every frame implements [`RotateToIcrf`], giving its rotation to and from
8//! ICRF. The blanket [`TryRotation`] impl uses these to rotate between any two
9//! frames.
10//!
11//! [`TryRotation`]: crate::rotations::TryRotation
12
13use lox_bodies::TryRotationalElements;
14use lox_time::{
15    Time,
16    offsets::TryOffset,
17    time_scales::{ContinuousTimeScale, Tdb, Tt, Ut1},
18};
19
20use crate::{
21    Frame,
22    frames::{Cirf, Iau, Icrf, Itrf, J2000, Mod, Pef, Teme, Tirf, Tod},
23    iers::{IersSystem, ReferenceSystem},
24    rotations::{Rotation, RotationError, RotationProvider, TryRotation},
25    traits::{FrameKey, ReferenceFrame, frame_key},
26};
27
28/// A frame that can produce its own rotation to and from ICRF from a provider's data.
29pub trait RotateToIcrf<T: ContinuousTimeScale, P> {
30    /// The error type returned when the rotation cannot be computed.
31    type Error;
32
33    /// Returns the rotation from this frame to ICRF at `time`.
34    fn rotation_to_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error>;
35
36    /// Returns the rotation from ICRF to this frame at `time`.
37    fn rotation_from_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error>;
38}
39
40/// Rotation from `origin` to `target`, composed through ICRF.
41pub fn rotation_via_icrf<T, P, O, Tg>(
42    provider: &P,
43    origin: O,
44    target: Tg,
45    time: Time<T>,
46) -> Result<Rotation, O::Error>
47where
48    T: ContinuousTimeScale + Copy,
49    O: RotateToIcrf<T, P>,
50    Tg: RotateToIcrf<T, P, Error = O::Error>,
51{
52    let origin_to_icrf = origin.rotation_to_icrf(provider, time)?;
53    let icrf_to_target = target.rotation_from_icrf(provider, time)?;
54    Ok(origin_to_icrf.compose(icrf_to_target))
55}
56
57/// Blanket rotation between any two frames that know their route to ICRF.
58impl<T, O, Tg, P> TryRotation<O, Tg, T> for P
59where
60    T: ContinuousTimeScale + Copy,
61    O: ReferenceFrame + RotateToIcrf<T, P, Error = RotationError>,
62    Tg: ReferenceFrame + RotateToIcrf<T, P, Error = RotationError>,
63{
64    type Error = RotationError;
65
66    fn try_rotation(&self, origin: O, target: Tg, time: Time<T>) -> Result<Rotation, Self::Error> {
67        // Skip work cheaply via frame keys: identical frames need no rotation,
68        // and when one endpoint is ICRF a single leg suffices (no composition).
69        let origin_key = frame_key(&origin);
70        let target_key = frame_key(&target);
71        if origin_key.is_some() && origin_key == target_key {
72            Ok(Rotation::IDENTITY)
73        } else if origin_key == Some(FrameKey::Icrf) {
74            target.rotation_from_icrf(self, time)
75        } else if target_key == Some(FrameKey::Icrf) {
76            origin.rotation_to_icrf(self, time)
77        } else {
78            rotation_via_icrf(self, origin, target, time)
79        }
80    }
81}
82
83// ---- the hub ---------------------------------------------------------------
84
85impl<T, P> RotateToIcrf<T, P> for Icrf
86where
87    T: ContinuousTimeScale + Copy,
88    P: RotationProvider<T>,
89{
90    type Error = RotationError;
91
92    fn rotation_to_icrf(&self, _provider: &P, _time: Time<T>) -> Result<Rotation, Self::Error> {
93        Ok(Rotation::IDENTITY)
94    }
95
96    fn rotation_from_icrf(&self, _provider: &P, _time: Time<T>) -> Result<Rotation, Self::Error> {
97        Ok(Rotation::IDENTITY)
98    }
99}
100
101// ---- quasi-inertial: frame bias / body-fixed -------------------------------
102
103impl<T, P> RotateToIcrf<T, P> for J2000
104where
105    T: ContinuousTimeScale + Copy,
106    P: RotationProvider<T>,
107{
108    type Error = RotationError;
109
110    fn rotation_to_icrf(&self, provider: &P, _time: Time<T>) -> Result<Rotation, Self::Error> {
111        Ok(provider.j2000_to_icrf())
112    }
113
114    fn rotation_from_icrf(&self, provider: &P, _time: Time<T>) -> Result<Rotation, Self::Error> {
115        Ok(provider.icrf_to_j2000())
116    }
117}
118
119impl<T, P, R> RotateToIcrf<T, P> for Iau<R>
120where
121    T: ContinuousTimeScale + Copy,
122    R: TryRotationalElements + Copy,
123    P: RotationProvider<T> + TryOffset<T, Tdb>,
124{
125    type Error = RotationError;
126
127    fn rotation_to_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
128        provider.iau_to_icrf(time, *self)
129    }
130
131    fn rotation_from_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
132        provider.icrf_to_iau(time, *self)
133    }
134}
135
136// ---- CIO branch: ICRF ← CIRF ← TIRF ← ITRF ---------------------------------
137
138impl<T, P> RotateToIcrf<T, P> for Cirf
139where
140    T: ContinuousTimeScale + Copy,
141    P: RotationProvider<T> + TryOffset<T, Tdb>,
142{
143    type Error = RotationError;
144
145    fn rotation_to_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
146        provider.cirf_to_icrf(time)
147    }
148
149    fn rotation_from_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
150        provider.icrf_to_cirf(time)
151    }
152}
153
154impl<T, P> RotateToIcrf<T, P> for Tirf
155where
156    T: ContinuousTimeScale + Copy,
157    P: RotationProvider<T> + TryOffset<T, Tdb> + TryOffset<T, Ut1>,
158{
159    type Error = RotationError;
160
161    fn rotation_to_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
162        Ok(provider
163            .tirf_to_cirf(time)?
164            .compose(provider.cirf_to_icrf(time)?))
165    }
166
167    fn rotation_from_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
168        Ok(provider
169            .icrf_to_cirf(time)?
170            .compose(provider.cirf_to_tirf(time)?))
171    }
172}
173
174impl<T, P> RotateToIcrf<T, P> for Itrf
175where
176    T: ContinuousTimeScale + Copy,
177    P: RotationProvider<T> + TryOffset<T, Tt> + TryOffset<T, Tdb> + TryOffset<T, Ut1>,
178{
179    type Error = RotationError;
180
181    fn rotation_to_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
182        provider.itrf_to_icrf(time)
183    }
184
185    fn rotation_from_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
186        provider.icrf_to_itrf(time)
187    }
188}
189
190// ---- equinox branch: ICRF ← MOD ← TOD ← PEF --------------------------------
191
192impl<T, P, C> RotateToIcrf<T, P> for Mod<C>
193where
194    T: ContinuousTimeScale + Copy,
195    C: IersSystem + Into<ReferenceSystem> + Copy,
196    P: RotationProvider<T> + TryOffset<T, Tt>,
197{
198    type Error = RotationError;
199
200    fn rotation_to_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
201        provider.mod_to_icrf(time, self.0.into())
202    }
203
204    fn rotation_from_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
205        provider.icrf_to_mod(time, self.0.into())
206    }
207}
208
209impl<T, P, C> RotateToIcrf<T, P> for Tod<C>
210where
211    T: ContinuousTimeScale + Copy,
212    C: IersSystem + Into<ReferenceSystem> + Copy,
213    P: RotationProvider<T> + TryOffset<T, Tt> + TryOffset<T, Tdb>,
214{
215    type Error = RotationError;
216
217    fn rotation_to_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
218        // The convention (and its nutation model) comes from the frame value,
219        // so `Tod(Iers2003(B))` genuinely computes the 2000B nutation.
220        let sys: ReferenceSystem = self.0.into();
221        Ok(provider
222            .tod_to_mod(time, sys)?
223            .compose(provider.mod_to_icrf(time, sys)?))
224    }
225
226    fn rotation_from_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
227        let sys: ReferenceSystem = self.0.into();
228        Ok(provider
229            .icrf_to_mod(time, sys)?
230            .compose(provider.mod_to_tod(time, sys)?))
231    }
232}
233
234impl<T, P, C> RotateToIcrf<T, P> for Pef<C>
235where
236    T: ContinuousTimeScale + Copy,
237    C: IersSystem + Into<ReferenceSystem> + Copy,
238    P: RotationProvider<T> + TryOffset<T, Tt> + TryOffset<T, Tdb> + TryOffset<T, Ut1>,
239{
240    type Error = RotationError;
241
242    fn rotation_to_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
243        let sys: ReferenceSystem = self.0.into();
244        Ok(provider
245            .pef_to_tod(time, sys)?
246            .compose(provider.tod_to_mod(time, sys)?)
247            .compose(provider.mod_to_icrf(time, sys)?))
248    }
249
250    fn rotation_from_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
251        let sys: ReferenceSystem = self.0.into();
252        Ok(provider
253            .icrf_to_mod(time, sys)?
254            .compose(provider.mod_to_tod(time, sys)?)
255            .compose(provider.tod_to_pef(time, sys)?))
256    }
257}
258
259// ---- TEME: tied to the IAU 1976/FK5 (IERS1996) equinox chain ---------------
260
261impl<T, P> RotateToIcrf<T, P> for Teme
262where
263    T: ContinuousTimeScale + Copy,
264    P: RotationProvider<T> + TryOffset<T, Tt> + TryOffset<T, Tdb>,
265{
266    type Error = RotationError;
267
268    fn rotation_to_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
269        provider.teme_to_icrf(time)
270    }
271
272    fn rotation_from_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
273        provider.icrf_to_teme(time)
274    }
275}
276
277// ---- dynamic dispatch ------------------------------------------------------
278
279impl<T, P> RotateToIcrf<T, P> for Frame
280where
281    T: ContinuousTimeScale + Copy,
282    P: RotationProvider<T> + TryOffset<T, Tt> + TryOffset<T, Tdb> + TryOffset<T, Ut1>,
283{
284    type Error = RotationError;
285
286    fn rotation_to_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
287        match *self {
288            Frame::Icrf => Icrf.rotation_to_icrf(provider, time),
289            Frame::J2000 => J2000.rotation_to_icrf(provider, time),
290            Frame::Cirf => Cirf.rotation_to_icrf(provider, time),
291            Frame::Tirf => Tirf.rotation_to_icrf(provider, time),
292            Frame::Itrf => Itrf.rotation_to_icrf(provider, time),
293            Frame::Iau(origin) => Iau::try_new(origin)?.rotation_to_icrf(provider, time),
294            Frame::Mod(sys) => Mod(sys).rotation_to_icrf(provider, time),
295            Frame::Tod(sys) => Tod(sys).rotation_to_icrf(provider, time),
296            Frame::Pef(sys) => Pef(sys).rotation_to_icrf(provider, time),
297            Frame::Teme => Teme.rotation_to_icrf(provider, time),
298        }
299    }
300
301    fn rotation_from_icrf(&self, provider: &P, time: Time<T>) -> Result<Rotation, Self::Error> {
302        match *self {
303            Frame::Icrf => Icrf.rotation_from_icrf(provider, time),
304            Frame::J2000 => J2000.rotation_from_icrf(provider, time),
305            Frame::Cirf => Cirf.rotation_from_icrf(provider, time),
306            Frame::Tirf => Tirf.rotation_from_icrf(provider, time),
307            Frame::Itrf => Itrf.rotation_from_icrf(provider, time),
308            Frame::Iau(origin) => Iau::try_new(origin)?.rotation_from_icrf(provider, time),
309            Frame::Mod(sys) => Mod(sys).rotation_from_icrf(provider, time),
310            Frame::Tod(sys) => Tod(sys).rotation_from_icrf(provider, time),
311            Frame::Pef(sys) => Pef(sys).rotation_from_icrf(provider, time),
312            Frame::Teme => Teme.rotation_from_icrf(provider, time),
313        }
314    }
315}
316
317#[cfg(test)]
318mod tests {
319    use lox_approx::assert_approx_eq;
320    use lox_core::glam::DMat3;
321
322    use lox_bodies::Origin;
323    use lox_time::time_scales::Tai;
324
325    use crate::iers::{Iau2000Model, Iers2003};
326    use crate::providers::DefaultRotationProvider;
327
328    use super::*;
329
330    fn epoch() -> Time<Tt> {
331        Time::from_two_part_julian_date(Tt, 2454195.5, 0.500754444444444)
332    }
333
334    fn max_abs_diff(a: DMat3, b: DMat3) -> f64 {
335        let d = a - b;
336        d.x_axis
337            .abs()
338            .max_element()
339            .max(d.y_axis.abs().max_element())
340            .max(d.z_axis.abs().max_element())
341    }
342
343    #[test]
344    fn roundtrip_icrf_itrf() {
345        let t = epoch();
346        let fwd = DefaultRotationProvider.try_rotation(Icrf, Itrf, t).unwrap();
347        let bwd = DefaultRotationProvider.try_rotation(Itrf, Icrf, t).unwrap();
348        assert_approx_eq!(fwd.m * bwd.m, DMat3::IDENTITY, atol <= 1e-14);
349    }
350
351    #[test]
352    fn rotates_between_two_non_icrf_frames() {
353        // Neither endpoint is ICRF, so the composition goes through both legs of
354        // `rotation_via_icrf`; the round-trip must return to identity.
355        let t = epoch();
356        let tod = Tod(Iers2003(Iau2000Model::A));
357        let fwd = DefaultRotationProvider.try_rotation(tod, Itrf, t).unwrap();
358        let bwd = DefaultRotationProvider.try_rotation(Itrf, tod, t).unwrap();
359        assert_approx_eq!(fwd.m * bwd.m, DMat3::IDENTITY, atol <= 1e-14);
360    }
361
362    #[test]
363    fn dynamic_rotates_between_two_non_icrf_frames() {
364        let t = epoch();
365        let tod = Frame::Tod(ReferenceSystem::Iers2003(Iau2000Model::A));
366        let fwd = DefaultRotationProvider
367            .try_rotation(tod, Frame::Itrf, t)
368            .unwrap();
369        let bwd = DefaultRotationProvider
370            .try_rotation(Frame::Itrf, tod, t)
371            .unwrap();
372        assert_approx_eq!(fwd.m * bwd.m, DMat3::IDENTITY, atol <= 1e-14);
373    }
374
375    #[test]
376    fn roundtrip_icrf_j2000() {
377        let t = epoch();
378        let fwd = DefaultRotationProvider
379            .try_rotation(Icrf, J2000, t)
380            .unwrap();
381        let bwd = DefaultRotationProvider
382            .try_rotation(J2000, Icrf, t)
383            .unwrap();
384        assert_approx_eq!(fwd.m * bwd.m, DMat3::IDENTITY, atol <= 1e-15);
385        // J2000 differs from ICRF only by the small frame bias.
386        assert!(!fwd.m.abs_diff_eq(DMat3::IDENTITY, 1e-9));
387        assert!(fwd.m.abs_diff_eq(DMat3::IDENTITY, 1e-6));
388    }
389
390    #[test]
391    fn threads_2000b_model() {
392        // The hub route reads the nutation model from the frame value, so 2000A
393        // and 2000B genuinely differ (mas-level) instead of collapsing to 2000A.
394        let t = epoch();
395        let tod_a = DefaultRotationProvider
396            .try_rotation(Icrf, Tod(Iers2003(Iau2000Model::A)), t)
397            .unwrap();
398        let tod_b = DefaultRotationProvider
399            .try_rotation(Icrf, Tod(Iers2003(Iau2000Model::B)), t)
400            .unwrap();
401        assert!(max_abs_diff(tod_a.m, tod_b.m) > 1e-9);
402    }
403
404    // ---- mixed concrete <-> Frame, served by the blanket impl -----------
405
406    fn tai_j2000() -> Time<Tai> {
407        Time::j2000(Tai)
408    }
409
410    #[test]
411    fn mixed_icrf_to_dynframe() {
412        let rot = DefaultRotationProvider
413            .try_rotation(Icrf, Frame::Icrf, tai_j2000())
414            .unwrap();
415        assert!(rot.m.abs_diff_eq(DMat3::IDENTITY, 1e-14));
416    }
417
418    #[test]
419    fn mixed_dynframe_to_icrf() {
420        let rot = DefaultRotationProvider
421            .try_rotation(Frame::Icrf, Icrf, tai_j2000())
422            .unwrap();
423        assert!(rot.m.abs_diff_eq(DMat3::IDENTITY, 1e-14));
424    }
425
426    #[test]
427    fn mixed_iau_dynorigin_and_dynframe() {
428        let iau_earth = Iau::try_new(Origin::Earth).unwrap();
429        let fwd = DefaultRotationProvider
430            .try_rotation(Icrf, Frame::Iau(Origin::Earth), tai_j2000())
431            .unwrap();
432        let bwd = DefaultRotationProvider
433            .try_rotation(iau_earth, Frame::Icrf, tai_j2000())
434            .unwrap();
435        // Non-trivial body-fixed rotation, with a clean round-trip across the
436        // concrete↔dynamic boundary.
437        assert!(!fwd.m.abs_diff_eq(DMat3::IDENTITY, 1e-6));
438        assert!((fwd.m * bwd.m).abs_diff_eq(DMat3::IDENTITY, 1e-14));
439    }
440}