celestial_coords/frames/
galactic.rs

1use crate::{
2    constants::GALACTIC_TO_ICRS, transforms::CoordinateFrame, CoordResult, Distance, ICRSPosition,
3};
4use celestial_core::Angle;
5use celestial_time::TT;
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 GalacticPosition {
13    l: Angle,
14    b: Angle,
15    distance: Option<Distance>,
16}
17
18impl GalacticPosition {
19    pub fn new(l: Angle, b: Angle) -> CoordResult<Self> {
20        let l = l.validate_longitude(true)?;
21        let b = b.validate_latitude()?;
22
23        Ok(Self {
24            l,
25            b,
26            distance: None,
27        })
28    }
29
30    pub fn with_distance(l: Angle, b: Angle, distance: Distance) -> CoordResult<Self> {
31        let mut pos = Self::new(l, b)?;
32        pos.distance = Some(distance);
33        Ok(pos)
34    }
35
36    pub fn from_degrees(l_deg: f64, b_deg: f64) -> CoordResult<Self> {
37        Self::new(Angle::from_degrees(l_deg), Angle::from_degrees(b_deg))
38    }
39
40    pub fn longitude(&self) -> Angle {
41        self.l
42    }
43
44    pub fn latitude(&self) -> Angle {
45        self.b
46    }
47
48    pub fn distance(&self) -> Option<Distance> {
49        self.distance
50    }
51
52    pub fn set_distance(&mut self, distance: Distance) {
53        self.distance = Some(distance);
54    }
55
56    pub fn galactic_center() -> Self {
57        Self {
58            l: Angle::ZERO,
59            b: Angle::ZERO,
60            distance: None,
61        }
62    }
63
64    pub fn galactic_anticenter() -> Self {
65        Self {
66            l: Angle::PI,
67            b: Angle::ZERO,
68            distance: None,
69        }
70    }
71
72    pub fn north_galactic_pole() -> Self {
73        Self {
74            l: Angle::ZERO,
75            b: Angle::HALF_PI,
76            distance: None,
77        }
78    }
79
80    pub fn south_galactic_pole() -> Self {
81        Self {
82            l: Angle::ZERO,
83            b: -Angle::HALF_PI,
84            distance: None,
85        }
86    }
87
88    pub fn is_near_galactic_plane(&self) -> bool {
89        self.b.abs().degrees() < 10.0
90    }
91
92    pub fn is_in_galactic_bulge(&self) -> bool {
93        self.b.abs().degrees() < 10.0 && (self.l.degrees() < 30.0 || self.l.degrees() > 330.0)
94    }
95
96    pub fn is_near_galactic_pole(&self) -> bool {
97        self.b.abs().degrees() > 80.0
98    }
99
100    pub fn angular_distance_from_gc(&self) -> Angle {
101        let gc = Self::galactic_center();
102        self.angular_separation(&gc)
103    }
104
105    pub fn angular_separation(&self, other: &Self) -> Angle {
106        let (sin_b1, cos_b1) = self.b.sin_cos();
107        let (sin_b2, cos_b2) = other.b.sin_cos();
108        let delta_l = (self.l - other.l).radians();
109
110        let angle_rad = celestial_core::math::vincenty_angular_separation(
111            sin_b1, cos_b1, sin_b2, cos_b2, delta_l,
112        );
113
114        Angle::from_radians(angle_rad)
115    }
116}
117
118impl CoordinateFrame for GalacticPosition {
119    fn to_icrs(&self, _epoch: &TT) -> CoordResult<ICRSPosition> {
120        let (sin_b, cos_b) = self.b.sin_cos();
121        let (sin_l, cos_l) = self.l.sin_cos();
122        let gal_cartesian = [cos_l * cos_b, sin_l * cos_b, sin_b];
123
124        // Matrix multiplication: icrs = M^T * gal (transpose because matrix is stored as columns)
125        // This works correctly because GALACTIC_TO_ICRS is orthonormal (rotation matrix).
126        let icrs_cartesian = [
127            GALACTIC_TO_ICRS[0][0] * gal_cartesian[0]
128                + GALACTIC_TO_ICRS[1][0] * gal_cartesian[1]
129                + GALACTIC_TO_ICRS[2][0] * gal_cartesian[2],
130            GALACTIC_TO_ICRS[0][1] * gal_cartesian[0]
131                + GALACTIC_TO_ICRS[1][1] * gal_cartesian[1]
132                + GALACTIC_TO_ICRS[2][1] * gal_cartesian[2],
133            GALACTIC_TO_ICRS[0][2] * gal_cartesian[0]
134                + GALACTIC_TO_ICRS[1][2] * gal_cartesian[1]
135                + GALACTIC_TO_ICRS[2][2] * gal_cartesian[2],
136        ];
137
138        let d2 = icrs_cartesian[0] * icrs_cartesian[0] + icrs_cartesian[1] * icrs_cartesian[1];
139        let ra = if d2 != 0.0 {
140            libm::atan2(icrs_cartesian[1], icrs_cartesian[0])
141        } else {
142            0.0
143        };
144        let dec = if d2 != 0.0 || icrs_cartesian[2] != 0.0 {
145            libm::atan2(icrs_cartesian[2], libm::sqrt(d2))
146        } else {
147            0.0
148        };
149
150        let mut icrs = ICRSPosition::new(Angle::from_radians(ra), Angle::from_radians(dec))?;
151
152        if let Some(distance) = self.distance {
153            icrs.set_distance(distance);
154        }
155
156        Ok(icrs)
157    }
158
159    fn from_icrs(icrs: &ICRSPosition, _epoch: &TT) -> CoordResult<Self> {
160        let (sin_dec, cos_dec) = icrs.dec().sin_cos();
161        let (sin_ra, cos_ra) = icrs.ra().sin_cos();
162        let icrs_cartesian = [cos_ra * cos_dec, sin_ra * cos_dec, sin_dec];
163
164        // Matrix multiplication: gal = M * icrs (standard row-major access)
165        // For orthonormal matrices, M^T = M^(-1), so this is the inverse of to_icrs.
166        let gal_cartesian = [
167            GALACTIC_TO_ICRS[0][0] * icrs_cartesian[0]
168                + GALACTIC_TO_ICRS[0][1] * icrs_cartesian[1]
169                + GALACTIC_TO_ICRS[0][2] * icrs_cartesian[2],
170            GALACTIC_TO_ICRS[1][0] * icrs_cartesian[0]
171                + GALACTIC_TO_ICRS[1][1] * icrs_cartesian[1]
172                + GALACTIC_TO_ICRS[1][2] * icrs_cartesian[2],
173            GALACTIC_TO_ICRS[2][0] * icrs_cartesian[0]
174                + GALACTIC_TO_ICRS[2][1] * icrs_cartesian[1]
175                + GALACTIC_TO_ICRS[2][2] * icrs_cartesian[2],
176        ];
177
178        let d2 = gal_cartesian[0] * gal_cartesian[0] + gal_cartesian[1] * gal_cartesian[1];
179        let l = if d2 != 0.0 {
180            libm::atan2(gal_cartesian[1], gal_cartesian[0])
181        } else {
182            0.0
183        };
184        let b = if d2 != 0.0 || gal_cartesian[2] != 0.0 {
185            libm::atan2(gal_cartesian[2], libm::sqrt(d2))
186        } else {
187            0.0
188        };
189
190        let mut galactic = Self::new(Angle::from_radians(l), Angle::from_radians(b))?;
191
192        if let Some(distance) = icrs.distance() {
193            galactic.set_distance(distance);
194        }
195
196        Ok(galactic)
197    }
198}
199
200impl std::fmt::Display for GalacticPosition {
201    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
202        write!(
203            f,
204            "Galactic(l={:.6}°, b={:.6}°",
205            self.l.degrees(),
206            self.b.degrees()
207        )?;
208
209        if let Some(distance) = self.distance {
210            write!(f, ", d={}", distance)?;
211        }
212
213        write!(f, ")")
214    }
215}
216
217#[cfg(test)]
218mod tests {
219    use super::*;
220
221    #[test]
222    fn test_galactic_creation() {
223        let pos = GalacticPosition::from_degrees(45.0, 30.0).unwrap();
224        assert!((pos.longitude().degrees() - 45.0).abs() < 1e-12);
225        assert!((pos.latitude().degrees() - 30.0).abs() < 1e-12);
226        assert!(pos.distance().is_none());
227    }
228
229    #[test]
230    fn test_galactic_validation() {
231        // Valid coordinates
232        assert!(GalacticPosition::from_degrees(0.0, 0.0).is_ok());
233        assert!(GalacticPosition::from_degrees(359.999, 89.999).is_ok());
234
235        // Longitude gets normalized
236        let pos = GalacticPosition::from_degrees(380.0, 45.0).unwrap();
237        assert!((pos.longitude().degrees() - 20.0).abs() < 1e-12);
238
239        // Invalid latitude
240        assert!(GalacticPosition::from_degrees(0.0, 95.0).is_err());
241        assert!(GalacticPosition::from_degrees(0.0, -95.0).is_err());
242    }
243
244    #[test]
245    fn test_special_positions() {
246        let gc = GalacticPosition::galactic_center();
247        assert_eq!(gc.longitude().degrees(), 0.0);
248        assert_eq!(gc.latitude().degrees(), 0.0);
249
250        let gac = GalacticPosition::galactic_anticenter();
251        assert!((gac.longitude().degrees() - 180.0).abs() < 1e-12);
252        assert_eq!(gac.latitude().degrees(), 0.0);
253
254        let ngp = GalacticPosition::north_galactic_pole();
255        assert!((ngp.latitude().degrees() - 90.0).abs() < 1e-12);
256
257        let sgp = GalacticPosition::south_galactic_pole();
258        assert!((sgp.latitude().degrees() - (-90.0)).abs() < 1e-12);
259    }
260
261    #[test]
262    fn test_galactic_regions() {
263        // Galactic plane
264        let plane_pos = GalacticPosition::from_degrees(45.0, 5.0).unwrap();
265        assert!(plane_pos.is_near_galactic_plane());
266        assert!(!plane_pos.is_near_galactic_pole());
267
268        // Galactic bulge
269        let bulge_pos = GalacticPosition::from_degrees(5.0, 5.0).unwrap();
270        assert!(bulge_pos.is_in_galactic_bulge());
271
272        // Galactic pole
273        let pole_pos = GalacticPosition::from_degrees(0.0, 85.0).unwrap();
274        assert!(pole_pos.is_near_galactic_pole());
275        assert!(!pole_pos.is_near_galactic_plane());
276    }
277
278    #[test]
279    fn test_angular_separation() {
280        let pos1 = GalacticPosition::from_degrees(0.0, 0.0).unwrap();
281        let pos2 = GalacticPosition::from_degrees(90.0, 0.0).unwrap();
282
283        let sep = pos1.angular_separation(&pos2);
284        // Should be approximately 90 degrees
285        assert!((sep.degrees() - 90.0).abs() < 1.0); // Allow 1° tolerance for approximation
286
287        // Distance from galactic center
288        let gc_dist = pos2.angular_distance_from_gc();
289        assert!((gc_dist.degrees() - 90.0).abs() < 1.0);
290    }
291
292    #[test]
293    fn test_coordinate_transformations() {
294        let epoch = TT::j2000();
295        let gal_pos = GalacticPosition::from_degrees(45.0, 30.0).unwrap();
296
297        // Test Galactic -> ICRS -> Galactic round trip
298        let icrs = gal_pos.to_icrs(&epoch).unwrap();
299        let gal_recovered = GalacticPosition::from_icrs(&icrs, &epoch).unwrap();
300
301        assert!(
302            (gal_recovered.longitude().degrees() - gal_pos.longitude().degrees()).abs() < 1e-12
303        );
304        assert!((gal_recovered.latitude().degrees() - gal_pos.latitude().degrees()).abs() < 1e-12);
305    }
306
307    #[test]
308    fn test_with_distance() {
309        let distance = Distance::from_parsecs(100.0).unwrap();
310        let pos = GalacticPosition::with_distance(
311            Angle::from_degrees(45.0),
312            Angle::from_degrees(30.0),
313            distance,
314        )
315        .unwrap();
316
317        assert_eq!(pos.distance().unwrap().parsecs(), 100.0);
318
319        // Test coordinate transformation preserves distance
320        let epoch = TT::j2000();
321        let icrs = pos.to_icrs(&epoch).unwrap();
322        assert_eq!(icrs.distance().unwrap().parsecs(), 100.0);
323    }
324}
celestial_coords/frames/galactic.rs

celestial_coords/frames/
galactic.rs
