celestial_time/transforms/
rotation.rs

1use crate::{JulianDate, TimeResult};
2use celestial_core::angle::wrap_0_2pi;
3use celestial_core::constants::{J2000_JD, TWOPI};
4use celestial_core::math::fmod;
5
6pub fn earth_rotation_angle(ut1_jd: &JulianDate) -> TimeResult<f64> {
7    let ut1_jd1 = ut1_jd.jd1();
8    let ut1_jd2 = ut1_jd.jd2();
9
10    let (d1, d2) = if ut1_jd1 < ut1_jd2 {
11        (ut1_jd1, ut1_jd2)
12    } else {
13        (ut1_jd2, ut1_jd1)
14    };
15
16    let t = d1 + (d2 - J2000_JD);
17
18    let f = fmod(d1, 1.0) + fmod(d2, 1.0);
19
20    let theta = wrap_0_2pi(TWOPI * (f + 0.7790572732640 + 0.00273781191135448 * t));
21
22    Ok(theta)
23}
24
25#[cfg(test)]
26mod tests {
27    use super::*;
28    use crate::UT1;
29    use celestial_core::constants::J2000_JD;
30
31    #[test]
32    fn test_era_j2000() {
33        let ut1 = UT1::j2000();
34        let era = earth_rotation_angle(&ut1.to_julian_date()).unwrap();
35
36        assert!((era - 4.894961212823757).abs() < 1e-12);
37    }
38
39    #[test]
40    fn test_era_precision() {
41        let test_cases = [
42            (J2000_JD, 0.0),
43            (2451545.5, 0.0),
44            (2440587.5, 0.0),
45            (J2000_JD, 0.5),
46        ];
47
48        for (jd1, jd2) in test_cases {
49            let ut1 = UT1::from_julian_date(JulianDate::new(jd1, jd2));
50            let era = earth_rotation_angle(&ut1.to_julian_date()).unwrap();
51
52            assert!((0.0..2.0 * celestial_core::constants::PI).contains(&era));
53        }
54    }
55}
celestial_time/transforms/rotation.rs

celestial_time/transforms/
rotation.rs
