Implement IAU 2006A nutation model (#30)

crates/lox_core/src/bodies/nutation.rs

@@ -3,13 +3,15 @@ use std::ops::Add;
 use crate::bodies::nutation::iau1980::nutation_iau1980;
 use crate::bodies::nutation::iau2000::nutation_iau2000a;
 use crate::bodies::nutation::iau2000::nutation_iau2000b;
+use crate::bodies::nutation::iau2006::nutation_iau2006a;
 use crate::math::RADIANS_IN_ARCSECOND;
 use crate::time::epochs::Epoch;
-use crate::time::intervals::{tdb_julian_centuries_since_j2000, TDBJulianCenturiesSinceJ2000};
+use crate::time::intervals::tdb_julian_centuries_since_j2000;
 use crate::types::Radians;
 
 mod iau1980;
 mod iau2000;
+mod iau2006;
 
 /// The supported IAU nutation models.
 pub enum Model {
@@ -61,10 +63,6 @@ pub fn nutation(model: Model, epoch: Epoch) -> Nutation {
     }
 }
 
-fn nutation_iau2006a(_t: TDBJulianCenturiesSinceJ2000) -> Nutation {
-    todo!()
-}
-
 const RADIANS_IN_POINT_ONE_MILLIARCSECOND: Radians = RADIANS_IN_ARCSECOND / 1e4;
 
 /// Units of 0.1 mas are returned by certain nutation calculations before being converted to
@@ -132,6 +130,18 @@ mod tests {
         assert_float_eq!(expected.obliquity, actual.obliquity, rel <= TOLERANCE);
     }
 
+    #[test]
+    fn test_nutation_iau2006a() {
+        let epoch = Epoch::j2000(TimeScale::TT);
+        let expected = Nutation {
+            longitude: -0.00006754425598969513,
+            obliquity: -0.00002797083119237414,
+        };
+        let actual = nutation(Model::IAU2006A, epoch);
+        assert_float_eq!(expected.longitude, actual.longitude, rel <= TOLERANCE);
+        assert_float_eq!(expected.obliquity, actual.obliquity, rel <= TOLERANCE);
+    }
+
     #[test]
     fn test_point1_milliarcsec_to_rad() {
         assert_float_eq!(point1_milliarcsec_to_rad(0.0), 0.0, abs <= TOLERANCE);

crates/lox_core/src/bodies/nutation/iau2006.rs

@@ -0,0 +1,67 @@
+/*
+ * Copyright (c) 2023. Helge Eichhorn and the LOX contributors
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, you can obtain one at https://mozilla.org/MPL/2.0/.
+ */
+
+use crate::bodies::nutation::iau2000::nutation_iau2000a;
+use crate::bodies::nutation::Nutation;
+use crate::bodies::Earth;
+use crate::time::intervals::TDBJulianCenturiesSinceJ2000;
+
+/// The IAU 2000A nutation model adjusted to match the IAU 2006 precession model per
+/// Wallace & Capitaine, 2006.
+pub fn nutation_iau2006a(t: TDBJulianCenturiesSinceJ2000) -> Nutation {
+    let mut nutation = nutation_iau2000a(t);
+    let j2_correction = Earth::j2_correction_factor(t);
+
+    nutation.longitude += nutation.longitude * (0.4697e-6 + j2_correction);
+    nutation.obliquity += nutation.obliquity * j2_correction;
+
+    nutation
+}
+
+impl Earth {
+    /// Factor correcting for secular variation of J₂.
+    #[inline]
+    fn j2_correction_factor(t: TDBJulianCenturiesSinceJ2000) -> f64 {
+        -2.7774e-6 * t
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use float_eq::assert_float_eq;
+
+    use crate::time::intervals::TDBJulianCenturiesSinceJ2000;
+
+    use super::nutation_iau2006a;
+
+    const TOLERANCE: f64 = 1e-11;
+
+    #[test]
+    fn test_nutation_iau2006a_jd0() {
+        let jd0: TDBJulianCenturiesSinceJ2000 = -67.11964407939767;
+        let actual = nutation_iau2006a(jd0);
+        assert_float_eq!(0.00000737285641780423, actual.longitude, rel <= TOLERANCE);
+        assert_float_eq!(0.00004132905772755788, actual.obliquity, rel <= TOLERANCE);
+    }
+
+    #[test]
+    fn test_nutation_iau2006a_j2000() {
+        let j2000: TDBJulianCenturiesSinceJ2000 = 0.0;
+        let actual = nutation_iau2006a(j2000);
+        assert_float_eq!(-0.00006754425598969513, actual.longitude, rel <= TOLERANCE);
+        assert_float_eq!(-0.00002797083119237414, actual.obliquity, rel <= TOLERANCE);
+    }
+
+    #[test]
+    fn test_nutation_iau2006a_j2100() {
+        let j2100: TDBJulianCenturiesSinceJ2000 = 1.0;
+        let actual = nutation_iau2006a(j2100);
+        assert_float_eq!(0.00001585983730501046, actual.longitude, rel <= TOLERANCE);
+        assert_float_eq!(0.00004162315218980551, actual.obliquity, rel <= TOLERANCE);
+    }
+}