Skip to content

Commit

Permalink
Move coordinate transformations to earth module
Browse files Browse the repository at this point in the history
  • Loading branch information
AngusGMorrison committed Dec 23, 2023
1 parent ac0a82b commit 37f7f67
Show file tree
Hide file tree
Showing 8 changed files with 62 additions and 35 deletions.
1 change: 1 addition & 0 deletions crates/lox_core/src/earth.rs
Original file line number Diff line number Diff line change
Expand Up @@ -11,4 +11,5 @@

mod cio;
mod cip;
mod coordinate_transformations;
mod nutation;
2 changes: 2 additions & 0 deletions crates/lox_core/src/earth/cio.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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;
11 changes: 7 additions & 4 deletions crates/lox_core/src/earth/cio/s06.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
3 changes: 3 additions & 0 deletions crates/lox_core/src/earth/cip.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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;
33 changes: 17 additions & 16 deletions crates/lox_core/src/earth/cip/xy06.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;

Expand All @@ -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);
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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)]
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand All @@ -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();
Expand All @@ -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;
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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,
Expand Down
11 changes: 11 additions & 0 deletions crates/lox_core/src/earth/nutation.rs
Original file line number Diff line number Diff line change
@@ -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;
Expand Down
1 change: 0 additions & 1 deletion crates/lox_core/src/frames.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down

0 comments on commit 37f7f67

Please sign in to comment.