Corrected floating point error in rhumb line end point and distance calculation. Ported solution from WWA and updated a few unit tests.

Author: Zach Glueckert

Diff for: src/gov/nasa/worldwind/geom/LatLon.java

@@ -24,6 +24,12 @@ public class LatLon
 {
     public static final LatLon ZERO = new LatLon(Angle.ZERO, Angle.ZERO);
 
+    /**
+     * A near zero threshold used in some of the rhumb line calculations where floating point calculations cause
+     * errors.
+     */
+    protected final static double NEAR_ZERO_THRESHOLD = 1e-15;
+
     public final Angle latitude;
     public final Angle longitude;
 
@@ -719,12 +725,18 @@ public static Angle rhumbDistance(LatLon p1, LatLon p2)
         // Taken from http://www.movable-type.co.uk/scripts/latlong.html
         double dLat = lat2 - lat1;
         double dLon = lon2 - lon1;
-        double dPhi = Math.log(Math.tan(lat2 / 2.0 + Math.PI / 4.0) / Math.tan(lat1 / 2.0 + Math.PI / 4.0));
-        double q = dLat / dPhi;
-        if (Double.isNaN(dPhi) || Double.isNaN(q))
+
+        double q;
+        if (Math.abs(dLat) < NEAR_ZERO_THRESHOLD)
         {
             q = Math.cos(lat1);
         }
+        else
+        {
+            double dPhi = Math.log(Math.tan(lat2 / 2.0 + Math.PI / 4.0) / Math.tan(lat1 / 2.0 + Math.PI / 4.0));
+            q = dLat / dPhi;
+        }
+
         // If lonChange over 180 take shorter rhumb across 180 meridian.
         if (Math.abs(dLon) > Math.PI)
         {
@@ -811,13 +823,19 @@ public static LatLon rhumbEndPosition(LatLon p, Angle rhumbAzimuth, Angle pathLe
             return p;
 
         // Taken from http://www.movable-type.co.uk/scripts/latlong.html
-        double lat2 = lat1 + distance * Math.cos(azimuth);
-        double dPhi = Math.log(Math.tan(lat2 / 2.0 + Math.PI / 4.0) / Math.tan(lat1 / 2.0 + Math.PI / 4.0));
-        double q = (lat2 - lat1) / dPhi;
-        if (Double.isNaN(dPhi) || Double.isNaN(q) || Double.isInfinite(q))
+        double dLat = distance * Math.cos(azimuth);
+        double lat2 = lat1 + dLat;
+        double q;
+        if (Math.abs(dLat) < NEAR_ZERO_THRESHOLD)
         {
             q = Math.cos(lat1);
         }
+        else
+        {
+            double dPhi = Math.log(Math.tan(lat2 / 2.0 + Math.PI / 4.0) / Math.tan(lat1 / 2.0 + Math.PI / 4.0));
+            q = (lat2 - lat1) / dPhi;
+        }
+
         double dLon = distance * Math.sin(azimuth) / q;
         // Handle latitude passing over either pole.
         if (Math.abs(lat2) > Math.PI / 2.0)

Diff for: test/gov/nasa/worldwind/geom/LatLonTest.java

@@ -653,8 +653,8 @@ public void testTrivialAzimuthB()
             double azimuthRadians = Math.toRadians(90.0);
             double distanceRadians = Math.toRadians(360.0);
             LatLon end = LatLon.rhumbEndPosition(begin, azimuthRadians, distanceRadians);
-            assertEquals("Trivial Azimuth B (lat)", 2.204291436880279e-14, end.getLatitude().degrees, THRESHOLD);
-            assertEquals("Trivial Azimuth B (lon)", -96.67061650574995, end.getLongitude().degrees,
+            assertEquals("Trivial Azimuth B (lat)", 0.0, end.getLatitude().degrees, THRESHOLD);
+            assertEquals("Trivial Azimuth B (lon)", 0.0, end.getLongitude().degrees,
                 1e-1); // Custom threshold
         }
 
@@ -681,6 +681,24 @@ public void testKnownPointsB()
             assertEquals("Known points B (lat)", 56.9679782407693, end.getLatitude().degrees, THRESHOLD);
             assertEquals("Known points B (lon)", 95.78434282105843, end.getLongitude().degrees, THRESHOLD);
         }
+
+        //////////////////////////////////////////////////////////
+        // Test problem points.
+        //////////////////////////////////////////////////////////
+
+        public void testProblemPointsA()
+        {
+            // This specific lat/lon and distance identified a floating point error
+            Angle initialLat = Angle.fromDegrees(4.076552742498428);
+            Angle initialLon = Angle.fromDegrees(-21.377644877408443);
+            Angle azimuth = Angle.fromDegrees(90.0);
+            Angle distance = Angle.fromDegrees(8.963656110719409);
+            LatLon begin = LatLon.fromRadians(initialLat.getRadians(), initialLon.getRadians());
+            LatLon end = LatLon.rhumbEndPosition(begin, azimuth, distance);
+            assertEquals("Problem points A (lat)", initialLat.getDegrees(), end.getLatitude().getDegrees(), THRESHOLD);
+            assertEquals("Problem points A (lon)", -12.391252821313167, end.getLongitude().getDegrees(), THRESHOLD);
+        }
+
     }
 
     public static class EllipsoidalDistanceTests extends TestCase