Skip to content
Snippets Groups Projects
Commit f4c3d770 authored by Luc Maisonobe's avatar Luc Maisonobe
Browse files

Improved tests.

parent f26bab6a
No related branches found
No related tags found
No related merge requests found
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
# 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 ""
File deleted
......@@ -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);
}
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment