celestial_core/obliquity.rs
1//! Mean obliquity of the ecliptic.
2//!
3//! The obliquity is the angle between Earth's equatorial plane and the ecliptic
4//! (the plane of Earth's orbit around the Sun). It's approximately 23.4° and
5//! decreases slowly due to gravitational perturbations from other planets.
6//!
7//! This module provides two IAU models:
8//!
9//! | Function | Model | J2000.0 Value | Polynomial Order |
10//! |----------|-------|---------------|------------------|
11//! | [`iau_2006_mean_obliquity`] | IAU 2006 | 84381.406″ | 5th order |
12//! | [`iau_1980_mean_obliquity`] | IAU 1980 | 84381.448″ | 3rd order |
13//!
14//! Both return the *mean* obliquity — the smoothly varying component without
15//! short-period nutation oscillations. For the *true* obliquity (mean + nutation
16//! in obliquity), add [`NutationResult::delta_eps`](crate::nutation::NutationResult::delta_eps).
17//!
18//! # Time Argument
19//!
20//! Both functions accept a two-part Julian Date in TDB. Split as `(jd1, jd2)`
21//! where typically `jd1 = 2451545.0` (J2000.0) and `jd2` is days from that epoch.
22//!
23//! # Example
24//!
25//! ```
26//! use celestial_core::obliquity::iau_2006_mean_obliquity;
27//! use celestial_core::constants::J2000_JD;
28//!
29//! // At J2000.0
30//! let eps = iau_2006_mean_obliquity(J2000_JD, 0.0);
31//! let eps_deg = eps.to_degrees();
32//! assert!((eps_deg - 23.4392794).abs() < 1e-6);
33//! ```
34
35use crate::constants::J2000_JD;
36
37/// Mean obliquity of the ecliptic using the IAU 2006 precession model.
38///
39/// Returns the mean obliquity in radians. This is a 5th-order polynomial
40/// valid for several centuries around J2000.0.
41///
42/// At J2000.0: ε₀ = 84381.406″ ≈ 23°26′21.406″
43pub fn iau_2006_mean_obliquity(date1: f64, date2: f64) -> f64 {
44 let t = ((date1 - J2000_JD) + date2) / crate::constants::DAYS_PER_JULIAN_CENTURY;
45
46 let obliquity_arcsec = 84381.406
47 + (-46.836769
48 + (-0.0001831 + (0.00200340 + (-0.000000576 + (-0.0000000434) * t) * t) * t) * t)
49 * t;
50
51 obliquity_arcsec * (crate::constants::PI / (180.0 * 3600.0))
52}
53
54/// Mean obliquity of the ecliptic using the IAU 1980 model.
55///
56/// Returns the mean obliquity in radians. This is a 3rd-order polynomial,
57/// less accurate than the IAU 2006 model but still used with IAU 1980 nutation.
58///
59/// At J2000.0: ε₀ = 84381.448″ ≈ 23°26′21.448″
60pub fn iau_1980_mean_obliquity(date1: f64, date2: f64) -> f64 {
61 let t = ((date1 - J2000_JD) + date2) / crate::constants::DAYS_PER_JULIAN_CENTURY;
62
63 let obliquity_arcsec = 84381.448 + (-46.8150 + (-0.00059 + (0.001813) * t) * t) * t;
64
65 obliquity_arcsec * (crate::constants::PI / (180.0 * 3600.0))
66}