From 7e2886a6846766a53fc04a7eca513d8349d105ad Mon Sep 17 00:00:00 2001 From: Dar Dahlen Date: Mon, 13 Jan 2025 15:29:36 -0800 Subject: [PATCH] Add Earth precession calculation (#169) * Add earth precession calculation * bump pprof version --- CHANGELOG.md | 1 + src/kete/conversion.py | 3 +- src/kete/rust/frame.rs | 39 ++++++++++++ src/kete/rust/lib.rs | 1 + src/kete_core/Cargo.toml | 4 +- src/kete_core/src/frames/definitions.rs | 82 +++++++++++++++++++++++++ 6 files changed, 127 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fb74b86..c0f4278 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Added support for SPICE kernels of type 9, this allows reading SOHO spice files. - Added support for SPICE kernels of type 18, this allows reading Rosetta spice files. +- Added Equatorial frame as observed calculation. ### Changed diff --git a/src/kete/conversion.py b/src/kete/conversion.py index 36ebd92..1f958d1 100644 --- a/src/kete/conversion.py +++ b/src/kete/conversion.py @@ -9,7 +9,7 @@ from numpy.typing import NDArray from . import _core from . import constants -from ._core import compute_obliquity +from ._core import compute_obliquity, earth_precession_rotation __all__ = [ "bin_data", @@ -25,6 +25,7 @@ "compute_tisserand", "dec_degrees_to_dms", "dec_dms_to_degrees", + "earth_precession_rotation", "flux_to_mag", "mag_to_flux", "ra_degrees_to_hms", diff --git a/src/kete/rust/frame.rs b/src/kete/rust/frame.rs index e38900d..8eac0cf 100644 --- a/src/kete/rust/frame.rs +++ b/src/kete/rust/frame.rs @@ -95,3 +95,42 @@ pub fn ecef_to_wgs_lat_lon(x: f64, y: f64, z: f64) -> (f64, f64, f64) { pub fn calc_obliquity_py(time: f64) -> f64 { calc_obliquity(time).to_degrees() } + +/// Calculate rotation matrix which transforms a vector from the J2000 Equatorial +/// frame to the desired epoch. +/// +/// Earth's north pole precesses at a rate of about 50 arcseconds per year. +/// This means there was an approximately 20 arcminute rotation of the Equatorial +/// axis from the year 2000 to 2025. +/// +/// This implementation is valid for around 200 years on either side of 2000 to +/// within sub micro-arcsecond accuracy. +/// +/// This function is an implementation equation (21) from this paper: +/// "Expressions for IAU 2000 precession quantities" +/// Capitaine, N. ; Wallace, P. T. ; Chapront, J. +/// Astronomy and Astrophysics, v.412, p.567-586 (2003) +/// +/// It is recommended to first look at the following paper, as it provides useful +/// discussion to help understand the above model. This defines the model used +/// by JPL Horizons: +/// "Precession matrix based on IAU (1976) system of astronomical constants." +/// Lieske, J. H. +/// Astronomy and Astrophysics, vol. 73, no. 3, Mar. 1979, p. 282-284. +/// +/// The IAU 2000 model paper improves accuracy by approximately ~300 mas/century over +/// the 1976 model. +/// +/// Parameters +/// ---------- +/// tdb_time: +/// Time in TDB scaled Julian Days. +#[pyfunction] +#[pyo3(name = "earth_precession_rotation")] +pub fn calc_earth_precession(time: f64) -> Vec> { + earth_precession_rotation(time) + .matrix() + .row_iter() + .map(|x| x.iter().cloned().collect()) + .collect() +} diff --git a/src/kete/rust/lib.rs b/src/kete/rust/lib.rs index 85fa97a..18a88eb 100644 --- a/src/kete/rust/lib.rs +++ b/src/kete/rust/lib.rs @@ -85,6 +85,7 @@ fn _core(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(frame::wgs_lat_lon_to_ecef, m)?)?; m.add_function(wrap_pyfunction!(frame::ecef_to_wgs_lat_lon, m)?)?; m.add_function(wrap_pyfunction!(frame::calc_obliquity_py, m)?)?; + m.add_function(wrap_pyfunction!(frame::calc_earth_precession, m)?)?; m.add_function(wrap_pyfunction!(kepler::compute_eccentric_anomaly_py, m)?)?; m.add_function(wrap_pyfunction!(kepler::propagation_kepler_py, m)?)?; diff --git a/src/kete_core/Cargo.toml b/src/kete_core/Cargo.toml index 1c57c53..b6350cc 100644 --- a/src/kete_core/Cargo.toml +++ b/src/kete_core/Cargo.toml @@ -35,11 +35,11 @@ criterion = { version = "*", features = ["html_reports"] } # pprof is used for flame graphs, this is failing on windows currently so only # linux and mac are supported here. [target.'cfg(target_os = "linux")'.dev-dependencies] -pprof = { version = "0.13", features = ["flamegraph", "criterion"] } +pprof = { version = "0.14", features = ["flamegraph", "criterion"] } # macos needs the frame-pointer flag, however this doesn't function on linux [target.'cfg(target_os = "macos")'.dev-dependencies] -pprof = { version = "0.13", features = ["flamegraph", "criterion", "frame-pointer"] } +pprof = { version = "0.14", features = ["flamegraph", "criterion", "frame-pointer"] } [[bench]] name = "propagation" diff --git a/src/kete_core/src/frames/definitions.rs b/src/kete_core/src/frames/definitions.rs index 7bbc7b4..6f25109 100644 --- a/src/kete_core/src/frames/definitions.rs +++ b/src/kete_core/src/frames/definitions.rs @@ -362,12 +362,94 @@ pub fn rotate_around( rot.transform_vector(vector) } +/// Rotation which transforms a vector from the J2000 Equatorial frame to the +/// desired epoch. +/// +/// Earth's north pole precesses at a rate of about 50 arcseconds per year. +/// This means there was an approximately 20 arcminute rotation of the Equatorial +/// axis from the year 2000 to 2025. +/// +/// This implementation is valid for around 200 years on either side of 2000 to +/// within sub micro-arcsecond accuracy. +/// +/// This function is an implementation equation (21) from this paper: +/// "Expressions for IAU 2000 precession quantities" +/// Capitaine, N. ; Wallace, P. T. ; Chapront, J. +/// Astronomy and Astrophysics, v.412, p.567-586 (2003) +/// +/// It is recommended to first look at the following paper, as it provides useful +/// discussion to help understand the above model. This defines the model used +/// by JPL Horizons: +/// "Precession matrix based on IAU (1976) system of astronomical constants." +/// Lieske, J. H. +/// Astronomy and Astrophysics, vol. 73, no. 3, Mar. 1979, p. 282-284. +/// +/// The IAU 2000 model paper improves accuracy by approximately ~300 mas/century over +/// the 1976 model. +/// +/// # Arguments +/// +/// * `tdb_time` - Time in TDB scaled Julian Days. +/// +#[inline(always)] +pub fn earth_precession_rotation(tdb_time: f64) -> Rotation3 { + // centuries since 2000 + let t = (tdb_time - 2451545.0) / 36525.0; + + // angles as defined in the cited paper, equations (21) + // Note that equation 45 is an even more developed model, which takes into + // account frame bias in addition to simple precession, however more clarity + // on the DE source and interpretation is probably required to take advantage + // of this increased precision. + let angle_c = -((2.5976176 + + (2306.0809506 + (0.3019015 + (0.0179663 + (-0.0000327 - 0.0000002 * t) * t) * t) * t) + * t) + / 3600.0) + .to_radians(); + let angle_a = -((-2.5976176 + + (2306.0803226 + (1.094779 + (0.0182273 + (0.000047 - 0.0000003 * t) * t) * t) * t) * t) + / 3600.0) + .to_radians(); + let angle_b = ((2004.1917476 + + (-0.4269353 + (-0.0418251 + (-0.0000601 - 0.0000001 * t) * t) * t) * t) + * t + / 3600.0) + .to_radians(); + let z_axis = Vector3::z_axis(); + Rotation3::from_axis_angle(&z_axis, angle_a) + * Rotation3::from_axis_angle(&Vector3::y_axis(), angle_b) + * Rotation3::from_axis_angle(&z_axis, angle_c) +} + #[cfg(test)] mod tests { use std::f64::consts::{FRAC_PI_2, TAU}; + use nalgebra::Matrix3; + use crate::frames::*; + #[test] + fn test_earth_precession() { + let rot = earth_precession_rotation(2451545.0); + assert!((rot.matrix() - Matrix3::identity()).norm() < 1e-16); + + let rot = earth_precession_rotation(2433282.42345905); + + let expected = Matrix3::new( + 0.9999257168056067, + -0.011178271729385899, + -0.004858715050078201, + 0.011178271443436222, + 0.9999375208015913, + -2.72158794420131e-05, + 0.004858715707952332, + -2.70981779365330e-05, + 0.9999881960040119 + ); + assert!((rot.matrix() - expected).norm() < 1e-16); + } + #[test] fn test_ecliptic_rot_roundtrip() { let vec = ecliptic_to_equatorial(&[1.0, 2.0, 3.0].into());