Implement IAU1980 nutation

crates/lox_core/src/bodies.rs

@@ -15,6 +15,7 @@ mod generated;
 pub use generated::*;
 
 pub mod fundamental;
+pub mod nutation;
 
 #[derive(Clone, Copy, Debug, Eq, PartialEq)]
 #[repr(transparent)]

crates/lox_core/src/bodies/fundamental.rs

@@ -9,12 +9,12 @@
 //! Functions for calculating fundamental astronomical parameters as specified by IERS Conventions
 //! (2003).
 
+use crate::math::{arcsec_to_rad, arcsec_to_rad_two_pi, normalize_two_pi};
+use crate::types::Radians;
 use std::f64::consts::TAU;
 
 use super::{Earth, Jupiter, Mars, Mercury, Moon, Neptune, Saturn, Sun, Uranus, Venus};
 
-type Radians = f64;
-
 /// Strictly TDB, TT is sufficient for most applications.
 type TDBJulianCenturiesSinceJ2000 = f64;
 
@@ -35,7 +35,7 @@ pub fn mean_moon_sun_elongation(t: TDBJulianCenturiesSinceJ2000) -> Radians {
             -0.00003169,
         ],
     );
-    arcsec_to_rad_wrapping(arcsec)
+    arcsec_to_rad_two_pi(arcsec)
 }
 
 pub trait MeanAnomaly {
@@ -54,7 +54,7 @@ impl MeanAnomaly for Sun {
                 -0.00001149,
             ],
         );
-        arcsec_to_rad_wrapping(arcsec)
+        arcsec_to_rad_two_pi(arcsec)
     }
 }
 
@@ -70,7 +70,7 @@ impl MeanAnomaly for Moon {
                 -0.00024470,
             ],
         );
-        arcsec_to_rad_wrapping(arcsec)
+        arcsec_to_rad_two_pi(arcsec)
     }
 }
 
@@ -90,15 +90,15 @@ impl Moon {
                 0.00000417,
             ],
         );
-        arcsec_to_rad_wrapping(arcsec)
+        arcsec_to_rad_two_pi(arcsec)
     }
 
     pub fn ascending_node_mean_longitude(&self, t: TDBJulianCenturiesSinceJ2000) -> Radians {
         let arcsec = fast_polynomial::poly_array(
             t,
             &[450160.398036, -6962890.5431, 7.4722, 0.007702, -0.00005939],
         );
-        arcsec_to_rad_wrapping(arcsec)
+        arcsec_to_rad_two_pi(arcsec)
     }
 }
 
@@ -162,20 +162,6 @@ impl MeanLongitude for Uranus {
     }
 }
 
-const ARCSECONDS_IN_CIRCLE: f64 = 360.0 * 60.0 * 60.0;
-
-const RADIANS_IN_ARCSECOND: Radians = TAU / ARCSECONDS_IN_CIRCLE;
-
-#[inline]
-fn arcsec_to_rad_wrapping(arcsec: f64) -> Radians {
-    arcsec_to_rad(arcsec % ARCSECONDS_IN_CIRCLE)
-}
-
-#[inline]
-fn arcsec_to_rad(arcsec: f64) -> Radians {
-    arcsec * RADIANS_IN_ARCSECOND
-}
-
 #[cfg(test)]
 #[allow(clippy::approx_constant)]
 mod tests {

crates/lox_core/src/bodies/nutation.rs

@@ -0,0 +1,74 @@
+use crate::bodies::nutation::iau1980::nutation_iau1980;
+use crate::time::epochs::Epoch;
+use crate::types::Radians;
+
+mod iau1980;
+
+pub enum Model {
+    IAU1980,
+    IAU2000A,
+    IAU2000B,
+    IAU2006A,
+}
+
+pub struct Nutation {
+    /// δψ
+    pub longitude: Radians,
+    /// δε
+    pub obliquity: Radians,
+}
+
+impl Default for Nutation {
+    fn default() -> Self {
+        Self {
+            longitude: 0.0,
+            obliquity: 0.0,
+        }
+    }
+}
+
+struct Coefficients {
+    /// `l`.
+    l: f64,
+    /// `l'`.
+    lp: f64,
+    /// `F`.
+    f: f64,
+    /// `D`.
+    d: f64,
+    /// `Ω`.
+    om: f64,
+    /// Longitude sine.
+    long_sin_1: f64,
+    long_sin_t: f64,
+    /// Obliquity cosine.
+    obl_cos_1: f64,
+    obl_cos_t: f64,
+}
+
+/// The interval between J2000 and a given Julian date.
+pub type JulianInterval = f64;
+
+pub fn nutation(model: Model, epoch: Epoch) -> Nutation {
+    // TODO: This call is placeholder. We need to ensure correct calculation of the Julian interval
+    // from the epoch.
+    let t: JulianInterval = epoch.j2000();
+    match model {
+        Model::IAU1980 => nutation_iau1980(t),
+        Model::IAU2000A => nutation_iau2000a(t),
+        Model::IAU2000B => nutation_iau2000b(t),
+        Model::IAU2006A => nutation_iau2006a(t),
+    }
+}
+
+fn nutation_iau2000a(t: JulianInterval) -> Nutation {
+    todo!()
+}
+
+fn nutation_iau2000b(t: JulianInterval) -> Nutation {
+    todo!()
+}
+
+fn nutation_iau2006a(t: JulianInterval) -> Nutation {
+    todo!()
+}

crates/lox_core/src/lib.rs

@@ -13,3 +13,6 @@ pub mod errors;
 pub mod frames;
 pub mod time;
 pub mod two_body;
+pub mod types;
+
+pub(crate) mod math;

crates/lox_core/src/math.rs

@@ -0,0 +1,31 @@
+use crate::types::{Arcsec, MilliArcsec, Radians};
+use std::f64::consts::{PI, TAU};
+
+/// Module math provides common mathematical functions shared by many parts of the library.
+
+/// Normalizes an angle `a` to the range [center-π, center+π].
+pub(crate) fn normalize_two_pi(a: Radians, center: Radians) -> Radians {
+    a - 2.0 * PI * ((a + PI - center) / (2.0 * PI)).floor()
+}
+
+pub(crate) const ARCSECONDS_IN_CIRCLE: f64 = 360.0 * 60.0 * 60.0;
+
+pub(crate) const RADIANS_IN_ARCSECOND: Radians = TAU / ARCSECONDS_IN_CIRCLE;
+
+pub(crate) const RADIANS_IN_MILLIARCSECOND: Radians = RADIANS_IN_ARCSECOND / 1000.0;
+
+/// Converts arcseconds to radians, modulo 2π.
+#[inline]
+pub(crate) fn arcsec_to_rad_two_pi(arcsec: Arcsec) -> Radians {
+    arcsec_to_rad(arcsec % ARCSECONDS_IN_CIRCLE)
+}
+
+#[inline]
+pub(crate) fn arcsec_to_rad(arcsec: Arcsec) -> Radians {
+    arcsec * RADIANS_IN_ARCSECOND
+}
+
+#[inline]
+pub(crate) fn milliarcsec_to_rad(mas: MilliArcsec) -> Radians {
+    mas * RADIANS_IN_MILLIARCSECOND
+}

crates/lox_core/src/two_body/elements.rs

@@ -6,6 +6,7 @@
  * file, you can obtain one at https://mozilla.org/MPL/2.0/.
  */
 
+use crate::math::normalize_two_pi;
 use float_eq::float_eq;
 use glam::{DMat3, DVec3};
 use std::f64::consts::PI;
@@ -113,10 +114,6 @@ pub fn cartesian_to_keplerian(grav_param: f64, pos: DVec3, vel: DVec3) -> Kepler
     )
 }
 
-fn normalize_two_pi(a: f64, center: f64) -> f64 {
-    a - 2.0 * PI * ((a + PI - center) / (2.0 * PI)).floor()
-}
-
 fn mod_two_pi(a: f64) -> f64 {
     let w = a % (2.0 * PI);
     if w < 0.0 {

crates/lox_core/src/types.rs

@@ -0,0 +1,7 @@
+///! Contains common types and aliases useful in many parts of the library.
+
+pub type MilliArcsec = f64;
+
+pub type Arcsec = f64;
+
+pub type Radians = f64;