celestial_coords/frames/
tirs.rs

1use crate::{frames::ITRSPosition, CoordResult};
2use celestial_core::constants::ARCSEC_TO_RAD;
3use celestial_core::Vector3;
4use celestial_time::{scales::conversions::ToUT1WithDeltaT, transforms::earth_rotation_angle, TT};
5
6#[cfg(feature = "serde")]
7use serde::{Deserialize, Serialize};
8
9#[derive(Debug, Clone, PartialEq)]
10#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
11pub struct TIRSPosition {
12    x: f64,
13    y: f64,
14    z: f64,
15    epoch: TT,
16}
17
18impl TIRSPosition {
19    pub fn new(x: f64, y: f64, z: f64, epoch: TT) -> Self {
20        Self { x, y, z, epoch }
21    }
22
23    pub fn x(&self) -> f64 {
24        self.x
25    }
26
27    pub fn y(&self) -> f64 {
28        self.y
29    }
30
31    pub fn z(&self) -> f64 {
32        self.z
33    }
34
35    pub fn epoch(&self) -> TT {
36        self.epoch
37    }
38
39    pub fn position_vector(&self) -> Vector3 {
40        Vector3::new(self.x, self.y, self.z)
41    }
42
43    pub fn from_position_vector(pos: Vector3, epoch: TT) -> Self {
44        Self::new(pos.x, pos.y, pos.z, epoch)
45    }
46
47    pub fn geocentric_distance(&self) -> f64 {
48        libm::sqrt(self.x * self.x + self.y * self.y + self.z * self.z)
49    }
50
51    pub fn distance_to(&self, other: &Self) -> f64 {
52        let dx = self.x - other.x;
53        let dy = self.y - other.y;
54        let dz = self.z - other.z;
55        libm::sqrt(dx * dx + dy * dy + dz * dz)
56    }
57
58    fn polar_motion_matrix(xp: f64, yp: f64, sp: f64) -> [[f64; 3]; 3] {
59        let mut matrix = [[0.0; 3]; 3];
60        matrix[0][0] = 1.0;
61        matrix[1][1] = 1.0;
62        matrix[2][2] = 1.0;
63
64        Self::apply_rotation_z(&mut matrix, sp);
65        Self::apply_rotation_y(&mut matrix, -xp);
66        Self::apply_rotation_x(&mut matrix, -yp);
67
68        matrix
69    }
70
71    fn apply_rotation_z(matrix: &mut [[f64; 3]; 3], angle: f64) {
72        let (s, c) = libm::sincos(angle);
73        let temp = *matrix;
74
75        for j in 0..3 {
76            matrix[0][j] = c * temp[0][j] + s * temp[1][j];
77            matrix[1][j] = -s * temp[0][j] + c * temp[1][j];
78        }
79    }
80
81    fn apply_rotation_y(matrix: &mut [[f64; 3]; 3], angle: f64) {
82        let (s, c) = libm::sincos(angle);
83        let temp = *matrix;
84
85        for j in 0..3 {
86            matrix[0][j] = c * temp[0][j] - s * temp[2][j];
87            matrix[2][j] = s * temp[0][j] + c * temp[2][j];
88        }
89    }
90
91    fn apply_rotation_x(matrix: &mut [[f64; 3]; 3], angle: f64) {
92        let (s, c) = libm::sincos(angle);
93        let temp = *matrix;
94
95        for j in 0..3 {
96            matrix[1][j] = c * temp[1][j] + s * temp[2][j];
97            matrix[2][j] = -s * temp[1][j] + c * temp[2][j];
98        }
99    }
100
101    /// Computes ΔT (TT - UT1) using EOP parameters.
102    ///
103    /// # EOP Freshness Check
104    /// Rejects EOP data >1 day from the target epoch. This ensures UT1-UTC is current, as
105    /// Earth's rotation is irregular. For sparse datasets, consider:
106    /// - Using an interpolator to fill gaps (see `EopManager`)
107    /// - Relaxing this check if lower precision is acceptable (modify this threshold)
108    /// - Pre-fetching/caching EOP data for your observation window
109    ///
110    /// Typical use cases handle this via interpolation, so sparse raw data should be rare.
111    pub fn compute_delta_t(epoch: &TT, eop: &crate::eop::EopParameters) -> CoordResult<f64> {
112        use celestial_time::scales::conversions::{ToTAI, ToUTC};
113
114        let epoch_jd = epoch.to_julian_date().to_f64();
115        let eop_jd = eop.mjd + celestial_core::constants::MJD_ZERO_POINT;
116
117        if (epoch_jd - eop_jd).abs() > 1.0 {
118            return Err(crate::CoordError::data_unavailable(
119                "EOP data is more than 1 day from epoch - Delta-T computation may be inaccurate",
120            ));
121        }
122
123        let tai = epoch
124            .to_tai()
125            .map_err(|e| crate::CoordError::external_library("TT to TAI", &e.to_string()))?;
126        let utc = tai
127            .to_utc()
128            .map_err(|e| crate::CoordError::external_library("TAI to UTC", &e.to_string()))?;
129
130        let utc_jd = utc.to_julian_date().to_f64();
131        let ut1_jd = utc_jd + eop.ut1_utc / celestial_core::constants::SECONDS_PER_DAY_F64;
132
133        let delta_t_days = epoch_jd - ut1_jd;
134        let delta_t_seconds = delta_t_days * celestial_core::constants::SECONDS_PER_DAY_F64;
135
136        Ok(delta_t_seconds)
137    }
138
139    pub fn to_itrs(
140        &self,
141        epoch: &TT,
142        eop: &crate::eop::EopParameters,
143    ) -> CoordResult<ITRSPosition> {
144        let delta_t_seconds = Self::compute_delta_t(epoch, eop)?;
145        let ut1 = epoch.to_ut1_with_delta_t(delta_t_seconds)?;
146        let era = earth_rotation_angle(&ut1.to_julian_date())?;
147
148        let (sin_era, cos_era) = libm::sincos(era);
149
150        let x_after_era = cos_era * self.x + sin_era * self.y;
151        let y_after_era = -sin_era * self.x + cos_era * self.y;
152        let z_after_era = self.z;
153
154        let xp_rad = eop.x_p * ARCSEC_TO_RAD;
155        let yp_rad = eop.y_p * ARCSEC_TO_RAD;
156
157        let polar_motion_matrix = Self::polar_motion_matrix(xp_rad, yp_rad, eop.s_prime);
158
159        let x_itrs = polar_motion_matrix[0][0] * x_after_era
160            + polar_motion_matrix[0][1] * y_after_era
161            + polar_motion_matrix[0][2] * z_after_era;
162        let y_itrs = polar_motion_matrix[1][0] * x_after_era
163            + polar_motion_matrix[1][1] * y_after_era
164            + polar_motion_matrix[1][2] * z_after_era;
165        let z_itrs = polar_motion_matrix[2][0] * x_after_era
166            + polar_motion_matrix[2][1] * y_after_era
167            + polar_motion_matrix[2][2] * z_after_era;
168
169        Ok(ITRSPosition::new(x_itrs, y_itrs, z_itrs, *epoch))
170    }
171
172    pub fn from_itrs(
173        itrs: &ITRSPosition,
174        epoch: &TT,
175        eop: &crate::eop::EopParameters,
176    ) -> CoordResult<Self> {
177        let xp_rad = eop.x_p * ARCSEC_TO_RAD;
178        let yp_rad = eop.y_p * ARCSEC_TO_RAD;
179
180        let polar_motion_matrix = Self::polar_motion_matrix(xp_rad, yp_rad, eop.s_prime);
181
182        let x_before_pm = polar_motion_matrix[0][0] * itrs.x()
183            + polar_motion_matrix[1][0] * itrs.y()
184            + polar_motion_matrix[2][0] * itrs.z();
185        let y_before_pm = polar_motion_matrix[0][1] * itrs.x()
186            + polar_motion_matrix[1][1] * itrs.y()
187            + polar_motion_matrix[2][1] * itrs.z();
188        let z_before_pm = polar_motion_matrix[0][2] * itrs.x()
189            + polar_motion_matrix[1][2] * itrs.y()
190            + polar_motion_matrix[2][2] * itrs.z();
191
192        let delta_t_seconds = Self::compute_delta_t(epoch, eop)?;
193        let ut1 = epoch.to_ut1_with_delta_t(delta_t_seconds)?;
194        let era = earth_rotation_angle(&ut1.to_julian_date())?;
195
196        let (sin_era, cos_era) = libm::sincos(era);
197
198        let x_tirs = cos_era * x_before_pm - sin_era * y_before_pm;
199        let y_tirs = sin_era * x_before_pm + cos_era * y_before_pm;
200        let z_tirs = z_before_pm;
201
202        Ok(Self::new(x_tirs, y_tirs, z_tirs, *epoch))
203    }
204
205    pub fn from_cirs(
206        cirs_vec: Vector3,
207        epoch: &TT,
208        eop: &crate::eop::EopParameters,
209    ) -> CoordResult<Self> {
210        let delta_t_seconds = Self::compute_delta_t(epoch, eop)?;
211        let ut1 = epoch.to_ut1_with_delta_t(delta_t_seconds)?;
212        let era = earth_rotation_angle(&ut1.to_julian_date())?;
213
214        let (sin_era, cos_era) = libm::sincos(era);
215
216        let x_tirs = cos_era * cirs_vec.x - sin_era * cirs_vec.y;
217        let y_tirs = sin_era * cirs_vec.x + cos_era * cirs_vec.y;
218        let z_tirs = cirs_vec.z;
219
220        Ok(Self::new(x_tirs, y_tirs, z_tirs, *epoch))
221    }
222
223    pub fn to_cirs(
224        &self,
225        eop: &crate::eop::EopParameters,
226    ) -> CoordResult<crate::frames::CIRSPosition> {
227        use crate::frames::CIRSPosition;
228
229        let delta_t_seconds = Self::compute_delta_t(&self.epoch, eop)?;
230        let ut1 = self.epoch.to_ut1_with_delta_t(delta_t_seconds)?;
231        let era = earth_rotation_angle(&ut1.to_julian_date())?;
232
233        let (sin_era, cos_era) = libm::sincos(era);
234
235        let x_cirs = cos_era * self.x + sin_era * self.y;
236        let y_cirs = -sin_era * self.x + cos_era * self.y;
237        let z_cirs = self.z;
238
239        let cirs_vec = Vector3::new(x_cirs, y_cirs, z_cirs);
240        CIRSPosition::from_unit_vector(cirs_vec, self.epoch)
241    }
242}
243
244impl std::fmt::Display for TIRSPosition {
245    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
246        write!(
247            f,
248            "TIRS(X={:.3}m, Y={:.3}m, Z={:.3}m, epoch=J{:.1})",
249            self.x,
250            self.y,
251            self.z,
252            self.epoch.julian_year()
253        )
254    }
255}
256
257#[cfg(test)]
258mod tests {
259    use super::*;
260    use celestial_core::constants::TWOPI;
261
262    #[test]
263    fn test_tirs_creation() {
264        let epoch = TT::j2000();
265        let pos = TIRSPosition::new(1000000.0, 2000000.0, 3000000.0, epoch);
266
267        assert_eq!(pos.x(), 1000000.0);
268        assert_eq!(pos.y(), 2000000.0);
269        assert_eq!(pos.z(), 3000000.0);
270        assert_eq!(pos.epoch(), epoch);
271    }
272
273    #[test]
274    fn test_vector_operations() {
275        let epoch = TT::j2000();
276        let original = TIRSPosition::new(1000.0, 2000.0, 3000.0, epoch);
277
278        let vec = original.position_vector();
279        assert_eq!(vec.x, 1000.0);
280        assert_eq!(vec.y, 2000.0);
281        assert_eq!(vec.z, 3000.0);
282
283        let recovered = TIRSPosition::from_position_vector(vec, epoch);
284        assert_eq!(recovered.x(), original.x());
285        assert_eq!(recovered.y(), original.y());
286        assert_eq!(recovered.z(), original.z());
287    }
288
289    #[test]
290    fn test_distance_calculations() {
291        let epoch = TT::j2000();
292
293        let pos1 = TIRSPosition::new(1000.0, 0.0, 0.0, epoch);
294        let pos2 = TIRSPosition::new(2000.0, 0.0, 0.0, epoch);
295
296        // Distance between positions
297        assert_eq!(pos1.distance_to(&pos2), 1000.0);
298        assert_eq!(pos2.distance_to(&pos1), 1000.0);
299
300        // Distance from origin
301        assert_eq!(pos1.geocentric_distance(), 1000.0);
302        assert_eq!(pos2.geocentric_distance(), 2000.0);
303    }
304
305    #[test]
306    fn test_itrs_transformation_roundtrip() {
307        use crate::eop::EopRecord;
308
309        let epoch = TT::j2000();
310        let original_tirs = TIRSPosition::new(4000000.0, 3000000.0, 5000000.0, epoch);
311
312        let eop = EopRecord::new(51544.5, 0.0, 0.0, 0.3, 0.0)
313            .unwrap()
314            .to_parameters();
315
316        // Transform to ITRS and back
317        let itrs = original_tirs.to_itrs(&epoch, &eop).unwrap();
318        let recovered_tirs = TIRSPosition::from_itrs(&itrs, &epoch, &eop).unwrap();
319
320        // Should roundtrip with floating-point precision (rotation with EOP-based Delta-T)
321        // Allow 2 ULP due to additional operations in Delta-T calculation
322        celestial_core::test_helpers::assert_ulp_le(
323            recovered_tirs.x(),
324            original_tirs.x(),
325            2,
326            "X roundtrip",
327        );
328        celestial_core::test_helpers::assert_ulp_le(
329            recovered_tirs.y(),
330            original_tirs.y(),
331            2,
332            "Y roundtrip",
333        );
334        celestial_core::test_helpers::assert_ulp_le(
335            recovered_tirs.z(),
336            original_tirs.z(),
337            2,
338            "Z roundtrip",
339        );
340    }
341
342    #[test]
343    fn test_itrs_transformation_properties() {
344        use crate::eop::EopRecord;
345
346        let epoch = TT::j2000();
347        let tirs = TIRSPosition::new(1000000.0, 0.0, 0.0, epoch);
348
349        let eop = EopRecord::new(51544.5, 0.0, 0.0, 0.3, 0.0)
350            .unwrap()
351            .to_parameters();
352
353        let itrs = tirs.to_itrs(&epoch, &eop).unwrap();
354
355        // Z coordinate should be unchanged (rotation around Z-axis)
356        assert_eq!(itrs.z(), tirs.z());
357
358        // Distance from origin should be preserved
359        celestial_core::test_helpers::assert_ulp_le(
360            itrs.geocentric_distance(),
361            tirs.geocentric_distance(),
362            1,
363            "Distance preservation",
364        );
365
366        // X and Y should change due to Earth rotation (unless ERA = 0)
367        let delta_t_seconds = 69.184;
368        let ut1 = epoch.to_ut1_with_delta_t(delta_t_seconds).unwrap();
369        let era = earth_rotation_angle(&ut1.to_julian_date()).unwrap();
370        if era != 0.0 {
371            // Only test if ERA is non-zero (it should be for J2000)
372            assert!(itrs.x() != tirs.x() || itrs.y() != tirs.y());
373        }
374    }
375
376    #[test]
377    fn test_display_formatting() {
378        let epoch = TT::j2000();
379        let pos = TIRSPosition::new(1234567.89, -987654.32, 555666.77, epoch);
380
381        let display = format!("{}", pos);
382        assert!(display.contains("TIRS"));
383        assert!(display.contains("1234567.890m"));
384        assert!(display.contains("-987654.320m"));
385        assert!(display.contains("555666.770m"));
386        assert!(display.contains("J2000.0"));
387    }
388
389    #[test]
390    fn test_itrs_consistency() {
391        use crate::eop::EopRecord;
392
393        let epoch = TT::j2000();
394        let itrs_original = ITRSPosition::new(2000000.0, 1000000.0, 6000000.0, epoch);
395
396        let eop = EopRecord::new(51544.5, 0.0, 0.0, 0.3, 0.0)
397            .unwrap()
398            .to_parameters();
399
400        // Transform ITRS -> TIRS -> ITRS
401        let tirs = TIRSPosition::from_itrs(&itrs_original, &epoch, &eop).unwrap();
402        let itrs_recovered = tirs.to_itrs(&epoch, &eop).unwrap();
403
404        // Should roundtrip with floating-point precision
405        celestial_core::test_helpers::assert_ulp_le(
406            itrs_recovered.x(),
407            itrs_original.x(),
408            1,
409            "ITRS X roundtrip",
410        );
411        celestial_core::test_helpers::assert_ulp_le(
412            itrs_recovered.y(),
413            itrs_original.y(),
414            1,
415            "ITRS Y roundtrip",
416        );
417        celestial_core::test_helpers::assert_ulp_le(
418            itrs_recovered.z(),
419            itrs_original.z(),
420            1,
421            "ITRS Z roundtrip",
422        );
423    }
424
425    #[test]
426    fn test_earth_rotation_angle_usage() {
427        use crate::eop::EopRecord;
428
429        let epoch = TT::j2000();
430
431        let eop = EopRecord::new(51544.5, 0.0, 0.0, 0.3, 0.0)
432            .unwrap()
433            .to_parameters();
434
435        let delta_t_seconds = TIRSPosition::compute_delta_t(&epoch, &eop).unwrap();
436        let ut1 = epoch.to_ut1_with_delta_t(delta_t_seconds).unwrap();
437        let era = earth_rotation_angle(&ut1.to_julian_date()).unwrap();
438
439        assert!(era >= 0.0);
440        assert!(era < TWOPI);
441
442        let tirs_x_axis = TIRSPosition::new(6378137.0, 0.0, 0.0, epoch);
443        let itrs = tirs_x_axis.to_itrs(&epoch, &eop).unwrap();
444
445        let expected_x = 6378137.0 * libm::cos(era);
446        let expected_y = -6378137.0 * libm::sin(era);
447
448        celestial_core::test_helpers::assert_ulp_le(
449            itrs.x(),
450            expected_x,
451            1,
452            "ERA X transformation",
453        );
454        celestial_core::test_helpers::assert_ulp_le(
455            itrs.y(),
456            expected_y,
457            1,
458            "ERA Y transformation",
459        );
460        assert_eq!(itrs.z(), 0.0);
461    }
462
463    #[test]
464    fn test_eop_timing_check() {
465        use crate::eop::EopRecord;
466
467        let epoch = TT::j2000();
468        let eop_old = EopRecord::new(51542.0, 0.0, 0.0, 0.3, 0.0)
469            .unwrap()
470            .to_parameters();
471
472        let result = TIRSPosition::compute_delta_t(&epoch, &eop_old);
473        assert!(result.is_err());
474    }
475
476    #[test]
477    fn test_cirs_transformation() {
478        use crate::eop::EopRecord;
479
480        let epoch = TT::j2000();
481        let eop = EopRecord::new(51544.5, 0.0, 0.0, 0.3, 0.0)
482            .unwrap()
483            .to_parameters();
484
485        let tirs = TIRSPosition::new(6378137.0, 0.0, 0.0, epoch);
486
487        let cirs = tirs.to_cirs(&eop).unwrap();
488
489        assert!(cirs.ra().degrees() >= 0.0);
490        assert!(cirs.ra().degrees() < 360.0);
491        assert!(cirs.dec().degrees() >= -90.0);
492        assert!(cirs.dec().degrees() <= 90.0);
493    }
494}
celestial_coords/frames/tirs.rs

celestial_coords/frames/
tirs.rs
