celestial_time/sidereal/
gmst.rs

1use super::angle::SiderealAngle;
2use crate::scales::{TT, UT1};
3use crate::JulianDate;
4use crate::TimeResult;
5use celestial_core::angle::wrap_0_2pi;
6use celestial_core::constants::{J2000_JD, TWOPI};
7use celestial_core::math::fmod;
8
9#[cfg(feature = "serde")]
10use serde::{Deserialize, Serialize};
11
12#[derive(Debug, Clone, Copy, PartialEq)]
13#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
14pub struct GMST(SiderealAngle);
15
16impl GMST {
17    pub fn from_ut1_and_tt(ut1: &UT1, tt: &TT) -> TimeResult<Self> {
18        let gmst_rad = calculate_gmst_iau2006(ut1, tt)?;
19        let angle = SiderealAngle::from_radians_exact(gmst_rad);
20        Ok(Self(angle))
21    }
22
23    pub fn from_hours(hours: f64) -> Self {
24        Self(SiderealAngle::from_hours(hours))
25    }
26
27    pub fn from_degrees(degrees: f64) -> Self {
28        Self(SiderealAngle::from_degrees(degrees))
29    }
30
31    pub fn from_radians(radians: f64) -> Self {
32        Self(SiderealAngle::from_radians(radians))
33    }
34
35    pub fn j2000() -> TimeResult<Self> {
36        let ut1 = UT1::j2000();
37        let tt = TT::j2000();
38        Self::from_ut1_and_tt(&ut1, &tt)
39    }
40
41    pub fn angle(&self) -> SiderealAngle {
42        self.0
43    }
44
45    pub fn hours(&self) -> f64 {
46        self.0.hours()
47    }
48
49    pub fn degrees(&self) -> f64 {
50        self.0.degrees()
51    }
52
53    pub fn radians(&self) -> f64 {
54        self.0.radians()
55    }
56
57    pub fn hour_angle_to_target(&self, target_ra_hours: f64) -> f64 {
58        self.0.hour_angle_to_target(target_ra_hours)
59    }
60
61    pub fn to_lmst(&self, location: &celestial_core::Location) -> crate::sidereal::LMST {
62        use super::angle::SiderealAngle;
63        use crate::sidereal::LMST;
64
65        let gmst_rad = self.radians();
66        let lmst_rad = gmst_rad + location.longitude;
67
68        let lmst_normalized = wrap_0_2pi(lmst_rad);
69        let angle = SiderealAngle::from_radians_exact(lmst_normalized);
70
71        LMST::from_radians(angle.radians(), location)
72    }
73}
74
75impl std::fmt::Display for GMST {
76    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
77        write!(f, "GMST {}", self.0)
78    }
79}
80
81fn calculate_gmst_iau2006(ut1: &UT1, tt: &TT) -> TimeResult<f64> {
82    let ut1_jd = ut1.to_julian_date();
83    let tt_jd = tt.to_julian_date();
84
85    let JulianDate {
86        jd1: ut1_jd1,
87        jd2: ut1_jd2,
88    } = ut1_jd;
89    let JulianDate {
90        jd1: tt_jd1,
91        jd2: tt_jd2,
92    } = tt_jd;
93
94    let t = ((tt_jd1 - J2000_JD) + tt_jd2) / celestial_core::constants::DAYS_PER_JULIAN_CENTURY;
95
96    let era = calculate_era00(ut1_jd1, ut1_jd2)?;
97
98    // IAU 2006 polynomial correction (Horner's method for precision)
99    let polynomial_arcsec = 0.014506
100        + t * (4612.156534
101            + t * (1.3915817 + t * (-0.00000044 + t * (-0.000029956 + t * (-0.0000000368)))));
102
103    let gmst = era + polynomial_arcsec * celestial_core::constants::ARCSEC_TO_RAD;
104
105    Ok(wrap_0_2pi(gmst))
106}
107
108fn calculate_era00(ut1_jd1: f64, ut1_jd2: f64) -> TimeResult<f64> {
109    let (d1, d2) = if ut1_jd1 < ut1_jd2 {
110        (ut1_jd1, ut1_jd2)
111    } else {
112        (ut1_jd2, ut1_jd1)
113    };
114
115    let t = d1 + (d2 - J2000_JD);
116
117    if t.is_infinite() || t.is_nan() || t.abs() > 1e12 {
118        return Err(crate::TimeError::CalculationError(format!(
119            "Time value out of valid range: {} days from J2000",
120            t
121        )));
122    }
123
124    let f = fmod(d1, 1.0) + fmod(d2, 1.0);
125
126    let rotation_term = 0.00273781191135448 * t;
127
128    if rotation_term.is_infinite() || rotation_term.is_nan() {
129        return Err(crate::TimeError::CalculationError(
130            "Earth rotation calculation overflow".to_string(),
131        ));
132    }
133
134    let theta = TWOPI * (f + 0.7790572732640 + rotation_term);
135
136    Ok(wrap_0_2pi(theta))
137}
138
139#[cfg(test)]
140mod tests {
141    use super::*;
142
143    #[test]
144    fn test_gmst_j2000() {
145        let gmst = GMST::j2000().unwrap();
146
147        let hours = gmst.hours();
148        assert!(
149            (0.0..24.0).contains(&hours),
150            "GMST should be in [0, 24) hours: {}",
151            hours
152        );
153        assert!(
154            hours > 18.0 && hours < 19.0,
155            "GMST at J2000.0 should be ~18.7 hours: {}",
156            hours
157        );
158    }
159
160    #[test]
161    fn test_gmst_from_ut1_and_tt() {
162        let ut1 = UT1::j2000();
163        let tt = TT::j2000();
164        let gmst = GMST::from_ut1_and_tt(&ut1, &tt).unwrap();
165
166        let gmst_j2000 = GMST::j2000().unwrap();
167        assert!((gmst.hours() - gmst_j2000.hours()).abs() < 1e-10);
168    }
169
170    #[test]
171    fn test_hour_angle_calculation() {
172        let gmst = GMST::from_hours(12.0);
173        let target_ra = 6.0;
174        let hour_angle = gmst.hour_angle_to_target(target_ra);
175        assert_eq!(hour_angle, 6.0);
176    }
177
178    #[test]
179    fn test_gmst_constructors_and_accessors() {
180        let gmst_deg = GMST::from_degrees(180.0);
181        assert_eq!(gmst_deg.degrees(), 180.0);
182        assert_eq!(gmst_deg.hours(), 12.0);
183
184        let gmst_rad = GMST::from_radians(celestial_core::constants::PI);
185        assert!((gmst_rad.radians() - celestial_core::constants::PI).abs() < 1e-15);
186        assert_eq!(gmst_rad.hours(), 12.0);
187
188        let angle = gmst_deg.angle();
189        assert_eq!(angle.degrees(), 180.0);
190
191        let degrees = gmst_deg.degrees();
192        assert_eq!(degrees, 180.0);
193
194        let display_str = format!("{}", gmst_deg);
195        assert!(display_str.contains("GMST"));
196        assert!(display_str.contains("12.000000h"));
197    }
198
199    #[test]
200    fn test_overflow_protection() {
201        use crate::scales::{TT, UT1};
202        use crate::JulianDate;
203
204        let extreme_jd = JulianDate::new(J2000_JD + 1e13, 0.0);
205        let extreme_ut1 = UT1::from_julian_date(extreme_jd);
206        let extreme_tt = TT::from_julian_date(extreme_jd);
207
208        let result = GMST::from_ut1_and_tt(&extreme_ut1, &extreme_tt);
209        assert!(result.is_err(), "Expected overflow protection to trigger");
210
211        if let Err(err) = result {
212            let error_message = format!("{}", err);
213            assert!(
214                error_message.contains("Time value out of valid range"),
215                "Expected overflow error message, got: {}",
216                error_message
217            );
218        }
219
220        let infinite_jd = JulianDate::new(f64::INFINITY, 0.0);
221        let infinite_ut1 = UT1::from_julian_date(infinite_jd);
222        let infinite_tt = TT::from_julian_date(infinite_jd);
223
224        let result = GMST::from_ut1_and_tt(&infinite_ut1, &infinite_tt);
225        assert!(
226            result.is_err(),
227            "Expected infinite value protection to trigger"
228        );
229
230        let nan_jd = JulianDate::new(f64::NAN, 0.0);
231        let nan_ut1 = UT1::from_julian_date(nan_jd);
232        let nan_tt = TT::from_julian_date(nan_jd);
233
234        let result = GMST::from_ut1_and_tt(&nan_ut1, &nan_tt);
235        assert!(result.is_err(), "Expected NaN value protection to trigger");
236    }
237}
celestial_time/sidereal/gmst.rs

celestial_time/sidereal/
gmst.rs
