diff --git a/crates/lox_core/src/earth.rs b/crates/lox_core/src/earth.rs index d6508771..a429d4a0 100644 --- a/crates/lox_core/src/earth.rs +++ b/crates/lox_core/src/earth.rs @@ -11,4 +11,5 @@ mod cio; mod cip; +mod coordinate_transformations; mod nutation; diff --git a/crates/lox_core/src/earth/cio.rs b/crates/lox_core/src/earth/cio.rs index d6834acd..c34a26b0 100644 --- a/crates/lox_core/src/earth/cio.rs +++ b/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; diff --git a/crates/lox_core/src/earth/cio/s06.rs b/crates/lox_core/src/earth/cio/s06.rs index b21bcae8..3a196311 100644 --- a/crates/lox_core/src/earth/cio/s06.rs +++ b/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); diff --git a/crates/lox_core/src/earth/cip.rs b/crates/lox_core/src/earth/cip.rs index 0373b517..9e46f0ef 100644 --- a/crates/lox_core/src/earth/cip.rs +++ b/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; diff --git a/crates/lox_core/src/earth/cip/xy06.rs b/crates/lox_core/src/earth/cip/xy06.rs index 6f60bf58..f77066a3 100644 --- a/crates/lox_core/src/earth/cip/xy06.rs +++ b/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)] diff --git a/crates/lox_core/src/frames/transformations.rs b/crates/lox_core/src/earth/coordinate_transformations.rs similarity index 83% rename from crates/lox_core/src/frames/transformations.rs rename to crates/lox_core/src/earth/coordinate_transformations.rs index 4eb16ecb..17cfafd7 100644 --- a/crates/lox_core/src/frames/transformations.rs +++ b/crates/lox_core/src/earth/coordinate_transformations.rs @@ -6,14 +6,12 @@ * file, you can obtain one at https://mozilla.org/MPL/2.0/. */ -use crate::types::Radians; -use glam::DMat3; +//! Module coordinate_transformations provides functions for transforming coordinates between +//! reference systems. + +use glam::{DMat3, DVec2}; -// TODO: Decide on correct home for this type once module structure discussed with -// Helge. (An array is preferable to a tuple because the CIP calculation requires the dynamic -// access an array provides, and we want to avoid excess allocations moving between steps of the -// pipeline.) -type XY = [f64; 2]; +use crate::types::Radians; /// The spherical angles E and d. struct SphericalAngles { @@ -22,7 +20,7 @@ struct SphericalAngles { } impl SphericalAngles { - fn new(cip: XY) -> Self { + 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(); @@ -36,8 +34,7 @@ impl SphericalAngles { /// /// Note that the signs of all angles are reversed relative to ERFA, which uses left-handed /// coordinates, whereas glam is right-handed. -#[allow(dead_code)] // TODO: Remove this once all module components are actively used. -pub fn celestial_to_intermediate_frame_of_date_matrix(cip: XY, s: Radians) -> DMat3 { +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; @@ -47,15 +44,19 @@ pub fn celestial_to_intermediate_frame_of_date_matrix(cip: XY, s: Radians) -> DM #[cfg(test)] mod tests { - use super::*; use float_eq::assert_float_eq; + use super::*; + // TODO: Is this sufficient? const TOLERANCE: f64 = 1e-9; #[test] fn test_celestial_to_intermediate_frame_of_date_matrix_jd0() { - let cip = [-0.4088355637476968, -0.38359667445777073]; + let cip = DVec2 { + x: -0.4088355637476968, + y: -0.38359667445777073, + }; let s = -0.0723985415686306; let expected = [ 0.899981235912944, @@ -74,7 +75,10 @@ mod tests { #[test] fn test_celestial_to_intermediate_frame_of_date_matrix_j2000() { - let cip = [-0.0000269463795685740, -0.00002800472282281282]; + let cip = DVec2 { + x: -0.0000269463795685740, + y: -0.00002800472282281282, + }; let s = -0.00000001013396519178; let expected = [ 0.999999999636946, @@ -93,7 +97,10 @@ mod tests { #[test] fn test_celestial_to_intermediate_frame_of_date_matrix_j2100() { - let cip = [0.00972070446172924, -0.0000673058699616719]; + let cip = DVec2 { + x: 0.00972070446172924, + y: -0.0000673058699616719, + }; let s = -0.00000000480511934533; let expected = [ 0.999952752836184, diff --git a/crates/lox_core/src/earth/nutation.rs b/crates/lox_core/src/earth/nutation.rs index 8981ac49..ae03eeeb 100644 --- a/crates/lox_core/src/earth/nutation.rs +++ b/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; diff --git a/crates/lox_core/src/frames.rs b/crates/lox_core/src/frames.rs index 4c5ee68d..3dc8cd24 100644 --- a/crates/lox_core/src/frames.rs +++ b/crates/lox_core/src/frames.rs @@ -11,7 +11,6 @@ use std::fmt::{Debug, Display, Formatter}; use glam::{DMat3, DVec3}; pub mod iau; -mod transformations; // TODO: Replace with proper `Epoch` type type Epoch = f64;