Calculate celestial to intermediate-frame-of-date matrix (#38)

Exposes a function to compute the celestial to intermediate-frame-of-date matrix given the CIP (X, Y) coordinates and the CIO locator, s. This matrix is the first step in transforming CRS to TRS coordinates.

crates/lox-space/src/prelude.rs

@@ -6,8 +6,6 @@
  * file, you can obtain one at https://mozilla.org/MPL/2.0/.
  */
 
-pub use lox_core::bodies::barycenters::*;
-pub use lox_core::bodies::planets::*;
 pub use lox_core::bodies::*;
 
 pub use lox_core::time::dates::*;

crates/lox_core/src/earth.rs

@@ -11,4 +11,5 @@
 
 mod cio;
 mod cip;
+mod coordinate_transformations;
 mod nutation;

crates/lox_core/src/earth/cio.rs

@@ -6,4 +6,6 @@
  * file, you can obtain one at https://mozilla.org/MPL/2.0/.
  */
 
+//! Module cio exposes functions for calculating the Celestial Intermediate Origin (CIO) locator, s.
+
 mod s06;

crates/lox_core/src/earth/cio/s06.rs

@@ -6,23 +6,26 @@
  * file, you can obtain one at https://mozilla.org/MPL/2.0/.
  */
 
+//! Module s06 exposes a function for calculating the Celestial Intermediate Origin (CIO) locator,
+//! s, using IAU 2006 precession and IAU 2000A nutation.
+
 mod terms;
 
 use crate::bodies::fundamental::iers03::{
     general_accum_precession_in_longitude_iers03, mean_moon_sun_elongation_iers03,
 };
 use crate::bodies::{Earth, Moon, Sun, Venus};
-use crate::earth::cip::xy06::XY;
 use crate::math::arcsec_to_rad;
 use crate::time::intervals::TDBJulianCenturiesSinceJ2000;
 use crate::types::Radians;
+use glam::DVec2;
 
 /// l, l', F, D, Ω, LVe, LE and pA.
 type FundamentalArgs = [Radians; 8];
 
-/// Computes the Celestial Intermediate Origin (CIO) locator s, in radians, given the (X, Y) coordinates of
-/// the Celestial Intermediate Pole (CIP). Based on IAU 2006 precession and IAU 2000A nutation.
-pub fn s(t: TDBJulianCenturiesSinceJ2000, xy: XY) -> Radians {
+/// Computes the Celestial Intermediate Origin (CIO) locator s, in radians, given the (X, Y)
+/// coordinates of the Celestial Intermediate Pole (CIP).
+pub fn s(t: TDBJulianCenturiesSinceJ2000, xy: DVec2) -> Radians {
     let fundamental_args = fundamental_args(t);
     let evaluated_terms = evaluate_terms(&fundamental_args);
     let arcsec = fast_polynomial::poly_array(t, &evaluated_terms);

crates/lox_core/src/earth/cip.rs

@@ -6,4 +6,7 @@
  * file, you can obtain one at https://mozilla.org/MPL/2.0/.
  */
 
+//! Module cip exposes functions for calculating the position of the
+//! Celestial Intermediate Pole (CIP).
+
 pub mod xy06;

crates/lox_core/src/earth/cip/xy06.rs

@@ -6,6 +6,9 @@
  * file, you can obtain one at https://mozilla.org/MPL/2.0/.
  */
 
+//! Module xy06 provides a function to calculate the (X, Y) position of the Celestial Intermediate
+//! Pole (CIP) using the IAU 2006 precession and IAU 2000A nutation models.
+
 mod amplitudes;
 mod luni_solar;
 mod planetary;
@@ -18,10 +21,7 @@ use crate::bodies::{Earth, Jupiter, Mars, Mercury, Moon, Neptune, Saturn, Sun, U
 use crate::math::arcsec_to_rad;
 use crate::time::intervals::TDBJulianCenturiesSinceJ2000;
 use crate::types::Radians;
-
-/// A convenient type for performing batch mathematical operations on X and Y components. This
-/// type may change or become unexported as the needs of upstream components become clearer.
-pub type XY = [f64; 2];
+use glam::DVec2;
 
 const MAX_POWER_OF_T: usize = 5;
 
@@ -33,13 +33,13 @@ type MicroArcsecond = f64;
 
 #[derive(Debug, Default)]
 struct NutationComponents {
-    planetary: XY,
-    luni_solar: XY,
+    planetary: DVec2,
+    luni_solar: DVec2,
 }
 
-/// (X, Y) coordinates of the Celestial Intermediate Pole (CIP) using the the IAU 2006 precession
-/// and IAU 2000A nutation models.
-pub fn xy(t: TDBJulianCenturiesSinceJ2000) -> XY {
+/// Calculates the (X, Y) coordinates of the Celestial Intermediate Pole (CIP) using the the IAU
+/// 2006 precession and IAU 2000A nutation models.
+pub fn xy(t: TDBJulianCenturiesSinceJ2000) -> DVec2 {
     let powers_of_t = powers_of_t(t);
     let fundamental_args = fundamental_args(t);
     let polynomial_components = polynomial_components(&powers_of_t);
@@ -78,8 +78,8 @@ fn fundamental_args(t: TDBJulianCenturiesSinceJ2000) -> FundamentalArgs {
     ]
 }
 
-fn polynomial_components(powers_of_t: &PowersOfT) -> XY {
-    let mut result = [0.0; 2];
+fn polynomial_components(powers_of_t: &PowersOfT) -> DVec2 {
+    let mut result = DVec2::default();
     for (i, power_of_t) in powers_of_t.iter().enumerate().rev() {
         result[0] += polynomial::COEFFICIENTS.x[i] * power_of_t;
         result[1] += polynomial::COEFFICIENTS.y[i] * power_of_t;
@@ -166,16 +166,17 @@ fn nutation_components(
 }
 
 fn calculate_cip_unit_vector(
-    polynomial_components: &XY,
+    polynomial_components: &DVec2,
     nutation_components: &NutationComponents,
-) -> XY {
+) -> DVec2 {
     let x_arcsec = polynomial_components[0]
         + (nutation_components.planetary[0] + nutation_components.luni_solar[0]) / 1e6;
     let y_arcsec = polynomial_components[1]
         + (nutation_components.planetary[1] + nutation_components.luni_solar[1]) / 1e6;
-    let x = arcsec_to_rad(x_arcsec);
-    let y = arcsec_to_rad(y_arcsec);
-    [x, y]
+    DVec2 {
+        x: arcsec_to_rad(x_arcsec),
+        y: arcsec_to_rad(y_arcsec),
+    }
 }
 
 #[cfg(test)]

crates/lox_core/src/earth/coordinate_transformations.rs

@@ -0,0 +1,134 @@
+/*
+ * 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/.
+ */
+
+//! Module coordinate_transformations provides functions for transforming coordinates between
+//! reference systems.
+
+use glam::{DMat3, DVec2};
+
+use crate::types::Radians;
+
+/// The spherical angles E and d.
+struct SphericalAngles {
+    e: Radians,
+    d: Radians,
+}
+
+impl SphericalAngles {
+    fn new(cip: DVec2) -> Self {
+        let r2 = cip[0] * cip[0] + cip[1] * cip[1];
+        let e = cip[1].atan2(cip[0]);
+        let d = (r2 / (1.0 - r2)).sqrt().atan();
+        Self { e, d }
+    }
+}
+
+/// Compute the celestial to intermediate-frame-of-date matrix given the CIP (X, Y) coordinates
+/// and the CIO locator, s. This matrix is the first step in transforming CRS to
+/// TRS coordinates.
+///
+/// Note that the signs of all angles are reversed relative to ERFA, which uses left-handed
+/// coordinates, whereas glam is right-handed.
+pub fn celestial_to_intermediate_frame_of_date_matrix(cip: DVec2, s: Radians) -> DMat3 {
+    let spherical_angles = SphericalAngles::new(cip);
+    let mut result = DMat3::default();
+    result = DMat3::from_rotation_z(-spherical_angles.e) * result;
+    result = DMat3::from_rotation_y(-spherical_angles.d) * result;
+    DMat3::from_rotation_z(spherical_angles.e + s) * result
+}
+
+#[cfg(test)]
+mod tests {
+    use float_eq::assert_float_eq;
+
+    use super::*;
+
+    const TOLERANCE: f64 = 1e-9;
+
+    #[test]
+    fn test_celestial_to_intermediate_frame_of_date_matrix_jd0() {
+        let cip = DVec2 {
+            x: -0.4088355637476968,
+            y: -0.38359667445777073,
+        };
+        let s = -0.0723985415686306;
+        let expected = [
+            0.899981235912944,
+            -0.151285348992267,
+            -0.408835563747697,
+            -0.019051024078611,
+            0.923304202214251,
+            -0.383596674457771,
+            0.435512150790498,
+            0.353018545339750,
+            0.8280743162060046,
+        ];
+        let actual = celestial_to_intermediate_frame_of_date_matrix(cip, s).to_cols_array();
+        assert_mat3_eq(&expected, &actual)
+    }
+
+    #[test]
+    fn test_celestial_to_intermediate_frame_of_date_matrix_j2000() {
+        let cip = DVec2 {
+            x: -0.0000269463795685740,
+            y: -0.00002800472282281282,
+        };
+        let s = -0.00000001013396519178;
+        let expected = [
+            0.999999999636946,
+            -0.00000001051127817488,
+            -0.000026946379569,
+            0.00000000975665225778,
+            0.999999999607868,
+            -0.000028004722823,
+            0.000026946379852,
+            0.000028004722550,
+            0.999999999244814,
+        ];
+        let actual = celestial_to_intermediate_frame_of_date_matrix(cip, s).to_cols_array();
+        assert_mat3_eq(&expected, &actual)
+    }
+
+    #[test]
+    fn test_celestial_to_intermediate_frame_of_date_matrix_j2100() {
+        let cip = DVec2 {
+            x: 0.00972070446172924,
+            y: -0.0000673058699616719,
+        };
+        let s = -0.00000000480511934533;
+        let expected = [
+            0.999952752836184,
+            0.00000032233307144280,
+            0.009720704461729,
+            0.00000033194308309287,
+            0.999999997734904,
+            -0.00006730586996167191,
+            -0.009720704461405,
+            0.00006730591667081695,
+            0.999952750571089,
+        ];
+        let actual = celestial_to_intermediate_frame_of_date_matrix(cip, s).to_cols_array();
+        assert_mat3_eq(&expected, &actual)
+    }
+
+    fn assert_mat3_eq(expected: &[f64; 9], actual: &[f64; 9]) {
+        for i in 0..9 {
+            assert_float_eq!(
+                expected[i],
+                actual[i],
+                abs <= TOLERANCE,
+                "actual matrix differed from expected matrix at col {}, row {}:\n\t\
+                    expected: {},\n\tactual: {}",
+                i / 3,
+                i % 3,
+                DMat3::from_cols_array(expected),
+                DMat3::from_cols_array(actual),
+            );
+        }
+    }
+}

crates/lox_core/src/earth/nutation.rs

@@ -1,3 +1,14 @@
+/*
+ * 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/.
+ */
+
+//! Module nutation exposes a function for calculating Earth nutation using a number of IAU nutation
+//! models.
+
 use std::ops::Add;
 
 use crate::earth::nutation::iau1980::nutation_iau1980;