Skip to main content

sphereql_core/
conversions.rs

1use crate::types::{CartesianPoint, GeoPoint, SphericalPoint};
2use std::f64::consts::TAU;
3
4/// Converts a spherical point to Cartesian coordinates.
5///
6/// ```
7/// use sphereql_core::{SphericalPoint, spherical_to_cartesian};
8///
9/// let p = SphericalPoint::new(1.0, 0.0, std::f64::consts::FRAC_PI_2).unwrap();
10/// let c = spherical_to_cartesian(&p);
11/// assert!((c.x - 1.0).abs() < 1e-10);
12/// ```
13pub fn spherical_to_cartesian(p: &SphericalPoint) -> CartesianPoint {
14    let x = p.r * p.phi.sin() * p.theta.cos();
15    let y = p.r * p.phi.sin() * p.theta.sin();
16    let z = p.r * p.phi.cos();
17    CartesianPoint::new(x, y, z)
18}
19
20/// Converts a Cartesian point to spherical coordinates.
21///
22/// ```
23/// use sphereql_core::*;
24///
25/// let original = SphericalPoint::new(1.0, 0.5, 0.7).unwrap();
26/// let roundtrip = cartesian_to_spherical(&spherical_to_cartesian(&original));
27/// assert!((roundtrip.r - original.r).abs() < 1e-10);
28/// ```
29pub fn cartesian_to_spherical(p: &CartesianPoint) -> SphericalPoint {
30    let r = (p.x * p.x + p.y * p.y + p.z * p.z).sqrt();
31    if r < f64::EPSILON {
32        return SphericalPoint::new_unchecked(0.0, 0.0, 0.0);
33    }
34    let theta = normalize_theta(p.y.atan2(p.x));
35    let phi = (p.z / r).clamp(-1.0, 1.0).acos();
36    SphericalPoint::new_unchecked(r, theta, phi)
37}
38
39/// Converts a spherical point to geographic coordinates (lat/lon/alt).
40///
41/// Assumes the unit sphere (r=1.0) represents the surface. Altitude is
42/// the excess radius above 1.0, clamped to zero for sub-unit radii.
43pub fn spherical_to_geo(p: &SphericalPoint) -> GeoPoint {
44    let lat = 90.0 - p.phi.to_degrees();
45    let lon = theta_to_longitude(p.theta);
46    let alt = (p.r - 1.0).max(0.0);
47    GeoPoint::new_unchecked(lat, lon, alt)
48}
49
50/// Converts geographic coordinates to spherical coordinates.
51///
52/// Latitude maps to polar angle phi, longitude to azimuthal angle theta,
53/// and altitude adds to the unit radius.
54pub fn geo_to_spherical(p: &GeoPoint) -> SphericalPoint {
55    let phi = (90.0 - p.lat).to_radians();
56    let theta = normalize_theta(p.lon.to_radians());
57    let r = 1.0 + p.alt;
58    SphericalPoint::new_unchecked(r, theta, phi)
59}
60
61/// Converts a Cartesian point to geographic coordinates via spherical.
62pub fn cartesian_to_geo(p: &CartesianPoint) -> GeoPoint {
63    spherical_to_geo(&cartesian_to_spherical(p))
64}
65
66/// Converts geographic coordinates to Cartesian via spherical.
67pub fn geo_to_cartesian(p: &GeoPoint) -> CartesianPoint {
68    spherical_to_cartesian(&geo_to_spherical(p))
69}
70
71fn normalize_theta(mut theta: f64) -> f64 {
72    theta %= TAU;
73    if theta < 0.0 {
74        theta += TAU;
75    }
76    theta
77}
78
79fn theta_to_longitude(theta: f64) -> f64 {
80    let deg = theta.to_degrees();
81    if deg > 180.0 { deg - 360.0 } else { deg }
82}
83
84impl From<&SphericalPoint> for CartesianPoint {
85    fn from(p: &SphericalPoint) -> Self {
86        spherical_to_cartesian(p)
87    }
88}
89
90impl From<&CartesianPoint> for SphericalPoint {
91    fn from(p: &CartesianPoint) -> Self {
92        cartesian_to_spherical(p)
93    }
94}
95
96impl From<&SphericalPoint> for GeoPoint {
97    fn from(p: &SphericalPoint) -> Self {
98        spherical_to_geo(p)
99    }
100}
101
102impl From<&GeoPoint> for SphericalPoint {
103    fn from(p: &GeoPoint) -> Self {
104        geo_to_spherical(p)
105    }
106}
107
108impl From<&CartesianPoint> for GeoPoint {
109    fn from(p: &CartesianPoint) -> Self {
110        cartesian_to_geo(p)
111    }
112}
113
114impl From<&GeoPoint> for CartesianPoint {
115    fn from(p: &GeoPoint) -> Self {
116        geo_to_cartesian(p)
117    }
118}
119
120#[cfg(test)]
121mod tests {
122    use super::*;
123    use approx::assert_relative_eq;
124    use std::f64::consts::{FRAC_PI_2, PI, TAU};
125
126    fn spherical(r: f64, theta: f64, phi: f64) -> SphericalPoint {
127        SphericalPoint::new_unchecked(r, theta, phi)
128    }
129
130    fn geo(lat: f64, lon: f64, alt: f64) -> GeoPoint {
131        GeoPoint::new_unchecked(lat, lon, alt)
132    }
133
134    // --- Roundtrip: spherical -> cartesian -> spherical ---
135
136    #[test]
137    fn roundtrip_spherical_cartesian_unit_sphere() {
138        let original = spherical(1.0, 1.0, 0.8);
139        let cart = spherical_to_cartesian(&original);
140        let back = cartesian_to_spherical(&cart);
141        assert_relative_eq!(back.r, original.r, epsilon = 1e-12);
142        assert_relative_eq!(back.theta, original.theta, epsilon = 1e-12);
143        assert_relative_eq!(back.phi, original.phi, epsilon = 1e-12);
144    }
145
146    #[test]
147    fn roundtrip_spherical_cartesian_large_radius() {
148        let original = spherical(42.0, 3.5, 1.2);
149        let cart = spherical_to_cartesian(&original);
150        let back = cartesian_to_spherical(&cart);
151        assert_relative_eq!(back.r, original.r, epsilon = 1e-10);
152        assert_relative_eq!(back.theta, original.theta, epsilon = 1e-10);
153        assert_relative_eq!(back.phi, original.phi, epsilon = 1e-10);
154    }
155
156    #[test]
157    fn roundtrip_spherical_cartesian_near_tau() {
158        let original = spherical(1.0, TAU - 0.001, FRAC_PI_2);
159        let cart = spherical_to_cartesian(&original);
160        let back = cartesian_to_spherical(&cart);
161        assert_relative_eq!(back.r, original.r, epsilon = 1e-12);
162        assert_relative_eq!(back.theta, original.theta, epsilon = 1e-6);
163        assert_relative_eq!(back.phi, original.phi, epsilon = 1e-12);
164    }
165
166    // --- Roundtrip: spherical -> geo -> spherical ---
167
168    #[test]
169    fn roundtrip_spherical_geo_equator() {
170        let original = spherical(1.0, 0.5, FRAC_PI_2);
171        let g = spherical_to_geo(&original);
172        let back = geo_to_spherical(&g);
173        assert_relative_eq!(back.r, original.r, epsilon = 1e-12);
174        assert_relative_eq!(back.theta, original.theta, epsilon = 1e-12);
175        assert_relative_eq!(back.phi, original.phi, epsilon = 1e-12);
176    }
177
178    #[test]
179    fn roundtrip_spherical_geo_with_altitude() {
180        let original = spherical(2.5, 1.0, 0.8);
181        let g = spherical_to_geo(&original);
182        let back = geo_to_spherical(&g);
183        assert_relative_eq!(back.r, original.r, epsilon = 1e-12);
184        assert_relative_eq!(back.theta, original.theta, epsilon = 1e-12);
185        assert_relative_eq!(back.phi, original.phi, epsilon = 1e-12);
186    }
187
188    #[test]
189    fn roundtrip_geo_spherical_negative_lon() {
190        let original = geo(45.0, -90.0, 0.0);
191        let s = geo_to_spherical(&original);
192        let back = spherical_to_geo(&s);
193        assert_relative_eq!(back.lat, original.lat, epsilon = 1e-12);
194        assert_relative_eq!(back.lon, original.lon, epsilon = 1e-12);
195        assert_relative_eq!(back.alt, original.alt, epsilon = 1e-12);
196    }
197
198    // --- Known values: origin ---
199
200    #[test]
201    fn origin_spherical_to_cartesian() {
202        let p = spherical(0.0, 0.0, 0.0);
203        let c = spherical_to_cartesian(&p);
204        assert_relative_eq!(c.x, 0.0, epsilon = 1e-15);
205        assert_relative_eq!(c.y, 0.0, epsilon = 1e-15);
206        assert_relative_eq!(c.z, 0.0, epsilon = 1e-15);
207    }
208
209    #[test]
210    fn origin_cartesian_to_spherical() {
211        let c = CartesianPoint::new(0.0, 0.0, 0.0);
212        let s = cartesian_to_spherical(&c);
213        assert_relative_eq!(s.r, 0.0, epsilon = 1e-15);
214        assert_relative_eq!(s.theta, 0.0, epsilon = 1e-15);
215        assert_relative_eq!(s.phi, 0.0, epsilon = 1e-15);
216    }
217
218    #[test]
219    fn near_zero_cartesian_to_spherical() {
220        let c = CartesianPoint::new(1e-20, 1e-20, 1e-20);
221        let s = cartesian_to_spherical(&c);
222        assert_relative_eq!(s.r, 0.0, epsilon = 1e-15);
223    }
224
225    // --- Known values: north pole (phi=0) ---
226
227    #[test]
228    fn north_pole_spherical_to_cartesian() {
229        let p = spherical(1.0, 0.0, 0.0);
230        let c = spherical_to_cartesian(&p);
231        assert_relative_eq!(c.x, 0.0, epsilon = 1e-15);
232        assert_relative_eq!(c.y, 0.0, epsilon = 1e-15);
233        assert_relative_eq!(c.z, 1.0, epsilon = 1e-15);
234    }
235
236    #[test]
237    fn north_pole_to_geo() {
238        let p = spherical(1.0, 0.0, 0.0);
239        let g = spherical_to_geo(&p);
240        assert_relative_eq!(g.lat, 90.0, epsilon = 1e-12);
241        assert_relative_eq!(g.alt, 0.0, epsilon = 1e-12);
242    }
243
244    // --- Known values: south pole (phi=π) ---
245
246    #[test]
247    fn south_pole_spherical_to_cartesian() {
248        let p = spherical(1.0, 0.0, PI);
249        let c = spherical_to_cartesian(&p);
250        assert_relative_eq!(c.x, 0.0, epsilon = 1e-15);
251        assert_relative_eq!(c.y, 0.0, epsilon = 1e-15);
252        assert_relative_eq!(c.z, -1.0, epsilon = 1e-15);
253    }
254
255    #[test]
256    fn south_pole_to_geo() {
257        let p = spherical(1.0, 0.0, PI);
258        let g = spherical_to_geo(&p);
259        assert_relative_eq!(g.lat, -90.0, epsilon = 1e-12);
260        assert_relative_eq!(g.alt, 0.0, epsilon = 1e-12);
261    }
262
263    // --- Known values: equator points ---
264
265    #[test]
266    fn equator_theta_zero() {
267        let p = spherical(1.0, 0.0, FRAC_PI_2);
268        let c = spherical_to_cartesian(&p);
269        assert_relative_eq!(c.x, 1.0, epsilon = 1e-15);
270        assert_relative_eq!(c.y, 0.0, epsilon = 1e-15);
271        assert_relative_eq!(c.z, 0.0, epsilon = 1e-15);
272    }
273
274    #[test]
275    fn equator_theta_half_pi() {
276        let p = spherical(1.0, FRAC_PI_2, FRAC_PI_2);
277        let c = spherical_to_cartesian(&p);
278        assert_relative_eq!(c.x, 0.0, epsilon = 1e-15);
279        assert_relative_eq!(c.y, 1.0, epsilon = 1e-15);
280        assert_relative_eq!(c.z, 0.0, epsilon = 1e-15);
281    }
282
283    #[test]
284    fn equator_theta_pi() {
285        let p = spherical(1.0, PI, FRAC_PI_2);
286        let c = spherical_to_cartesian(&p);
287        assert_relative_eq!(c.x, -1.0, epsilon = 1e-15);
288        assert_relative_eq!(c.y, 0.0, epsilon = 1e-12);
289        assert_relative_eq!(c.z, 0.0, epsilon = 1e-15);
290    }
291
292    #[test]
293    fn equator_to_geo_gives_zero_lat() {
294        let p = spherical(1.0, 0.5, FRAC_PI_2);
295        let g = spherical_to_geo(&p);
296        assert_relative_eq!(g.lat, 0.0, epsilon = 1e-12);
297    }
298
299    // --- Edge cases: theta boundaries ---
300
301    #[test]
302    fn theta_zero_roundtrip() {
303        let original = spherical(1.0, 0.0, 1.0);
304        let back = cartesian_to_spherical(&spherical_to_cartesian(&original));
305        assert_relative_eq!(back.theta, 0.0, epsilon = 1e-12);
306    }
307
308    #[test]
309    fn theta_near_tau_roundtrip() {
310        let original = spherical(1.0, TAU - 1e-10, 1.0);
311        let back = cartesian_to_spherical(&spherical_to_cartesian(&original));
312        assert_relative_eq!(back.theta, original.theta, epsilon = 1e-6);
313    }
314
315    // --- Edge cases: phi boundaries ---
316
317    #[test]
318    fn phi_zero_roundtrip() {
319        let p = spherical(1.0, 0.0, 0.0);
320        let c = spherical_to_cartesian(&p);
321        let back = cartesian_to_spherical(&c);
322        assert_relative_eq!(back.r, 1.0, epsilon = 1e-12);
323        assert_relative_eq!(back.phi, 0.0, epsilon = 1e-12);
324    }
325
326    #[test]
327    fn phi_pi_roundtrip() {
328        let p = spherical(1.0, 0.0, PI);
329        let c = spherical_to_cartesian(&p);
330        let back = cartesian_to_spherical(&c);
331        assert_relative_eq!(back.r, 1.0, epsilon = 1e-12);
332        assert_relative_eq!(back.phi, PI, epsilon = 1e-12);
333    }
334
335    // --- Geo edge cases ---
336
337    #[test]
338    fn geo_longitude_180_boundary() {
339        let p = spherical(1.0, PI, FRAC_PI_2);
340        let g = spherical_to_geo(&p);
341        assert_relative_eq!(g.lon, 180.0, epsilon = 1e-12);
342    }
343
344    #[test]
345    fn geo_negative_longitude() {
346        let p = spherical(1.0, PI + 0.5, FRAC_PI_2);
347        let g = spherical_to_geo(&p);
348        assert!(g.lon < 0.0);
349        let back = geo_to_spherical(&g);
350        assert_relative_eq!(back.theta, PI + 0.5, epsilon = 1e-12);
351    }
352
353    #[test]
354    fn geo_altitude_clamped_to_zero() {
355        let p = spherical(0.5, 0.0, FRAC_PI_2);
356        let g = spherical_to_geo(&p);
357        assert_relative_eq!(g.alt, 0.0, epsilon = 1e-15);
358    }
359
360    // --- From trait tests ---
361
362    #[test]
363    fn from_trait_spherical_to_cartesian() {
364        let s = spherical(1.0, 0.0, FRAC_PI_2);
365        let c: CartesianPoint = (&s).into();
366        assert_relative_eq!(c.x, 1.0, epsilon = 1e-15);
367    }
368
369    #[test]
370    fn from_trait_cartesian_to_spherical() {
371        let c = CartesianPoint::new(1.0, 0.0, 0.0);
372        let s: SphericalPoint = (&c).into();
373        assert_relative_eq!(s.r, 1.0, epsilon = 1e-12);
374        assert_relative_eq!(s.theta, 0.0, epsilon = 1e-12);
375        assert_relative_eq!(s.phi, FRAC_PI_2, epsilon = 1e-12);
376    }
377
378    #[test]
379    fn from_trait_geo_to_cartesian() {
380        let g = geo(0.0, 0.0, 0.0);
381        let c: CartesianPoint = (&g).into();
382        assert_relative_eq!(c.x, 1.0, epsilon = 1e-12);
383        assert_relative_eq!(c.y, 0.0, epsilon = 1e-12);
384        assert_relative_eq!(c.z, 0.0, epsilon = 1e-12);
385    }
386
387    #[test]
388    fn from_trait_cartesian_to_geo() {
389        let c = CartesianPoint::new(0.0, 0.0, 1.0);
390        let g: GeoPoint = (&c).into();
391        assert_relative_eq!(g.lat, 90.0, epsilon = 1e-12);
392        assert_relative_eq!(g.alt, 0.0, epsilon = 1e-12);
393    }
394}