refactor(lox-earth): refactor nutation functions

helgee · helgee · commit 72b632653e01 · 2025-11-09T11:10:34.000+01:00
diff --git a/crates/lox-earth/src/nutation.rs b/crates/lox-earth/src/nutation.rs
@@ -8,28 +8,20 @@
 
 use std::ops::Add;
 
+use glam::DMat3;
 use lox_test_utils::ApproxEq;
-use lox_time::Time;
-use lox_time::julian_dates::JulianDate;
-use lox_time::time_scales::Tdb;
+use lox_time::{
+    Time,
+    time_scales::{Tdb, Tt},
+};
 use lox_units::Angle;
 
-use crate::nutation::iau2000::nutation_iau2000a;
-use crate::nutation::iau2000::nutation_iau2000b;
-use crate::nutation::iau2006::nutation_iau2006a;
+use crate::ecliptic::MeanObliquity;
 
 mod iau1980;
 mod iau2000;
 mod iau2006;
 
-/// The supported IAU nutation models.
-pub enum Model {
-    IAU1980,
-    IAU2000A,
-    IAU2000B,
-    IAU2006A,
-}
-
 /// Nutation components with respect to some ecliptic of date.
 #[derive(Debug, Default, Clone, PartialEq, ApproxEq)]
 pub struct Nutation {
@@ -39,6 +31,22 @@ pub struct Nutation {
     pub obliquity: Angle,
 }
 
+impl Nutation {
+    fn nutation_matrix(&self, mean_obliquity: MeanObliquity) -> DMat3 {
+        let rot1 = mean_obliquity.0.rotation_x();
+        let rot2 = (-self.longitude).rotation_z();
+        let rot3 = (-(self.obliquity + mean_obliquity.0)).rotation_x();
+        rot3 * rot2 * rot1
+    }
+
+    pub fn nutation_matrix_iau1980(tdb: Time<Tdb>) -> DMat3 {
+        let tt = tdb.with_scale(Tt);
+        let nut = Self::iau1980(tdb);
+        let obl = MeanObliquity::iau1980(tt);
+        nut.nutation_matrix(obl)
+    }
+}
+
 impl Add for Nutation {
     type Output = Self;
 
@@ -61,67 +69,52 @@ impl Add<&Self> for Nutation {
     }
 }
 
-/// Calculate nutation coefficients at `time` using the given [Model].
-pub fn nutation(model: Model, time: Time<Tdb>) -> Nutation {
-    let t = time.centuries_since_j2000();
-    match model {
-        Model::IAU1980 => Nutation::iau1980(t),
-        Model::IAU2000A => nutation_iau2000a(t),
-        Model::IAU2000B => nutation_iau2000b(t),
-        Model::IAU2006A => nutation_iau2006a(t),
-    }
-}
-
 #[cfg(test)]
 mod tests {
     use lox_test_utils::assert_approx_eq;
     use lox_units::AngleUnits;
 
     use super::*;
 
-    const TOLERANCE: f64 = 1e-12;
-
-    #[test]
-    fn test_nutation_iau1980() {
-        let time = Time::j2000(Tdb);
-        let expected = Nutation {
-            longitude: -0.00006750247617532478.rad(),
-            obliquity: -0.00002799221238377013.rad(),
-        };
-        let actual = nutation(Model::IAU1980, time);
-        assert_approx_eq!(expected, actual, rtol <= TOLERANCE);
-    }
-
     #[test]
-    fn test_nutation_iau2000a() {
-        let time = Time::j2000(Tdb);
-        let expected = Nutation {
-            longitude: -0.00006754422426417299.rad(),
-            obliquity: -0.00002797083119237414.rad(),
+    fn test_nutation_matrix() {
+        let nut = Nutation {
+            longitude: -9.630_909_107_115_582e-6.rad(),
+            obliquity: 4.063_239_174_001_679e-5.rad(),
         };
-        let actual = nutation(Model::IAU2000A, time);
-        assert_approx_eq!(expected, actual, rtol <= TOLERANCE);
+        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);
     }
 
     #[test]
-    fn test_nutation_iau2000b() {
-        let time = Time::j2000(Tdb);
-        let expected = Nutation {
-            longitude: -0.00006754261253992235.rad(),
-            obliquity: -0.00002797092331098565.rad(),
-        };
-        let actual = nutation(Model::IAU2000B, time);
-        assert_approx_eq!(expected, actual, rtol <= TOLERANCE);
-    }
-
-    #[test]
-    fn test_nutation_iau2006a() {
-        let time = Time::j2000(Tdb);
-        let expected = Nutation {
-            longitude: -0.00006754425598969513.rad(),
-            obliquity: -0.00002797083119237414.rad(),
-        };
-        let actual = nutation(Model::IAU2006A, time);
-        assert_approx_eq!(expected, actual, rtol <= TOLERANCE);
+    fn test_nutation_matrix_iau1980() {
+        let time = Time::from_two_part_julian_date(Tdb, 2400000.5, 53736.0);
+        let exp = DMat3::from_cols_array(&[
+            0.9999999999534999268,
+            0.8847935789636432161e-5,
+            0.3835906502164019142e-5,
+            -0.8847780042583435924e-5,
+            0.9999999991366569963,
+            -0.4060052702727130809e-4,
+            -0.3836265729708478796e-5,
+            0.4060049308612638555e-4,
+            0.9999999991684415129,
+        ])
+        .transpose();
+        let act = Nutation::nutation_matrix_iau1980(time);
+        assert_approx_eq!(act, exp, rtol <= 1e-12);
     }
 }
diff --git a/crates/lox-earth/src/nutation/iau1980.rs b/crates/lox-earth/src/nutation/iau1980.rs
@@ -6,6 +6,7 @@
 use std::f64::consts::TAU;
 
 use lox_core::types::units::{Arcseconds, JulianCenturies};
+use lox_time::{Time, julian_dates::JulianDate, time_scales::Tdb};
 use lox_units::{Angle, AngleUnits};
 
 use crate::nutation::Nutation;
@@ -28,12 +29,13 @@ struct Coefficients {
 }
 
 impl Nutation {
-    pub fn iau1980(centuries_since_j2000_tdb: JulianCenturies) -> Self {
-        let l = l(centuries_since_j2000_tdb);
-        let lp = lp(centuries_since_j2000_tdb);
-        let f = f(centuries_since_j2000_tdb);
-        let d = d(centuries_since_j2000_tdb);
-        let om = omega(centuries_since_j2000_tdb);
+    pub fn iau1980(time: Time<Tdb>) -> Self {
+        let t = time.centuries_since_j2000();
+        let l = l(t);
+        let lp = lp(t);
+        let f = f(t);
+        let d = d(t);
+        let om = omega(t);
 
         let (dpsi, deps) = COEFFICIENTS
             .iter()
@@ -45,8 +47,8 @@ impl Nutation {
                 let arg = coeff.l * l + coeff.lp * lp + coeff.f * f + coeff.d * d + coeff.om * om;
 
                 // Accumulate current term.
-                let sin = coeff.sin_psi + coeff.sin_psi_t * centuries_since_j2000_tdb;
-                let cos = coeff.cos_eps + coeff.cos_eps_t * centuries_since_j2000_tdb;
+                let sin = coeff.sin_psi + coeff.sin_psi_t * t;
+                let cos = coeff.cos_eps + coeff.cos_eps_t * t;
                 dpsi += sin * arg.sin();
                 deps += cos * arg.cos();
 
@@ -227,46 +229,19 @@ const COEFFICIENTS: [Coefficients; 106] = [
 // @formatter:on
 
 #[cfg(test)]
-/// All fixtures and assertion values were generated using the ERFA C library unless otherwise
-/// stated.
 mod tests {
-    use lox_core::types::units::JulianCenturies;
     use lox_test_utils::assert_approx_eq;
 
     use super::*;
 
-    const TOLERANCE: f64 = 1e-12;
-
-    #[test]
-    fn test_nutation_iau1980_jd0() {
-        let jd0: JulianCenturies = -67.11964407939767;
-        let actual = Nutation::iau1980(jd0);
-        let expected = Nutation {
-            longitude: 0.00000693404778664026.rad(),
-            obliquity: 0.00004131255061383108.rad(),
-        };
-        assert_approx_eq!(expected, actual, rtol <= TOLERANCE);
-    }
-
-    #[test]
-    fn test_nutation_iau1980_j2000() {
-        let j2000: JulianCenturies = 0.0;
-        let actual = Nutation::iau1980(j2000);
-        let expected = Nutation {
-            longitude: -0.00006750247617532478.rad(),
-            obliquity: -0.00002799221238377013.rad(),
-        };
-        assert_approx_eq!(expected, actual, rtol <= TOLERANCE);
-    }
-
     #[test]
-    fn test_nutation_iau1980_j2100() {
-        let j2100: JulianCenturies = 1.0;
-        let actual = Nutation::iau1980(j2100);
+    fn test_nutation_iau1980() {
+        let time = Time::from_two_part_julian_date(Tdb, 2400000.5, 53736.0);
+        let actual = Nutation::iau1980(time);
         let expected = Nutation {
-            longitude: 0.00001584138015187132.rad(),
-            obliquity: 0.00004158958379918889.rad(),
+            longitude: -9.643_658_353_226_563e-6.rad(),
+            obliquity: 4.060_051_006_879_713e-5.rad(),
         };
-        assert_approx_eq!(expected, actual, rtol <= TOLERANCE);
+        assert_approx_eq!(expected, actual, rtol <= 1e-13);
     }
 }
diff --git a/crates/lox-earth/src/nutation/iau2000.rs b/crates/lox-earth/src/nutation/iau2000.rs
@@ -3,8 +3,6 @@
 //
 // SPDX-License-Identifier: MPL-2.0
 
-pub(crate) use iau2000a::nutation_iau2000a;
-pub(crate) use iau2000b::nutation_iau2000b;
 use lox_core::types::units::JulianCenturies;
 use lox_units::{Angle, AngleUnits};
 
diff --git a/crates/lox-earth/src/nutation/iau2000/iau2000a.rs b/crates/lox-earth/src/nutation/iau2000/iau2000a.rs
@@ -9,6 +9,9 @@ use lox_bodies::fundamental::mhb2000::{
 };
 use lox_bodies::*;
 use lox_core::types::units::JulianCenturies;
+use lox_time::Time;
+use lox_time::julian_dates::JulianDate;
+use lox_time::time_scales::Tdb;
 use lox_units::AngleUnits;
 
 use crate::nutation::Nutation;
@@ -46,30 +49,28 @@ struct PlanetaryCoefficients {
     cos_eps: f64,
 }
 
-pub(crate) fn nutation_iau2000a(centuries_since_j2000_tdb: JulianCenturies) -> Nutation {
-    let luni_solar_args = DelaunayArguments {
-        l: Moon.mean_anomaly_iers03(centuries_since_j2000_tdb),
-        lp: Sun.mean_anomaly_mhb2000(centuries_since_j2000_tdb),
-        f: Moon
-            .mean_longitude_minus_ascending_node_mean_longitude_iers03(centuries_since_j2000_tdb),
-        d: mean_moon_sun_elongation_mhb2000_luni_solar(centuries_since_j2000_tdb),
-        om: Moon.ascending_node_mean_longitude_iers03(centuries_since_j2000_tdb),
-    };
-
-    let planetary_args = DelaunayArguments {
-        l: Moon.mean_anomaly_mhb2000(centuries_since_j2000_tdb),
-        lp: 0.0.rad(), // unused
-        f: Moon
-            .mean_longitude_minus_ascending_node_mean_longitude_mhb2000(centuries_since_j2000_tdb),
-        d: mean_moon_sun_elongation_mhb2000_planetary(centuries_since_j2000_tdb),
-        om: Moon.ascending_node_mean_longitude_mhb2000(centuries_since_j2000_tdb),
-    };
-
-    luni_solar_nutation(
-        centuries_since_j2000_tdb,
-        &luni_solar_args,
-        &luni_solar::COEFFICIENTS,
-    ) + planetary_nutation(centuries_since_j2000_tdb, planetary_args)
+impl Nutation {
+    pub fn iau2000a(time: Time<Tdb>) -> Nutation {
+        let t = time.centuries_since_j2000();
+        let luni_solar_args = DelaunayArguments {
+            l: Moon.mean_anomaly_iers03(t),
+            lp: Sun.mean_anomaly_mhb2000(t),
+            f: Moon.mean_longitude_minus_ascending_node_mean_longitude_iers03(t),
+            d: mean_moon_sun_elongation_mhb2000_luni_solar(t),
+            om: Moon.ascending_node_mean_longitude_iers03(t),
+        };
+
+        let planetary_args = DelaunayArguments {
+            l: Moon.mean_anomaly_mhb2000(t),
+            lp: 0.0.rad(), // unused
+            f: Moon.mean_longitude_minus_ascending_node_mean_longitude_mhb2000(t),
+            d: mean_moon_sun_elongation_mhb2000_planetary(t),
+            om: Moon.ascending_node_mean_longitude_mhb2000(t),
+        };
+
+        luni_solar_nutation(t, &luni_solar_args, &luni_solar::COEFFICIENTS)
+            + planetary_nutation(t, planetary_args)
+    }
 }
 
 fn planetary_nutation(
@@ -118,46 +119,20 @@ fn planetary_nutation(
 /// All fixtures and assertion values were generated using the ERFA C library unless otherwise
 /// stated.
 mod tests {
-    use lox_core::types::units::JulianCenturies;
     use lox_test_utils::assert_approx_eq;
+    use lox_time::{Time, time_scales::Tdb};
     use lox_units::AngleUnits;
 
     use crate::nutation::Nutation;
 
-    use super::nutation_iau2000a;
-
-    const TOLERANCE: f64 = 1e-11;
-
-    #[test]
-    fn test_nutation_iau2000a_jd0() {
-        let jd0: JulianCenturies = -67.11964407939767;
-        let expected = Nutation {
-            longitude: 0.00000737147877835653.rad(),
-            obliquity: 0.00004132135467915123.rad(),
-        };
-        let actual = nutation_iau2000a(jd0);
-        assert_approx_eq!(expected, actual, rtol <= TOLERANCE);
-    }
-
-    #[test]
-    fn test_nutation_iau2000a_j2000() {
-        let j2000: JulianCenturies = 0.0;
-        let expected = Nutation {
-            longitude: -0.00006754422426417299.rad(),
-            obliquity: -0.00002797083119237414.rad(),
-        };
-        let actual = nutation_iau2000a(j2000);
-        assert_approx_eq!(expected, actual, rtol <= TOLERANCE);
-    }
-
     #[test]
-    fn test_nutation_iau2000a_j2100() {
-        let j2100: JulianCenturies = 1.0;
+    fn test_nutation_iau2000a() {
+        let time = Time::from_two_part_julian_date(Tdb, 2400000.5, 53736.0);
         let expected = Nutation {
-            longitude: 0.00001585987390484147.rad(),
-            obliquity: 0.00004162326779426948.rad(),
+            longitude: -9.630_909_107_115_518e-6.rad(),
+            obliquity: 4.063_239_174_001_679e-5.rad(),
         };
-        let actual = nutation_iau2000a(j2100);
-        assert_approx_eq!(expected, actual, rtol <= TOLERANCE);
+        let actual = Nutation::iau2000a(time);
+        assert_approx_eq!(expected, actual, rtol <= 1e-13);
     }
 }
diff --git a/crates/lox-earth/src/nutation/iau2000/iau2000b.rs b/crates/lox-earth/src/nutation/iau2000/iau2000b.rs
diff --git a/crates/lox-earth/src/nutation/iau2006.rs b/crates/lox-earth/src/nutation/iau2006.rs
