celestial_core/location/
geodesy.rs

1//! Geodetic to geocentric coordinate conversions for Earth-based observers.
2//!
3//! # Geodetic vs Geocentric Coordinates
4//!
5//! **Geodetic coordinates** (what GPS gives you) define position relative to the WGS84
6//! reference ellipsoid:
7//! - Latitude: angle between the equatorial plane and the ellipsoid surface normal
8//! - Longitude: angle from the prime meridian
9//! - Height: distance above the ellipsoid surface
10//!
11//! **Geocentric coordinates** define position relative to Earth's center of mass:
12//! - The Earth is modeled as an oblate spheroid (equatorial bulge)
13//! - At mid-latitudes, geodetic and geocentric latitude differ by up to ~11 arcminutes
14//!
15//! Topocentric corrections (parallax, aberration, refraction) require knowing the
16//! observer's true position in space, not their position on the reference ellipsoid.
17//! The geocentric coordinates returned here are the cylindrical components needed for:
18//!
19//! - **Diurnal parallax**: Moon position shifts by up to 1° depending on observer
20//! - **Stellar parallax**: precise baseline for nearby star distances
21//! - **Satellite tracking**: ground station positions in Earth-centered frame
22//!
23//! # WGS84 Ellipsoid Parameters
24//!
25//! This module uses the WGS84 reference ellipsoid:
26//! - Semi-major axis (equatorial radius): 6,378,137.0 m (or 6378.137 km)
27//! - Flattening: 1/298.257223563
28//! - First eccentricity squared: ~0.00669438
29//!
30//! # Output Format
31//!
32//! Both conversion methods return `(u, v)` where:
33//! - `u`: distance from Earth's rotation axis (equatorial component)
34//! - `v`: distance from equatorial plane (polar component)
35//!
36//! These are cylindrical coordinates centered on Earth's center of mass.
37//! To get Cartesian XYZ, you'd combine with longitude: `x = u*cos(lon)`, `y = u*sin(lon)`, `z = v`.
38
39use crate::constants::{
40    WGS84_ECCENTRICITY_SQUARED, WGS84_SEMI_MAJOR_AXIS, WGS84_SEMI_MAJOR_AXIS_KM,
41};
42use crate::errors::{AstroError, AstroResult, MathErrorKind};
43
44use super::Location;
45
46impl Location {
47    /// Converts geodetic coordinates to geocentric cylindrical coordinates in kilometers.
48    ///
49    /// Uses the WGS84 ellipsoid to compute the observer's position relative to
50    /// Earth's center of mass. The result accounts for Earth's equatorial bulge.
51    ///
52    /// # Returns
53    ///
54    /// `(u, v)` in kilometers where:
55    /// - `u`: perpendicular distance from Earth's rotation axis
56    /// - `v`: distance from the equatorial plane (positive north)
57    ///
58    /// # Example
59    ///
60    /// ```
61    /// use celestial_core::Location;
62    ///
63    /// let obs = Location::from_degrees(45.0, 0.0, 0.0)?;
64    /// let (u, v) = obs.to_geocentric_km()?;
65    ///
66    /// // At 45 degrees, u and v are similar but u > v due to Earth's shape
67    /// assert!(u > 4500.0 && u < 4600.0);
68    /// assert!(v > 4400.0 && v < 4500.0);
69    /// # Ok::<(), celestial_core::AstroError>(())
70    /// ```
71    ///
72    /// # Errors
73    ///
74    /// Returns an error for degenerate latitude values that would cause division by zero.
75    /// In practice, this is unlikely with validated `Location` values.
76    pub fn to_geocentric_km(&self) -> AstroResult<(f64, f64)> {
77        let lat = self.latitude;
78        let height_km = self.height / 1000.0;
79
80        let (sin_lat, cos_lat) = libm::sincos(lat);
81
82        let denominator = 1.0 - WGS84_ECCENTRICITY_SQUARED * sin_lat * sin_lat;
83        if denominator <= f64::EPSILON {
84            return Err(AstroError::math_error(
85                "geocentric_conversion",
86                MathErrorKind::DivisionByZero,
87                "Latitude too close to critical value causing division by zero",
88            ));
89        }
90
91        let n = WGS84_SEMI_MAJOR_AXIS_KM / libm::sqrt(denominator);
92
93        let u = (n + height_km) * cos_lat;
94
95        let v = (n * (1.0 - WGS84_ECCENTRICITY_SQUARED) + height_km) * sin_lat;
96
97        Ok((u, v))
98    }
99
100    /// Converts geodetic coordinates to geocentric cylindrical coordinates in meters.
101    ///
102    /// Same algorithm as [`to_geocentric_km`](Self::to_geocentric_km) but returns
103    /// results in meters for applications requiring higher precision or SI units.
104    ///
105    /// This method uses a slightly different formulation internally (computing the
106    /// flattening ratio explicitly) but produces equivalent results to the km version.
107    ///
108    /// # Returns
109    ///
110    /// `(u, v)` in meters where:
111    /// - `u`: perpendicular distance from Earth's rotation axis
112    /// - `v`: distance from the equatorial plane (positive north)
113    ///
114    /// # Example
115    ///
116    /// ```
117    /// use celestial_core::Location;
118    ///
119    /// // Equator at sea level
120    /// let equator = Location::from_degrees(0.0, 0.0, 0.0)?;
121    /// let (u, v) = equator.to_geocentric_meters()?;
122    ///
123    /// // At equator: u equals semi-major axis, v is zero
124    /// assert!((u - 6_378_137.0).abs() < 1.0);
125    /// assert!(v.abs() < 1e-9);
126    /// # Ok::<(), celestial_core::AstroError>(())
127    /// ```
128    pub fn to_geocentric_meters(&self) -> AstroResult<(f64, f64)> {
129        let height_m = self.height;
130
131        let (phi_sin, phi_cos) = libm::sincos(self.latitude);
132
133        let wgs_flattened = 1.0 / 298.257223563;
134        let axis_ratio = 1.0 - wgs_flattened;
135        let axis_ratio_sq = axis_ratio * axis_ratio;
136
137        let norm_sq = phi_cos * phi_cos + axis_ratio_sq * phi_sin * phi_sin;
138        if norm_sq <= f64::EPSILON {
139            return Err(AstroError::math_error(
140                "geocentric_conversion",
141                MathErrorKind::DivisionByZero,
142                "Latitude too close to critical value causing division by zero",
143            ));
144        }
145
146        let prime_vertical_radius = WGS84_SEMI_MAJOR_AXIS / libm::sqrt(norm_sq);
147        let as_val = axis_ratio_sq * prime_vertical_radius;
148
149        let equatorial_radius = (prime_vertical_radius + height_m) * phi_cos;
150        let z_coordinate = (as_val + height_m) * phi_sin;
151
152        Ok((equatorial_radius, z_coordinate))
153    }
154}
155
156#[cfg(test)]
157mod tests {
158    use super::*;
159
160    #[test]
161    fn test_geocentric_at_equator() {
162        let loc = Location::from_degrees(0.0, 0.0, 0.0).unwrap();
163        let (u, v) = loc.to_geocentric_km().unwrap();
164
165        assert_eq!(
166            u,
167            crate::constants::WGS84_SEMI_MAJOR_AXIS_KM,
168            "u = {} km, expected ~6378.137 km",
169            u
170        );
171        assert_eq!(v, 0.0, "v = {} km, expected ~0 km", v);
172    }
173
174    #[test]
175    fn test_geocentric_at_north_pole() {
176        let loc = Location::from_degrees(90.0, 0.0, 0.0).unwrap();
177        let (u, v) = loc.to_geocentric_km().unwrap();
178
179        assert!(u.abs() < 1e-10, "u = {} km, expected very close to 0 km", u);
180        let expected_polar_radius =
181            crate::constants::WGS84_SEMI_MAJOR_AXIS_KM * (1.0 - 1.0 / 298.257223563);
182        assert_eq!(
183            v, expected_polar_radius,
184            "v = {} km, expected ~{} km",
185            v, expected_polar_radius
186        );
187    }
188
189    #[test]
190    fn test_geocentric_km_rejects_degenerate_latitude() {
191        use crate::constants::{HALF_PI, PI};
192        let critical_lat = HALF_PI - 1e-10;
193        let critical_lat_deg = critical_lat * 180.0 / PI;
194
195        let mut test_lat = critical_lat_deg;
196        while test_lat < 90.0 {
197            let loc = Location::from_degrees(test_lat, 0.0, 0.0).unwrap();
198            if loc.to_geocentric_km().is_err() {
199                return;
200            }
201            test_lat += 1e-12;
202        }
203    }
204
205    #[test]
206    fn test_geocentric_meters_rejects_degenerate_latitude() {
207        use crate::constants::{HALF_PI, PI};
208        let critical_lat = HALF_PI - 1e-10;
209        let critical_lat_deg = critical_lat * 180.0 / PI;
210
211        let mut test_lat = critical_lat_deg;
212        while test_lat < 90.0 {
213            let loc = Location::from_degrees(test_lat, 0.0, 0.0).unwrap();
214            if loc.to_geocentric_meters().is_err() {
215                return;
216            }
217            test_lat += 1e-12;
218        }
219    }
220
221    #[test]
222    fn test_geocentric_meters_handles_equator() {
223        let loc = Location::from_degrees(0.0, 0.0, 0.0).unwrap();
224        let (u, v) = loc.to_geocentric_meters().unwrap();
225        crate::test_helpers::assert_float_eq(u, WGS84_SEMI_MAJOR_AXIS, 1);
226        crate::test_helpers::assert_float_eq(v, 0.0, 1);
227    }
228
229    #[test]
230    fn test_geocentric_meters_handles_north_pole() {
231        let loc = Location::from_degrees(90.0, 0.0, 0.0).unwrap();
232        let (u, v) = loc.to_geocentric_meters().unwrap();
233        assert!(u.abs() < 1e-9);
234        crate::test_helpers::assert_ulp_le(v, 6356752.314245179, 1, "v at pole");
235    }
236
237    #[test]
238    fn test_geocentric_at_45_degrees() {
239        let loc = Location::from_degrees(45.0, 0.0, 0.0).unwrap();
240        let (u, v) = loc.to_geocentric_km().unwrap();
241
242        assert!(u > 4000.0 && u < 5000.0, "u = {} km, expected ~4500 km", u);
243        assert!(v > 4000.0 && v < 5000.0, "v = {} km, expected ~4500 km", v);
244
245        assert!(
246            u > v,
247            "u should be larger than v due to Earth's oblate shape: u={}, v={}",
248            u,
249            v
250        );
251        assert!(
252            (u - v).abs() < 100.0,
253            "At 45°, u and v should be similar: u={}, v={}",
254            u,
255            v
256        );
257    }
258
259    #[test]
260    fn test_geocentric_with_height() {
261        let loc_sea_level = Location::from_degrees(0.0, 0.0, 0.0).unwrap();
262        let loc_elevated = Location::from_degrees(0.0, 0.0, 1000.0).unwrap();
263
264        let (u1, v1) = loc_sea_level.to_geocentric_km().unwrap();
265        let (u2, v2) = loc_elevated.to_geocentric_km().unwrap();
266
267        assert!(
268            (u2 - u1 - 1.0).abs() < 0.001,
269            "1km elevation should increase u by ~1km"
270        );
271        assert!(
272            (v2 - v1).abs() < 0.001,
273            "At equator, elevation shouldn't affect v much"
274        );
275    }
276
277    #[test]
278    fn test_negative_latitude() {
279        let loc = Location::from_degrees(-45.0, 0.0, 0.0).unwrap();
280        let (u, v) = loc.to_geocentric_km().unwrap();
281
282        assert!(u > 0.0, "u should be positive: {}", u);
283        assert!(
284            v < 0.0,
285            "v should be negative in southern hemisphere: {}",
286            v
287        );
288    }
289
290    #[test]
291    fn test_geocentric_division_by_zero() {
292        let north_pole = Location {
293            latitude: crate::constants::HALF_PI,
294            longitude: 0.0,
295            height: 0.0,
296        };
297
298        let result = north_pole.to_geocentric_km();
299        assert!(result.is_ok());
300    }
301}
celestial_core/location/geodesy.rs

celestial_core/location/
geodesy.rs
