1use crate::{lunar, transforms::CoordinateFrame, CoordResult, Distance, ICRSPosition};
2use celestial_core::constants::HALF_PI;
3use celestial_core::matrix::RotationMatrix3;
4use celestial_core::utils::normalize_angle_to_positive;
5use celestial_core::Angle;
6use celestial_time::TT;
7
8#[cfg(feature = "serde")]
9use serde::{Deserialize, Serialize};
10
11#[derive(Debug, Clone, PartialEq)]
12#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
13pub struct SelenographicPosition {
14 latitude: Angle,
15 longitude: Angle,
16 radius: Option<Distance>,
17}
18
19impl SelenographicPosition {
20 pub fn new(latitude: Angle, longitude: Angle) -> CoordResult<Self> {
21 let latitude = latitude.validate_latitude()?;
22 let longitude = longitude.validate_longitude(true)?;
23
24 Ok(Self {
25 latitude,
26 longitude,
27 radius: None,
28 })
29 }
30
31 pub fn with_radius(latitude: Angle, longitude: Angle, radius: Distance) -> CoordResult<Self> {
32 let mut pos = Self::new(latitude, longitude)?;
33 pos.radius = Some(radius);
34 Ok(pos)
35 }
36
37 pub fn from_degrees(lat_deg: f64, lon_deg: f64) -> CoordResult<Self> {
38 Self::new(Angle::from_degrees(lat_deg), Angle::from_degrees(lon_deg))
39 }
40
41 pub fn latitude(&self) -> Angle {
42 self.latitude
43 }
44
45 pub fn longitude(&self) -> Angle {
46 self.longitude
47 }
48
49 pub fn radius(&self) -> Option<Distance> {
50 self.radius
51 }
52
53 pub fn set_radius(&mut self, radius: Distance) {
54 self.radius = Some(radius);
55 }
56
57 pub fn sub_earth_point(epoch: &TT) -> CoordResult<Self> {
58 let (lon, lat) = lunar::compute_sub_earth_point(epoch);
59 Self::new(lat, lon)
60 }
61
62 pub fn nearside_center() -> Self {
63 Self {
64 latitude: Angle::ZERO,
65 longitude: Angle::ZERO,
66 radius: None,
67 }
68 }
69
70 pub fn farside_center() -> Self {
71 Self {
72 latitude: Angle::ZERO,
73 longitude: Angle::PI,
74 radius: None,
75 }
76 }
77
78 pub fn north_pole() -> Self {
79 Self {
80 latitude: Angle::HALF_PI,
81 longitude: Angle::ZERO,
82 radius: None,
83 }
84 }
85
86 pub fn south_pole() -> Self {
87 Self {
88 latitude: -Angle::HALF_PI,
89 longitude: Angle::ZERO,
90 radius: None,
91 }
92 }
93
94 pub fn angular_separation(&self, other: &Self) -> Angle {
95 let (sin_lat1, cos_lat1) = self.latitude.sin_cos();
96 let (sin_lat2, cos_lat2) = other.latitude.sin_cos();
97 let delta_lon = (self.longitude - other.longitude).radians();
98
99 let angle_rad = celestial_core::math::vincenty_angular_separation(
100 sin_lat1, cos_lat1, sin_lat2, cos_lat2, delta_lon,
101 );
102
103 Angle::from_radians(angle_rad)
104 }
105
106 pub fn is_visible_from_earth(&self, epoch: &TT) -> bool {
107 let sub_earth = Self::sub_earth_point(epoch).unwrap_or_else(|_| Self::nearside_center());
108 let separation = self.angular_separation(&sub_earth);
109 separation.degrees() < 90.0
110 }
111}
112
113fn selenographic_to_icrs_matrix(epoch: &TT) -> CoordResult<RotationMatrix3> {
114 let orientation = lunar::compute_lunar_orientation(epoch);
115 let lib_lon = orientation.optical_libration.longitude.radians();
116 let lib_lat = orientation.optical_libration.latitude.radians();
117 let c = orientation.position_angle.radians();
118
119 let moon_icrs = lunar::get_moon_icrs(epoch)?;
120 let moon_ra = moon_icrs.ra().radians();
121 let moon_dec = moon_icrs.dec().radians();
122
123 let mut m = RotationMatrix3::identity();
124 m.rotate_z(-lib_lon);
125 m.rotate_y(-lib_lat);
126 m.rotate_z(c);
127 m.rotate_y(moon_dec - HALF_PI);
128 m.rotate_z(-moon_ra);
129 Ok(m)
130}
131
132impl CoordinateFrame for SelenographicPosition {
133 fn to_icrs(&self, epoch: &TT) -> CoordResult<ICRSPosition> {
134 let m = selenographic_to_icrs_matrix(epoch)?;
135 let (ra, dec) = m
136 .transpose()
137 .transform_spherical(self.longitude.radians(), self.latitude.radians());
138
139 let mut icrs = ICRSPosition::new(
140 Angle::from_radians(normalize_angle_to_positive(ra)),
141 Angle::from_radians(dec),
142 )?;
143
144 if let Some(radius) = self.radius {
145 icrs.set_distance(radius);
146 }
147 Ok(icrs)
148 }
149
150 fn from_icrs(icrs: &ICRSPosition, epoch: &TT) -> CoordResult<Self> {
151 let m = selenographic_to_icrs_matrix(epoch)?;
152 let (lon, lat) = m.transform_spherical(icrs.ra().radians(), icrs.dec().radians());
153
154 let mut pos = Self::new(
155 Angle::from_radians(lat),
156 Angle::from_radians(normalize_angle_to_positive(lon)),
157 )?;
158
159 if let Some(dist) = icrs.distance() {
160 pos.set_radius(dist);
161 }
162 Ok(pos)
163 }
164}
165
166impl std::fmt::Display for SelenographicPosition {
167 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
168 write!(
169 f,
170 "Selenographic(lat={:.6}°, lon={:.6}°",
171 self.latitude.degrees(),
172 self.longitude.degrees()
173 )?;
174
175 if let Some(radius) = self.radius {
176 write!(f, ", r={}", radius)?;
177 }
178
179 write!(f, ")")
180 }
181}
182
183#[cfg(test)]
184mod tests {
185 use super::*;
186
187 #[test]
188 fn test_selenographic_creation() {
189 let pos = SelenographicPosition::from_degrees(45.0, 30.0).unwrap();
190 assert!((pos.latitude().degrees() - 45.0).abs() < 1e-12);
191 assert!((pos.longitude().degrees() - 30.0).abs() < 1e-12);
192 assert!(pos.radius().is_none());
193 }
194
195 #[test]
196 fn test_selenographic_validation() {
197 assert!(SelenographicPosition::from_degrees(0.0, 0.0).is_ok());
198 assert!(SelenographicPosition::from_degrees(90.0, 180.0).is_ok());
199 assert!(SelenographicPosition::from_degrees(-90.0, 359.0).is_ok());
200
201 assert!(SelenographicPosition::from_degrees(95.0, 0.0).is_err());
202 assert!(SelenographicPosition::from_degrees(-95.0, 0.0).is_err());
203 }
204
205 #[test]
206 fn test_special_positions() {
207 let nearside = SelenographicPosition::nearside_center();
208 assert_eq!(nearside.latitude().degrees(), 0.0);
209 assert_eq!(nearside.longitude().degrees(), 0.0);
210
211 let farside = SelenographicPosition::farside_center();
212 assert_eq!(farside.latitude().degrees(), 0.0);
213 assert_eq!(farside.longitude().degrees(), 180.0);
214
215 let north_pole = SelenographicPosition::north_pole();
216 assert_eq!(north_pole.latitude().degrees(), 90.0);
217
218 let south_pole = SelenographicPosition::south_pole();
219 assert_eq!(south_pole.latitude().degrees(), -90.0);
220 }
221
222 #[test]
223 fn test_angular_separation() {
224 let nearside = SelenographicPosition::nearside_center();
225 let farside = SelenographicPosition::farside_center();
226
227 let sep = nearside.angular_separation(&farside);
228 assert!((sep.degrees() - 180.0).abs() < 1e-10);
229
230 let north = SelenographicPosition::north_pole();
231 let sep_to_north = nearside.angular_separation(&north);
232 assert!((sep_to_north.degrees() - 90.0).abs() < 1e-10);
233 }
234
235 #[test]
236 fn test_visibility_from_earth_farside() {
237 let epoch = TT::j2000();
238
239 let farside = SelenographicPosition::farside_center();
240 assert!(!farside.is_visible_from_earth(&epoch));
241 }
242
243 #[test]
244 fn test_sub_earth_point() {
245 let epoch = TT::j2000();
246 let sub_earth = SelenographicPosition::sub_earth_point(&epoch).unwrap();
247
248 assert!(
249 sub_earth.latitude().degrees().abs() <= 7.5,
250 "Sub-earth latitude = {}",
251 sub_earth.latitude().degrees()
252 );
253 assert!(
254 sub_earth.longitude().degrees() >= 0.0 && sub_earth.longitude().degrees() < 360.0,
255 "Sub-earth longitude = {}",
256 sub_earth.longitude().degrees()
257 );
258 }
259
260 #[test]
261 fn test_coordinate_frame_to_icrs() {
262 let epoch = TT::j2000();
263 let original = SelenographicPosition::from_degrees(0.0, 0.0).unwrap();
264
265 let icrs = original.to_icrs(&epoch).unwrap();
266
267 assert!(icrs.ra().degrees() >= 0.0 && icrs.ra().degrees() < 360.0);
268 assert!(icrs.dec().degrees() >= -90.0 && icrs.dec().degrees() <= 90.0);
269 }
270
271 #[test]
272 fn test_coordinate_frame_roundtrip() {
273 let epoch = TT::j2000();
274 let test_cases = [
275 (0.0, 0.0),
276 (5.0, 30.0),
277 (-5.0, 90.0),
278 (3.0, 180.0),
279 (-3.0, 270.0),
280 ];
281
282 for (lat, lon) in test_cases {
283 let original = SelenographicPosition::from_degrees(lat, lon).unwrap();
284 let icrs = original.to_icrs(&epoch).unwrap();
285 let recovered = SelenographicPosition::from_icrs(&icrs, &epoch).unwrap();
286
287 let lat_err = (original.latitude().degrees() - recovered.latitude().degrees()).abs();
288 let lon_diff = (original.longitude().radians() - recovered.longitude().radians()).abs();
289 let lon_err = if lon_diff > std::f64::consts::PI {
290 std::f64::consts::TAU - lon_diff
291 } else {
292 lon_diff
293 } * celestial_core::constants::RAD_TO_DEG;
294
295 assert!(
296 lat_err < 1.0 / 3600.0,
297 "({}, {}): Latitude error {:.6} arcsec",
298 lat,
299 lon,
300 lat_err * 3600.0,
301 );
302 assert!(
303 lon_err < 1.0 / 3600.0,
304 "({}, {}): Longitude error {:.6} arcsec",
305 lat,
306 lon,
307 lon_err * 3600.0,
308 );
309 }
310 }
311
312 #[test]
313 fn test_with_radius() {
314 let radius = Distance::from_au(0.00257).unwrap();
315 let pos = SelenographicPosition::with_radius(
316 Angle::from_degrees(0.0),
317 Angle::from_degrees(0.0),
318 radius,
319 )
320 .unwrap();
321
322 assert!(pos.radius().is_some());
323 assert_eq!(pos.radius().unwrap(), radius);
324 }
325
326 #[test]
327 fn test_display_formatting() {
328 let pos = SelenographicPosition::from_degrees(45.123456, 30.654321).unwrap();
329 let display = format!("{}", pos);
330 assert!(display.contains("45.123456"));
331 assert!(display.contains("30.654321"));
332 assert!(display.contains("Selenographic"));
333 }
334}