From 9c2ad6e1e5e060cb63e96352b45d492a6b08e9d9 Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@orekit.org>
Date: Tue, 20 Sep 2016 19:37:50 +0200
Subject: [PATCH] Improved computation accuracy.
---
.../rugged/refraction/MultiLayerModel.java | 69 +++++++++++++------
.../refraction/MultiLayerModelTest.java | 4 +-
2 files changed, 51 insertions(+), 22 deletions(-)
diff --git a/src/main/java/org/orekit/rugged/refraction/MultiLayerModel.java b/src/main/java/org/orekit/rugged/refraction/MultiLayerModel.java
index b6c3bfc6..7b110ff2 100644
--- a/src/main/java/org/orekit/rugged/refraction/MultiLayerModel.java
+++ b/src/main/java/org/orekit/rugged/refraction/MultiLayerModel.java
@@ -16,6 +16,10 @@
*/
package org.orekit.rugged.refraction;
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.List;
+
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
@@ -26,10 +30,6 @@ import org.orekit.rugged.intersection.IntersectionAlgorithm;
import org.orekit.rugged.utils.ExtendedEllipsoid;
import org.orekit.rugged.utils.NormalizedGeodeticPoint;
-import java.util.ArrayList;
-import java.util.Collections;
-import java.util.List;
-
/**
* Atmospheric refraction model based on multiple layers with associated refractive index.
* @author Sergio Esteves, Luc Maisonobe
@@ -102,7 +102,7 @@ public class MultiLayerModel implements AtmosphericRefraction {
Vector3D pos = satPos;
Vector3D los = satLos.normalize();
- double previousRefractiveIndex = -1;
+ double n1 = -1;
GeodeticPoint gp = ellipsoid.transform(satPos, ellipsoid.getBodyFrame(), null);
for (ConstantRefractionLayer refractionLayer : refractionLayers) {
@@ -111,7 +111,9 @@ public class MultiLayerModel implements AtmosphericRefraction {
continue;
}
- if (previousRefractiveIndex > 0) {
+ final double n2 = refractionLayer.getRefractiveIndex();
+
+ if (n1 > 0) {
// when we get here, we have already performed one iteration in the loop
// so gp is the los intersection with the layers interface (it was a
@@ -146,24 +148,51 @@ public class MultiLayerModel implements AtmosphericRefraction {
// = sin² θ₂ / sin² θâ‚
// as α is positive, and both θ₠and θ₂ are between 0 and π/2, we finally get
// α = sin θ₂ / sin θâ‚
- // (5) α = n₠/ n₂
+ // (5) α = nâ‚/nâ‚‚
// the α coefficient is independent from the incidence angle,
// it depends only on the ratio of refractive indices!
//
// back to equation (4) and using again the fact θ₂ is between 0 and π/2, we can now write
// β = α cos θ₠- cos θ₂
- // = n₠/ n₂ cos θ₠- cos θ₂
- // = n₠/ n₂ cos θ₠- √(1 - sin² θ₂)
- // = nâ‚ / nâ‚‚ cos θ₠- √(1 - (nâ‚ / nâ‚‚)² sin² θâ‚)
- // = nâ‚ / nâ‚‚ cos θ₠- √(1 - (nâ‚ / nâ‚‚)² (1 - cos² θâ‚))
- // = n₠/ n₂ cos θ₠- √(1 + (n₠/ n₂)² cos² θ₠- (n₠/ n₂)²)
- // (6) β = -k - √(1 + k² - (n₠/ n₂)²)
- // where k = n₠/ n₂ u·z
- final double n1On2 = previousRefractiveIndex / refractionLayer.getRefractiveIndex();
- final double k = n1On2 * Vector3D.dotProduct(los, gp.getZenith());
- los = new Vector3D(n1On2, los,
- -k - FastMath.sqrt(1 + k * k - n1On2 * n1On2), gp.getZenith());
-
+ // = nâ‚/nâ‚‚ cos θ₠- cos θ₂
+ // = nâ‚/nâ‚‚ cos θ₠- √(1 - sin² θ₂)
+ // = nâ‚/nâ‚‚ cos θ₠- √(1 - (nâ‚/nâ‚‚)² sin² θâ‚)
+ // = nâ‚/nâ‚‚ cos θ₠- √(1 - (nâ‚/nâ‚‚)² (1 - cos² θâ‚))
+ // = nâ‚/nâ‚‚ cos θ₠- √(1 - (nâ‚/nâ‚‚)² + (nâ‚/nâ‚‚)² cos² θâ‚)
+ // (6) β = -k - √(k² - ζ)
+ // where ζ = (nâ‚/nâ‚‚)² - 1 and k = nâ‚/nâ‚‚ u·z, which is negative, and close to -1 for
+ // nadir observations. As we expect atmosphere models to have small transitions between
+ // layers, we have to handle accurately the case where nâ‚/nâ‚‚ ≈ 1 so ζ ≈ 0. In this case,
+ // a cancellation occurs inside the square root: √(k² - ζ) ≈ √k² ≈ -k (because k is negative).
+ // So β ≈ -k + k ≈ 0 and another cancellation occurs, outside of the square root.
+ // This means that despite equation (6) is mathematically correct, it is prone to numerical
+ // errors when consecutive layers have close refractive indices. A better equivalent
+ // expression is needed. The fact β is close to 0 in this case was expected because
+ // equation (1) reads v = α u + β z, and α = nâ‚/nâ‚‚, so when nâ‚/nâ‚‚ ≈ 1, we have
+ // α ≈ 1 and β ≈ 0, so v ≈ u: when two layers have similar refractive indices, the
+ // propagation direction is almost unchanged.
+ //
+ // The first step for the accurate computation of β is to compute ζ = (nâ‚/nâ‚‚)² - 1
+ // accurately and avoid a cancellation just after a division (which is the least accurate
+ // of the four operations) and a squaring. We will simply use:
+ // ζ = (nâ‚/nâ‚‚)² - 1
+ // = (n₠- n₂) (n₠+ n₂) / n₂²
+ // The cancellation is still there, but it occurs in the subtraction nâ‚ - nâ‚‚, which are
+ // the most accurate values we can get.
+ // The second step for the accurate computation of β is to rewrite equation (6)
+ // by both multiplying and dividing by the dual factor -k + √(k² - ζ):
+ // β = -k - √(k² - ζ)
+ // = (-k - √(k² - ζ)) * (-k + √(k² - ζ)) / (-k + √(k² - ζ))
+ // = (k² - (k² - ζ)) / (-k + √(k² - ζ))
+ // (7) β = ζ / (-k + √(k² - ζ))
+ // expression (7) is more stable numerically than expression (6), because when ζ ≈ 0
+ // its denominator is about -2k, there are no cancellation anymore after the square root.
+ // β is computed with the same accuracy as ζ
+ final double alpha = n1 / n2;
+ final double k = alpha * Vector3D.dotProduct(los, gp.getZenith());
+ final double zeta = (n1 - n2) * (n1 + n2) / (n2 * n2);
+ final double beta = zeta / (FastMath.sqrt(k * k - zeta) - k);
+ los = new Vector3D(alpha, los, beta, gp.getZenith());
}
if (rawIntersection.getAltitude() > refractionLayer.getLowestAltitude()) {
@@ -174,7 +203,7 @@ public class MultiLayerModel implements AtmosphericRefraction {
pos = ellipsoid.pointAtAltitude(pos, los, refractionLayer.getLowestAltitude());
gp = ellipsoid.transform(pos, ellipsoid.getBodyFrame(), null);
- previousRefractiveIndex = refractionLayer.getRefractiveIndex();
+ n1 = n2;
}
diff --git a/src/test/java/org/orekit/rugged/refraction/MultiLayerModelTest.java b/src/test/java/org/orekit/rugged/refraction/MultiLayerModelTest.java
index 27c595f6..6df709c1 100644
--- a/src/test/java/org/orekit/rugged/refraction/MultiLayerModelTest.java
+++ b/src/test/java/org/orekit/rugged/refraction/MultiLayerModelTest.java
@@ -56,7 +56,7 @@ public class MultiLayerModelTest extends AbstractAlgorithmTest {
// this is almost a Nadir observation (LOS deviates between 1.4 and 1.6 degrees from vertical)
// so the refraction correction is small
- Assert.assertEquals(0.0553797, distance, 1.0e-6);
+ Assert.assertEquals(0.0553796, distance, 1.0e-6);
}
@@ -83,7 +83,7 @@ public class MultiLayerModelTest extends AbstractAlgorithmTest {
model = new MultiLayerModel(earth, refractionLayers);
correctedIntersection = model.applyCorrection(position, los, rawIntersection, algorithm);
distance = Vector3D.distance(earth.transform(rawIntersection), earth.transform(correctedIntersection));
- Assert.assertEquals(0.0, distance, 1.0e-9);
+ Assert.assertEquals(0.0, distance, 1.0e-20);
}
--
GitLab