diff --git a/atmospheric-deviation.csv b/atmospheric-deviation.csv
deleted file mode 100644
index 2807245d5cfeb5959c8dc80a22722c2b6f87614a..0000000000000000000000000000000000000000
--- a/atmospheric-deviation.csv
+++ /dev/null
@@ -1,41 +0,0 @@
-angle,correction
-0.0,0
-0.01,0
-0.02,0
-0.03,0
-0.04,0
-0.05,0
-0.060000000000000005,0
-0.07,0
-0.08,0
-0.09,0
-0.09999999999999999,0
-0.10999999999999999,0
-0.11999999999999998,0
-0.12999999999999998,0
-0.13999999999999999,0
-0.15,0
-0.16,1
-0.17,1
-0.18000000000000002,1
-0.19000000000000003,1
-0.20000000000000004,1
-0.21000000000000005,1
-0.22000000000000006,1
-0.23000000000000007,1
-0.24000000000000007,1
-0.25000000000000006,1
-0.26000000000000006,1
-0.2700000000000001,1
-0.2800000000000001,1
-0.2900000000000001,1
-0.3000000000000001,1
-0.3100000000000001,1
-0.3200000000000001,1
-0.3300000000000001,1
-0.34000000000000014,1
-0.35000000000000014,1
-0.36000000000000015,1
-0.37000000000000016,1
-0.38000000000000017,2
-0.3900000000000002,2
diff --git a/atmospheric-deviation.gp b/atmospheric-deviation.gp
deleted file mode 100644
index 88716c6e85efa2af4d61f2ddb0aaf927e17c5fee..0000000000000000000000000000000000000000
--- a/atmospheric-deviation.gp
+++ /dev/null
@@ -1,13 +0,0 @@
-# gnuplot script file for plotting bandwidth over time
-#!/usr/bin/gnuplot
-
-
-set terminal pdf size 3,2
-set output 'atmospheric-deviation.pdf'
-set grid
-
-set title "Atmospheric Correction"
-set xlabel "angle (radians)"
-set ylabel "correction (meters)"
-
-plot "atmospheric-deviation.csv" using 1:0 with lines title ""
diff --git a/atmospheric-deviation.pdf b/atmospheric-deviation.pdf
deleted file mode 100644
index 1ce2e9bf36894482b4061c28bf544f2d59192897..0000000000000000000000000000000000000000
Binary files a/atmospheric-deviation.pdf and /dev/null differ
diff --git a/src/test/java/org/orekit/rugged/refraction/MultiLayerModelTest.java b/src/test/java/org/orekit/rugged/refraction/MultiLayerModelTest.java
index 409a6a7180c494ebcc8968d9f27976da4523689b..27c595f6ff9329e7659ee35df7bbec4901ec26a0 100644
--- a/src/test/java/org/orekit/rugged/refraction/MultiLayerModelTest.java
+++ b/src/test/java/org/orekit/rugged/refraction/MultiLayerModelTest.java
@@ -16,46 +16,63 @@
*/
package org.orekit.rugged.refraction;
+import org.hipparchus.analysis.UnivariateFunction;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.RotationConvention;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
-import org.orekit.bodies.GeodeticPoint;
import org.orekit.errors.OrekitException;
import org.orekit.rugged.errors.RuggedException;
+import org.orekit.rugged.errors.RuggedMessages;
import org.orekit.rugged.intersection.AbstractAlgorithmTest;
import org.orekit.rugged.intersection.IntersectionAlgorithm;
import org.orekit.rugged.intersection.duvenhage.DuvenhageAlgorithm;
import org.orekit.rugged.raster.TileUpdater;
import org.orekit.rugged.utils.NormalizedGeodeticPoint;
-import java.io.FileNotFoundException;
-import java.io.PrintWriter;
+import java.lang.reflect.Field;
import java.util.ArrayList;
+import java.util.Collections;
import java.util.List;
public class MultiLayerModelTest extends AbstractAlgorithmTest {
@Test
- public void testApplyCorrection() throws OrekitException, RuggedException, FileNotFoundException {
+ public void testAlmostNadir() throws OrekitException, RuggedException {
setUpMayonVolcanoContext();
final IntersectionAlgorithm algorithm = createAlgorithm(updater, 8);
final Vector3D position = new Vector3D(-3787079.6453602533, 5856784.405679551, 1655869.0582939098);
final Vector3D los = new Vector3D( 0.5127552821932051, -0.8254313129088879, -0.2361041470463311);
- final NormalizedGeodeticPoint rawIntersection = algorithm.refineIntersection(earth, position, los,
- algorithm.intersection(earth, position, los));
+ final NormalizedGeodeticPoint rawIntersection =
+ algorithm.refineIntersection(earth, position, los,
+ algorithm.intersection(earth, position, los));
MultiLayerModel model = new MultiLayerModel(earth);
- GeodeticPoint correctedIntersection = model.applyCorrection(position, los, rawIntersection, algorithm);
+ NormalizedGeodeticPoint correctedIntersection = model.applyCorrection(position, los, rawIntersection, algorithm);
double distance = Vector3D.distance(earth.transform(rawIntersection), earth.transform(correctedIntersection));
// 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);
+ }
+
+ @Test
+ public void testNoOpRefraction() throws OrekitException, RuggedException {
+
+ setUpMayonVolcanoContext();
+ final IntersectionAlgorithm algorithm = createAlgorithm(updater, 8);
+ final Vector3D position = new Vector3D(-3787079.6453602533, 5856784.405679551, 1655869.0582939098);
+ final Vector3D los = los(position, FastMath.toRadians(50.0));
+ final NormalizedGeodeticPoint rawIntersection = algorithm.refineIntersection(earth, position, los,
+ algorithm.intersection(earth, position, los));
+
+ MultiLayerModel model = new MultiLayerModel(earth);
+ NormalizedGeodeticPoint correctedIntersection = model.applyCorrection(position, los, rawIntersection, algorithm);
+ double distance = Vector3D.distance(earth.transform(rawIntersection), earth.transform(correctedIntersection));
// a test with indices all set to 1.0 - correction must be zero
final int numberOfLayers = 16;
@@ -66,78 +83,185 @@ 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, 0.001);
+ Assert.assertEquals(0.0, distance, 1.0e-9);
+ }
+
+ @Test
+ public void testReversedAtmosphere()
+ throws OrekitException, RuggedException, IllegalArgumentException,
+ IllegalAccessException, NoSuchFieldException, SecurityException {
+
+ setUpMayonVolcanoContext();
+ final IntersectionAlgorithm algorithm = createAlgorithm(updater, 8);
+ final Vector3D position = new Vector3D(-3787079.6453602533, 5856784.405679551, 1655869.0582939098);
+ final Vector3D los = los(position, FastMath.toRadians(50.0));
+ final NormalizedGeodeticPoint rawIntersection = algorithm.refineIntersection(earth, position, los,
+ algorithm.intersection(earth, position, los));
+
+ MultiLayerModel baseModel = new MultiLayerModel(earth);
+ NormalizedGeodeticPoint correctedIntersection = baseModel.applyCorrection(position, los, rawIntersection, algorithm);
// an intentionally flawed atmosphere with refractive indices decreasing with altitude,
// that should exhibit a LOS bending upwards
- refractionLayers = new ArrayList<ConstantRefractionLayer>(numberOfLayers);
- for(int i = numberOfLayers - 1; i >= 0; i--) {
- refractionLayers.add(new ConstantRefractionLayer(i * 1.0e4, 1.0 + i*i*1e-6));
+ Field refractionLayersField = MultiLayerModel.class.getDeclaredField("refractionLayers");
+ refractionLayersField.setAccessible(true);
+ @SuppressWarnings("unchecked")
+ List<ConstantRefractionLayer> baseRefractionLayers =
+ (List<ConstantRefractionLayer>) refractionLayersField.get(baseModel);
+ List<ConstantRefractionLayer> denserRefractionLayers = new ArrayList<>();
+ for (final ConstantRefractionLayer layer : baseRefractionLayers) {
+ denserRefractionLayers.add(new ConstantRefractionLayer(layer.getLowestAltitude(),
+ 1.0 / layer.getRefractiveIndex()));
}
- model = new MultiLayerModel(earth, refractionLayers);
- correctedIntersection = model.applyCorrection(position, los, rawIntersection, algorithm);
- double anglePosRawIntersection = Vector3D.angle(position, earth.transform(rawIntersection));
+ MultiLayerModel reversedModel = new MultiLayerModel(earth, denserRefractionLayers);
+ NormalizedGeodeticPoint reversedIntersection = reversedModel.applyCorrection(position, los, rawIntersection, algorithm);
+ double anglePosRawIntersection = Vector3D.angle(position, earth.transform(rawIntersection));
double anglePosCorrectedIntersection = Vector3D.angle(position, earth.transform(correctedIntersection));
- Assert.assertTrue(anglePosRawIntersection < anglePosCorrectedIntersection);
+ double anglePosReversedIntersection = Vector3D.angle(position, earth.transform(reversedIntersection));
+
+ // with regular atmosphere, the ray bends downwards,
+ // so the ground point is closer to the sub-satellite point than the raw intersection
+ Assert.assertTrue(anglePosCorrectedIntersection < anglePosRawIntersection);
+
+ // with reversed atmosphere, the ray bends upwards,
+ // so the ground point is farther from the sub-satellite point than the raw intersection
+ Assert.assertTrue(anglePosReversedIntersection > anglePosRawIntersection);
+ // the points are almost aligned (for distances around 20m, Earth curvature is small enough)
+ double dRawCorrected = Vector3D.distance(earth.transform(rawIntersection), earth.transform(correctedIntersection));
+ double dRawReversed = Vector3D.distance(earth.transform(rawIntersection), earth.transform(reversedIntersection));
+ double dReversedCorrected = Vector3D.distance(earth.transform(reversedIntersection), earth.transform(correctedIntersection));
+ Assert.assertEquals(dRawCorrected + dRawReversed, dReversedCorrected, 1.0e-12 * dReversedCorrected);
+
+ }
+
+ @Test
+ public void testTwoAtmospheres()
+ throws OrekitException, RuggedException, IllegalArgumentException,
+ IllegalAccessException, NoSuchFieldException, SecurityException {
+
+ setUpMayonVolcanoContext();
+ final IntersectionAlgorithm algorithm = createAlgorithm(updater, 8);
+ final Vector3D position = new Vector3D(-3787079.6453602533, 5856784.405679551, 1655869.0582939098);
+ final Vector3D los = los(position, FastMath.toRadians(50.0));
+ final NormalizedGeodeticPoint rawIntersection = algorithm.refineIntersection(earth, position, los,
+ algorithm.intersection(earth, position, los));
// a comparison between two atmospheres, one more dense than the other and showing correction
// is more important with high indices
- List<ConstantRefractionLayer> baseRefracLayers = new ArrayList<ConstantRefractionLayer>(numberOfLayers);
- List<ConstantRefractionLayer> denserRefracLayers = new ArrayList<ConstantRefractionLayer>(numberOfLayers);
- for(int i = numberOfLayers - 1; i >= 0; i--) {
- double baseRefractiveIndex = FastMath.pow(numberOfLayers - i + 1, 2) * 1e-6;
- baseRefracLayers.add(new ConstantRefractionLayer(i * 1.0e4, baseRefractiveIndex));
- denserRefracLayers.add(new ConstantRefractionLayer(i * 1.0e4, baseRefractiveIndex + 1e-3));
+ MultiLayerModel baseModel = new MultiLayerModel(earth);
+ Field refractionLayersField = MultiLayerModel.class.getDeclaredField("refractionLayers");
+ refractionLayersField.setAccessible(true);
+ @SuppressWarnings("unchecked")
+ List<ConstantRefractionLayer> baseRefractionLayers =
+ (List<ConstantRefractionLayer>) refractionLayersField.get(baseModel);
+ List<ConstantRefractionLayer> denserRefractionLayers = new ArrayList<>();
+ double previousBaseN = 1.0;
+ double previousDenserN = 1.0;
+ double factor = 1.00001;
+ for (final ConstantRefractionLayer layer : baseRefractionLayers) {
+ final double currentBaseN = layer.getRefractiveIndex();
+ final double baseRatio = currentBaseN / previousBaseN;
+ final double currentDenserN = previousDenserN * factor * baseRatio;
+ denserRefractionLayers.add(new ConstantRefractionLayer(layer.getLowestAltitude(),
+ currentDenserN));
+ previousBaseN = currentBaseN;
+ previousDenserN = currentDenserN;
}
- MultiLayerModel baseModel = new MultiLayerModel(earth, baseRefracLayers);
- MultiLayerModel denserModel = new MultiLayerModel(earth, denserRefracLayers);
- GeodeticPoint baseIntersection = baseModel.applyCorrection(position, los, rawIntersection, algorithm);
- GeodeticPoint denserIntersection = denserModel.applyCorrection(position, los, rawIntersection, algorithm);
- double baseDistance = Vector3D.distance(earth.transform(rawIntersection), earth.transform(baseIntersection));
+ MultiLayerModel denserModel = new MultiLayerModel(earth, denserRefractionLayers);
+ NormalizedGeodeticPoint baseIntersection = baseModel.applyCorrection(position, los, rawIntersection, algorithm);
+ NormalizedGeodeticPoint denserIntersection = denserModel.applyCorrection(position, los, rawIntersection, algorithm);
+ double baseDistance = Vector3D.distance(earth.transform(rawIntersection),
+ earth.transform(baseIntersection));
double denserDistance = Vector3D.distance(earth.transform(rawIntersection),
- earth.transform(denserIntersection));
- // denserDistance: 291.6042252928431, baseDistance: 2710.1036961651967
- // Assert.assertTrue(denserDistance > baseDistance);
+ earth.transform(denserIntersection));
+ Assert.assertTrue(denserDistance > baseDistance);
+ }
- // a test with a single refraction layer
- refractionLayers = new ArrayList<ConstantRefractionLayer>(numberOfLayers);
- refractionLayers.add(new ConstantRefractionLayer(0, 1.2));
- 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, 0.0);
+ @Test
+ public void testMissingLayers() throws OrekitException, RuggedException {
+ setUpMayonVolcanoContext();
+ final IntersectionAlgorithm algorithm = createAlgorithm(updater, 8);
+ final Vector3D position = new Vector3D(-3787079.6453602533, 5856784.405679551, 1655869.0582939098);
+ final Vector3D los = los(position, FastMath.toRadians(50.0));
+ final NormalizedGeodeticPoint rawIntersection = algorithm.refineIntersection(earth, position, los,
+ algorithm.intersection(earth, position, los));
+ final double h = rawIntersection.getAltitude();
- // deviation should increase as the angle between los and zenith increases
- PrintWriter writer = new PrintWriter("atmospheric-deviation.csv");
- writer.println("angle,correction");
+ MultiLayerModel model = new MultiLayerModel(earth,
+ Collections.singletonList(new ConstantRefractionLayer(h + 100.0,
+ 1.5)));
+ try {
+ model.applyCorrection(position, los, rawIntersection, algorithm);
+ Assert.fail("an exception should have been thrown");
+ } catch (RuggedException re) {
+ Assert.assertEquals(RuggedMessages.NO_LAYER_DATA, re.getSpecifier());
+ Assert.assertEquals(h, ((Double) re.getParts()[0]).doubleValue(), 1.0e-6);
+ Assert.assertEquals(h + 100.0, ((Double) re.getParts()[1]).doubleValue(), 1.0e-6);
+ }
- GeodeticPoint satGP = earth.transform(position, earth.getBodyFrame(), null);
- Vector3D nadir = satGP.getNadir();
- Vector3D horizontal = nadir.orthogonal();
- for (double alpha = 0; alpha < 0.4; alpha += 0.01) {
- Vector3D rotatingLos = new Rotation(horizontal, alpha, RotationConvention.VECTOR_OPERATOR).applyTo(nadir);
- NormalizedGeodeticPoint uncorrectedIntersection = algorithm.refineIntersection(earth, position, rotatingLos,
- algorithm.intersection(earth, position, rotatingLos));
+ }
- model = new MultiLayerModel(earth);
- correctedIntersection = model.applyCorrection(position, rotatingLos, uncorrectedIntersection, algorithm);
- distance = Vector3D.distance(earth.transform(uncorrectedIntersection),
- earth.transform(correctedIntersection));
+ @Test
+ public void testLayersBelowDEM() throws OrekitException, RuggedException {
- writer.println(alpha + "," + FastMath.round(distance));
- }
- writer.close();
+ setUpMayonVolcanoContext();
+ final IntersectionAlgorithm algorithm = createAlgorithm(updater, 8);
+ final Vector3D position = new Vector3D(-3787079.6453602533, 5856784.405679551, 1655869.0582939098);
+ final Vector3D los = los(position, FastMath.toRadians(50.0));
+ final NormalizedGeodeticPoint rawIntersection = algorithm.refineIntersection(earth, position, los,
+ algorithm.intersection(earth, position, los));
+ MultiLayerModel model = new MultiLayerModel(earth,
+ Collections.singletonList(new ConstantRefractionLayer(rawIntersection.getAltitude() - 100.0,
+ 1.5)));
+ NormalizedGeodeticPoint correctedIntersection = model.applyCorrection(position, los, rawIntersection, algorithm);
+ double distance = Vector3D.distance(earth.transform(rawIntersection), earth.transform(correctedIntersection));
+ Assert.assertEquals(0.0, distance, 1.0e-20);
+ }
+
+ @Test
+ public void testDivingAngleChange() throws OrekitException, RuggedException {
+
+ setUpMayonVolcanoContext();
+ final IntersectionAlgorithm algorithm = createAlgorithm(updater, 8);
+ final Vector3D position = new Vector3D(-3787079.6453602533, 5856784.405679551, 1655869.0582939098);
+ AtmosphericRefraction model = new MultiLayerModel(earth);
+
+ // deviation should increase from 0 to about 17m
+ // as the angle between los and nadir increases from 0 to 50 degrees
+ // the reference model below has been obtained by fitting the test results themselves
+ // it is NOT considered a full featured model, it's just a helper function for this specific test
+ UnivariateFunction reference = alpha -> 1.17936 * FastMath.tan((2.94613 - 1.40162 * alpha) * alpha);
+
+ for (double alpha = 0; alpha < FastMath.toRadians(50.0); alpha += 0.1) {
+ final Vector3D rotatingLos = los(position, alpha);
+ final NormalizedGeodeticPoint rawIntersection = algorithm.refineIntersection(earth, position, rotatingLos,
+ algorithm.intersection(earth, position, rotatingLos));
+
+ final NormalizedGeodeticPoint correctedIntersection = model.applyCorrection(position, rotatingLos, rawIntersection, algorithm);
+ final double distance = Vector3D.distance(earth.transform(rawIntersection),
+ earth.transform(correctedIntersection));
+ Assert.assertEquals(reference.value(alpha), distance, 0.12);
+ }
}
+ private Vector3D los(final Vector3D position, final double angleFromNadir)
+ throws OrekitException {
+ final Vector3D nadir = earth.transform(position, earth.getBodyFrame(), null).getNadir();
+ final Rotation losRotation = new Rotation(nadir.orthogonal(), angleFromNadir,
+ RotationConvention.VECTOR_OPERATOR);
+ return losRotation.applyTo(nadir);
+ }
+
@Override
protected IntersectionAlgorithm createAlgorithm(TileUpdater updater, int maxCachedTiles) {
return new DuvenhageAlgorithm(updater, maxCachedTiles, false);
}
+
}