celestial_coords/eop/
record.rs

1use crate::{CoordError, CoordResult};
2use celestial_core::constants::{
3    ARCSEC_TO_RAD, DAYS_PER_JULIAN_CENTURY, J2000_JD, MILLIARCSEC_TO_RAD, MJD_ZERO_POINT,
4};
5use celestial_time::{transforms::earth_rotation_angle, JulianDate};
6
7#[cfg(feature = "serde")]
8use serde::{Deserialize, Serialize};
9
10#[derive(Debug, Clone, PartialEq)]
11#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
12pub struct EopRecord {
13    pub mjd: f64,
14
15    pub x_p_encoded: i32,
16
17    pub y_p_encoded: i32,
18
19    pub ut1_utc_encoded: i32,
20
21    pub lod_encoded: i32,
22
23    pub dx_encoded: Option<i16>,
24
25    pub dy_encoded: Option<i16>,
26
27    pub xrt_encoded: Option<i32>,
28
29    pub yrt_encoded: Option<i32>,
30
31    pub flags: EopFlags,
32}
33
34#[derive(Debug, Clone, Copy, PartialEq, Eq)]
35#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
36pub struct EopFlags {
37    pub source: EopSource,
38
39    pub quality: EopQuality,
40
41    pub has_polar_motion: bool,
42
43    pub has_ut1_utc: bool,
44
45    pub has_cip_offsets: bool,
46
47    pub has_pole_rates: bool,
48}
49
50#[derive(Debug, Clone, Copy, PartialEq, Eq)]
51#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
52pub enum EopSource {
53    IersC04,
54
55    IersFinals,
56
57    IersPrediction,
58
59    UserData,
60
61    Interpolated,
62}
63
64#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
65#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
66pub enum EopQuality {
67    HighPrecision,
68
69    Standard,
70
71    LowPrecision,
72
73    Predicted,
74}
75
76impl EopRecord {
77    const ARCSEC_TO_UNITS: f64 = 10_000_000.0;
78    const SEC_TO_UNITS: f64 = 10_000_000.0;
79    const MILLIARCSEC_TO_UNITS: f64 = 10_000.0;
80
81    pub fn new(
82        mjd: f64,
83        x_p_arcsec: f64,
84        y_p_arcsec: f64,
85        ut1_utc_sec: f64,
86        lod_sec: f64,
87    ) -> CoordResult<Self> {
88        if x_p_arcsec.abs() > 6.0 {
89            return Err(CoordError::invalid_coordinate(format!(
90                "X polar motion out of range: {} arcsec",
91                x_p_arcsec
92            )));
93        }
94
95        if y_p_arcsec.abs() > 6.0 {
96            return Err(CoordError::invalid_coordinate(format!(
97                "Y polar motion out of range: {} arcsec",
98                y_p_arcsec
99            )));
100        }
101
102        if ut1_utc_sec.abs() > 1.0 {
103            return Err(CoordError::invalid_coordinate(format!(
104                "UT1-UTC out of range: {} sec",
105                ut1_utc_sec
106            )));
107        }
108
109        if lod_sec.abs() > 0.01 {
110            return Err(CoordError::invalid_coordinate(format!(
111                "LOD out of range: {} sec",
112                lod_sec
113            )));
114        }
115
116        Ok(Self {
117            mjd,
118            x_p_encoded: libm::round(x_p_arcsec * Self::ARCSEC_TO_UNITS) as i32,
119            y_p_encoded: libm::round(y_p_arcsec * Self::ARCSEC_TO_UNITS) as i32,
120            ut1_utc_encoded: libm::round(ut1_utc_sec * Self::SEC_TO_UNITS) as i32,
121            lod_encoded: libm::round(lod_sec * Self::SEC_TO_UNITS) as i32,
122            dx_encoded: None,
123            dy_encoded: None,
124            xrt_encoded: None,
125            yrt_encoded: None,
126            flags: EopFlags::default(),
127        })
128    }
129
130    pub fn with_cip_offsets(
131        mut self,
132        dx_milliarcsec: f64,
133        dy_milliarcsec: f64,
134    ) -> CoordResult<Self> {
135        if dx_milliarcsec.abs() > 1000.0 {
136            return Err(CoordError::invalid_coordinate(format!(
137                "CIP dX out of range: {} mas",
138                dx_milliarcsec
139            )));
140        }
141
142        if dy_milliarcsec.abs() > 1000.0 {
143            return Err(CoordError::invalid_coordinate(format!(
144                "CIP dY out of range: {} mas",
145                dy_milliarcsec
146            )));
147        }
148
149        self.dx_encoded = Some(libm::round(dx_milliarcsec * Self::MILLIARCSEC_TO_UNITS) as i16);
150        self.dy_encoded = Some(libm::round(dy_milliarcsec * Self::MILLIARCSEC_TO_UNITS) as i16);
151        self.flags.has_cip_offsets = true;
152
153        Ok(self)
154    }
155
156    pub fn with_pole_rates(
157        mut self,
158        xrt_arcsec_per_day: f64,
159        yrt_arcsec_per_day: f64,
160    ) -> CoordResult<Self> {
161        if xrt_arcsec_per_day.abs() > 1.0 {
162            return Err(CoordError::invalid_coordinate(format!(
163                "Pole rate xrt out of range: {} arcsec/day",
164                xrt_arcsec_per_day
165            )));
166        }
167        if yrt_arcsec_per_day.abs() > 1.0 {
168            return Err(CoordError::invalid_coordinate(format!(
169                "Pole rate yrt out of range: {} arcsec/day",
170                yrt_arcsec_per_day
171            )));
172        }
173        self.xrt_encoded = Some(libm::round(xrt_arcsec_per_day * Self::ARCSEC_TO_UNITS) as i32);
174        self.yrt_encoded = Some(libm::round(yrt_arcsec_per_day * Self::ARCSEC_TO_UNITS) as i32);
175        self.flags.has_pole_rates = true;
176        Ok(self)
177    }
178
179    pub fn with_flags(mut self, flags: EopFlags) -> Self {
180        self.flags = flags;
181        self
182    }
183
184    pub fn to_parameters(&self) -> EopParameters {
185        EopParameters {
186            mjd: self.mjd,
187            x_p: self.x_p_encoded as f64 / Self::ARCSEC_TO_UNITS,
188            y_p: self.y_p_encoded as f64 / Self::ARCSEC_TO_UNITS,
189            ut1_utc: self.ut1_utc_encoded as f64 / Self::SEC_TO_UNITS,
190            lod: self.lod_encoded as f64 / Self::SEC_TO_UNITS,
191            dx: self
192                .dx_encoded
193                .map(|v| v as f64 / Self::MILLIARCSEC_TO_UNITS),
194            dy: self
195                .dy_encoded
196                .map(|v| v as f64 / Self::MILLIARCSEC_TO_UNITS),
197            xrt: self.xrt_encoded.map(|v| v as f64 / Self::ARCSEC_TO_UNITS),
198            yrt: self.yrt_encoded.map(|v| v as f64 / Self::ARCSEC_TO_UNITS),
199            s_prime: 0.0,
200            flags: self.flags,
201        }
202    }
203}
204
205#[derive(Debug, Clone, PartialEq)]
206pub struct EopParameters {
207    pub mjd: f64,
208
209    pub x_p: f64,
210
211    pub y_p: f64,
212
213    pub ut1_utc: f64,
214
215    pub lod: f64,
216
217    pub dx: Option<f64>,
218
219    pub dy: Option<f64>,
220
221    pub xrt: Option<f64>,
222
223    pub yrt: Option<f64>,
224
225    pub s_prime: f64,
226
227    pub flags: EopFlags,
228}
229
230impl EopParameters {
231    /// Computes the TIO locator s' using the IAU 2000 linear approximation.
232    ///
233    /// s' ≈ -47 µas/century × t, where t is Julian centuries from J2000.0.
234    /// Result is in radians (matching IAU SOFA convention).
235    ///
236    /// Note: For maximum precision, use `compute_s_prime_jd()` with 2-part JD.
237    pub fn compute_s_prime(&mut self) {
238        self.compute_s_prime_jd(self.mjd + MJD_ZERO_POINT, 0.0);
239    }
240
241    /// Computes s' from 2-part TT Julian Date for maximum precision.
242    ///
243    /// The 2-part JD allows for precision-preserving arithmetic when
244    /// computing time intervals from J2000.0.
245    pub fn compute_s_prime_jd(&mut self, tt1: f64, tt2: f64) {
246        let t = ((tt1 - J2000_JD) + tt2) / DAYS_PER_JULIAN_CENTURY;
247        self.s_prime = -47e-6 * t * ARCSEC_TO_RAD;
248    }
249
250    pub fn compute_era(&self) -> CoordResult<f64> {
251        let ut1_jd = self.mjd
252            + MJD_ZERO_POINT
253            + self.ut1_utc / celestial_core::constants::SECONDS_PER_DAY_F64;
254
255        let ut1_jd1 = libm::floor(ut1_jd);
256        let ut1_jd2 = ut1_jd - ut1_jd1;
257
258        let jd = JulianDate::new(ut1_jd1, ut1_jd2);
259        earth_rotation_angle(&jd)
260            .map_err(|e| CoordError::external_library("Earth Rotation Angle", &e.to_string()))
261    }
262
263    /// Returns CIP X coordinate corrected by dX offset (if available), in radians.
264    pub fn corrected_cip_x(&self, x_iau: f64) -> f64 {
265        x_iau + self.dx.unwrap_or(0.0) * MILLIARCSEC_TO_RAD
266    }
267
268    /// Returns CIP Y coordinate corrected by dY offset (if available), in radians.
269    pub fn corrected_cip_y(&self, y_iau: f64) -> f64 {
270        y_iau + self.dy.unwrap_or(0.0) * MILLIARCSEC_TO_RAD
271    }
272}
273
274impl Default for EopFlags {
275    fn default() -> Self {
276        Self {
277            source: EopSource::UserData,
278            quality: EopQuality::Standard,
279            has_polar_motion: true,
280            has_ut1_utc: true,
281            has_cip_offsets: false,
282            has_pole_rates: false,
283        }
284    }
285}
286
287impl std::fmt::Display for EopParameters {
288    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
289        write!(
290            f,
291            "EOP(MJD={:.1}, xp={:.6}\", yp={:.6}\", UT1-UTC={:.7}s)",
292            self.mjd, self.x_p, self.y_p, self.ut1_utc
293        )
294    }
295}
296
297#[cfg(test)]
298mod tests {
299    use super::*;
300
301    #[test]
302    fn test_eop_record_encoding() {
303        let record = EopRecord::new(
304            59945.0,   // MJD 2023-01-01
305            0.123456,  // x_p arcsec
306            0.234567,  // y_p arcsec
307            0.0123456, // UT1-UTC sec
308            0.0012345, // LOD sec
309        )
310        .unwrap();
311
312        let params = record.to_parameters();
313
314        // Check precision preservation (should be within 0.1 µas for angles)
315        assert!((params.x_p - 0.123456).abs() == 0.0);
316        assert!((params.y_p - 0.234567).abs() == 0.0);
317        assert!((params.ut1_utc - 0.0123456).abs() == 0.0);
318        assert!((params.lod - 0.0012345).abs() == 0.0);
319        assert_eq!(params.mjd, 59945.0);
320    }
321
322    #[test]
323    fn test_cip_offsets() {
324        let record = EopRecord::new(59945.0, 0.1, 0.2, 0.01, 0.001)
325            .unwrap()
326            .with_cip_offsets(0.5, -0.3) // mas
327            .unwrap();
328
329        let params = record.to_parameters();
330
331        assert_eq!(params.dx, Some(0.5));
332        assert_eq!(params.dy, Some(-0.3));
333        assert!(params.flags.has_cip_offsets);
334    }
335
336    #[test]
337    fn test_parameter_display() {
338        let mut params = EopParameters {
339            mjd: 59945.0,
340            x_p: 0.123456,
341            y_p: 0.234567,
342            ut1_utc: 0.0123456,
343            lod: 0.001,
344            dx: None,
345            dy: None,
346            xrt: None,
347            yrt: None,
348            s_prime: 0.0,
349            flags: EopFlags::default(),
350        };
351
352        params.compute_s_prime();
353
354        let display = format!("{}", params);
355        assert!(display.contains("MJD=59945.0"));
356        assert!(display.contains("xp=0.123456"));
357        assert!(display.contains("yp=0.234567"));
358    }
359
360    #[test]
361    fn test_validation_x_polar_motion_out_of_range() {
362        let result = EopRecord::new(59945.0, 6.1, 0.2, 0.01, 0.001);
363        assert!(result.is_err());
364        let err = result.unwrap_err();
365        assert!(err.to_string().contains("X polar motion out of range"));
366    }
367
368    #[test]
369    fn test_validation_y_polar_motion_out_of_range() {
370        let result = EopRecord::new(59945.0, 0.1, -6.1, 0.01, 0.001);
371        assert!(result.is_err());
372        let err = result.unwrap_err();
373        assert!(err.to_string().contains("Y polar motion out of range"));
374    }
375
376    #[test]
377    fn test_validation_ut1_utc_out_of_range() {
378        let result = EopRecord::new(59945.0, 0.1, 0.2, 1.1, 0.001);
379        assert!(result.is_err());
380        let err = result.unwrap_err();
381        assert!(err.to_string().contains("UT1-UTC out of range"));
382    }
383
384    #[test]
385    fn test_validation_lod_out_of_range() {
386        let result = EopRecord::new(59945.0, 0.1, 0.2, 0.01, 0.011);
387        assert!(result.is_err());
388        let err = result.unwrap_err();
389        assert!(err.to_string().contains("LOD out of range"));
390    }
391
392    #[test]
393    fn test_validation_cip_dx_out_of_range() {
394        let result = EopRecord::new(59945.0, 0.1, 0.2, 0.01, 0.001)
395            .unwrap()
396            .with_cip_offsets(1001.0, 0.0);
397        assert!(result.is_err());
398        let err = result.unwrap_err();
399        assert!(err.to_string().contains("CIP dX out of range"));
400    }
401
402    #[test]
403    fn test_validation_cip_dy_out_of_range() {
404        let result = EopRecord::new(59945.0, 0.1, 0.2, 0.01, 0.001)
405            .unwrap()
406            .with_cip_offsets(0.0, -1001.0);
407        assert!(result.is_err());
408        let err = result.unwrap_err();
409        assert!(err.to_string().contains("CIP dY out of range"));
410    }
411}
celestial_coords/eop/record.rs

celestial_coords/eop/
record.rs
