1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.rugged.refraction;
18
19 import java.lang.reflect.Field;
20 import java.util.ArrayList;
21 import java.util.Collections;
22 import java.util.List;
23
24 import org.hipparchus.analysis.UnivariateFunction;
25 import org.hipparchus.geometry.euclidean.threed.Rotation;
26 import org.hipparchus.geometry.euclidean.threed.RotationConvention;
27 import org.hipparchus.geometry.euclidean.threed.Vector3D;
28 import org.hipparchus.util.FastMath;
29 import org.junit.jupiter.api.Assertions;
30 import org.junit.jupiter.api.Test;
31 import org.orekit.rugged.errors.RuggedException;
32 import org.orekit.rugged.errors.RuggedMessages;
33 import org.orekit.rugged.intersection.AbstractAlgorithmTest;
34 import org.orekit.rugged.intersection.IntersectionAlgorithm;
35 import org.orekit.rugged.intersection.duvenhage.DuvenhageAlgorithm;
36 import org.orekit.rugged.raster.TileUpdater;
37 import org.orekit.rugged.utils.NormalizedGeodeticPoint;
38
39 public class MultiLayerModelTest extends AbstractAlgorithmTest {
40
41 @Test
42 public void testAlmostNadir() {
43
44 setUpMayonVolcanoContext();
45
46 final IntersectionAlgorithm algorithm = createAlgorithm(updater, 8, true);
47 final Vector3D position = new Vector3D(-3787079.6453602533, 5856784.405679551, 1655869.0582939098);
48 final Vector3D los = new Vector3D( 0.5127552821932051, -0.8254313129088879, -0.2361041470463311);
49 final NormalizedGeodeticPoint rawIntersection =
50 algorithm.refineIntersection(earth, position, los,
51 algorithm.intersection(earth, position, los));
52
53 MultiLayerModel model = new MultiLayerModel(earth);
54 NormalizedGeodeticPoint correctedIntersection = model.applyCorrection(position, los, rawIntersection, algorithm);
55 double distance = Vector3D.distance(earth.transform(rawIntersection), earth.transform(correctedIntersection));
56
57
58
59 Assertions.assertEquals(0.0553796, distance, 1.0e-6);
60
61 }
62
63 @Test
64 public void testNoOpRefraction() {
65
66 setUpMayonVolcanoContext();
67
68 final IntersectionAlgorithm algorithm = createAlgorithm(updater, 8, true);
69 final Vector3D position = new Vector3D(-3787079.6453602533, 5856784.405679551, 1655869.0582939098);
70 final Vector3D los = los(position, FastMath.toRadians(50.0));
71 final NormalizedGeodeticPoint rawIntersection = algorithm.refineIntersection(earth, position, los,
72 algorithm.intersection(earth, position, los));
73
74 MultiLayerModel model = new MultiLayerModel(earth);
75 NormalizedGeodeticPoint correctedIntersection = model.applyCorrection(position, los, rawIntersection, algorithm);
76 double distance = Vector3D.distance(earth.transform(rawIntersection), earth.transform(correctedIntersection));
77
78
79 final int numberOfLayers = 16;
80 List<ConstantRefractionLayer> refractionLayers = new ArrayList<ConstantRefractionLayer>(numberOfLayers);
81 for(int i = numberOfLayers - 1; i >= 0; i--) {
82 refractionLayers.add(new ConstantRefractionLayer(i * 1.0e4, 1.0));
83 }
84 model = new MultiLayerModel(earth, refractionLayers);
85 correctedIntersection = model.applyCorrection(position, los, rawIntersection, algorithm);
86 distance = Vector3D.distance(earth.transform(rawIntersection), earth.transform(correctedIntersection));
87 Assertions.assertEquals(0.0, distance, 1.7e-9);
88
89 }
90
91 @Test
92 public void testReversedAtmosphere()
93 throws IllegalArgumentException, IllegalAccessException, NoSuchFieldException, SecurityException {
94
95 setUpMayonVolcanoContext();
96
97 final IntersectionAlgorithm algorithm = createAlgorithm(updater, 8, true);
98 final Vector3D position = new Vector3D(-3787079.6453602533, 5856784.405679551, 1655869.0582939098);
99 final Vector3D los = los(position, FastMath.toRadians(50.0));
100 final NormalizedGeodeticPoint rawIntersection = algorithm.refineIntersection(earth, position, los,
101 algorithm.intersection(earth, position, los));
102
103 MultiLayerModel baseModel = new MultiLayerModel(earth);
104 NormalizedGeodeticPoint correctedIntersection = baseModel.applyCorrection(position, los, rawIntersection, algorithm);
105
106
107
108 Field refractionLayersField = MultiLayerModel.class.getDeclaredField("refractionLayers");
109 refractionLayersField.setAccessible(true);
110 @SuppressWarnings("unchecked")
111 List<ConstantRefractionLayer> baseRefractionLayers =
112 (List<ConstantRefractionLayer>) refractionLayersField.get(baseModel);
113 List<ConstantRefractionLayer> denserRefractionLayers = new ArrayList<>();
114 for (final ConstantRefractionLayer layer : baseRefractionLayers) {
115 denserRefractionLayers.add(new ConstantRefractionLayer(layer.getLowestAltitude(),
116 1.0 / layer.getRefractiveIndex()));
117 }
118 MultiLayerModel reversedModel = new MultiLayerModel(earth, denserRefractionLayers);
119 NormalizedGeodeticPoint reversedIntersection = reversedModel.applyCorrection(position, los, rawIntersection, algorithm);
120 double anglePosRawIntersection = Vector3D.angle(position, earth.transform(rawIntersection));
121 double anglePosCorrectedIntersection = Vector3D.angle(position, earth.transform(correctedIntersection));
122 double anglePosReversedIntersection = Vector3D.angle(position, earth.transform(reversedIntersection));
123
124
125
126 Assertions.assertTrue(anglePosCorrectedIntersection < anglePosRawIntersection);
127
128
129
130 Assertions.assertTrue(anglePosReversedIntersection > anglePosRawIntersection);
131
132
133 double dRawCorrected = Vector3D.distance(earth.transform(rawIntersection), earth.transform(correctedIntersection));
134 double dRawReversed = Vector3D.distance(earth.transform(rawIntersection), earth.transform(reversedIntersection));
135 double dReversedCorrected = Vector3D.distance(earth.transform(reversedIntersection), earth.transform(correctedIntersection));
136 Assertions.assertEquals(dRawCorrected + dRawReversed, dReversedCorrected, 1.0e-12 * dReversedCorrected);
137
138 }
139
140 @Test
141 public void testTwoAtmospheres()
142 throws IllegalArgumentException, IllegalAccessException, NoSuchFieldException, SecurityException {
143
144 setUpMayonVolcanoContext();
145
146 final IntersectionAlgorithm algorithm = createAlgorithm(updater, 8, true);
147 final Vector3D position = new Vector3D(-3787079.6453602533, 5856784.405679551, 1655869.0582939098);
148 final Vector3D los = los(position, FastMath.toRadians(50.0));
149 final NormalizedGeodeticPoint rawIntersection = algorithm.refineIntersection(earth, position, los,
150 algorithm.intersection(earth, position, los));
151
152
153
154 MultiLayerModel baseModel = new MultiLayerModel(earth);
155 Field refractionLayersField = MultiLayerModel.class.getDeclaredField("refractionLayers");
156 refractionLayersField.setAccessible(true);
157 @SuppressWarnings("unchecked")
158 List<ConstantRefractionLayer> baseRefractionLayers =
159 (List<ConstantRefractionLayer>) refractionLayersField.get(baseModel);
160 List<ConstantRefractionLayer> denserRefractionLayers = new ArrayList<>();
161 double previousBaseN = 1.0;
162 double previousDenserN = 1.0;
163 double factor = 1.00001;
164 for (final ConstantRefractionLayer layer : baseRefractionLayers) {
165 final double currentBaseN = layer.getRefractiveIndex();
166 final double baseRatio = currentBaseN / previousBaseN;
167 final double currentDenserN = previousDenserN * factor * baseRatio;
168 denserRefractionLayers.add(new ConstantRefractionLayer(layer.getLowestAltitude(),
169 currentDenserN));
170 previousBaseN = currentBaseN;
171 previousDenserN = currentDenserN;
172 }
173 MultiLayerModel denserModel = new MultiLayerModel(earth, denserRefractionLayers);
174 NormalizedGeodeticPoint baseIntersection = baseModel.applyCorrection(position, los, rawIntersection, algorithm);
175 NormalizedGeodeticPoint denserIntersection = denserModel.applyCorrection(position, los, rawIntersection, algorithm);
176 double baseDistance = Vector3D.distance(earth.transform(rawIntersection),
177 earth.transform(baseIntersection));
178 double denserDistance = Vector3D.distance(earth.transform(rawIntersection),
179 earth.transform(denserIntersection));
180 Assertions.assertTrue(denserDistance > baseDistance);
181
182 }
183
184 @Test
185 public void testMissingLayers() {
186
187 setUpMayonVolcanoContext();
188
189 final IntersectionAlgorithm algorithm = createAlgorithm(updater, 8, true);
190 final Vector3D position = new Vector3D(-3787079.6453602533, 5856784.405679551, 1655869.0582939098);
191 final Vector3D los = los(position, FastMath.toRadians(50.0));
192 final NormalizedGeodeticPoint rawIntersection = algorithm.refineIntersection(earth, position, los,
193 algorithm.intersection(earth, position, los));
194 final double h = rawIntersection.getAltitude();
195
196 MultiLayerModel model = new MultiLayerModel(earth,
197 Collections.singletonList(new ConstantRefractionLayer(h + 100.0,
198 1.5)));
199 try {
200 model.applyCorrection(position, los, rawIntersection, algorithm);
201 Assertions.fail("an exception should have been thrown");
202 } catch (RuggedException re) {
203 Assertions.assertEquals(RuggedMessages.NO_LAYER_DATA, re.getSpecifier());
204 Assertions.assertEquals(h, ((Double) re.getParts()[0]).doubleValue(), 1.0e-6);
205 Assertions.assertEquals(h + 100.0, ((Double) re.getParts()[1]).doubleValue(), 1.0e-6);
206 }
207
208 }
209
210 @Test
211 public void testLayersBelowDEM() {
212
213 setUpMayonVolcanoContext();
214
215 final IntersectionAlgorithm algorithm = createAlgorithm(updater, 8, true);
216 final Vector3D position = new Vector3D(-3787079.6453602533, 5856784.405679551, 1655869.0582939098);
217 final Vector3D los = los(position, FastMath.toRadians(50.0));
218 final NormalizedGeodeticPoint rawIntersection = algorithm.refineIntersection(earth, position, los,
219 algorithm.intersection(earth, position, los));
220
221 MultiLayerModel model = new MultiLayerModel(earth,
222 Collections.singletonList(new ConstantRefractionLayer(rawIntersection.getAltitude() - 100.0,
223 1.5)));
224 NormalizedGeodeticPoint correctedIntersection = model.applyCorrection(position, los, rawIntersection, algorithm);
225 double distance = Vector3D.distance(earth.transform(rawIntersection), earth.transform(correctedIntersection));
226 Assertions.assertEquals(0.0, distance, 1.3e-9);
227
228 }
229
230 @Test
231 public void testDivingAngleChange() {
232
233 setUpMayonVolcanoContext();
234
235 final IntersectionAlgorithm algorithm = createAlgorithm(updater, 8, true);
236 final Vector3D position = new Vector3D(-3787079.6453602533, 5856784.405679551, 1655869.0582939098);
237 AtmosphericRefraction model = new MultiLayerModel(earth);
238
239
240
241
242
243 UnivariateFunction reference = alpha -> 1.17936 * FastMath.tan((2.94613 - 1.40162 * alpha) * alpha);
244
245 for (double alpha = 0; alpha < FastMath.toRadians(50.0); alpha += 0.1) {
246 final Vector3D rotatingLos = los(position, alpha);
247 final NormalizedGeodeticPoint rawIntersection = algorithm.refineIntersection(earth, position, rotatingLos,
248 algorithm.intersection(earth, position, rotatingLos));
249
250 final NormalizedGeodeticPoint correctedIntersection = model.applyCorrection(position, rotatingLos, rawIntersection, algorithm);
251 final double distance = Vector3D.distance(earth.transform(rawIntersection),
252 earth.transform(correctedIntersection));
253 Assertions.assertEquals(reference.value(alpha), distance, 0.12);
254 }
255
256 }
257
258 private Vector3D los(final Vector3D position, final double angleFromNadir) {
259
260 final Vector3D nadir = earth.transform(position, earth.getBodyFrame(), null).getNadir();
261 final Rotation losRotation = new Rotation(nadir.orthogonal(), angleFromNadir,
262 RotationConvention.VECTOR_OPERATOR);
263 return losRotation.applyTo(nadir);
264 }
265
266 @Override
267 protected IntersectionAlgorithm createAlgorithm(TileUpdater updater, int maxCachedTiles, boolean isOverlappingTiles) {
268 return new DuvenhageAlgorithm(updater, maxCachedTiles, false, isOverlappingTiles);
269 }
270 }