From 542bc4f684beb46809402d94c16a7b6837bfe479 Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@orekit.org>
Date: Mon, 13 Oct 2014 14:04:53 +0200
Subject: [PATCH] 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.
---
.../duvenhage/DuvenhageAlgorithm.java | 40 ++++++++++++-------
.../rugged/utils/ExtendedEllipsoidTest.java | 19 +++++++++
src/site/xdoc/changes.xml | 4 ++
3 files changed, 49 insertions(+), 14 deletions(-)
diff --git a/core/src/main/java/org/orekit/rugged/intersection/duvenhage/DuvenhageAlgorithm.java b/core/src/main/java/org/orekit/rugged/intersection/duvenhage/DuvenhageAlgorithm.java
index ac359e18..d390dab1 100644
--- a/core/src/main/java/org/orekit/rugged/intersection/duvenhage/DuvenhageAlgorithm.java
+++ b/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());
diff --git a/core/src/test/java/org/orekit/rugged/utils/ExtendedEllipsoidTest.java b/core/src/test/java/org/orekit/rugged/utils/ExtendedEllipsoidTest.java
index 7d015237..c098c4ca 100644
--- a/core/src/test/java/org/orekit/rugged/utils/ExtendedEllipsoidTest.java
+++ b/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 {
diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml
index b9b33c67..10a25a72 100644
--- a/src/site/xdoc/changes.xml
+++ b/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>
--
GitLab