globe_rs/
geographic.rs

1use crate::CartesianPoint;
2use std::f64::consts::{FRAC_PI_2, PI};
3use std::ops::{Add, Mul, Neg, Sub};
4use wasm_bindgen::prelude::wasm_bindgen;
5
6/// Represents a point using the geographic system of coordinates.
7#[wasm_bindgen]
8#[derive(Debug, Default, Clone, Copy, PartialEq)]
9pub struct GeographicPoint {
10    longitude: f64,
11    latitude: f64,
12    altitude: f64,
13}
14
15impl From<CartesianPoint> for GeographicPoint {
16    fn from(value: CartesianPoint) -> Self {
17        GeographicPoint::from_cartesian(&value)
18    }
19}
20
21#[wasm_bindgen]
22impl GeographicPoint {
23    pub fn new(longitude: f64, latitude: f64, altitude: f64) -> Self {
24        Self {
25            longitude,
26            latitude,
27            altitude,
28        }
29    }
30
31    /// Returns the equivalent [`GeographicPoint`] of the given [`CartesianPoint`]
32    pub fn from_cartesian(point: &CartesianPoint) -> Self {
33        // see: https://en.wikipedia.org/wiki/Spherical_coordinate_system
34        let radius = f64::sqrt(point.x().powi(2) + point.y().powi(2) + point.z().powi(2));
35
36        let theta = match (point.x(), point.y(), point.z()) {
37            (x, y, z) if z > 0. => f64::atan(f64::sqrt(x.powi(2) + y.powi(2)) / z),
38            (x, y, z) if z < 0. => PI + f64::atan(f64::sqrt(x.powi(2) + y.powi(2)) / z),
39            (x, y, z) if z == 0. && x * y != 0. => FRAC_PI_2,
40            (x, y, z) if x == y && y == z => FRAC_PI_2, // fallback value
41
42            _ => FRAC_PI_2, // fallback value
43        };
44
45        let phi = match (point.x(), point.y()) {
46            (x, y) if x > 0. => f64::atan(y / x),
47            (x, y) if x < 0. && y >= 0. => f64::atan(y / x) + PI,
48            (x, y) if x < 0. && y < 0. => f64::atan(y / x) - PI,
49            (x, y) if x == 0. && y > 0. => FRAC_PI_2,
50            (x, y) if x == 0. && y < 0. => -FRAC_PI_2,
51            (x, y) if x == 0. && y == 0. => 0., // fallback value
52
53            _ => 0., // fallback value
54        };
55
56        Self::default()
57            .with_longitude(phi)
58            .with_latitude(FRAC_PI_2 - theta)
59            .with_altitude(radius)
60    }
61
62    /// Calls set_longitude on self and returns it.
63    pub fn with_longitude(mut self, value: f64) -> Self {
64        self.set_longitude(value);
65        self
66    }
67
68    /// Calls set_latitude on self and returns it.
69    pub fn with_latitude(mut self, value: f64) -> Self {
70        self.set_latitude(value);
71        self
72    }
73
74    /// Calls set_altitude on self and returns it.
75    pub fn with_altitude(mut self, value: f64) -> Self {
76        self.set_altitude(value);
77        self
78    }
79
80    /// Sets the given longitude (in radiants) to the point.
81    ///
82    /// ## Definition
83    /// Since the longitude of a point on a sphere is the angle east (positive) or
84    /// west (negative) in reference of the maridian zero, the longitude value must
85    /// be in the range __[-π, +π)__. Any other value will be recomputed in order
86    /// to set its equivalent inside the range.
87    ///
88    /// ### Longitude adjustment
89    /// Both boundaries of the longitude range are consecutive, which means that
90    /// overflowing one is the same as continuing from the other in the same
91    /// direction.
92    ///
93    /// ## Example
94    /// ```
95    /// use globe_rs::GeographicPoint;
96    /// use std::f64::consts::PI;
97    /// use float_cmp::approx_eq;
98    ///
99    /// let mut point = GeographicPoint::default();
100    /// point.set_longitude(PI + 1_f64);
101    ///
102    /// assert!(approx_eq!(f64, point.longitude(), -PI + 1_f64, ulps = 2));
103    /// ```
104    pub fn set_longitude(&mut self, value: f64) {
105        self.longitude = (-PI..=PI)
106            .contains(&value)
107            .then_some(value)
108            .unwrap_or_else(|| {
109                // Both boundaries of the range are consecutive, which means that
110                // overflowing one is the same as continuing from the other in the
111                // same direction.
112                value.add(PI).rem_euclid(2_f64.mul(PI)).sub(PI)
113            })
114    }
115
116    /// Sets the given latitude (in radiants) to the point.
117    ///
118    /// ## Definition
119    /// Since the latitude of a point on a sphere is the angle between the
120    /// equatorial plane and the straight line that passes through that point and
121    /// through the center of the sphere, the latitude value must be in the range
122    /// __\[-π/2, +π/2\]__. Any other value will be recomputed in order to set its
123    /// equivalent inside the range. Notice that this action may recompute the
124    /// longitude as well.
125    ///
126    /// ### Latitude adjustment
127    /// Overflowing any of both boundaries of the latitude range behaves like
128    /// moving away from that point and getting closer to the oposite one.
129    ///
130    /// ### Longitude adjustment
131    /// Geometrically speaking, meridians are half of a circle going from the north
132    /// pole to the south one. The position of each meridian in the perimeter of
133    /// the sphere (horizontal axis) is set by the longitude itself. However, this
134    /// value may change when the latitude overflows its normalized range. This
135    /// happen since exceeding any of its established limits means moving from one
136    /// to the other half of the circle on which the meridian is drawn. And
137    /// therefore, the longitude gets increased by exactly `π` radiants.
138    ///
139    /// Of course, this mutation on the longitude only applies when the overflow of
140    /// the latitude is not enough to complete a full lap. If it is, the longitude
141    /// does not change at all.
142    ///
143    /// ## Example
144    /// ```
145    /// use globe_rs::GeographicPoint;
146    /// use std::f64::consts::PI;
147    /// use float_cmp::approx_eq;
148    ///
149    /// let mut point = GeographicPoint::default();
150    /// point.set_latitude(-5. * PI / 4.);
151    ///
152    /// assert!(approx_eq!(f64, point.latitude(), PI / 4., ulps = 2));
153    /// assert!(approx_eq!(f64, point.longitude(), -PI, ulps = 2));
154    /// ```
155    pub fn set_latitude(&mut self, value: f64) {
156        self.latitude = (-FRAC_PI_2..=FRAC_PI_2)
157            .contains(&value)
158            .then_some(value)
159            .unwrap_or_else(|| {
160                // The derivative of sin(x) is cos(x), and so, cos(x) determines if
161                // the sign of the longitude of the point must change.
162                if value.cos().is_sign_negative() {
163                    // Increasing the longitude of the point by ±π radiants (180º)
164                    // ensures the sign is changed while maintaining it in the same
165                    // pair of complementary meridians.
166                    let direction = self.longitude.cos().signum().neg(); // invert direction
167                    let rotation = PI * direction;
168                    self.set_longitude(self.longitude + rotation);
169                }
170
171                value.sin().asin()
172            });
173    }
174
175    /// Sets the given altitude to the point.
176    pub fn set_altitude(&mut self, value: f64) {
177        self.altitude = value;
178    }
179
180    /// Returns the longitude (in radiants) of the point.
181    pub fn longitude(&self) -> f64 {
182        self.longitude
183    }
184
185    /// Returns the latitude (in radiants) of the point.
186    pub fn latitude(&self) -> f64 {
187        self.latitude
188    }
189
190    /// Returns the altitude (in radiants) of the point.
191    pub fn altitude(&self) -> f64 {
192        self.altitude
193    }
194
195    /// Returns the result of dividing `π` to the longitude of the point, resulting
196    /// in a value in the range __[-1.0, 1.0)__
197    pub fn long_ratio(&self) -> f64 {
198        self.longitude / PI
199    }
200
201    /// Returns the result of dividing `π/2` to the latitude of the point, resulting
202    /// in a value in the range __\[-1.0, 1.0\]__
203    pub fn lat_ratio(&self) -> f64 {
204        self.latitude / FRAC_PI_2
205    }
206
207    /// Computes the [great-circle distance](https://en.wikipedia.org/wiki/Great-circle_distance) from self to the given point (in radiants).
208    pub fn distance(&self, other: &GeographicPoint) -> f64 {
209        (self.latitude().sin() * other.latitude().sin()
210            + self.latitude().cos()
211                * other.latitude().cos()
212                * (self.longitude() - other.longitude()).abs().cos())
213        .acos()
214    }
215}
216
217#[cfg(test)]
218mod tests {
219    use super::*;
220    use float_cmp::approx_eq;
221
222    const ULPS: i64 = 2;
223
224    #[test]
225    fn longitude_must_not_exceed_boundaries() {
226        struct TestCase {
227            name: &'static str,
228            input: f64,
229            longitude: f64,
230        }
231
232        vec![
233            TestCase {
234                name: "positive longitude value must not change",
235                input: 1.,
236                longitude: 1.,
237            },
238            TestCase {
239                name: "negative longitude value must not change",
240                input: -3.,
241                longitude: -3.,
242            },
243            TestCase {
244                name: "positive overflowing longitude must change",
245                input: PI + 1.,
246                longitude: -PI + 1.,
247            },
248            TestCase {
249                name: "negative overflowing longitude must change",
250                input: -PI - 1.,
251                longitude: PI - 1.,
252            },
253        ]
254        .into_iter()
255        .for_each(|test_case| {
256            let point = GeographicPoint::default().with_longitude(test_case.input);
257            assert!(
258                approx_eq!(f64, point.longitude, test_case.longitude, ulps = ULPS),
259                "{}: {} ±ε = {}",
260                test_case.name,
261                point.longitude,
262                test_case.longitude
263            );
264        });
265    }
266
267    #[test]
268    fn latitude_must_not_exceed_boundaries() {
269        struct TestCase {
270            name: &'static str,
271            input: f64,
272            latitude: f64,
273        }
274
275        vec![
276            TestCase {
277                name: "positive latitude value must not change",
278                input: 1.,
279                latitude: 1.,
280            },
281            TestCase {
282                name: "negative latitude value must not change",
283                input: -1.,
284                latitude: -1.,
285            },
286            TestCase {
287                name: "positive overflowing latitude must change",
288                input: 7. * PI / 4.,
289                latitude: -PI / 4.,
290            },
291            TestCase {
292                name: "negative overflowing latidude must change",
293                input: -7. * PI / 4.,
294                latitude: PI / 4.,
295            },
296        ]
297        .into_iter()
298        .for_each(|test_case| {
299            let point = GeographicPoint::default().with_latitude(test_case.input);
300            assert!(
301                approx_eq!(f64, point.latitude, test_case.latitude, ulps = ULPS),
302                "{}: {} ±ε = {}",
303                test_case.name,
304                point.latitude,
305                test_case.latitude
306            );
307        });
308    }
309
310    #[test]
311    fn longitude_may_change_based_on_latitude() {
312        struct TestCase {
313            name: &'static str,
314            input: f64,
315            longitude: f64,
316        }
317
318        vec![
319            TestCase {
320                name: "positive latitude value must not change longitude",
321                input: 1.,
322                longitude: 0.,
323            },
324            TestCase {
325                name: "negative latitude value must not change longitude",
326                input: -1.,
327                longitude: 0.,
328            },
329            TestCase {
330                name: "positive overflowing latitude must not change longitude",
331                input: 7. * PI / 4.,
332                longitude: 0.,
333            },
334            TestCase {
335                name: "negative overflowing latidude must not change longitude",
336                input: -7. * PI / 4.,
337                longitude: 0.,
338            },
339            TestCase {
340                name: "positive overflowing latitude must change longitude",
341                input: 5. * PI / 4.,
342                longitude: -PI,
343            },
344            TestCase {
345                name: "negative overflowing latidude must change longitude",
346                input: -PI,
347                longitude: -PI,
348            },
349        ]
350        .into_iter()
351        .for_each(|test_case| {
352            let point = GeographicPoint::default().with_latitude(test_case.input);
353            assert!(
354                approx_eq!(f64, point.longitude, test_case.longitude, ulps = ULPS),
355                "{}: {} ±ε = {}",
356                test_case.name,
357                point.longitude,
358                test_case.longitude
359            );
360        });
361    }
362
363    #[test]
364    fn long_ratio_must_not_exceed_boundaries() {
365        struct TestCase {
366            name: &'static str,
367            longitude: f64,
368            ratio: f64,
369        }
370
371        vec![
372            TestCase {
373                name: "zero longitude must be equals to zero",
374                longitude: 0.,
375                ratio: 0.,
376            },
377            TestCase {
378                name: "positive longitude boundary must equals to positive one",
379                longitude: PI,
380                ratio: 1.,
381            },
382            TestCase {
383                name: "arbitrary positive longitude must be positive",
384                longitude: FRAC_PI_2,
385                ratio: 0.5,
386            },
387            TestCase {
388                name: "negative longitude boundary must equals to negative one",
389                longitude: -PI,
390                ratio: -1.,
391            },
392            TestCase {
393                name: "arbitrary negative longitude must be negative",
394                longitude: -FRAC_PI_2,
395                ratio: -0.5,
396            },
397        ]
398        .into_iter()
399        .for_each(|test_case| {
400            let point = GeographicPoint::default().with_longitude(test_case.longitude);
401            assert!(
402                approx_eq!(f64, point.long_ratio(), test_case.ratio, ulps = ULPS),
403                "{}: {} ±ε = {}",
404                test_case.name,
405                point.long_ratio(),
406                test_case.ratio
407            );
408        });
409    }
410
411    #[test]
412    fn lat_ratio_must_not_exceed_boundaries() {
413        struct TestCase {
414            name: &'static str,
415            latitude: f64,
416            ratio: f64,
417        }
418
419        vec![
420            TestCase {
421                name: "zero latitude must be equals to zero",
422                latitude: 0.,
423                ratio: 0.,
424            },
425            TestCase {
426                name: "positive latitude boundary must equals to one",
427                latitude: FRAC_PI_2,
428                ratio: 1.,
429            },
430            TestCase {
431                name: "arbitrary positive latitude must be positive",
432                latitude: FRAC_PI_2 / 2.,
433                ratio: 0.5,
434            },
435            TestCase {
436                name: "negative latitude boundary must equals to negative one",
437                latitude: -FRAC_PI_2,
438                ratio: -1.,
439            },
440            TestCase {
441                name: "arbitrary negative latitude must be negative",
442                latitude: -FRAC_PI_2 / 2.,
443                ratio: -0.5,
444            },
445        ]
446        .into_iter()
447        .for_each(|test_case| {
448            let point = GeographicPoint::default().with_latitude(test_case.latitude);
449            assert!(
450                approx_eq!(f64, point.lat_ratio(), test_case.ratio, ulps = ULPS),
451                "{}: {} ±ε = {}",
452                test_case.name,
453                point.lat_ratio(),
454                test_case.ratio
455            );
456        });
457    }
458
459    #[test]
460    fn basic_operations_must_be_consistent() {
461        let mut point = GeographicPoint::default()
462            .with_longitude(-FRAC_PI_2)
463            .with_latitude(FRAC_PI_2 / 2.);
464
465        point.set_latitude(point.latitude().add(PI));
466
467        assert!(
468            approx_eq!(f64, point.longitude(), FRAC_PI_2, ulps = ULPS),
469            "longitude must switch to positive: {} ±ε = {}",
470            point.longitude(),
471            FRAC_PI_2
472        );
473
474        assert!(
475            approx_eq!(f64, point.latitude(), -FRAC_PI_2 / 2., ulps = ULPS),
476            "latitude must switch to negative: {} ±ε = {}",
477            point.latitude(),
478            -FRAC_PI_2 / 2.
479        );
480    }
481
482    #[test]
483    fn geographic_from_cartesian_must_not_fail() {
484        struct TestCase {
485            name: &'static str,
486            geographic: GeographicPoint,
487            cartesian: CartesianPoint,
488        }
489
490        vec![
491            TestCase {
492                name: "north point",
493                geographic: GeographicPoint::default()
494                    .with_latitude(FRAC_PI_2)
495                    .with_altitude(1.),
496                cartesian: CartesianPoint::new(0., 0., 1.),
497            },
498            TestCase {
499                name: "south point",
500                geographic: GeographicPoint::default()
501                    .with_latitude(-FRAC_PI_2)
502                    .with_altitude(1.),
503                cartesian: CartesianPoint::new(0., 0., -1.),
504            },
505            TestCase {
506                name: "east point",
507                geographic: GeographicPoint::default()
508                    .with_longitude(FRAC_PI_2)
509                    .with_altitude(1.),
510                cartesian: CartesianPoint::new(0., 1., 0.),
511            },
512            TestCase {
513                name: "weast point",
514                geographic: GeographicPoint::default()
515                    .with_longitude(-FRAC_PI_2)
516                    .with_altitude(1.),
517                cartesian: CartesianPoint::new(0., -1., 0.),
518            },
519            TestCase {
520                name: "front point",
521                geographic: GeographicPoint::default().with_altitude(1.),
522                cartesian: CartesianPoint::new(1., 0., 0.),
523            },
524            TestCase {
525                name: "back point",
526                geographic: GeographicPoint::default()
527                    .with_longitude(PI)
528                    .with_altitude(1.),
529                cartesian: CartesianPoint::new(-1., 0., 0.),
530            },
531        ]
532        .into_iter()
533        .for_each(|test_case| {
534            let point = GeographicPoint::from(test_case.cartesian);
535            assert!(
536                approx_eq!(
537                    f64,
538                    point.longitude(),
539                    test_case.geographic.longitude(),
540                    ulps = ULPS
541                ),
542                "{}: longitude {:#?} ±ε == {:#?}",
543                test_case.name,
544                point.longitude(),
545                test_case.geographic.longitude(),
546            );
547
548            assert!(
549                approx_eq!(
550                    f64,
551                    point.latitude(),
552                    test_case.geographic.latitude(),
553                    ulps = ULPS
554                ),
555                "{}: latitude {:#?} ±ε == {:#?}",
556                test_case.name,
557                point.latitude(),
558                test_case.geographic.latitude(),
559            );
560
561            assert!(
562                approx_eq!(
563                    f64,
564                    point.altitude(),
565                    test_case.geographic.altitude(),
566                    ulps = ULPS
567                ),
568                "{}: altitude {:#?} ±ε == {:#?}",
569                test_case.name,
570                point.altitude(),
571                test_case.geographic.altitude(),
572            );
573        });
574    }
575
576    #[test]
577    fn distance_must_not_fail() {
578        struct TestCase<'a> {
579            name: &'a str,
580            from: GeographicPoint,
581            to: GeographicPoint,
582            distance: f64,
583        }
584
585        vec![
586            TestCase {
587                name: "Same point must be zero",
588                from: GeographicPoint::default(),
589                to: GeographicPoint::default(),
590                distance: 0.,
591            },
592            TestCase {
593                name: "Oposite points in the horizontal",
594                from: GeographicPoint::default(),
595                to: GeographicPoint::default().with_longitude(-PI),
596                distance: PI,
597            },
598            TestCase {
599                name: "Oposite points in the vertical",
600                from: GeographicPoint::default().with_latitude(FRAC_PI_2),
601                to: GeographicPoint::default().with_latitude(-FRAC_PI_2),
602                distance: PI,
603            },
604        ]
605        .into_iter()
606        .for_each(|test_case| {
607            let got = test_case.from.distance(&test_case.to);
608
609            assert!(
610                approx_eq!(f64, got, test_case.distance, ulps = ULPS),
611                "{}: distance {:#?} ±ε == {:#?}",
612                test_case.name,
613                got,
614                test_case.distance,
615            )
616        });
617    }
618}
globe_rs/geographic.rs

globe_rs/
geographic.rs
