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

lox_frames/
rotations.rs
