diff --git a/crates/lox_core/src/bodies.rs b/crates/lox_core/src/bodies.rs index 59925935..2de7d7fd 100644 --- a/crates/lox_core/src/bodies.rs +++ b/crates/lox_core/src/bodies.rs @@ -14,6 +14,7 @@ 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; diff --git a/crates/lox_core/src/bodies/cio.rs b/crates/lox_core/src/bodies/cio.rs new file mode 100644 index 00000000..d233b26e --- /dev/null +++ b/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; diff --git a/crates/lox_core/src/bodies/cio/s06.rs b/crates/lox_core/src/bodies/cio/s06.rs new file mode 100644 index 00000000..c37a81ab --- /dev/null +++ b/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); + }); + } +} diff --git a/crates/lox_core/src/bodies/cio/s06/terms.rs b/crates/lox_core/src/bodies/cio/s06/terms.rs new file mode 100644 index 00000000..c1227f07 --- /dev/null +++ b/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 } +];