astrodyn_frames 0.1.1

//! J2000 Precession matrix computation.
//!
//! Faithful port of JEOD's `precession_j2000.cc` (IAU-76/FK5).
//!
//! Reference: Mulcahy & Bond, "The RNP Routine for the Standard Epoch J2000",
//! NASA JSC-24574, September 1990.

// Constants ported verbatim from JEOD — suppress excessive precision warnings.
#![allow(clippy::excessive_precision)]

use glam::{DMat3, DVec3};
use std::f64::consts::PI;

const DEG_TO_RAD: f64 = PI / 180.0;
const ARCSEC_TO_RAD: f64 = DEG_TO_RAD / 3600.0;

/// Compute the precession matrix.
///
/// Port of JEOD `precession_j2000.cc`.
///
/// # Arguments
/// * `time` — Julian centuries since J2000.0 TT
///
/// # Returns
/// 3x3 precession rotation matrix (stored as JEOD convention: row-major in the
/// `rotation` array, representing the transformation from J2000 mean equator
/// to mean equator of date).
pub fn precession_matrix(time: f64) -> DMat3 {
    let time2 = time * time;
    let time3 = time2 * time;

    // Precession parameters in arcseconds (Mulcahy & Bond, JSC-24574),
    // then converted to radians in place (matching JEOD convention).
    let mut zeta = 2306.2181 * time + 0.30188 * time2 + 0.017998 * time3;
    let mut theta = 2004.3109 * time - 0.42665 * time2 - 0.041833 * time3;
    let mut z = 2306.2181 * time + 1.09468 * time2 + 0.018203 * time3;

    zeta *= ARCSEC_TO_RAD;
    theta *= ARCSEC_TO_RAD;
    z *= ARCSEC_TO_RAD;

    let (s_zeta, c_zeta) = zeta.sin_cos();
    let (s_theta, c_theta) = theta.sin_cos();
    let (s_z, c_z) = z.sin_cos();

    // Precession matrix: rot_z(zeta) * rot_y(-theta) * rot_z(z)
    // Stored in JEOD convention (see precession_j2000.cc)
    DMat3::from_cols(
        DVec3::new(
            c_theta * c_z * c_zeta - s_z * s_zeta,
            -s_zeta * c_theta * c_z - c_zeta * s_z,
            -s_theta * c_z,
        ),
        DVec3::new(
            s_z * c_theta * c_zeta + s_zeta * c_z,
            -s_z * s_zeta * c_theta + c_z * c_zeta,
            -s_theta * s_z,
        ),
        DVec3::new(c_zeta * s_theta, -s_zeta * s_theta, c_theta),
    )
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn precession_at_j2000_is_identity() {
        let p = precession_matrix(0.0);
        for i in 0..3 {
            for j in 0..3 {
                let expected = if i == j { 1.0 } else { 0.0 };
                assert!(
                    (p.col(j)[i] - expected).abs() < 1e-15,
                    "precession[{}][{}] = {}, expected {}",
                    i,
                    j,
                    p.col(j)[i],
                    expected
                );
            }
        }
    }
}