lox_frames/
rotations.rs

1// SPDX-FileCopyrightText: 2025 Helge Eichhorn <git@helgeeichhorn.de>
2//
3// SPDX-License-Identifier: MPL-2.0
4
5use std::{fmt::Display, ops::Mul};
6
7use glam::{DMat3, DVec3};
8use lox_bodies::{TryRotationalElements, UndefinedOriginPropertyError};
9use lox_core::{coords::Cartesian, f64::consts::ROTATION_RATE_EARTH};
10use lox_test_utils::ApproxEq;
11use lox_time::{
12    Time,
13    julian_dates::JulianDate,
14    offsets::{OffsetProvider, TryOffset},
15    time_scales::{Tdb, TimeScale, Tt, Ut1},
16};
17use thiserror::Error;
18
19use crate::{
20    Iau, ReferenceFrame,
21    iau::icrf_to_iau,
22    iers::{
23        Corrections, ReferenceSystem,
24        cio::CioLocator,
25        cip::CipCoords,
26        earth_rotation::{EarthRotationAngle, EquationOfTheEquinoxes},
27        polar_motion::PoleCoords,
28        precession::frame_bias,
29    },
30};
31
32mod impls;
33
34pub trait TryRotation<Origin, Target, T>
35where
36    Origin: ReferenceFrame,
37    Target: ReferenceFrame,
38    T: TimeScale,
39{
40    type Error: std::error::Error + Send + Sync + 'static;
41
42    fn try_rotation(
43        &self,
44        origin: Origin,
45        target: Target,
46        time: Time<T>,
47    ) -> Result<Rotation, Self::Error>;
48}
49
50pub trait TryComposedRotation<T, P>
51where
52    T: TimeScale + Copy,
53{
54    fn try_composed_rotation(&self, provider: &P, time: Time<T>)
55    -> Result<Rotation, RotationError>;
56}
57
58impl<T, P, R1, R2, R3> TryComposedRotation<T, P> for (R1, R2, R3)
59where
60    T: TimeScale + Copy,
61    R1: ReferenceFrame + Copy,
62    R2: ReferenceFrame + Copy,
63    R3: ReferenceFrame + Copy,
64    P: RotationProvider<T>
65        + TryRotation<R1, R2, T, Error = RotationError>
66        + TryRotation<R2, R3, T, Error = RotationError>,
67{
68    fn try_composed_rotation(
69        &self,
70        provider: &P,
71        time: Time<T>,
72    ) -> Result<Rotation, RotationError> {
73        Ok(provider
74            .try_rotation(self.0, self.1, time)?
75            .compose(provider.try_rotation(self.1, self.2, time)?))
76    }
77}
78
79impl<T, P, R1, R2, R3, R4> TryComposedRotation<T, P> for (R1, R2, R3, R4)
80where
81    T: TimeScale + Copy,
82    R1: ReferenceFrame + Copy,
83    R2: ReferenceFrame + Copy,
84    R3: ReferenceFrame + Copy,
85    R4: ReferenceFrame + Copy,
86    P: RotationProvider<T>
87        + TryRotation<R1, R2, T, Error = RotationError>
88        + TryRotation<R2, R3, T, Error = RotationError>
89        + TryRotation<R3, R4, T, Error = RotationError>,
90{
91    fn try_composed_rotation(
92        &self,
93        provider: &P,
94        time: Time<T>,
95    ) -> Result<Rotation, RotationError> {
96        Ok(provider
97            .try_rotation(self.0, self.1, time)?
98            .compose(provider.try_rotation(self.1, self.2, time)?)
99            .compose(provider.try_rotation(self.2, self.3, time)?))
100    }
101}
102
103impl<T, P, R1, R2, R3, R4, R5> TryComposedRotation<T, P> for (R1, R2, R3, R4, R5)
104where
105    T: TimeScale + Copy,
106    R1: ReferenceFrame + Copy,
107    R2: ReferenceFrame + Copy,
108    R3: ReferenceFrame + Copy,
109    R4: ReferenceFrame + Copy,
110    R5: ReferenceFrame + Copy,
111    P: RotationProvider<T>
112        + TryRotation<R1, R2, T, Error = RotationError>
113        + TryRotation<R2, R3, T, Error = RotationError>
114        + TryRotation<R3, R4, T, Error = RotationError>
115        + TryRotation<R4, R5, T, Error = RotationError>,
116{
117    fn try_composed_rotation(
118        &self,
119        provider: &P,
120        time: Time<T>,
121    ) -> Result<Rotation, RotationError> {
122        Ok(provider
123            .try_rotation(self.0, self.1, time)?
124            .compose(provider.try_rotation(self.1, self.2, time)?)
125            .compose(provider.try_rotation(self.2, self.3, time)?)
126            .compose(provider.try_rotation(self.3, self.4, time)?))
127    }
128}
129
130impl<T, P, R1, R2, R3, R4, R5, R6> TryComposedRotation<T, P> for (R1, R2, R3, R4, R5, R6)
131where
132    T: TimeScale + Copy,
133    R1: ReferenceFrame + Copy,
134    R2: ReferenceFrame + Copy,
135    R3: ReferenceFrame + Copy,
136    R4: ReferenceFrame + Copy,
137    R5: ReferenceFrame + Copy,
138    R6: ReferenceFrame + Copy,
139    P: RotationProvider<T>
140        + TryRotation<R1, R2, T, Error = RotationError>
141        + TryRotation<R2, R3, T, Error = RotationError>
142        + TryRotation<R3, R4, T, Error = RotationError>
143        + TryRotation<R4, R5, T, Error = RotationError>
144        + TryRotation<R5, R6, T, Error = RotationError>,
145{
146    fn try_composed_rotation(
147        &self,
148        provider: &P,
149        time: Time<T>,
150    ) -> Result<Rotation, RotationError> {
151        Ok(provider
152            .try_rotation(self.0, self.1, time)?
153            .compose(provider.try_rotation(self.1, self.2, time)?)
154            .compose(provider.try_rotation(self.2, self.3, time)?)
155            .compose(provider.try_rotation(self.3, self.4, time)?)
156            .compose(provider.try_rotation(self.4, self.5, time)?))
157    }
158}
159
160#[derive(Debug)]
161pub enum RotationErrorKind {
162    Offset,
163    Eop,
164}
165
166impl Display for RotationErrorKind {
167    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
168        match self {
169            RotationErrorKind::Offset => "offset error".fmt(f),
170            RotationErrorKind::Eop => "EOP error".fmt(f),
171        }
172    }
173}
174
175#[derive(Debug, Error)]
176#[error("{kind}: {error}")]
177pub struct RotationError {
178    kind: RotationErrorKind,
179    error: Box<dyn std::error::Error + Send + Sync + 'static>,
180}
181
182impl RotationError {
183    pub fn offset(err: impl std::error::Error + Send + Sync + 'static) -> Self {
184        RotationError {
185            kind: RotationErrorKind::Offset,
186            error: Box::new(err),
187        }
188    }
189
190    pub fn eop(err: impl std::error::Error + Send + Sync + 'static) -> Self {
191        RotationError {
192            kind: RotationErrorKind::Eop,
193            error: Box::new(err),
194        }
195    }
196}
197
198#[derive(Debug, Error)]
199pub enum DynRotationError {
200    #[error(transparent)]
201    Offset(#[from] RotationError),
202    #[error("incompatible reference systems")]
203    IncompatibleReferenceSystems,
204    #[error(transparent)]
205    UndefinedProperty(#[from] UndefinedOriginPropertyError),
206}
207
208pub trait RotationProvider<T: TimeScale>: OffsetProvider {
209    type EopError: std::error::Error + Send + Sync + 'static;
210
211    fn corrections(
212        &self,
213        time: Time<T>,
214        sys: ReferenceSystem,
215    ) -> Result<Corrections, Self::EopError>;
216    fn pole_coords(&self, time: Time<T>) -> Result<PoleCoords, Self::EopError>;
217
218    fn icrf_to_iau<R>(&self, time: Time<T>, frame: Iau<R>) -> Result<Rotation, RotationError>
219    where
220        T: TimeScale + Copy,
221        R: TryRotationalElements,
222        Self: TryOffset<T, Tdb>,
223    {
224        let seconds = time
225            .try_to_scale(Tdb, self)
226            .map_err(RotationError::offset)?
227            .seconds_since_j2000();
228        let angles = frame.rotational_elements(seconds);
229        let rates = frame.rotational_element_rates(seconds);
230
231        Ok(icrf_to_iau(angles, rates))
232    }
233    fn iau_to_icrf<R>(&self, time: Time<T>, frame: Iau<R>) -> Result<Rotation, RotationError>
234    where
235        T: TimeScale + Copy,
236        R: TryRotationalElements,
237        Self: TryOffset<T, Tdb>,
238    {
239        Ok(self.icrf_to_iau(time, frame)?.transpose())
240    }
241
242    // TODO: Support other IERS conventions
243    fn icrf_to_itrf(&self, time: Time<T>) -> Result<Rotation, RotationError>
244    where
245        T: TimeScale + Copy,
246        Self: TryOffset<T, Tdb> + TryOffset<T, Tt> + TryOffset<T, Ut1>,
247    {
248        Ok(self
249            .icrf_to_cirf(time)?
250            .compose(self.cirf_to_tirf(time)?)
251            .compose(self.tirf_to_itrf(time)?))
252    }
253    fn itrf_to_icrf(&self, time: Time<T>) -> Result<Rotation, RotationError>
254    where
255        T: TimeScale + Copy,
256        Self: TryOffset<T, Tdb> + TryOffset<T, Tt> + TryOffset<T, Ut1>,
257    {
258        Ok(self.icrf_to_itrf(time)?.transpose())
259    }
260
261    fn icrf_to_j2000(&self) -> Rotation {
262        Rotation::new(frame_bias())
263    }
264
265    fn j2000_to_icrf(&self) -> Rotation {
266        Rotation::new(frame_bias().transpose())
267    }
268
269    fn j2000_to_mod(&self, time: Time<T>, sys: ReferenceSystem) -> Result<Rotation, RotationError>
270    where
271        T: TimeScale + Copy,
272        Self: TryOffset<T, Tt>,
273    {
274        let time = time.try_to_scale(Tt, self).map_err(RotationError::offset)?;
275        Ok(sys.precession_matrix(time).into())
276    }
277
278    fn mod_to_j2000(&self, time: Time<T>, sys: ReferenceSystem) -> Result<Rotation, RotationError>
279    where
280        T: TimeScale + Copy,
281        Self: TryOffset<T, Tt>,
282    {
283        Ok(self.j2000_to_mod(time, sys)?.transpose())
284    }
285
286    fn icrf_to_mod(&self, time: Time<T>, sys: ReferenceSystem) -> Result<Rotation, RotationError>
287    where
288        T: TimeScale + Copy,
289        Self: TryOffset<T, Tt>,
290    {
291        let time = time.try_to_scale(Tt, self).map_err(RotationError::offset)?;
292        Ok(sys.bias_precession_matrix(time).into())
293    }
294    fn mod_to_icrf(&self, time: Time<T>, sys: ReferenceSystem) -> Result<Rotation, RotationError>
295    where
296        T: TimeScale + Copy,
297        Self: TryOffset<T, Tt>,
298    {
299        Ok(self.icrf_to_mod(time, sys)?.transpose())
300    }
301
302    fn mod_to_tod(&self, time: Time<T>, sys: ReferenceSystem) -> Result<Rotation, RotationError>
303    where
304        T: TimeScale + Copy,
305        Self: TryOffset<T, Tdb>,
306    {
307        let tdb = time
308            .try_to_scale(Tdb, self)
309            .map_err(RotationError::offset)?;
310        let corr = self.corrections(time, sys).map_err(RotationError::eop)?;
311        Ok(sys.nutation_matrix(tdb, corr).into())
312    }
313    fn tod_to_mod(&self, time: Time<T>, sys: ReferenceSystem) -> Result<Rotation, RotationError>
314    where
315        T: TimeScale + Copy,
316        Self: TryOffset<T, Tdb>,
317    {
318        Ok(self.mod_to_tod(time, sys)?.transpose())
319    }
320
321    fn tod_to_pef(&self, time: Time<T>, sys: ReferenceSystem) -> Result<Rotation, RotationError>
322    where
323        T: TimeScale + Copy,
324        Self: TryOffset<T, Tt> + TryOffset<T, Ut1>,
325    {
326        let tt = time.try_to_scale(Tt, self).map_err(RotationError::offset)?;
327        let ut1 = time
328            .try_to_scale(Ut1, self)
329            .map_err(RotationError::offset)?;
330        let corr = self.corrections(time, sys).map_err(RotationError::eop)?;
331        Ok(
332            Rotation::new(sys.earth_rotation(tt, ut1, corr)).with_angular_velocity(DVec3::new(
333                0.0,
334                0.0,
335                ROTATION_RATE_EARTH,
336            )),
337        )
338    }
339    fn pef_to_tod(&self, time: Time<T>, sys: ReferenceSystem) -> Result<Rotation, RotationError>
340    where
341        T: TimeScale + Copy,
342        Self: TryOffset<T, Tt> + TryOffset<T, Ut1>,
343    {
344        Ok(self.tod_to_pef(time, sys)?.transpose())
345    }
346
347    fn pef_to_itrf(&self, time: Time<T>, sys: ReferenceSystem) -> Result<Rotation, RotationError>
348    where
349        T: TimeScale + Copy,
350        Self: TryOffset<T, Tt>,
351    {
352        let tt = time.try_to_scale(Tt, self).map_err(RotationError::offset)?;
353        let pole_coords = self.pole_coords(time).map_err(RotationError::eop)?;
354        Ok(sys.polar_motion_matrix(tt, pole_coords).into())
355    }
356
357    fn itrf_to_pef(&self, time: Time<T>, sys: ReferenceSystem) -> Result<Rotation, RotationError>
358    where
359        T: TimeScale + Copy,
360        Self: TryOffset<T, Tt>,
361    {
362        Ok(self.pef_to_itrf(time, sys)?.transpose())
363    }
364
365    fn tod_to_teme(&self, time: Time<T>) -> Result<Rotation, RotationError>
366    where
367        T: TimeScale + Copy,
368        Self: TryOffset<T, Tdb>,
369    {
370        let tdb = time
371            .try_to_scale(Tdb, self)
372            .map_err(RotationError::offset)?;
373
374        let eoe = EquationOfTheEquinoxes::iau1994(tdb);
375
376        // TOD to TEME is R_z(EoE)
377        Ok(Rotation::new(eoe.0.rotation_z()))
378    }
379
380    fn teme_to_tod(&self, time: Time<T>) -> Result<Rotation, RotationError>
381    where
382        T: TimeScale + Copy,
383        Self: TryOffset<T, Tdb>,
384    {
385        Ok(self.tod_to_teme(time)?.transpose())
386    }
387
388    fn icrf_to_cirf(&self, time: Time<T>) -> Result<Rotation, RotationError>
389    where
390        T: TimeScale + Copy,
391        Self: TryOffset<T, Tdb>,
392    {
393        let tdb = time
394            .try_to_scale(Tdb, self)
395            .map_err(RotationError::offset)?;
396        let mut xy = CipCoords::iau2006(tdb);
397        let s = CioLocator::iau2006(tdb, xy);
398
399        xy += self
400            .corrections(time, ReferenceSystem::Iers2010)
401            .unwrap_or_default();
402
403        Ok(Rotation::new(xy.celestial_to_intermediate_matrix(s)))
404    }
405    fn cirf_to_icrf(&self, time: Time<T>) -> Result<Rotation, RotationError>
406    where
407        T: TimeScale + Copy,
408        Self: TryOffset<T, Tdb>,
409    {
410        Ok(self.icrf_to_cirf(time)?.transpose())
411    }
412
413    fn cirf_to_tirf(&self, time: Time<T>) -> Result<Rotation, RotationError>
414    where
415        T: TimeScale + Copy,
416        Self: TryOffset<T, Ut1>,
417    {
418        let time = time
419            .try_to_scale(Ut1, self)
420            .map_err(RotationError::offset)?;
421        let era = EarthRotationAngle::iau2000(time);
422        Ok(
423            Rotation::new(era.0.rotation_z()).with_angular_velocity(DVec3::new(
424                0.0,
425                0.0,
426                ROTATION_RATE_EARTH,
427            )),
428        )
429    }
430    fn tirf_to_cirf(&self, time: Time<T>) -> Result<Rotation, RotationError>
431    where
432        T: TimeScale + Copy,
433        Self: TryOffset<T, Ut1>,
434    {
435        Ok(self.cirf_to_tirf(time)?.transpose())
436    }
437
438    fn tirf_to_itrf(&self, time: Time<T>) -> Result<Rotation, RotationError>
439    where
440        T: TimeScale + Copy,
441        Self: TryOffset<T, Tt>,
442    {
443        let tt = time.try_to_scale(Tt, self).map_err(RotationError::offset)?;
444        let pole_coords = self.pole_coords(time).map_err(RotationError::eop)?;
445        Ok(ReferenceSystem::Iers2010
446            .polar_motion_matrix(tt, pole_coords)
447            .into())
448    }
449
450    fn itrf_to_tirf(&self, time: Time<T>) -> Result<Rotation, RotationError>
451    where
452        T: TimeScale + Copy,
453        Self: TryOffset<T, Tt>,
454    {
455        Ok(self.tirf_to_itrf(time)?.transpose())
456    }
457}
458
459fn rotation_matrix_derivative(m: DMat3, v: DVec3) -> DMat3 {
460    let sx = DVec3::new(0.0, v.z, v.y);
461    let sy = DVec3::new(-v.z, 0.0, v.x);
462    let sz = DVec3::new(v.y, -v.x, 0.0);
463    let s = DMat3::from_cols(sx, sy, sz);
464    -s * m
465}
466
467#[derive(Debug, Clone, Copy, PartialEq, ApproxEq)]
468#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
469pub struct Rotation {
470    /// Rotation matrix
471    pub m: DMat3,
472    /// Time derivative of the rotation matrix
473    pub dm: DMat3,
474}
475
476impl Rotation {
477    pub const IDENTITY: Self = Self {
478        m: DMat3::IDENTITY,
479        dm: DMat3::ZERO,
480    };
481
482    pub fn new(m: DMat3) -> Self {
483        Self { m, dm: DMat3::ZERO }
484    }
485
486    pub fn with_derivative(mut self, dm: DMat3) -> Self {
487        self.dm = dm;
488        self
489    }
490
491    pub fn with_angular_velocity(mut self, v: DVec3) -> Self {
492        self.dm = rotation_matrix_derivative(self.m, v);
493        self
494    }
495
496    pub fn position_matrix(&self) -> DMat3 {
497        self.m
498    }
499
500    pub fn velocity_matrix(&self) -> DMat3 {
501        self.dm
502    }
503
504    pub fn compose(self, other: Self) -> Self {
505        Self {
506            m: other.m * self.m,
507            dm: other.dm * self.m + other.m * self.dm,
508        }
509    }
510
511    pub fn transpose(&self) -> Self {
512        let m = self.m.transpose();
513        let dm = self.dm.transpose();
514        Self { m, dm }
515    }
516
517    pub fn rotate_position(&self, pos: DVec3) -> DVec3 {
518        self.m * pos
519    }
520
521    pub fn rotate_velocity(&self, pos: DVec3, vel: DVec3) -> DVec3 {
522        self.dm * pos + self.m * vel
523    }
524
525    pub fn rotate_state(&self, pos: DVec3, vel: DVec3) -> (DVec3, DVec3) {
526        (self.rotate_position(pos), self.rotate_velocity(pos, vel))
527    }
528}
529
530impl Default for Rotation {
531    fn default() -> Self {
532        Self {
533            m: DMat3::IDENTITY,
534            dm: DMat3::ZERO,
535        }
536    }
537}
538
539impl Mul<DVec3> for Rotation {
540    type Output = DVec3;
541
542    fn mul(self, rhs: DVec3) -> Self::Output {
543        self.m * rhs
544    }
545}
546
547impl Mul<Cartesian> for Rotation {
548    type Output = Cartesian;
549
550    fn mul(self, rhs: Cartesian) -> Self::Output {
551        let pos = self.m * rhs.position();
552        let vel = self.dm * rhs.position() + self.m * rhs.velocity();
553        Cartesian::from_vecs(pos, vel)
554    }
555}
556
557impl From<DMat3> for Rotation {
558    fn from(matrix: DMat3) -> Self {
559        Rotation::new(matrix)
560    }
561}
562
563#[cfg(test)]
564mod tests {
565    use std::convert::Infallible;
566
567    use lox_test_utils::assert_approx_eq;
568    use lox_time::{deltas::TimeDelta, offsets::OffsetProvider};
569    use lox_units::AngleUnits;
570
571    use crate::iers::Iau2000Model;
572
573    use super::*;
574
575    #[derive(Debug)]
576    struct TestRotationProvider;
577
578    impl OffsetProvider for TestRotationProvider {
579        type Error = Infallible;
580
581        fn tai_to_ut1(&self, _delta: TimeDelta) -> Result<TimeDelta, Self::Error> {
582            Ok(TimeDelta::from_seconds_f64(-33.072073684954375))
583        }
584
585        fn ut1_to_tai(&self, _delta: TimeDelta) -> Result<TimeDelta, Self::Error> {
586            unreachable!()
587        }
588    }
589
590    impl<T> RotationProvider<T> for TestRotationProvider
591    where
592        T: TimeScale,
593    {
594        type EopError = Infallible;
595
596        fn corrections(
597            &self,
598            _time: Time<T>,
599            sys: ReferenceSystem,
600        ) -> Result<Corrections, Infallible> {
601            match sys {
602                ReferenceSystem::Iers1996 => {
603                    Ok(Corrections(-55.0655e-3.arcsec(), -6.3580e-3.arcsec()))
604                }
605                ReferenceSystem::Iers2003(_) => {
606                    Ok(Corrections(0.1725e-3.arcsec(), -0.2650e-3.arcsec()))
607                }
608                ReferenceSystem::Iers2010 => {
609                    Ok(Corrections(0.1750e-3.arcsec(), -0.2259e-3.arcsec()))
610                }
611            }
612        }
613
614        fn pole_coords(&self, _time: Time<T>) -> Result<PoleCoords, Infallible> {
615            Ok(PoleCoords {
616                xp: 0.0349282.arcsec(),
617                yp: 0.4833163.arcsec(),
618            })
619        }
620    }
621
622    #[test]
623    fn test_celestial_to_terrestrial_iers1996() {
624        let tt = Time::from_two_part_julian_date(Tt, 2454195.5, 0.500754444444444);
625        let sys = ReferenceSystem::Iers1996;
626
627        let npb_exp = DMat3::from_cols_array(&[
628            0.999998403176203,
629            -0.001639032970562,
630            -0.000712190961847,
631            0.001639000942243,
632            0.999998655799521,
633            -0.000045552846624,
634            0.000712264667137,
635            0.000044385492226,
636            0.999999745354454,
637        ])
638        .transpose();
639        let c2t_exp = DMat3::from_cols_array(&[
640            0.973104317592265,
641            0.230363826166883,
642            -0.000703332813776,
643            -0.230363798723533,
644            0.973104570754697,
645            0.000120888299841,
646            0.000712264667137,
647            0.000044385492226,
648            0.999999745354454,
649        ])
650        .transpose();
651        let c2t_pm_exp = DMat3::from_cols_array(&[
652            0.973104317712772,
653            0.230363826174782,
654            -0.000703163477127,
655            -0.230363800391868,
656            0.973104570648022,
657            0.000118545116892,
658            0.000711560100206,
659            0.000046626645796,
660            0.999999745754058,
661        ])
662        .transpose();
663
664        let npb_act = TestRotationProvider
665            .j2000_to_mod(tt.with_scale(Tdb), sys)
666            .unwrap()
667            .compose(
668                TestRotationProvider
669                    .mod_to_tod(tt.with_scale(Tdb), sys)
670                    .unwrap(),
671            );
672        assert_approx_eq!(npb_act.m, npb_exp, atol <= 1e-12);
673
674        let c2t_act = npb_act.compose(TestRotationProvider.tod_to_pef(tt, sys).unwrap());
675        assert_approx_eq!(c2t_act.m, c2t_exp, atol <= 1e-12);
676
677        let c2t_pm_act = c2t_act.compose(TestRotationProvider.pef_to_itrf(tt, sys).unwrap());
678        assert_approx_eq!(c2t_pm_act.m, c2t_pm_exp, atol <= 1e-12);
679    }
680
681    #[test]
682    fn test_celestial_to_terrestrial_iers2003() {
683        let tt = Time::from_two_part_julian_date(Tt, 2454195.5, 0.500754444444444);
684        let sys = ReferenceSystem::Iers2003(Iau2000Model::A);
685
686        let npb_exp = DMat3::from_cols_array(&[
687            0.999998402755640,
688            -0.001639289519579,
689            -0.000712191013215,
690            0.001639257491365,
691            0.999998655379006,
692            -0.000045552787478,
693            0.000712264729795,
694            0.000044385250265,
695            0.999999745354420,
696        ])
697        .transpose();
698        let c2t_exp = DMat3::from_cols_array(&[
699            0.973104317573209,
700            0.230363826247361,
701            -0.000703332818999,
702            -0.230363798803834,
703            0.973104570735656,
704            0.000120888549787,
705            0.000712264729795,
706            0.000044385250265,
707            0.999999745354420,
708        ])
709        .transpose();
710        let c2t_pm_exp = DMat3::from_cols_array(&[
711            0.973104317697618,
712            0.230363826238780,
713            -0.000703163482352,
714            -0.230363800455689,
715            0.973104570632883,
716            0.000118545366826,
717            0.000711560162864,
718            0.000046626403835,
719            0.999999745754024,
720        ])
721        .transpose();
722
723        let npb_act = TestRotationProvider
724            .icrf_to_mod(tt, sys)
725            .unwrap()
726            .compose(TestRotationProvider.mod_to_tod(tt, sys).unwrap());
727        assert_approx_eq!(npb_act.m, npb_exp, atol <= 1e-12);
728
729        let c2t_act = npb_act.compose(TestRotationProvider.tod_to_pef(tt, sys).unwrap());
730        assert_approx_eq!(c2t_act.m, c2t_exp, atol <= 1e-12);
731
732        let c2t_pm_act = c2t_act.compose(TestRotationProvider.pef_to_itrf(tt, sys).unwrap());
733        assert_approx_eq!(c2t_pm_act.m, c2t_pm_exp, atol <= 1e-12);
734    }
735
736    #[test]
737    fn test_celestial_to_terrestrial_iau2006() {
738        let tt = Time::from_two_part_julian_date(Tt, 2454195.5, 0.500754444444444);
739
740        let npb_exp = DMat3::from_cols_array(&[
741            0.999999746339445,
742            -0.000000005138822,
743            -0.000712264730072,
744            -0.000000026475227,
745            0.999999999014975,
746            -0.000044385242827,
747            0.000712264729599,
748            0.000044385250426,
749            0.999999745354420,
750        ])
751        .transpose();
752        let c2t_exp = DMat3::from_cols_array(&[
753            0.973104317573127,
754            0.230363826247709,
755            -0.000703332818845,
756            -0.230363798804182,
757            0.973104570735574,
758            0.000120888549586,
759            0.000712264729599,
760            0.000044385250426,
761            0.999999745354420,
762        ])
763        .transpose();
764        let c2t_pm_exp = DMat3::from_cols_array(&[
765            0.973104317697535,
766            0.230363826239128,
767            -0.000703163482198,
768            -0.230363800456037,
769            0.973104570632801,
770            0.000118545366625,
771            0.000711560162668,
772            0.000046626403995,
773            0.999999745754024,
774        ])
775        .transpose();
776
777        let npb_act = TestRotationProvider.icrf_to_cirf(tt).unwrap();
778        assert_approx_eq!(npb_act.m, npb_exp, atol <= 1e-11);
779
780        let c2t_act = npb_act.compose(TestRotationProvider.cirf_to_tirf(tt).unwrap());
781        assert_approx_eq!(c2t_act.m, c2t_exp, atol <= 1e-11);
782
783        let c2t_pm_act = c2t_act.compose(TestRotationProvider.tirf_to_itrf(tt).unwrap());
784        assert_approx_eq!(c2t_pm_act.m, c2t_pm_exp, atol <= 1e-11);
785    }
786
787    #[test]
788    fn test_tod_to_teme() {
789        // Use the same time as the EoE test to verify against known reference value
790        // EoE at this time = 5.357_758_254_609_257e-5 radians (from test_equation_of_the_equinoxes_iau1994)
791        let tdb = Time::from_two_part_julian_date(Tdb, 2400000.5, 41234.0);
792        let eoe: f64 = 5.357_758_254_609_257e-5; // radians
793
794        let rotation = TestRotationProvider.tod_to_teme(tdb).unwrap();
795
796        // TOD to TEME is R_z(EoE), so the rotation matrix should be:
797        // [cos(EoE)  -sin(EoE)  0]
798        // [sin(EoE)   cos(EoE)  0]
799        // [0          0         1]
800        let (sin_eoe, cos_eoe) = eoe.sin_cos();
801        let expected =
802            DMat3::from_cols_array(&[cos_eoe, sin_eoe, 0.0, -sin_eoe, cos_eoe, 0.0, 0.0, 0.0, 1.0])
803                .transpose();
804
805        assert_approx_eq!(rotation.m, expected, atol <= 1e-15);
806
807        // Verify round-trip
808        let roundtrip = rotation.compose(TestRotationProvider.teme_to_tod(tdb).unwrap());
809        assert_approx_eq!(roundtrip.m, DMat3::IDENTITY, atol <= 1e-15);
810    }
811
812    #[test]
813    fn test_teme_icrf_roundtrip() {
814        // Test the full TEME <-> ICRF transformation chain
815        // Path: ICRF -> MOD -> TOD -> TEME and back
816        let tt = Time::from_two_part_julian_date(Tt, 2454195.5, 0.500754444444444);
817        let sys = ReferenceSystem::Iers1996;
818
819        // Build the full ICRF to TEME transformation
820        let icrf_to_mod = TestRotationProvider.icrf_to_mod(tt, sys).unwrap();
821        let mod_to_tod = TestRotationProvider.mod_to_tod(tt, sys).unwrap();
822        let tod_to_teme = TestRotationProvider.tod_to_teme(tt).unwrap();
823
824        let icrf_to_teme = icrf_to_mod.compose(mod_to_tod).compose(tod_to_teme);
825
826        // Build the reverse transformation
827        let teme_to_tod = TestRotationProvider.teme_to_tod(tt).unwrap();
828        let tod_to_mod = TestRotationProvider.tod_to_mod(tt, sys).unwrap();
829        let mod_to_icrf = TestRotationProvider.mod_to_icrf(tt, sys).unwrap();
830
831        let teme_to_icrf = teme_to_tod.compose(tod_to_mod).compose(mod_to_icrf);
832
833        // Round-trip should give identity
834        let roundtrip = icrf_to_teme.compose(teme_to_icrf);
835        assert_approx_eq!(roundtrip.m, DMat3::IDENTITY, atol <= 1e-14);
836    }
837
838    #[test]
839    fn test_icrf_j2000_roundtrip() {
840        let fwd =
841            <TestRotationProvider as RotationProvider<Tt>>::icrf_to_j2000(&TestRotationProvider);
842        let rev =
843            <TestRotationProvider as RotationProvider<Tt>>::j2000_to_icrf(&TestRotationProvider);
844        let roundtrip = fwd.compose(rev);
845        assert_approx_eq!(roundtrip.m, DMat3::IDENTITY, atol <= 1e-15);
846    }
847
848    #[test]
849    fn test_j2000_mod_equivalence_iers1996() {
850        let tt = Time::from_two_part_julian_date(Tt, 2454195.5, 0.500754444444444);
851        let sys = ReferenceSystem::Iers1996;
852
853        // ICRF -> J2000 -> MOD should equal ICRF -> MOD
854        let via_j2000 =
855            <TestRotationProvider as RotationProvider<Tt>>::icrf_to_j2000(&TestRotationProvider)
856                .compose(TestRotationProvider.j2000_to_mod(tt, sys).unwrap());
857        let direct = TestRotationProvider.icrf_to_mod(tt, sys).unwrap();
858
859        assert_approx_eq!(via_j2000.m, direct.m, atol <= 1e-14);
860    }
861
862    #[test]
863    fn test_j2000_mod_equivalence_iers2003() {
864        let tt = Time::from_two_part_julian_date(Tt, 2454195.5, 0.500754444444444);
865        let sys = ReferenceSystem::Iers2003(Iau2000Model::A);
866
867        let via_j2000 =
868            <TestRotationProvider as RotationProvider<Tt>>::icrf_to_j2000(&TestRotationProvider)
869                .compose(TestRotationProvider.j2000_to_mod(tt, sys).unwrap());
870        let direct = TestRotationProvider.icrf_to_mod(tt, sys).unwrap();
871
872        assert_approx_eq!(via_j2000.m, direct.m, atol <= 1e-14);
873    }
874
875    #[test]
876    fn test_j2000_mod_equivalence_iers2010() {
877        let tt = Time::from_two_part_julian_date(Tt, 2454195.5, 0.500754444444444);
878        let sys = ReferenceSystem::Iers2010;
879
880        let via_j2000 =
881            <TestRotationProvider as RotationProvider<Tt>>::icrf_to_j2000(&TestRotationProvider)
882                .compose(TestRotationProvider.j2000_to_mod(tt, sys).unwrap());
883        let direct = TestRotationProvider.icrf_to_mod(tt, sys).unwrap();
884
885        assert_approx_eq!(via_j2000.m, direct.m, atol <= 1e-14);
886    }
887
888    #[test]
889    fn test_j2000_full_chain_equivalence() {
890        // ICRF -> J2000 -> MOD -> TOD -> PEF -> ITRF == ICRF -> MOD -> TOD -> PEF -> ITRF
891        let tt = Time::from_two_part_julian_date(Tt, 2454195.5, 0.500754444444444);
892        let sys = ReferenceSystem::Iers2003(Iau2000Model::A);
893
894        let via_j2000 =
895            <TestRotationProvider as RotationProvider<Tt>>::icrf_to_j2000(&TestRotationProvider)
896                .compose(TestRotationProvider.j2000_to_mod(tt, sys).unwrap())
897                .compose(TestRotationProvider.mod_to_tod(tt, sys).unwrap())
898                .compose(TestRotationProvider.tod_to_pef(tt, sys).unwrap())
899                .compose(TestRotationProvider.pef_to_itrf(tt, sys).unwrap());
900
901        let direct = TestRotationProvider
902            .icrf_to_mod(tt, sys)
903            .unwrap()
904            .compose(TestRotationProvider.mod_to_tod(tt, sys).unwrap())
905            .compose(TestRotationProvider.tod_to_pef(tt, sys).unwrap())
906            .compose(TestRotationProvider.pef_to_itrf(tt, sys).unwrap());
907
908        assert_approx_eq!(via_j2000.m, direct.m, atol <= 1e-14);
909    }
910}
lox_frames/rotations.rs

lox_frames/
rotations.rs
