celestial_core/cio/
mod.rs

1//! Celestial Intermediate Origin (CIO) based celestial-to-terrestrial transformations.
2//!
3//! This module implements the IAU 2000/2006 CIO-based transformation from GCRS (Geocentric
4//! Celestial Reference System) to CIRS (Celestial Intermediate Reference System). The CIO
5//! approach is the modern replacement for the classical equinox-based method.
6//!
7//! # Components
8//!
9//! - [`CipCoordinates`]: X/Y coordinates of the Celestial Intermediate Pole
10//! - [`CioLocator`]: The CIO locator `s`, positioning the origin on the CIP equator
11//! - [`EquationOfOrigins`]: Relates CIO-based and equinox-based right ascension
12//! - [`CioSolution`]: Bundles all CIO quantities for a given epoch
13//!
14//! # Usage
15//!
16//! For most use cases, compute a [`CioSolution`] from the NPB (nutation-precession-bias) matrix:
17//!
18//! ```ignore
19//! let solution = CioSolution::calculate(&npb_matrix, tt_centuries)?;
20//! let cirs_matrix = gcrs_to_cirs_matrix(solution.cip.x, solution.cip.y, solution.s);
21//! ```
22
23pub mod coordinates;
24pub mod locator;
25pub mod origins;
26
27pub use coordinates::CipCoordinates;
28pub use locator::CioLocator;
29pub use origins::EquationOfOrigins;
30
31use crate::errors::AstroResult;
32use crate::matrix::RotationMatrix3;
33
34/// Builds the GCRS-to-CIRS rotation matrix from CIP coordinates and CIO locator.
35///
36/// This implements the IAU 2006 CIO-based transformation using three rotations:
37/// R₃(E) · R₂(d) · R₃(-(E+s)) where E = atan2(Y, X) and d = atan(sqrt(X²+Y²/(1-X²-Y²))).
38pub fn gcrs_to_cirs_matrix(x: f64, y: f64, s: f64) -> RotationMatrix3 {
39    let r2 = x * x + y * y;
40    let e = if r2 > 0.0 { libm::atan2(y, x) } else { 0.0 };
41    let d = libm::atan(libm::sqrt(r2 / (1.0 - r2)));
42
43    let mut matrix = RotationMatrix3::identity();
44    matrix.rotate_z(e);
45    matrix.rotate_y(d);
46    matrix.rotate_z(-(e + s));
47
48    matrix
49}
50
51/// All CIO-based quantities for a given epoch.
52///
53/// Bundles CIP coordinates, CIO locator, and equation of origins — everything needed
54/// for the GCRS↔CIRS transformation.
55#[derive(Debug, Clone, PartialEq)]
56pub struct CioSolution {
57    /// CIP X/Y coordinates (radians)
58    pub cip: CipCoordinates,
59    /// CIO locator s (radians)
60    pub s: f64,
61    /// Equation of origins (radians) — difference between CIO-based and equinox-based RA
62    pub equation_of_origins: f64,
63}
64
65impl CioSolution {
66    /// Computes all CIO quantities from an NPB matrix and TT centuries since J2000.
67    pub fn calculate(
68        npb_matrix: &crate::matrix::RotationMatrix3,
69        tt_centuries: f64,
70    ) -> AstroResult<Self> {
71        let cip = CipCoordinates::from_npb_matrix(npb_matrix)?;
72
73        let locator = CioLocator::iau2006a(tt_centuries);
74        let s = locator.calculate(cip.x, cip.y)?;
75
76        let equation_of_origins = EquationOfOrigins::from_npb_and_locator(npb_matrix, s)?;
77
78        Ok(Self {
79            cip,
80            s,
81            equation_of_origins,
82        })
83    }
84}
85
86#[cfg(test)]
87mod tests {
88    use super::*;
89
90    #[test]
91    fn cio_identity_matrix_returns_zero_components() {
92        let identity = crate::matrix::RotationMatrix3::identity();
93        let solution = CioSolution::calculate(&identity, 0.0).unwrap();
94
95        assert!(solution.cip.x.abs() < 1e-15);
96        assert!(solution.cip.y.abs() < 1e-15);
97        assert!(solution.s.abs() < 1e-8);
98        assert!(solution.equation_of_origins.abs() < 1e-8);
99    }
100
101    #[test]
102    fn gcrs_to_cirs_matrix_with_zero_inputs_returns_identity() {
103        let matrix = gcrs_to_cirs_matrix(0.0, 0.0, 0.0);
104        let identity = RotationMatrix3::identity();
105        assert!(matrix.max_difference(&identity) < 1e-15);
106    }
107
108    #[test]
109    fn gcrs_to_cirs_matrix_is_rotation_matrix() {
110        let matrix = gcrs_to_cirs_matrix(1e-6, 2e-6, 5e-9);
111        assert!(matrix.is_rotation_matrix(1e-14));
112    }
113
114    #[test]
115    fn gcrs_to_cirs_matrix_small_cip_produces_near_identity() {
116        let x = 1e-6;
117        let y = 1e-6;
118        let s = 1e-9;
119        let matrix = gcrs_to_cirs_matrix(x, y, s);
120        let identity = RotationMatrix3::identity();
121        assert!(matrix.max_difference(&identity) < 1e-5);
122    }
123}
celestial_core/cio/mod.rs

celestial_core/cio/
mod.rs
