lox_frames/iers/
nutation.rs

1// SPDX-FileCopyrightText: 2023 Angus Morrison <github@angus-morrison.com>
2// SPDX-FileCopyrightText: 2024 Helge Eichhorn <git@helgeeichhorn.de>
3//
4// SPDX-License-Identifier: MPL-2.0
5
6//! Module nutation exposes a function for calculating Earth nutation using a number of IAU nutation
7//! models.
8
9use std::ops::Add;
10
11use glam::DMat3;
12use lox_core::units::Angle;
13use lox_test_utils::ApproxEq;
14use lox_time::{
15    Time,
16    time_scales::{Tdb, Tt},
17};
18
19use crate::iers::{Corrections, Iau2000Model, ReferenceSystem, ecliptic::MeanObliquity};
20
21mod iau1980;
22mod iau2000;
23mod iau2006;
24
25impl ReferenceSystem {
26    /// Computes nutation for the active IERS convention.
27    pub fn nutation(&self, time: Time<Tdb>) -> Nutation {
28        match self {
29            ReferenceSystem::Iers1996 => Nutation::iau1980(time),
30            ReferenceSystem::Iers2003(model) => match model {
31                Iau2000Model::A => Nutation::iau2000a(time),
32                Iau2000Model::B => Nutation::iau2000b(time),
33            },
34            ReferenceSystem::Iers2010 => Nutation::iau2006a(time),
35        }
36    }
37
38    /// Computes the nutation rotation matrix, applying the given corrections.
39    pub fn nutation_matrix(&self, time: Time<Tdb>, corr: Corrections) -> DMat3 {
40        match self {
41            ReferenceSystem::Iers1996 => {
42                let mut nut = self.nutation(time);
43                nut += corr;
44                nut.nutation_matrix(self.mean_obliquity(time.with_scale(Tt)))
45            }
46            ReferenceSystem::Iers2003(_) | ReferenceSystem::Iers2010 => {
47                let tt = time.with_scale(Tt);
48                let epsa = self.mean_obliquity(tt);
49                let mut nut = self.nutation(time);
50                let rpb = self.bias_precession_matrix(tt);
51
52                nut += self.ecliptic_corrections(corr, nut, epsa, rpb);
53
54                nut.nutation_matrix(epsa)
55            }
56        }
57    }
58}
59
60/// Nutation components with respect to some ecliptic of date.
61#[derive(Debug, Default, Clone, Copy, PartialEq, ApproxEq)]
62pub struct Nutation {
63    /// δψ
64    pub dpsi: Angle,
65    /// δε
66    pub deps: Angle,
67}
68
69impl Nutation {
70    /// Builds the nutation rotation matrix from the mean obliquity of the ecliptic.
71    pub fn nutation_matrix(&self, epsa: MeanObliquity) -> DMat3 {
72        let epsa = epsa.0;
73        let rot1 = epsa.rotation_x();
74        let rot2 = (-self.dpsi).rotation_z();
75        let rot3 = (-(self.deps + epsa)).rotation_x();
76        rot3 * rot2 * rot1
77    }
78}
79
80impl Add for Nutation {
81    type Output = Self;
82
83    fn add(self, rhs: Self) -> Self::Output {
84        Self {
85            dpsi: self.dpsi + rhs.dpsi,
86            deps: self.deps + rhs.deps,
87        }
88    }
89}
90
91impl Add<&Self> for Nutation {
92    type Output = Self;
93
94    fn add(self, rhs: &Self) -> Self::Output {
95        Nutation {
96            dpsi: self.dpsi + rhs.dpsi,
97            deps: self.deps + rhs.deps,
98        }
99    }
100}
101
102#[cfg(test)]
103mod tests {
104    use lox_core::units::AngleUnits;
105    use lox_test_utils::assert_approx_eq;
106
107    use super::*;
108
109    #[test]
110    fn test_nutation_matrix() {
111        let nut = Nutation {
112            dpsi: -9.630_909_107_115_582e-6.rad(),
113            deps: 4.063_239_174_001_679e-5.rad(),
114        };
115        let obl = MeanObliquity(Angle::new(0.409_078_976_335_651));
116        let exp = DMat3::from_cols_array(&[
117            0.999_999_999_953_622_8,
118            8.836_239_320_236_25e-6,
119            3.830_833_447_458_252e-6,
120            -8.836_083_657_016_69e-6,
121            0.999_999_999_135_465_4,
122            -4.063_240_865_361_857_4e-5,
123            -3.831_192_481_833_385_5e-6,
124            4.063_237_480_216_934e-5,
125            0.999_999_999_167_166_1,
126        ])
127        .transpose();
128        let act = nut.nutation_matrix(obl);
129        assert_approx_eq!(act, exp, rtol <= 1e-12);
130    }
131}
lox_frames/iers/nutation.rs

lox_frames/iers/
nutation.rs
