Attempted to avoid some extremely rare numerical problems.

The problems were encountered in direct localization. Some line-of-sight
was almost parallel to a latitude boundary, at 1.0e-14 level. This
triggered an exception in the Duvenhage algorithm.

core/src/main/java/org/orekit/rugged/intersection/duvenhage/DuvenhageAlgorithm.java

@@ -257,8 +257,18 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
                 if (longitude >= FastMath.min(entry.getLongitude(), exit.getLongitude()) &&
                     longitude <= FastMath.max(entry.getLongitude(), exit.getLongitude())) {
 
-                    final GeodeticPoint crossingGP;
-                    if (flatBody) {
+                    GeodeticPoint crossingGP = null;
+                    if (!flatBody) {
+                        try {
+                            // full computation of crossing point
+                            final Vector3D crossingP = ellipsoid.pointAtLongitude(position, los, longitude);
+                            crossingGP = ellipsoid.transform(crossingP, ellipsoid.getBodyFrame(), null);
+                        } catch  (RuggedException re) {
+                            // in some very rare cases of numerical noise, we miss the crossing point
+                            crossingGP = null;
+                        }
+                    }
+                    if (crossingGP == null) {
                         // linear approximation of crossing point
                         final double d  = exit.getLongitude() - entry.getLongitude();
                         final double cN = (exit.getLongitude() - longitude) / d;
@@ -266,10 +276,6 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
                         crossingGP = new GeodeticPoint(cN * entry.getLatitude() + cX * exit.getLatitude(),
                                                        longitude,
                                                        cN * entry.getAltitude() + cX * exit.getAltitude());
-                    } else {
-                        // full computation of crossing point
-                        final Vector3D crossingP = ellipsoid.pointAtLongitude(position, los, longitude);
-                        crossingGP = ellipsoid.transform(crossingP, ellipsoid.getBodyFrame(), null);
                     }
                     final int crossingLat = tile.getLatitudeIndex(crossingGP.getLatitude());
 
@@ -305,8 +311,20 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
                 if (latitude >= FastMath.min(entry.getLatitude(), exit.getLatitude()) &&
                     latitude <= FastMath.max(entry.getLatitude(), exit.getLatitude())) {
 
-                    final GeodeticPoint crossingGP;
-                    if (flatBody) {
+                    GeodeticPoint crossingGP = null;
+                    if (!flatBody) {
+                        // full computation of crossing point
+                        try {
+                            final Vector3D crossingP = ellipsoid.pointAtLatitude(position, los,
+                                                                                 tile.getLatitudeAtIndex(crossingLat),
+                                                                                 ellipsoid.transform(entry));
+                            crossingGP = ellipsoid.transform(crossingP, ellipsoid.getBodyFrame(), null);
+                        } catch  (RuggedException re) {
+                            // in some very rare cases of numerical noise, we miss the crossing point
+                            crossingGP = null;
+                        }
+                    }
+                    if (crossingGP == null) {
                         // linear approximation of crossing point
                         final double d  = exit.getLatitude() - entry.getLatitude();
                         final double cN = (exit.getLatitude() - latitude) / d;
@@ -314,12 +332,6 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
                         crossingGP = new GeodeticPoint(latitude,
                                                        cN * entry.getLongitude() + cX * exit.getLongitude(),
                                                        cN * entry.getAltitude()  + cX * exit.getAltitude());
-                    } else {
-                        // full computation of crossing point
-                        final Vector3D crossingP = ellipsoid.pointAtLatitude(position, los,
-                                                                             tile.getLatitudeAtIndex(crossingLat),
-                                                                             ellipsoid.transform(entry));
-                        crossingGP = ellipsoid.transform(crossingP, ellipsoid.getBodyFrame(), null);
                     }
                     final int crossingLon = tile.getLongitudeIndex(crossingGP.getLongitude());
 

core/src/test/java/org/orekit/rugged/utils/ExtendedEllipsoidTest.java

@@ -237,6 +237,25 @@ public class ExtendedEllipsoidTest {
 
     }
 
+    @Test
+    public void testPointAtLatitudeError() throws RuggedException, OrekitException {
+
+        Vector3D p = new Vector3D(-3052690.88784496, 6481300.309857268, 25258.7478104745);
+
+        Vector3D d = new Vector3D(0.3063261631703422, -0.951393802931811, -0.0318451487715828);
+        double latitude = 1.3959699281945737E-14;
+        Vector3D c = new Vector3D(-2809972.5765414005, 5727461.020250551, 26.163518446261833);
+
+        try {
+            ellipsoid.pointAtLatitude(p, d, latitude, c);
+            Assert.fail("an error should have been triggered");
+        } catch (RuggedException re) {
+            Assert.assertEquals(RuggedMessages.LINE_OF_SIGHT_NEVER_CROSSES_LATITUDE, re.getSpecifier());
+            Assert.assertEquals(FastMath.toDegrees(latitude), re.getParts()[0]);
+        }
+
+    }
+
     @Before
     public void setUp() {
         try {

src/site/xdoc/changes.xml

@@ -22,6 +22,10 @@
   <body>
     <release version="1.0" date="TBD"
              description="TBD">
+      <action dev="luc" type="fix" >
+        Added a protection agains some extremely rare numerical problems
+        in Duvenhage algorithm.
+      </action>
       <action dev="luc" type="add" due-to="Silvia Ríos Bergantiños and Beatriz Salazar">
         Added Spanish and Galician translations of error messages.
       </action>