Skip to main content

celestial_time/sidereal/
gast.rs

1use super::angle::SiderealAngle;
2use crate::scales::{TT, UT1};
3use crate::sidereal::LAST;
4use crate::transforms::earth_rotation_angle;
5use crate::transforms::nutation::NutationCalculator;
6use crate::TimeResult;
7use celestial_core::angle::wrap_0_2pi;
8use celestial_core::cio::CioSolution;
9
10#[cfg(feature = "serde")]
11use serde::{Deserialize, Serialize};
12
13#[derive(Debug, Clone, Copy, PartialEq)]
14#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
15pub struct GAST(SiderealAngle);
16
17impl GAST {
18    pub fn from_ut1_and_tt(ut1: &UT1, tt: &TT) -> TimeResult<Self> {
19        let gast_rad = calculate_gast_iau2006a(ut1, tt)?;
20        let angle = SiderealAngle::from_radians_exact(gast_rad);
21        Ok(Self(angle))
22    }
23
24    pub fn from_hours(hours: f64) -> Self {
25        Self(SiderealAngle::from_hours(hours))
26    }
27
28    pub fn from_degrees(degrees: f64) -> Self {
29        Self(SiderealAngle::from_degrees(degrees))
30    }
31
32    pub fn from_radians(radians: f64) -> Self {
33        Self(SiderealAngle::from_radians(radians))
34    }
35
36    pub fn j2000() -> TimeResult<Self> {
37        let ut1 = UT1::j2000();
38        let tt = TT::j2000();
39        Self::from_ut1_and_tt(&ut1, &tt)
40    }
41
42    pub fn angle(&self) -> SiderealAngle {
43        self.0
44    }
45
46    pub fn hours(&self) -> f64 {
47        self.0.hours()
48    }
49
50    pub fn degrees(&self) -> f64 {
51        self.0.degrees()
52    }
53
54    pub fn radians(&self) -> f64 {
55        self.0.radians()
56    }
57
58    pub fn hour_angle_to_target(&self, target_ra_hours: f64) -> f64 {
59        self.0.hour_angle_to_target(target_ra_hours)
60    }
61
62    pub fn to_last(&self, location: &celestial_core::Location) -> crate::sidereal::LAST {
63        let gast_rad = self.radians();
64        let last_rad = gast_rad + location.longitude;
65
66        let last_normalized = wrap_0_2pi(last_rad);
67        let angle = SiderealAngle::from_radians_exact(last_normalized);
68
69        LAST::from_radians(angle.radians(), location)
70    }
71}
72
73impl std::fmt::Display for GAST {
74    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
75        write!(f, "GAST {}", self.0)
76    }
77}
78
79fn calculate_gast_iau2006a(ut1: &UT1, tt: &TT) -> TimeResult<f64> {
80    let era = earth_rotation_angle(&ut1.to_julian_date())?;
81
82    let tt_centuries = tt_to_centuries(tt)?;
83    let npb_matrix = calculate_npb_matrix_iau2006a(tt)?;
84
85    let cio_solution = CioSolution::calculate(&npb_matrix, tt_centuries).map_err(|e| {
86        crate::TimeError::CalculationError(format!("CIO calculation failed: {}", e))
87    })?;
88
89    let gast = era - cio_solution.equation_of_origins;
90
91    Ok(wrap_0_2pi(gast))
92}
93
94fn calculate_npb_matrix_iau2006a(tt: &TT) -> TimeResult<celestial_core::RotationMatrix3> {
95    let tt_jd = tt.to_julian_date();
96    let t = celestial_core::utils::jd_to_centuries(tt_jd.jd1(), tt_jd.jd2());
97
98    let nutation_result = tt.nutation_iau2006a()?;
99    let dpsi = nutation_result.nutation_longitude();
100    let deps = nutation_result.nutation_obliquity();
101
102    let precession_calc = celestial_core::precession::PrecessionIAU2006::new();
103    let npb_matrix = precession_calc.npb_matrix_iau2006a(t, dpsi, deps);
104
105    Ok(npb_matrix)
106}
107
108fn tt_to_centuries(tt: &TT) -> TimeResult<f64> {
109    crate::transforms::nutation::tt_to_centuries(tt)
110}
111
112#[cfg(test)]
113mod tests {
114    use super::*;
115
116    #[test]
117    fn test_gast_j2000() {
118        let gast = GAST::j2000().unwrap();
119
120        let hours = gast.hours();
121        assert!(
122            (0.0..24.0).contains(&hours),
123            "GAST should be in [0, 24) hours: {}",
124            hours
125        );
126        assert!(
127            hours > 18.0 && hours < 19.0,
128            "GAST at J2000.0 should be ~18.7 hours: {}",
129            hours
130        );
131    }
132
133    #[test]
134    fn test_gast_from_ut1_and_tt() {
135        let ut1 = UT1::j2000();
136        let tt = TT::j2000();
137        let gast = GAST::from_ut1_and_tt(&ut1, &tt).unwrap();
138
139        let gast_j2000 = GAST::j2000().unwrap();
140        assert!((gast.hours() - gast_j2000.hours()).abs() < 1e-10);
141    }
142
143    #[test]
144    fn test_hour_angle_calculation() {
145        let gast = GAST::from_hours(12.0);
146        let target_ra = 6.0;
147        let hour_angle = gast.hour_angle_to_target(target_ra);
148        assert_eq!(hour_angle, 6.0);
149    }
150
151    #[test]
152    fn test_gast_constructors_and_accessors() {
153        let gast_deg = GAST::from_degrees(270.0);
154        assert_eq!(gast_deg.degrees(), 270.0);
155        assert_eq!(gast_deg.hours(), 18.0);
156
157        let gast_rad = GAST::from_radians(celestial_core::constants::HALF_PI);
158        assert!((gast_rad.radians() - celestial_core::constants::HALF_PI).abs() < 1e-15);
159        assert_eq!(gast_rad.hours(), 6.0);
160
161        let angle = gast_deg.angle();
162        assert_eq!(angle.degrees(), 270.0);
163
164        let degrees = gast_deg.degrees();
165        assert_eq!(degrees, 270.0);
166
167        let display_str = format!("{}", gast_deg);
168        assert!(display_str.contains("GAST"));
169        assert!(display_str.contains("18.000000h"));
170    }
171}