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    pub fn nutation(&self, time: Time<Tdb>) -> Nutation {
27        match self {
28            ReferenceSystem::Iers1996 => Nutation::iau1980(time),
29            ReferenceSystem::Iers2003(model) => match model {
30                Iau2000Model::A => Nutation::iau2000a(time),
31                Iau2000Model::B => Nutation::iau2000b(time),
32            },
33            ReferenceSystem::Iers2010 => Nutation::iau2006a(time),
34        }
35    }
36
37    pub fn nutation_matrix(&self, time: Time<Tdb>, corr: Corrections) -> DMat3 {
38        match self {
39            ReferenceSystem::Iers1996 => {
40                let mut nut = self.nutation(time);
41                nut += corr;
42                nut.nutation_matrix(self.mean_obliquity(time.with_scale(Tt)))
43            }
44            ReferenceSystem::Iers2003(_) | ReferenceSystem::Iers2010 => {
45                let tt = time.with_scale(Tt);
46                let epsa = self.mean_obliquity(tt);
47                let mut nut = self.nutation(time);
48                let rpb = self.bias_precession_matrix(tt);
49
50                nut += self.ecliptic_corrections(corr, nut, epsa, rpb);
51
52                nut.nutation_matrix(epsa)
53            }
54        }
55    }
56}
57
58/// Nutation components with respect to some ecliptic of date.
59#[derive(Debug, Default, Clone, Copy, PartialEq, ApproxEq)]
60pub struct Nutation {
61    /// δψ
62    pub dpsi: Angle,
63    /// δε
64    pub deps: Angle,
65}
66
67impl Nutation {
68    pub fn nutation_matrix(&self, epsa: MeanObliquity) -> DMat3 {
69        let epsa = epsa.0;
70        let rot1 = epsa.rotation_x();
71        let rot2 = (-self.dpsi).rotation_z();
72        let rot3 = (-(self.deps + epsa)).rotation_x();
73        rot3 * rot2 * rot1
74    }
75}
76
77impl Add for Nutation {
78    type Output = Self;
79
80    fn add(self, rhs: Self) -> Self::Output {
81        Self {
82            dpsi: self.dpsi + rhs.dpsi,
83            deps: self.deps + rhs.deps,
84        }
85    }
86}
87
88impl Add<&Self> for Nutation {
89    type Output = Self;
90
91    fn add(self, rhs: &Self) -> Self::Output {
92        Nutation {
93            dpsi: self.dpsi + rhs.dpsi,
94            deps: self.deps + rhs.deps,
95        }
96    }
97}
98
99#[cfg(test)]
100mod tests {
101    use lox_core::units::AngleUnits;
102    use lox_test_utils::assert_approx_eq;
103
104    use super::*;
105
106    #[test]
107    fn test_nutation_matrix() {
108        let nut = Nutation {
109            dpsi: -9.630_909_107_115_582e-6.rad(),
110            deps: 4.063_239_174_001_679e-5.rad(),
111        };
112        let obl = MeanObliquity(Angle::new(0.409_078_976_335_651));
113        let exp = DMat3::from_cols_array(&[
114            0.999_999_999_953_622_8,
115            8.836_239_320_236_25e-6,
116            3.830_833_447_458_252e-6,
117            -8.836_083_657_016_69e-6,
118            0.999_999_999_135_465_4,
119            -4.063_240_865_361_857_4e-5,
120            -3.831_192_481_833_385_5e-6,
121            4.063_237_480_216_934e-5,
122            0.999_999_999_167_166_1,
123        ])
124        .transpose();
125        let act = nut.nutation_matrix(obl);
126        assert_approx_eq!(act, exp, rtol <= 1e-12);
127    }
128}
lox_frames/iers/nutation.rs

lox_frames/iers/
nutation.rs
