Implement nutation model IAU 2000B (#29)

crates/lox_core/src/bodies/fundamental.rs

@@ -10,3 +10,4 @@
 
 pub mod iers03;
 pub mod mhb2000;
+pub mod simon1994;

crates/lox_core/src/bodies/fundamental/simon1994.rs

@@ -0,0 +1,154 @@
+//! Functions for calculating fundamental astronomical parameters as proposed by Simon et al.
+//! (1994).
+
+use crate::bodies::{Moon, Sun};
+use crate::math::arcsec_to_rad_two_pi;
+use crate::time::intervals::TDBJulianCenturiesSinceJ2000;
+use crate::types::{Arcsec, Radians};
+
+pub fn mean_moon_sun_elongation_simon1994(t: TDBJulianCenturiesSinceJ2000) -> Radians {
+    let arcsec: Arcsec = fast_polynomial::poly_array(t, &[1072260.70369, 1602961601.2090]);
+    arcsec_to_rad_two_pi(arcsec)
+}
+
+impl Sun {
+    pub fn mean_anomaly_simon1994(&self, t: TDBJulianCenturiesSinceJ2000) -> Radians {
+        let arcsec: Arcsec = fast_polynomial::poly_array(t, &[1287104.79305, 129596581.0481]);
+        arcsec_to_rad_two_pi(arcsec)
+    }
+}
+
+impl Moon {
+    pub fn mean_anomaly_simon1994(&self, t: TDBJulianCenturiesSinceJ2000) -> Radians {
+        let arcsec: Arcsec = fast_polynomial::poly_array(t, &[485868.249036, 1717915923.2178]);
+        arcsec_to_rad_two_pi(arcsec)
+    }
+
+    pub fn mean_argument_of_latitude_simon1994(&self, t: TDBJulianCenturiesSinceJ2000) -> Radians {
+        let arcsec: Arcsec = fast_polynomial::poly_array(t, &[335779.526232, 1739527262.8478]);
+        arcsec_to_rad_two_pi(arcsec)
+    }
+
+    pub fn ascending_node_mean_longitude_simon1994(
+        &self,
+        t: TDBJulianCenturiesSinceJ2000,
+    ) -> Radians {
+        let arcsec: Arcsec = fast_polynomial::poly_array(t, &[450160.398036, -6962890.5431]);
+        arcsec_to_rad_two_pi(arcsec)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use float_eq::assert_float_eq;
+
+    use super::*;
+
+    // Note that all expected values are outputs from the equivalent ERFA functions.
+
+    // Relative error tolerance for float_eq assertions.
+    // This is somewhat loose, being based on observations of how closely our implementations
+    // match ERFA outputs rather than any target tolerance.
+    // See https://github.com/lox-space/lox/pull/23#discussion_r1398485509
+    const TOLERANCE: f64 = 1e-12;
+
+    // Test cases for t.
+    const T_ZERO: TDBJulianCenturiesSinceJ2000 = 0.0;
+    const T_POSITIVE: TDBJulianCenturiesSinceJ2000 = 1.23456789;
+    const T_NEGATIVE: TDBJulianCenturiesSinceJ2000 = -1.23456789;
+
+    #[test]
+    fn test_mean_moon_sun_elongation_simon1994() {
+        assert_float_eq!(
+            mean_moon_sun_elongation_simon1994(T_ZERO),
+            5.198466588650503,
+            rel <= TOLERANCE
+        );
+        assert_float_eq!(
+            mean_moon_sun_elongation_simon1994(T_POSITIVE),
+            5.067187555274916,
+            rel <= TOLERANCE
+        );
+        assert_float_eq!(
+            mean_moon_sun_elongation_simon1994(T_NEGATIVE),
+            -0.953439685154148,
+            rel <= TOLERANCE
+        );
+    }
+
+    #[test]
+    fn test_sun_mean_anomaly_simon1994() {
+        assert_float_eq!(
+            Sun.mean_anomaly_simon1994(T_ZERO),
+            6.24006012692298,
+            rel <= TOLERANCE
+        );
+        assert_float_eq!(
+            Sun.mean_anomaly_simon1994(T_POSITIVE),
+            2.806501115480207,
+            rel <= TOLERANCE
+        );
+        assert_float_eq!(
+            Sun.mean_anomaly_simon1994(T_NEGATIVE),
+            -2.892751475993361,
+            rel <= TOLERANCE
+        );
+    }
+
+    #[test]
+    fn test_moon_mean_anomaly_simon1994() {
+        assert_float_eq!(
+            Moon.mean_anomaly_simon1994(T_ZERO),
+            2.355555743493879,
+            rel <= TOLERANCE
+        );
+        assert_float_eq!(
+            Moon.mean_anomaly_simon1994(T_POSITIVE),
+            5.399393108792649,
+            rel <= TOLERANCE
+        );
+        assert_float_eq!(
+            Moon.mean_anomaly_simon1994(T_NEGATIVE),
+            -0.688281621805333,
+            rel <= TOLERANCE
+        );
+    }
+
+    #[test]
+    fn test_moon_mean_argument_of_latitude_simon1994() {
+        assert_float_eq!(
+            Moon.mean_argument_of_latitude_simon1994(T_ZERO),
+            1.627905081537519,
+            rel <= TOLERANCE
+        );
+        assert_float_eq!(
+            Moon.mean_argument_of_latitude_simon1994(T_POSITIVE),
+            2.076369815616488,
+            rel <= TOLERANCE
+        );
+        assert_float_eq!(
+            Moon.mean_argument_of_latitude_simon1994(T_NEGATIVE),
+            -5.103744959722151,
+            rel <= TOLERANCE
+        );
+    }
+
+    #[test]
+    fn test_moon_ascending_node_mean_longitude_simon1994() {
+        assert_float_eq!(
+            Moon.ascending_node_mean_longitude_simon1994(T_ZERO),
+            2.182439196615671,
+            rel <= TOLERANCE
+        );
+        assert_float_eq!(
+            Moon.ascending_node_mean_longitude_simon1994(T_POSITIVE),
+            -1.793813955913912,
+            rel <= TOLERANCE
+        );
+        assert_float_eq!(
+            Moon.ascending_node_mean_longitude_simon1994(T_NEGATIVE),
+            6.158692349145257,
+            rel <= TOLERANCE
+        );
+    }
+}

crates/lox_core/src/bodies/nutation.rs

@@ -1,14 +1,15 @@
 use std::ops::Add;
 
 use crate::bodies::nutation::iau1980::nutation_iau1980;
-use crate::bodies::nutation::iau2000a::nutation_iau2000a;
+use crate::bodies::nutation::iau2000::nutation_iau2000a;
+use crate::bodies::nutation::iau2000::nutation_iau2000b;
 use crate::math::RADIANS_IN_ARCSECOND;
 use crate::time::epochs::Epoch;
 use crate::time::intervals::{tdb_julian_centuries_since_j2000, TDBJulianCenturiesSinceJ2000};
 use crate::types::Radians;
 
 mod iau1980;
-mod iau2000a;
+mod iau2000;
 
 /// The supported IAU nutation models.
 pub enum Model {
@@ -38,6 +39,17 @@ impl Add for Nutation {
     }
 }
 
+impl Add<&Self> for Nutation {
+    type Output = Self;
+
+    fn add(self, rhs: &Self) -> Self::Output {
+        Nutation {
+            longitude: self.longitude + rhs.longitude,
+            obliquity: self.obliquity + rhs.obliquity,
+        }
+    }
+}
+
 /// Calculate nutation coefficients at `epoch` using the given [Model].
 pub fn nutation(model: Model, epoch: Epoch) -> Nutation {
     let t = tdb_julian_centuries_since_j2000(epoch);
@@ -49,10 +61,6 @@ pub fn nutation(model: Model, epoch: Epoch) -> Nutation {
     }
 }
 
-fn nutation_iau2000b(_t: TDBJulianCenturiesSinceJ2000) -> Nutation {
-    todo!()
-}
-
 fn nutation_iau2006a(_t: TDBJulianCenturiesSinceJ2000) -> Nutation {
     todo!()
 }
@@ -112,6 +120,18 @@ mod tests {
         assert_float_eq!(expected.obliquity, actual.obliquity, rel <= TOLERANCE);
     }
 
+    #[test]
+    fn test_nutation_iau2000b() {
+        let epoch = Epoch::j2000(TimeScale::TT);
+        let expected = Nutation {
+            longitude: -0.00006754261253992235,
+            obliquity: -0.00002797092331098565,
+        };
+        let actual = nutation(Model::IAU2000B, 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/iau2000.rs

@@ -0,0 +1,74 @@
+mod iau2000a;
+mod iau2000b;
+
+use crate::bodies::nutation::{point1_microarcsec_to_rad, Nutation};
+use crate::time::intervals::TDBJulianCenturiesSinceJ2000;
+pub(super) use iau2000a::nutation_iau2000a;
+pub(super) use iau2000b::nutation_iau2000b;
+use std::f64::consts::TAU;
+
+/// IAU 2000A and 2000B use the same structure for luni-solar coefficients.
+struct LuniSolarCoefficients {
+    /// Coefficients of l, l', F, D and Ω.
+    l: f64,
+    lp: f64,
+    f: f64,
+    d: f64,
+    om: f64,
+
+    /// Longitude coefficients.
+    sin_psi: f64,
+    sin_psi_t: f64,
+    cos_psi: f64,
+
+    /// Obliquity coefficients.
+    cos_eps: f64,
+    cos_eps_t: f64,
+    sin_eps: f64,
+}
+
+struct DelaunayArguments {
+    l: f64,
+    lp: f64,
+    f: f64,
+    d: f64,
+    om: f64,
+}
+
+/// Calculate the luni-solar nutation for `t` given `args` and coefficients for either models A or
+/// B.
+fn luni_solar_nutation(
+    t: TDBJulianCenturiesSinceJ2000,
+    args: &DelaunayArguments,
+    coeffs: &[LuniSolarCoefficients],
+) -> Nutation {
+    let mut nutation = coeffs
+        .iter()
+        // The coefficients are given by descending magnitude but folded by ascending
+        // magnitude to minimise floating-point error.
+        .rev()
+        .fold(Nutation::default(), |mut nut, coeff| {
+            // Form argument for current term.
+            let arg = (coeff.l * args.l
+                + coeff.lp * args.lp
+                + coeff.f * args.f
+                + coeff.d * args.d
+                + coeff.om * args.om)
+                % TAU;
+
+            // Accumulate current term.
+            let sin_arg = arg.sin();
+            let cos_arg = arg.cos();
+            nut.longitude +=
+                (coeff.sin_psi + coeff.sin_psi_t * t) * sin_arg + coeff.cos_psi * cos_arg;
+            nut.obliquity +=
+                (coeff.cos_eps + coeff.cos_eps_t * t) * cos_arg + coeff.sin_eps * sin_arg;
+
+            nut
+        });
+
+    nutation.longitude = point1_microarcsec_to_rad(nutation.longitude);
+    nutation.obliquity = point1_microarcsec_to_rad(nutation.obliquity);
+
+    nutation
+}