Rhumb no longer goes past pole

Author: joe-saronic

geo/src/algorithm/rhumb/bearing.rs

@@ -83,7 +83,7 @@ mod test {
     #[test]
     fn consistent_with_destination() {
         let p_1 = point!(x: 9.177789688110352f64, y: 48.776781529534965);
-        let p_2 = p_1.rhumb_destination(45., 10000.);
+        let p_2 = p_1.rhumb_destination(45., 10000.).unwrap();
 
         let b_1 = p_1.rhumb_bearing(p_2);
         assert_relative_eq!(b_1, 45., epsilon = 1.0e-6);

geo/src/algorithm/rhumb/destination.rs

@@ -12,6 +12,9 @@ pub trait RhumbDestination<T: CoordFloat> {
     /// Returns the destination Point having travelled the given distance along a [rhumb line]
     /// from the origin Point with the given bearing
     ///
+    /// A rhumb line has a finite length, so if the path distance exceeds the location of the pole,
+    /// the result is None.
+    ///
     /// # Units
     ///
     /// - `bearing`: degrees, zero degrees is north
@@ -24,18 +27,18 @@ pub trait RhumbDestination<T: CoordFloat> {
     /// use geo::Point;
     ///
     /// let p_1 = Point::new(9.177789688110352, 48.776781529534965);
-    /// let p_2 = p_1.rhumb_destination(45., 10000.);
+    /// let p_2 = p_1.rhumb_destination(45., 10000.).unwrap();
     /// assert_eq!(p_2, Point::new(9.274348757829898, 48.84037308229984))
     /// ```
     /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line
-    fn rhumb_destination(&self, bearing: T, distance: T) -> Point<T>;
+    fn rhumb_destination(&self, bearing: T, distance: T) -> Option<Point<T>>;
 }
 
 impl<T> RhumbDestination<T> for Point<T>
 where
     T: CoordFloat + FromPrimitive,
 {
-    fn rhumb_destination(&self, bearing: T, distance: T) -> Point<T> {
+    fn rhumb_destination(&self, bearing: T, distance: T) -> Option<Point<T>> {
         let delta = distance / T::from(MEAN_EARTH_RADIUS).unwrap(); // angular distance in radians
         let lambda1 = self.x().to_radians();
         let phi1 = self.y().to_radians();
@@ -49,32 +52,38 @@ where
 mod test {
     use super::*;
     use crate::RhumbDistance;
-    use num_traits::pow;
 
     #[test]
     fn returns_a_new_point() {
         let p_1 = Point::new(9.177789688110352, 48.776781529534965);
-        let p_2 = p_1.rhumb_destination(45., 10000.);
+        let p_2 = p_1.rhumb_destination(45., 10000.).unwrap();
         assert_eq!(p_2, Point::new(9.274348757829898, 48.84037308229984));
         let distance = p_1.rhumb_distance(&p_2);
-        assert_relative_eq!(distance, 10000., epsilon = 1.0e-6)
+        assert_relative_eq!(distance, 10000., epsilon = 1.0e-6);
     }
 
     #[test]
     fn direct_and_indirect_destinations_are_close() {
         let p_1 = Point::new(9.177789688110352, 48.776781529534965);
-        let p_2 = p_1.rhumb_destination(45., 10000.);
-        let square_edge = { pow(10000., 2) / 2f64 }.sqrt();
-        let p_3 = p_1.rhumb_destination(0., square_edge);
-        let p_4 = p_3.rhumb_destination(90., square_edge);
+        let p_2 = p_1.rhumb_destination(45., 10000.).unwrap();
+        let square_edge = 10000.0 * 0.5f64.sqrt();
+        let p_3 = p_1.rhumb_destination(0., square_edge).unwrap();
+        let p_4 = p_3.rhumb_destination(90., square_edge).unwrap();
         assert_relative_eq!(p_4, p_2, epsilon = 1.0e-3);
     }
 
     #[test]
     fn bearing_zero_is_north() {
         let p_1 = Point::new(9.177789688110352, 48.776781529534965);
-        let p_2 = p_1.rhumb_destination(0., 1000.);
+        let p_2 = p_1.rhumb_destination(0., 1000.).unwrap();
         assert_relative_eq!(p_1.x(), p_2.x(), epsilon = 1.0e-6);
-        assert!(p_2.y() > p_1.y())
+        assert!(p_2.y() > p_1.y());
+    }
+
+    #[test]
+    fn past_the_pole() {
+        let p = Point::new(9.177789688110352, 48.776781529534965);
+        let result = p.rhumb_destination(0., 1e8);
+        assert!(result.is_none());
     }
 }

geo/src/algorithm/rhumb/mod.rs

@@ -53,7 +53,7 @@ impl<T: CoordFloat + FromPrimitive> RhumbCalculations<T> {
         let delta_psi = ((phi2 / two + pi / four).tan() / (phi1 / two + pi / four).tan()).ln();
         let delta_phi = phi2 - phi1;
 
-        RhumbCalculations {
+        Self {
             from: *from,
             to: *to,
             phi1,
@@ -82,7 +82,7 @@ impl<T: CoordFloat + FromPrimitive> RhumbCalculations<T> {
         let delta = fraction * self.delta();
         let theta = self.theta();
         let lambda1 = self.from.x().to_radians();
-        calculate_destination(delta, lambda1, self.phi1, theta)
+        calculate_destination(delta, lambda1, self.phi1, theta).unwrap()
     }
 
     fn intermediate_fill(&self, max_delta: T, include_ends: bool) -> Vec<Point<T>> {
@@ -112,7 +112,7 @@ impl<T: CoordFloat + FromPrimitive> RhumbCalculations<T> {
         while current_step < T::one() {
             let delta = total_delta * current_step;
             let point = calculate_destination(delta, lambda1, self.phi1, theta);
-            points.push(point);
+            points.push(point.unwrap());
             current_step = current_step + interval;
         }
 
@@ -129,22 +129,18 @@ fn calculate_destination<T: CoordFloat + FromPrimitive>(
     lambda1: T,
     phi1: T,
     theta: T,
-) -> Point<T> {
+) -> Option<Point<T>> {
     let pi = T::from(std::f64::consts::PI).unwrap();
     let two = T::one() + T::one();
     let four = two + two;
     let threshold = T::from(10.0e-12).unwrap();
 
     let delta_phi = delta * theta.cos();
-    let mut phi2 = phi1 + delta_phi;
+    let phi2 = phi1 + delta_phi;
 
     // check for some daft bugger going past the pole, normalise latitude if so
     if phi2.abs() > pi / two {
-        phi2 = if phi2 > T::zero() {
-            pi - phi2
-        } else {
-            -pi - phi2
-        };
+        return None;
     }
 
     let delta_psi = ((phi2 / two + pi / four).tan() / (phi1 / two + pi / four).tan()).ln();
@@ -158,8 +154,8 @@ fn calculate_destination<T: CoordFloat + FromPrimitive>(
     let delta_lambda = (delta * theta.sin()) / q;
     let lambda2 = lambda1 + delta_lambda;
 
-    point! {
+    Some(point! {
         x: normalize_longitude(lambda2.to_degrees()),
         y: phi2.to_degrees(),
-    }
+    })
 }

geo/src/lib.rs

@@ -242,7 +242,7 @@ extern crate log;
 /// https://link.springer.com/article/10.1007%2Fs001900050278
 /// https://sci-hub.se/https://doi.org/10.1007/s001900050278
 /// https://en.wikipedia.org/wiki/Earth_radius#Mean_radius
-const MEAN_EARTH_RADIUS: f64 = 6371008.8;
+const MEAN_EARTH_RADIUS: f64 = 6_371_008.8;
 
 // Radius of Earth at the equator in meters (derived from the WGS-84 ellipsoid)
 const EQUATORIAL_EARTH_RADIUS: f64 = 6_378_137.0;