lox-frames 0.1.0-alpha.16

Reference frame transformations for the Lox ecosystem
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

// SPDX-FileCopyrightText: 2023 Angus Morrison <github@angus-morrison.com>
// SPDX-FileCopyrightText: 2024 Helge Eichhorn <git@helgeeichhorn.de>
//
// SPDX-License-Identifier: MPL-2.0

//! Module nutation exposes a function for calculating Earth nutation using a number of IAU nutation
//! models.

use std::ops::Add;

use lox_core::glam::DMat3;
use lox_core::units::Angle;
use lox_test_utils::ApproxEq;
use lox_time::{
    Time,
    time_scales::{Tdb, Tt},
};

use crate::iers::{Corrections, Iau2000Model, ReferenceSystem, ecliptic::MeanObliquity};

mod iau1980;
mod iau2000;
mod iau2006;

impl ReferenceSystem {
    /// Computes nutation for the active IERS convention.
    pub fn nutation(&self, time: Time<Tdb>) -> Nutation {
        match self {
            ReferenceSystem::Iers1996 => Nutation::iau1980(time),
            ReferenceSystem::Iers2003(model) => match model {
                Iau2000Model::A => Nutation::iau2000a(time),
                Iau2000Model::B => Nutation::iau2000b(time),
            },
            ReferenceSystem::Iers2010 => Nutation::iau2006a(time),
        }
    }

    /// Computes the nutation rotation matrix, applying the given corrections.
    pub fn nutation_matrix(&self, time: Time<Tdb>, corr: Corrections) -> DMat3 {
        match self {
            ReferenceSystem::Iers1996 => {
                let mut nut = self.nutation(time);
                nut += corr;
                nut.nutation_matrix(self.mean_obliquity(time.with_scale(Tt)))
            }
            ReferenceSystem::Iers2003(_) | ReferenceSystem::Iers2010 => {
                let tt = time.with_scale(Tt);
                let epsa = self.mean_obliquity(tt);
                let mut nut = self.nutation(time);
                let rpb = self.bias_precession_matrix(tt);

                nut += self.ecliptic_corrections(corr, nut, epsa, rpb);

                nut.nutation_matrix(epsa)
            }
        }
    }
}

/// Nutation components with respect to some ecliptic of date.
#[derive(Debug, Default, Clone, Copy, PartialEq, ApproxEq)]
pub struct Nutation {
    /// δψ
    pub dpsi: Angle,
    /// δε
    pub deps: Angle,
}

impl Nutation {
    /// Builds the nutation rotation matrix from the mean obliquity of the ecliptic.
    pub fn nutation_matrix(&self, epsa: MeanObliquity) -> DMat3 {
        let epsa = epsa.0;
        let rot1 = epsa.rotation_x();
        let rot2 = (-self.dpsi).rotation_z();
        let rot3 = (-(self.deps + epsa)).rotation_x();
        rot3 * rot2 * rot1
    }
}

impl Add for Nutation {
    type Output = Self;

    fn add(self, rhs: Self) -> Self::Output {
        Self {
            dpsi: self.dpsi + rhs.dpsi,
            deps: self.deps + rhs.deps,
        }
    }
}

impl Add<&Self> for Nutation {
    type Output = Self;

    fn add(self, rhs: &Self) -> Self::Output {
        Nutation {
            dpsi: self.dpsi + rhs.dpsi,
            deps: self.deps + rhs.deps,
        }
    }
}

#[cfg(test)]
mod tests {
    use lox_core::units::AngleUnits;
    use lox_test_utils::assert_approx_eq;

    use super::*;

    #[test]
    fn test_nutation_matrix() {
        let nut = Nutation {
            dpsi: -9.630_909_107_115_582e-6.rad(),
            deps: 4.063_239_174_001_679e-5.rad(),
        };
        let obl = MeanObliquity(Angle::new(0.409_078_976_335_651));
        let exp = DMat3::from_cols_array(&[
            0.999_999_999_953_622_8,
            8.836_239_320_236_25e-6,
            3.830_833_447_458_252e-6,
            -8.836_083_657_016_69e-6,
            0.999_999_999_135_465_4,
            -4.063_240_865_361_857_4e-5,
            -3.831_192_481_833_385_5e-6,
            4.063_237_480_216_934e-5,
            0.999_999_999_167_166_1,
        ])
        .transpose();
        let act = nut.nutation_matrix(obl);
        assert_approx_eq!(act, exp, rtol <= 1e-12);
    }
}