celestial_coords/frames/
itrs.rs

1use crate::CoordResult;
2use celestial_core::{Angle, Vector3};
3use celestial_time::TT;
4
5#[cfg(feature = "serde")]
6use serde::{Deserialize, Serialize};
7
8#[derive(Debug, Clone, PartialEq)]
9#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
10pub struct ITRSPosition {
11    x: f64,
12    y: f64,
13    z: f64,
14    epoch: TT,
15}
16
17impl ITRSPosition {
18    pub fn new(x: f64, y: f64, z: f64, epoch: TT) -> Self {
19        Self { x, y, z, epoch }
20    }
21
22    pub fn from_geodetic(
23        longitude: Angle,
24        latitude: Angle,
25        height: f64,
26        epoch: TT,
27    ) -> CoordResult<Self> {
28        const A: f64 = 6378137.0;
29        const F: f64 = 1.0 / 298.257223563;
30
31        let (sin_lat, cos_lat) = latitude.sin_cos();
32        let (sin_lon, cos_lon) = longitude.sin_cos();
33
34        let w = 1.0 - F;
35        let w2 = w * w;
36        let d = cos_lat * cos_lat + w2 * sin_lat * sin_lat;
37        let ac = A / libm::sqrt(d);
38        let a_s = w2 * ac;
39
40        let r = (ac + height) * cos_lat;
41        let x = r * cos_lon;
42        let y = r * sin_lon;
43        let z = (a_s + height) * sin_lat;
44
45        Ok(Self::new(x, y, z, epoch))
46    }
47
48    pub fn x(&self) -> f64 {
49        self.x
50    }
51
52    pub fn y(&self) -> f64 {
53        self.y
54    }
55
56    pub fn z(&self) -> f64 {
57        self.z
58    }
59
60    pub fn epoch(&self) -> TT {
61        self.epoch
62    }
63
64    pub fn position_vector(&self) -> Vector3 {
65        Vector3::new(self.x, self.y, self.z)
66    }
67
68    pub fn from_position_vector(pos: Vector3, epoch: TT) -> Self {
69        Self::new(pos.x, pos.y, pos.z, epoch)
70    }
71
72    pub fn to_geodetic(&self) -> CoordResult<(Angle, Angle, f64)> {
73        const A: f64 = 6378137.0;
74        const F: f64 = 1.0 / 298.257223563;
75        const B: f64 = A * (1.0 - F);
76        const E2: f64 = F * (2.0 - F);
77        let p = libm::sqrt(self.x * self.x + self.y * self.y);
78        let longitude = libm::atan2(self.y, self.x);
79
80        let theta = libm::atan2(self.z * A, p * B);
81        let (sin_theta, cos_theta) = libm::sincos(theta);
82        let ep2 = E2 / (1.0 - E2);
83        let mut latitude = libm::atan2(
84            self.z + ep2 * B * sin_theta.powi(3),
85            p - E2 * A * cos_theta.powi(3),
86        );
87        let mut height = 0.0;
88
89        for _ in 0..5 {
90            let (sin_lat, cos_lat) = libm::sincos(latitude);
91            let n = A / libm::sqrt(1.0 - E2 * sin_lat * sin_lat);
92            height = p / cos_lat - n;
93            latitude = libm::atan2(self.z, p * (1.0 - E2 * n / (n + height)));
94        }
95
96        Ok((
97            Angle::from_radians(longitude),
98            Angle::from_radians(latitude),
99            height,
100        ))
101    }
102
103    pub fn geocentric_distance(&self) -> f64 {
104        libm::sqrt(self.x * self.x + self.y * self.y + self.z * self.z)
105    }
106
107    pub fn distance_to(&self, other: &Self) -> f64 {
108        let dx = self.x - other.x;
109        let dy = self.y - other.y;
110        let dz = self.z - other.z;
111        libm::sqrt(dx * dx + dy * dy + dz * dz)
112    }
113
114    pub fn to_tirs(
115        &self,
116        epoch: &TT,
117        eop: &crate::eop::EopParameters,
118    ) -> CoordResult<crate::frames::TIRSPosition> {
119        use crate::frames::TIRSPosition;
120        TIRSPosition::from_itrs(self, epoch, eop)
121    }
122}
123
124impl std::fmt::Display for ITRSPosition {
125    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
126        write!(
127            f,
128            "ITRS(X={:.3}m, Y={:.3}m, Z={:.3}m, epoch=J{:.1})",
129            self.x,
130            self.y,
131            self.z,
132            self.epoch.julian_year()
133        )
134    }
135}
136
137#[cfg(test)]
138mod tests {
139    use super::*;
140
141    #[test]
142    fn test_itrs_creation() {
143        let epoch = TT::j2000();
144        let pos = ITRSPosition::new(1000000.0, 2000000.0, 3000000.0, epoch);
145
146        assert_eq!(pos.x(), 1000000.0);
147        assert_eq!(pos.y(), 2000000.0);
148        assert_eq!(pos.z(), 3000000.0);
149        assert_eq!(pos.epoch(), epoch);
150    }
151
152    #[test]
153    fn test_vector_operations() {
154        let epoch = TT::j2000();
155        let original = ITRSPosition::new(1000.0, 2000.0, 3000.0, epoch);
156
157        let vec = original.position_vector();
158        assert_eq!(vec.x, 1000.0);
159        assert_eq!(vec.y, 2000.0);
160        assert_eq!(vec.z, 3000.0);
161
162        let recovered = ITRSPosition::from_position_vector(vec, epoch);
163        assert_eq!(recovered.x(), original.x());
164        assert_eq!(recovered.y(), original.y());
165        assert_eq!(recovered.z(), original.z());
166    }
167
168    #[test]
169    fn test_geodetic_conversion_roundtrip() {
170        let epoch = TT::j2000();
171
172        // Test known location: Greenwich Observatory
173        let greenwich_lon = Angle::from_degrees(0.0);
174        let greenwich_lat = Angle::from_degrees(51.4769);
175        let greenwich_height = 47.0; // meters
176
177        let itrs =
178            ITRSPosition::from_geodetic(greenwich_lon, greenwich_lat, greenwich_height, epoch)
179                .unwrap();
180
181        let (lon, lat, height) = itrs.to_geodetic().unwrap();
182
183        // Test roundtrip accuracy (should be exact for this conversion)
184        assert_eq!(lon.degrees(), greenwich_lon.degrees());
185        assert_eq!(lat.degrees(), greenwich_lat.degrees());
186        assert_eq!(height, greenwich_height);
187    }
188
189    #[test]
190    fn test_geodetic_conversion_equator() {
191        let epoch = TT::j2000();
192
193        // Test equatorial position
194        let pos = ITRSPosition::from_geodetic(
195            Angle::from_degrees(0.0),
196            Angle::from_degrees(0.0),
197            0.0,
198            epoch,
199        )
200        .unwrap();
201
202        // Should be exactly at Earth's equatorial radius
203        const A: f64 = 6378137.0;
204        assert_eq!(pos.x(), A);
205        assert_eq!(pos.y(), 0.0);
206        assert_eq!(pos.z(), 0.0);
207    }
208
209    #[test]
210    fn test_distance_calculations() {
211        let epoch = TT::j2000();
212
213        let pos1 = ITRSPosition::new(1000.0, 0.0, 0.0, epoch);
214        let pos2 = ITRSPosition::new(2000.0, 0.0, 0.0, epoch);
215
216        assert_eq!(pos1.distance_to(&pos2), 1000.0);
217        assert_eq!(pos2.distance_to(&pos1), 1000.0);
218
219        assert_eq!(pos1.geocentric_distance(), 1000.0);
220        assert_eq!(pos2.geocentric_distance(), 2000.0);
221    }
222
223    #[test]
224    fn test_display_formatting() {
225        let epoch = TT::j2000();
226        let pos = ITRSPosition::new(1234567.89, -987654.32, 555666.77, epoch);
227
228        let display = format!("{}", pos);
229        assert!(display.contains("ITRS"));
230        assert!(display.contains("1234567.890m"));
231        assert!(display.contains("-987654.320m"));
232        assert!(display.contains("555666.770m"));
233        assert!(display.contains("J2000.0"));
234    }
235}
celestial_coords/frames/itrs.rs

celestial_coords/frames/
itrs.rs
