Implement 06 CIO locator calculation

crates/lox_core/src/bodies.rs

@@ -14,6 +14,8 @@ use crate::time::constants::f64::{SECONDS_PER_DAY, SECONDS_PER_JULIAN_CENTURY};
 mod generated;
 pub use generated::*;
 
+mod cio;
+mod cip;
 pub mod fundamental;
 pub mod nutation;
 

crates/lox_core/src/bodies/cio.rs

@@ -0,0 +1,12 @@
+/*
+ * 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/.
+ */
+
+// TODO: Remove this once all module components are actively used.
+#![allow(dead_code)]
+
+mod s06;

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

@@ -0,0 +1,123 @@
+/*
+ * 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/.
+ */
+
+mod terms;
+
+use crate::bodies::cip::xy06::XY;
+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::math::arcsec_to_rad;
+use crate::time::intervals::TDBJulianCenturiesSinceJ2000;
+use crate::types::Radians;
+
+/// 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 {
+    let fundamental_args = fundamental_args(t);
+    let evaluated_terms = evaluate_terms(&fundamental_args);
+    let arcsec = fast_polynomial::poly_array(t, &evaluated_terms);
+    let radians = arcsec_to_rad(arcsec);
+    radians - xy[0] * xy[1] / 2.0
+}
+
+fn fundamental_args(t: TDBJulianCenturiesSinceJ2000) -> FundamentalArgs {
+    // The output of the CIO calculation is dependent on the ordering of these arguments. DO NOT
+    // EDIT.
+    [
+        Moon.mean_anomaly_iers03(t),
+        Sun.mean_anomaly_iers03(t),
+        Moon.mean_longitude_minus_ascending_node_mean_longitude_iers03(t),
+        mean_moon_sun_elongation_iers03(t),
+        Moon.ascending_node_mean_longitude_iers03(t),
+        Venus.mean_longitude_iers03(t),
+        Earth.mean_longitude_iers03(t),
+        general_accum_precession_in_longitude_iers03(t),
+    ]
+}
+
+fn evaluate_terms(args: &FundamentalArgs) -> [f64; 6] {
+    [
+        evaluate_single_order_terms(args, terms::COEFFICIENTS[0], &terms::ZERO_ORDER),
+        evaluate_single_order_terms(args, terms::COEFFICIENTS[1], &terms::FIRST_ORDER),
+        evaluate_single_order_terms(args, terms::COEFFICIENTS[2], &terms::SECOND_ORDER),
+        evaluate_single_order_terms(args, terms::COEFFICIENTS[3], &terms::THIRD_ORDER),
+        evaluate_single_order_terms(args, terms::COEFFICIENTS[4], &terms::FOURTH_ORDER),
+        terms::COEFFICIENTS[5],
+    ]
+}
+
+fn evaluate_single_order_terms(
+    args: &FundamentalArgs,
+    term_coefficient: f64,
+    terms: &[terms::Term],
+) -> f64 {
+    terms.iter().rev().fold(term_coefficient, |acc, term| {
+        let a = term
+            .fundamental_arg_coeffs
+            .iter()
+            .zip(args)
+            .fold(0.0, |acc, (coeff, arg)| acc + coeff * arg);
+
+        acc + term.sin_coeff * a.sin() + term.cos_coeff * a.cos()
+    })
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::bodies::cip::xy06::xy;
+    use float_eq::assert_float_eq;
+
+    const TOLERANCE: f64 = 1e-11;
+
+    #[test]
+    fn test_s_jd0() {
+        let jd0: TDBJulianCenturiesSinceJ2000 = -67.11964407939767;
+        let xy = xy(jd0);
+        assert_float_eq!(s(jd0, xy), -0.0723985415686306, rel <= TOLERANCE);
+    }
+
+    #[test]
+    fn test_s_j2000() {
+        let j2000: TDBJulianCenturiesSinceJ2000 = 0.0;
+        let xy = xy(j2000);
+        assert_float_eq!(s(j2000, xy), -0.00000001013396519178, rel <= TOLERANCE);
+    }
+
+    #[test]
+    fn test_s_j2100() {
+        let j2100: TDBJulianCenturiesSinceJ2000 = 1.0;
+        let xy = xy(j2100);
+        assert_float_eq!(s(j2100, xy), -0.00000000480511934533, rel <= TOLERANCE);
+    }
+
+    #[test]
+    fn test_fundamental_args_ordering() {
+        let j2000: TDBJulianCenturiesSinceJ2000 = 0.0;
+        let actual = fundamental_args(j2000);
+        let expected = [
+            Moon.mean_anomaly_iers03(j2000),
+            Sun.mean_anomaly_iers03(j2000),
+            Moon.mean_longitude_minus_ascending_node_mean_longitude_iers03(j2000),
+            mean_moon_sun_elongation_iers03(j2000),
+            Moon.ascending_node_mean_longitude_iers03(j2000),
+            Venus.mean_longitude_iers03(j2000),
+            Earth.mean_longitude_iers03(j2000),
+            general_accum_precession_in_longitude_iers03(j2000),
+        ];
+
+        expected.iter().enumerate().for_each(|(i, expected)| {
+            assert_float_eq!(*expected, actual[i], rel <= TOLERANCE);
+        });
+    }
+}

crates/lox_core/src/bodies/cio/s06/terms.rs

@@ -0,0 +1,118 @@
+/*
+ * 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::types::Radians;
+
+pub(super) const COEFFICIENTS: [f64; 6] = [
+    94.00e-6,
+    3808.65e-6,
+    -122.68e-6,
+    -72574.11e-6,
+    27.98e-6,
+    15.62e-6,
+];
+
+/// Coefficients of l, l', F, D, Ω, LVe, LE and pA.
+pub(super) type FundamentalArgCoeffs = [Radians; 8];
+
+pub(super) struct Term {
+    pub fundamental_arg_coeffs: FundamentalArgCoeffs,
+    pub sin_coeff: f64,
+    pub cos_coeff: f64,
+}
+
+#[rustfmt::skip]
+// @formatter:off (sometimes RustRover ignores the rustfmt skip)
+pub(super) const ZERO_ORDER: [Term; 33] = [
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  0.0], sin_coeff: -2640.73e-6, cos_coeff:  0.39e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  0.0,  0.0,  2.0,  0.0,  0.0,  0.0], sin_coeff:   -63.53e-6, cos_coeff:  0.02e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0, -2.0,  3.0,  0.0,  0.0,  0.0], sin_coeff:   -11.75e-6, cos_coeff: -0.01e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0, -2.0,  1.0,  0.0,  0.0,  0.0], sin_coeff:   -11.21e-6, cos_coeff: -0.01e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0, -2.0,  2.0,  0.0,  0.0,  0.0], sin_coeff:     4.57e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0,  0.0,  3.0,  0.0,  0.0,  0.0], sin_coeff:    -2.02e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0,  0.0,  1.0,  0.0,  0.0,  0.0], sin_coeff:    -1.98e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  0.0,  0.0,  3.0,  0.0,  0.0,  0.0], sin_coeff:     1.72e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  1.0,  0.0,  0.0,  1.0,  0.0,  0.0,  0.0], sin_coeff:     1.41e-6, cos_coeff:  0.01e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  1.0,  0.0,  0.0, -1.0,  0.0,  0.0,  0.0], sin_coeff:     1.26e-6, cos_coeff:  0.01e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0,  0.0,  0.0, -1.0,  0.0,  0.0,  0.0], sin_coeff:     0.63e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  0.0], sin_coeff:     0.63e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  1.0,  2.0, -2.0,  3.0,  0.0,  0.0,  0.0], sin_coeff:    -0.46e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  1.0,  2.0, -2.0,  1.0,  0.0,  0.0,  0.0], sin_coeff:    -0.45e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  4.0, -4.0,  4.0,  0.0,  0.0,  0.0], sin_coeff:    -0.36e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  1.0, -1.0,  1.0, -8.0, 12.0,  0.0], sin_coeff:     0.24e-6, cos_coeff:  0.12e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0,  0.0,  0.0,  0.0,  0.0,  0.0], sin_coeff:    -0.32e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0,  0.0,  2.0,  0.0,  0.0,  0.0], sin_coeff:    -0.28e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0,  2.0,  0.0,  3.0,  0.0,  0.0,  0.0], sin_coeff:    -0.27e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0,  2.0,  0.0,  1.0,  0.0,  0.0,  0.0], sin_coeff:    -0.26e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0, -2.0,  0.0,  0.0,  0.0,  0.0], sin_coeff:     0.21e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  1.0, -2.0,  2.0, -3.0,  0.0,  0.0,  0.0], sin_coeff:    -0.19e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  1.0, -2.0,  2.0, -1.0,  0.0,  0.0,  0.0], sin_coeff:    -0.18e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  0.0,  0.0,  0.0,  8.0,-13.0, -1.0], sin_coeff:     0.10e-6, cos_coeff: -0.05e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  0.0,  2.0,  0.0,  0.0,  0.0,  0.0], sin_coeff:    -0.15e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [2.0,  0.0, -2.0,  0.0, -1.0,  0.0,  0.0,  0.0], sin_coeff:     0.14e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  1.0,  2.0, -2.0,  2.0,  0.0,  0.0,  0.0], sin_coeff:     0.14e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0,  0.0, -2.0,  1.0,  0.0,  0.0,  0.0], sin_coeff:    -0.14e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0,  0.0, -2.0, -1.0,  0.0,  0.0,  0.0], sin_coeff:    -0.14e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  4.0, -2.0,  4.0,  0.0,  0.0,  0.0], sin_coeff:    -0.13e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0, -2.0,  4.0,  0.0,  0.0,  0.0], sin_coeff:     0.11e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0, -2.0,  0.0, -3.0,  0.0,  0.0,  0.0], sin_coeff:    -0.11e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0, -2.0,  0.0, -1.0,  0.0,  0.0,  0.0], sin_coeff:    -0.11e-6, cos_coeff:  0.00e-6 },
+];
+
+#[rustfmt::skip]
+// @formatter:off (sometimes RustRover ignores the rustfmt skip)
+pub(super) const FIRST_ORDER: [Term; 3] = [
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  0.0,  0.0,  2.0,  0.0,  0.0,  0.0], sin_coeff: -0.07e-6, cos_coeff:  3.57e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  0.0], sin_coeff:  1.73e-6, cos_coeff: -0.03e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0, -2.0,  3.0,  0.0,  0.0,  0.0], sin_coeff:  0.00e-6, cos_coeff:  0.48e-6 },
+];
+
+#[rustfmt::skip]
+// @formatter:off (sometimes RustRover ignores the rustfmt skip)
+pub(super) const SECOND_ORDER: [Term; 25] = [
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  0.0], sin_coeff: 743.52e-6, cos_coeff: -0.17e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0, -2.0,  2.0,  0.0,  0.0,  0.0], sin_coeff:  56.91e-6, cos_coeff:  0.06e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0,  0.0,  2.0,  0.0,  0.0,  0.0], sin_coeff:   9.84e-6, cos_coeff: -0.01e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  0.0,  0.0,  2.0,  0.0,  0.0,  0.0], sin_coeff:  -8.85e-6, cos_coeff:  0.01e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  1.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0], sin_coeff:  -6.38e-6, cos_coeff: -0.05e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0], sin_coeff:  -3.07e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  1.0,  2.0, -2.0,  2.0,  0.0,  0.0,  0.0], sin_coeff:   2.23e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0,  0.0,  1.0,  0.0,  0.0,  0.0], sin_coeff:   1.67e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0,  2.0,  0.0,  2.0,  0.0,  0.0,  0.0], sin_coeff:   1.30e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  1.0, -2.0,  2.0, -2.0,  0.0,  0.0,  0.0], sin_coeff:   0.93e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0,  0.0, -2.0,  0.0,  0.0,  0.0,  0.0], sin_coeff:   0.68e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0, -2.0,  1.0,  0.0,  0.0,  0.0], sin_coeff:  -0.55e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0, -2.0,  0.0, -2.0,  0.0,  0.0,  0.0], sin_coeff:   0.53e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  0.0,  2.0,  0.0,  0.0,  0.0,  0.0], sin_coeff:  -0.27e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  0.0], sin_coeff:  -0.27e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0, -2.0, -2.0, -2.0,  0.0,  0.0,  0.0], sin_coeff:  -0.26e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0,  0.0,  0.0, -1.0,  0.0,  0.0,  0.0], sin_coeff:  -0.25e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0,  2.0,  0.0,  1.0,  0.0,  0.0,  0.0], sin_coeff:   0.22e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [2.0,  0.0,  0.0, -2.0,  0.0,  0.0,  0.0,  0.0], sin_coeff:  -0.21e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [2.0,  0.0, -2.0,  0.0, -1.0,  0.0,  0.0,  0.0], sin_coeff:   0.20e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0,  2.0,  2.0,  0.0,  0.0,  0.0], sin_coeff:   0.17e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [2.0,  0.0,  2.0,  0.0,  2.0,  0.0,  0.0,  0.0], sin_coeff:   0.13e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [2.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0], sin_coeff:  -0.13e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [1.0,  0.0,  2.0, -2.0,  2.0,  0.0,  0.0,  0.0], sin_coeff:  -0.12e-6, cos_coeff:  0.00e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0,  0.0,  0.0,  0.0,  0.0,  0.0], sin_coeff:  -0.11e-6, cos_coeff:  0.00e-6 },
+];
+
+#[rustfmt::skip]
+// @formatter:off (sometimes RustRover ignores the rustfmt skip)
+pub(super) const THIRD_ORDER: [Term; 4] = [
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  0.0], sin_coeff:  0.30e-6, cos_coeff: -23.42e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0, -2.0,  2.0,  0.0,  0.0,  0.0], sin_coeff: -0.03e-6, cos_coeff:  -1.46e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  2.0,  0.0,  2.0,  0.0,  0.0,  0.0], sin_coeff: -0.01e-6, cos_coeff:  -0.25e-6 },
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  0.0,  0.0,  2.0,  0.0,  0.0,  0.0], sin_coeff:  0.00e-6, cos_coeff:   0.23e-6 },
+];
+
+#[rustfmt::skip]
+// @formatter:off (sometimes RustRover ignores the rustfmt skip)
+pub(super) const FOURTH_ORDER: [Term; 1] = [
+    Term{ fundamental_arg_coeffs: [0.0,  0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  0.0], sin_coeff: -0.26e-6, cos_coeff: -0.01e-6 }
+];

crates/lox_core/src/bodies/cip.rs

@@ -0,0 +1,12 @@
+/*
+ * 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/.
+ */
+
+// TODO: Remove this once all module components are actively used.
+#![allow(dead_code)]
+
+pub mod xy06;