1use crate::{solar, 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 HeliographicStonyhurst {
14 latitude: Angle,
15 longitude: Angle,
16 radius: Option<Distance>,
17}
18
19impl HeliographicStonyhurst {
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 to_carrington(&self, epoch: &TT) -> CoordResult<HeliographicCarrington> {
58 let l0 = solar::compute_l0(epoch);
59 let carrington_lon = self.longitude + l0;
60 let normalized_lon =
61 Angle::from_radians(normalize_angle_to_positive(carrington_lon.radians()));
62
63 let mut carr = HeliographicCarrington::new(self.latitude, normalized_lon)?;
64 if let Some(r) = self.radius {
65 carr.set_radius(r);
66 }
67 Ok(carr)
68 }
69
70 pub fn disk_center(epoch: &TT) -> Self {
71 let orientation = solar::compute_solar_orientation(epoch);
72 Self {
73 latitude: orientation.b0,
74 longitude: Angle::ZERO,
75 radius: None,
76 }
77 }
78}
79
80#[derive(Debug, Clone, PartialEq)]
81#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
82pub struct HeliographicCarrington {
83 latitude: Angle,
84 longitude: Angle,
85 radius: Option<Distance>,
86}
87
88impl HeliographicCarrington {
89 pub fn new(latitude: Angle, longitude: Angle) -> CoordResult<Self> {
90 let latitude = latitude.validate_latitude()?;
91 let longitude = longitude.validate_longitude(true)?;
92
93 Ok(Self {
94 latitude,
95 longitude,
96 radius: None,
97 })
98 }
99
100 pub fn with_radius(latitude: Angle, longitude: Angle, radius: Distance) -> CoordResult<Self> {
101 let mut pos = Self::new(latitude, longitude)?;
102 pos.radius = Some(radius);
103 Ok(pos)
104 }
105
106 pub fn from_degrees(lat_deg: f64, lon_deg: f64) -> CoordResult<Self> {
107 Self::new(Angle::from_degrees(lat_deg), Angle::from_degrees(lon_deg))
108 }
109
110 pub fn latitude(&self) -> Angle {
111 self.latitude
112 }
113
114 pub fn longitude(&self) -> Angle {
115 self.longitude
116 }
117
118 pub fn radius(&self) -> Option<Distance> {
119 self.radius
120 }
121
122 pub fn set_radius(&mut self, radius: Distance) {
123 self.radius = Some(radius);
124 }
125
126 pub fn to_stonyhurst(&self, epoch: &TT) -> CoordResult<HeliographicStonyhurst> {
127 let l0 = solar::compute_l0(epoch);
128 let stonyhurst_lon = self.longitude - l0;
129 let normalized_lon =
130 Angle::from_radians(normalize_angle_to_positive(stonyhurst_lon.radians()));
131
132 let mut stony = HeliographicStonyhurst::new(self.latitude, normalized_lon)?;
133 if let Some(r) = self.radius {
134 stony.set_radius(r);
135 }
136 Ok(stony)
137 }
138
139 pub fn carrington_rotation_number(epoch: &TT) -> f64 {
140 const CARRINGTON_EPOCH_JD: f64 = 2398220.0;
141 const CARRINGTON_PERIOD_DAYS: f64 = 25.38;
142
143 let jd = epoch.to_julian_date();
144 let d = jd.jd1() + jd.jd2() - CARRINGTON_EPOCH_JD;
145 d / CARRINGTON_PERIOD_DAYS
146 }
147}
148
149fn heliographic_to_icrs_matrix(epoch: &TT) -> CoordResult<RotationMatrix3> {
150 let orientation = solar::compute_solar_orientation(epoch);
151 let b0 = orientation.b0.radians();
152 let p = orientation.p.radians();
153
154 let sun_icrs = solar::get_sun_icrs(epoch)?;
155 let sun_ra = sun_icrs.ra().radians();
156 let sun_dec = sun_icrs.dec().radians();
157
158 let mut m = RotationMatrix3::identity();
159 m.rotate_y(-b0);
160 m.rotate_z(p);
161 m.rotate_y(sun_dec - HALF_PI);
162 m.rotate_z(-sun_ra);
163 Ok(m)
164}
165
166impl CoordinateFrame for HeliographicStonyhurst {
167 fn to_icrs(&self, epoch: &TT) -> CoordResult<ICRSPosition> {
168 let m = heliographic_to_icrs_matrix(epoch)?;
169 let (ra, dec) = m
170 .transpose()
171 .transform_spherical(self.longitude.radians(), self.latitude.radians());
172
173 let mut icrs = ICRSPosition::new(
174 Angle::from_radians(normalize_angle_to_positive(ra)),
175 Angle::from_radians(dec),
176 )?;
177
178 if let Some(radius) = self.radius {
179 icrs.set_distance(radius);
180 }
181 Ok(icrs)
182 }
183
184 fn from_icrs(icrs: &ICRSPosition, epoch: &TT) -> CoordResult<Self> {
185 let m = heliographic_to_icrs_matrix(epoch)?;
186 let (lon, lat) = m.transform_spherical(icrs.ra().radians(), icrs.dec().radians());
187
188 let mut pos = Self::new(
189 Angle::from_radians(lat),
190 Angle::from_radians(normalize_angle_to_positive(lon)),
191 )?;
192
193 if let Some(dist) = icrs.distance() {
194 pos.set_radius(dist);
195 }
196 Ok(pos)
197 }
198}
199
200impl CoordinateFrame for HeliographicCarrington {
201 fn to_icrs(&self, epoch: &TT) -> CoordResult<ICRSPosition> {
202 let stonyhurst = self.to_stonyhurst(epoch)?;
203 stonyhurst.to_icrs(epoch)
204 }
205
206 fn from_icrs(icrs: &ICRSPosition, epoch: &TT) -> CoordResult<Self> {
207 let stonyhurst = HeliographicStonyhurst::from_icrs(icrs, epoch)?;
208 stonyhurst.to_carrington(epoch)
209 }
210}
211
212impl std::fmt::Display for HeliographicStonyhurst {
213 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
214 write!(
215 f,
216 "HeliographicStonyhurst(lat={:.6}°, lon={:.6}°",
217 self.latitude.degrees(),
218 self.longitude.degrees()
219 )?;
220
221 if let Some(radius) = self.radius {
222 write!(f, ", r={}", radius)?;
223 }
224
225 write!(f, ")")
226 }
227}
228
229impl std::fmt::Display for HeliographicCarrington {
230 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
231 write!(
232 f,
233 "HeliographicCarrington(lat={:.6}°, lon={:.6}°",
234 self.latitude.degrees(),
235 self.longitude.degrees()
236 )?;
237
238 if let Some(radius) = self.radius {
239 write!(f, ", r={}", radius)?;
240 }
241
242 write!(f, ")")
243 }
244}
245
246#[cfg(test)]
247mod tests {
248 use super::*;
249
250 #[test]
251 fn test_stonyhurst_creation() {
252 let pos = HeliographicStonyhurst::from_degrees(45.0, 30.0).unwrap();
253 assert!((pos.latitude().degrees() - 45.0).abs() < 1e-12);
254 assert!((pos.longitude().degrees() - 30.0).abs() < 1e-12);
255 assert!(pos.radius().is_none());
256 }
257
258 #[test]
259 fn test_carrington_creation() {
260 let pos = HeliographicCarrington::from_degrees(-30.0, 180.0).unwrap();
261 assert!((pos.latitude().degrees() - (-30.0)).abs() < 1e-12);
262 assert!((pos.longitude().degrees() - 180.0).abs() < 1e-12);
263 assert!(pos.radius().is_none());
264 }
265
266 #[test]
267 fn test_stonyhurst_validation() {
268 assert!(HeliographicStonyhurst::from_degrees(0.0, 0.0).is_ok());
269 assert!(HeliographicStonyhurst::from_degrees(90.0, 180.0).is_ok());
270 assert!(HeliographicStonyhurst::from_degrees(-90.0, 359.0).is_ok());
271
272 assert!(HeliographicStonyhurst::from_degrees(95.0, 0.0).is_err());
273 assert!(HeliographicStonyhurst::from_degrees(-95.0, 0.0).is_err());
274 }
275
276 #[test]
277 fn test_stonyhurst_to_carrington_differs_by_l0() {
278 let epoch = TT::j2000();
279 let stonyhurst = HeliographicStonyhurst::from_degrees(15.0, 45.0).unwrap();
280 let carrington = stonyhurst.to_carrington(&epoch).unwrap();
281
282 assert_eq!(
283 stonyhurst.latitude().degrees(),
284 carrington.latitude().degrees()
285 );
286
287 let l0 = solar::compute_l0(&epoch);
288 let expected_carr_lon =
289 normalize_angle_to_positive((stonyhurst.longitude() + l0).radians())
290 * celestial_core::constants::RAD_TO_DEG;
291
292 assert!((carrington.longitude().degrees() - expected_carr_lon).abs() < 1e-10);
293 }
294
295 #[test]
296 fn test_carrington_to_stonyhurst_roundtrip() {
297 let epoch = TT::j2000();
298 let original = HeliographicCarrington::from_degrees(30.0, 120.0).unwrap();
299 let stonyhurst = original.to_stonyhurst(&epoch).unwrap();
300 let roundtrip = stonyhurst.to_carrington(&epoch).unwrap();
301
302 assert!((original.latitude().degrees() - roundtrip.latitude().degrees()).abs() < 1e-10);
303 assert!((original.longitude().degrees() - roundtrip.longitude().degrees()).abs() < 1e-10);
304 }
305
306 #[test]
307 fn test_disk_center() {
308 let epoch = TT::j2000();
309 let center = HeliographicStonyhurst::disk_center(&epoch);
310
311 let b0 = solar::compute_b0(&epoch);
312 assert!((center.latitude().degrees() - b0.degrees()).abs() < 1e-12);
313 assert_eq!(center.longitude().degrees(), 0.0);
314 }
315
316 #[test]
317 fn test_carrington_rotation_number() {
318 let epoch = TT::j2000();
319 let rotation = HeliographicCarrington::carrington_rotation_number(&epoch);
320
321 assert!(
322 rotation > 1900.0 && rotation < 2200.0,
323 "Carrington rotation number at J2000 = {} should be reasonable",
324 rotation
325 );
326 }
327
328 #[test]
329 fn test_coordinate_frame_roundtrip() {
330 let epoch = TT::j2000();
331 let test_cases = [
332 (20.0, 30.0),
333 (0.0, 0.0),
334 (45.0, 90.0),
335 (-7.0, 180.0),
336 (7.0, 270.0),
337 ];
338
339 for (lat, lon) in test_cases {
340 let original = HeliographicStonyhurst::from_degrees(lat, lon).unwrap();
341 let icrs = original.to_icrs(&epoch).unwrap();
342 let recovered = HeliographicStonyhurst::from_icrs(&icrs, &epoch).unwrap();
343
344 let lat_err = (original.latitude().degrees() - recovered.latitude().degrees()).abs();
345 let lon_diff = (original.longitude().radians() - recovered.longitude().radians()).abs();
346 let lon_err = if lon_diff > std::f64::consts::PI {
347 std::f64::consts::TAU - lon_diff
348 } else {
349 lon_diff
350 } * celestial_core::constants::RAD_TO_DEG;
351
352 assert!(
353 lat_err < 1.0 / 3600.0,
354 "({}, {}): Latitude error {:.6} arcsec",
355 lat,
356 lon,
357 lat_err * 3600.0,
358 );
359 assert!(
360 lon_err < 1.0 / 3600.0,
361 "({}, {}): Longitude error {:.6} arcsec",
362 lat,
363 lon,
364 lon_err * 3600.0,
365 );
366 }
367 }
368
369 #[test]
370 fn test_with_radius() {
371 let radius = Distance::from_au(0.00465047).unwrap();
372 let pos = HeliographicStonyhurst::with_radius(
373 Angle::from_degrees(0.0),
374 Angle::from_degrees(0.0),
375 radius,
376 )
377 .unwrap();
378
379 assert!(pos.radius().is_some());
380 assert_eq!(pos.radius().unwrap(), radius);
381 }
382
383 #[test]
384 fn test_display_formatting() {
385 let pos = HeliographicStonyhurst::from_degrees(45.123456, 30.654321).unwrap();
386 let display = format!("{}", pos);
387 assert!(display.contains("45.123456"));
388 assert!(display.contains("30.654321"));
389 assert!(display.contains("HeliographicStonyhurst"));
390 }
391}