lox_frames/iers/nutation/iau2000/
iau2000a.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
6use lox_core::types::units::JulianCenturies;
7use lox_core::units::AngleUnits;
8use lox_time::Time;
9use lox_time::julian_dates::JulianDate;
10use lox_time::time_scales::Tdb;
11
12use crate::iers::fundamental::iers03::{
13    earth_l_iers03, f_iers03, jupiter_l_iers03, l_iers03, mars_l_iers03, mercury_l_iers03,
14    omega_iers03, pa_iers03, saturn_l_iers03, uranus_l_iers03, venus_l_iers03,
15};
16use crate::iers::fundamental::mhb2000::{
17    d_mhb2000_luni_solar, d_mhb2000_planetary, f_mhb2000, l_mhb2000, lp_mhb2000, neptune_l_mhb2000,
18    omega_mhb2000,
19};
20use crate::iers::nutation::Nutation;
21use crate::iers::nutation::iau2000::{DelaunayArguments, luni_solar_nutation};
22
23mod luni_solar;
24mod planetary;
25
26struct PlanetaryCoefficients {
27    /// Coefficients of l, F, D and Ω.
28    l: f64,
29    f: f64,
30    d: f64,
31    om: f64,
32
33    /// Planetary longitude coefficients.
34    mercury: f64,
35    venus: f64,
36    earth: f64,
37    mars: f64,
38    jupiter: f64,
39    saturn: f64,
40    uranus: f64,
41    neptune: f64,
42
43    /// Coefficient of general precession.
44    pa: f64,
45
46    /// Longitude coefficients.
47    sin_psi: f64,
48    cos_psi: f64,
49
50    /// Obliquity coefficients.
51    sin_eps: f64,
52    cos_eps: f64,
53}
54
55impl Nutation {
56    pub fn iau2000a(time: Time<Tdb>) -> Nutation {
57        let t = time.centuries_since_j2000();
58        let luni_solar_args = DelaunayArguments {
59            l: l_iers03(t),
60            lp: lp_mhb2000(t),
61            f: f_iers03(t),
62            d: d_mhb2000_luni_solar(t),
63            om: omega_iers03(t),
64        };
65
66        let planetary_args = DelaunayArguments {
67            l: l_mhb2000(t),
68            lp: 0.0.rad(), // unused
69            f: f_mhb2000(t),
70            d: d_mhb2000_planetary(t),
71            om: omega_mhb2000(t),
72        };
73
74        luni_solar_nutation(t, &luni_solar_args, &luni_solar::COEFFICIENTS)
75            + planetary_nutation(t, planetary_args)
76    }
77}
78
79fn planetary_nutation(
80    centuries_since_j2000_tdb: JulianCenturies,
81    args: DelaunayArguments,
82) -> Nutation {
83    let (dpsi, deps) = planetary::COEFFICIENTS
84        .iter()
85        // The coefficients are given by descending magnitude but folded by ascending
86        // magnitude to minimise floating-point error.
87        .rev()
88        .fold((0.0, 0.0), |(mut dpsi, mut deps), coeff| {
89            // Form argument for current term.
90            let arg = (coeff.l * args.l
91                + coeff.f * args.f
92                + coeff.d * args.d
93                + coeff.om * args.om
94                + coeff.mercury * mercury_l_iers03(centuries_since_j2000_tdb)
95                + coeff.venus * venus_l_iers03(centuries_since_j2000_tdb)
96                + coeff.earth * earth_l_iers03(centuries_since_j2000_tdb)
97                + coeff.mars * mars_l_iers03(centuries_since_j2000_tdb)
98                + coeff.jupiter * jupiter_l_iers03(centuries_since_j2000_tdb)
99                + coeff.saturn * saturn_l_iers03(centuries_since_j2000_tdb)
100                + coeff.uranus * uranus_l_iers03(centuries_since_j2000_tdb)
101                + coeff.neptune * neptune_l_mhb2000(centuries_since_j2000_tdb)
102                + coeff.pa * pa_iers03(centuries_since_j2000_tdb))
103            .mod_two_pi_signed();
104
105            // Accumulate current term.
106            let sin_arg = arg.sin();
107            let cos_arg = arg.cos();
108            dpsi += coeff.sin_psi * sin_arg + coeff.cos_psi * cos_arg;
109            deps += coeff.sin_eps * sin_arg + coeff.cos_eps * cos_arg;
110
111            (dpsi, deps)
112        });
113
114    Nutation {
115        dpsi: (dpsi * 1e-1).uas(),
116        deps: (deps * 1e-1).uas(),
117    }
118}
119
120#[cfg(test)]
121/// All fixtures and assertion values were generated using the ERFA C library unless otherwise
122/// stated.
123mod tests {
124    use lox_core::units::AngleUnits;
125    use lox_test_utils::assert_approx_eq;
126    use lox_time::{Time, time_scales::Tdb};
127
128    use crate::iers::nutation::Nutation;
129
130    #[test]
131    fn test_nutation_iau2000a() {
132        let time = Time::from_two_part_julian_date(Tdb, 2400000.5, 53736.0);
133        let expected = Nutation {
134            dpsi: -9.630_909_107_115_518e-6.rad(),
135            deps: 4.063_239_174_001_679e-5.rad(),
136        };
137        let actual = Nutation::iau2000a(time);
138        assert_approx_eq!(expected, actual, rtol <= 1e-13);
139    }
140}
lox_frames/iers/nutation/iau2000/iau2000a.rs

lox_frames/iers/nutation/iau2000/
iau2000a.rs
