lox_frames/
rotations.rs

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

lox_frames/
rotations.rs
