1use crate::types::{CartesianPoint, GeoPoint, SphericalPoint};
2use std::f64::consts::TAU;
3
4pub 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
20pub 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
39pub 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
50pub 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
61pub fn cartesian_to_geo(p: &CartesianPoint) -> GeoPoint {
63 spherical_to_geo(&cartesian_to_spherical(p))
64}
65
66pub 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 #[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 #[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 #[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 #[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 #[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 #[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 #[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 #[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 #[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 #[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}