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    /// Computes nutation using the IAU 2000A model.
57    pub fn iau2000a(time: Time<Tdb>) -> Nutation {
58        let t = time.centuries_since_j2000();
59        let luni_solar_args = DelaunayArguments {
60            l: l_iers03(t),
61            lp: lp_mhb2000(t),
62            f: f_iers03(t),
63            d: d_mhb2000_luni_solar(t),
64            om: omega_iers03(t),
65        };
66
67        let planetary_args = DelaunayArguments {
68            l: l_mhb2000(t),
69            lp: 0.0.rad(), // unused
70            f: f_mhb2000(t),
71            d: d_mhb2000_planetary(t),
72            om: omega_mhb2000(t),
73        };
74
75        luni_solar_nutation(t, &luni_solar_args, &luni_solar::COEFFICIENTS)
76            + planetary_nutation(t, planetary_args)
77    }
78}
79
80fn planetary_nutation(
81    centuries_since_j2000_tdb: JulianCenturies,
82    args: DelaunayArguments,
83) -> Nutation {
84    let (dpsi, deps) = planetary::COEFFICIENTS
85        .iter()
86        // The coefficients are given by descending magnitude but folded by ascending
87        // magnitude to minimise floating-point error.
88        .rev()
89        .fold((0.0, 0.0), |(mut dpsi, mut deps), coeff| {
90            // Form argument for current term.
91            let arg = (coeff.l * args.l
92                + coeff.f * args.f
93                + coeff.d * args.d
94                + coeff.om * args.om
95                + coeff.mercury * mercury_l_iers03(centuries_since_j2000_tdb)
96                + coeff.venus * venus_l_iers03(centuries_since_j2000_tdb)
97                + coeff.earth * earth_l_iers03(centuries_since_j2000_tdb)
98                + coeff.mars * mars_l_iers03(centuries_since_j2000_tdb)
99                + coeff.jupiter * jupiter_l_iers03(centuries_since_j2000_tdb)
100                + coeff.saturn * saturn_l_iers03(centuries_since_j2000_tdb)
101                + coeff.uranus * uranus_l_iers03(centuries_since_j2000_tdb)
102                + coeff.neptune * neptune_l_mhb2000(centuries_since_j2000_tdb)
103                + coeff.pa * pa_iers03(centuries_since_j2000_tdb))
104            .mod_two_pi_signed();
105
106            // Accumulate current term.
107            let sin_arg = arg.sin();
108            let cos_arg = arg.cos();
109            dpsi += coeff.sin_psi * sin_arg + coeff.cos_psi * cos_arg;
110            deps += coeff.sin_eps * sin_arg + coeff.cos_eps * cos_arg;
111
112            (dpsi, deps)
113        });
114
115    Nutation {
116        dpsi: (dpsi * 1e-1).uas(),
117        deps: (deps * 1e-1).uas(),
118    }
119}
120
121#[cfg(test)]
122/// All fixtures and assertion values were generated using the ERFA C library unless otherwise
123/// stated.
124mod tests {
125    use lox_core::units::AngleUnits;
126    use lox_test_utils::assert_approx_eq;
127    use lox_time::{Time, time_scales::Tdb};
128
129    use crate::iers::nutation::Nutation;
130
131    #[test]
132    fn test_nutation_iau2000a() {
133        let time = Time::from_two_part_julian_date(Tdb, 2400000.5, 53736.0);
134        let expected = Nutation {
135            dpsi: -9.630_909_107_115_518e-6.rad(),
136            deps: 4.063_239_174_001_679e-5.rad(),
137        };
138        let actual = Nutation::iau2000a(time);
139        assert_approx_eq!(expected, actual, rtol <= 1e-13);
140    }
141}
lox_frames/iers/nutation/iau2000/iau2000a.rs

lox_frames/iers/nutation/iau2000/
iau2000a.rs
