1 package org.orekit.rugged.refraction;
2
3 import static org.junit.jupiter.api.Assertions.assertEquals;
4 import static org.junit.jupiter.api.Assertions.assertFalse;
5 import static org.junit.jupiter.api.Assertions.assertTrue;
6 import static org.junit.jupiter.api.Assertions.fail;
7
8 import java.io.File;
9 import java.net.URISyntaxException;
10 import java.util.ArrayList;
11 import java.util.Arrays;
12 import java.util.List;
13 import java.util.stream.Stream;
14
15 import org.hipparchus.analysis.differentiation.Derivative;
16 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
17 import org.hipparchus.geometry.euclidean.threed.Rotation;
18 import org.hipparchus.geometry.euclidean.threed.RotationConvention;
19 import org.hipparchus.geometry.euclidean.threed.Vector3D;
20 import org.hipparchus.util.FastMath;
21 import org.junit.jupiter.api.Assertions;
22 import org.junit.jupiter.api.Test;
23 import org.orekit.bodies.BodyShape;
24 import org.orekit.bodies.GeodeticPoint;
25 import org.orekit.data.DataContext;
26 import org.orekit.data.DirectoryCrawler;
27 import org.orekit.frames.Frame;
28 import org.orekit.frames.FramesFactory;
29 import org.orekit.models.earth.ReferenceEllipsoid;
30 import org.orekit.orbits.Orbit;
31 import org.orekit.rugged.TestUtils;
32 import org.orekit.rugged.api.AlgorithmId;
33 import org.orekit.rugged.api.BodyRotatingFrameId;
34 import org.orekit.rugged.api.EllipsoidId;
35 import org.orekit.rugged.api.InertialFrameId;
36 import org.orekit.rugged.api.Rugged;
37 import org.orekit.rugged.api.RuggedBuilder;
38 import org.orekit.rugged.errors.RuggedException;
39 import org.orekit.rugged.errors.RuggedMessages;
40 import org.orekit.rugged.linesensor.LineDatation;
41 import org.orekit.rugged.linesensor.LineSensor;
42 import org.orekit.rugged.linesensor.LinearLineDatation;
43 import org.orekit.rugged.linesensor.SensorPixel;
44 import org.orekit.rugged.los.LOSBuilder;
45 import org.orekit.rugged.los.TimeDependentLOS;
46 import org.orekit.rugged.raster.RandomLandscapeUpdater;
47 import org.orekit.rugged.raster.TileUpdater;
48 import org.orekit.rugged.utils.DerivativeGenerator;
49 import org.orekit.rugged.utils.GeodeticUtilities;
50 import org.orekit.time.AbsoluteDate;
51 import org.orekit.time.TimeScalesFactory;
52 import org.orekit.utils.AngularDerivativesFilter;
53 import org.orekit.utils.CartesianDerivativesFilter;
54 import org.orekit.utils.Constants;
55 import org.orekit.utils.IERSConventions;
56 import org.orekit.utils.ParameterDriver;
57 import org.orekit.utils.TimeStampedAngularCoordinates;
58 import org.orekit.utils.TimeStampedPVCoordinates;
59
60 public class AtmosphericRefractionTest {
61
62
63 @Test
64 public void testAtmosphericRefractionCorrection() throws URISyntaxException {
65
66 String sensorName = "line";
67 int dimension = 4000;
68
69 RuggedBuilder builder = initRuggedForAtmosphericTests(dimension, sensorName);
70
71
72 Rugged ruggedWithout = builder.build();
73
74
75 AtmosphericRefraction atmosphericRefraction = new MultiLayerModel(ruggedWithout.getEllipsoid());
76 int pixelStep = 100;
77 int lineStep = 100;
78 atmosphericRefraction.setGridSteps(pixelStep, lineStep);
79
80
81 builder.setRefractionCorrection(atmosphericRefraction);
82 Rugged ruggedWith = builder.build();
83
84
85 assertTrue(ruggedWith.getRefractionCorrection().getClass().isInstance(new MultiLayerModel(ruggedWith.getEllipsoid())));
86
87
88 LineSensor lineSensor = ruggedWithout.getLineSensor(sensorName);
89 int minLine = (int) FastMath.floor(lineSensor.getLine(ruggedWithout.getMinDate()));
90 int maxLine = (int) FastMath.ceil(lineSensor.getLine(ruggedWithout.getMaxDate()));
91
92 final double pixelThreshold = 1.e-3;
93 final double lineThreshold = 1.e-2;
94
95 final double epsilonPixel = pixelThreshold;
96 final double epsilonLine = lineThreshold;
97 final double earthRadius = ruggedWithout.getEllipsoid().getEquatorialRadius();
98
99
100
101 double chosenLine = 200.;
102 GeodeticPoint[] gpWithoutAtmosphericRefractionCorrection = ruggedWithout.directLocation(sensorName, chosenLine);
103 GeodeticPoint[] gpWithAtmosphericRefractionCorrection = ruggedWith.directLocation(sensorName, chosenLine);
104
105
106 for (int i = 0; i < gpWithAtmosphericRefractionCorrection.length; i++) {
107 double currentRadius = earthRadius + (gpWithAtmosphericRefractionCorrection[i].getAltitude()+ gpWithoutAtmosphericRefractionCorrection[i].getAltitude())/2.;
108 double distance = GeodeticUtilities.computeDistanceInMeter(currentRadius, gpWithAtmosphericRefractionCorrection[i], gpWithoutAtmosphericRefractionCorrection[i]);
109
110 Assertions.assertTrue(distance > 0.0);
111 Assertions.assertTrue(distance < 2.);
112 }
113
114
115
116 for (int i = 0; i < gpWithAtmosphericRefractionCorrection.length; i++) {
117
118
119 GeodeticPoint gpWith = gpWithAtmosphericRefractionCorrection[i];
120 SensorPixel sensorPixelReverseWith = ruggedWith.inverseLocation(sensorName, gpWith, minLine, maxLine);
121
122 if (sensorPixelReverseWith != null) {
123 assertEquals(i, sensorPixelReverseWith.getPixelNumber(), epsilonPixel);
124 assertEquals(chosenLine, sensorPixelReverseWith.getLineNumber(), epsilonLine);
125 } else {
126 fail("Inverse location failed for pixel " + i + " with atmospheric refraction correction for geodetic point computed with" );
127 }
128 }
129
130
131
132 double dummyLat = gpWithAtmosphericRefractionCorrection[0].getLatitude() + FastMath.PI/4.;
133 double dummyLon = gpWithAtmosphericRefractionCorrection[0].getLongitude() - FastMath.PI/4.;
134 GeodeticPoint dummyGP = new GeodeticPoint(dummyLat, dummyLon, 0.);
135 try {
136 ruggedWith.inverseLocation(sensorName, dummyGP, minLine, maxLine);
137 Assertions.fail("an exeption should have been thrown");
138 } catch (RuggedException re) {
139 Assertions.assertEquals(RuggedMessages.SENSOR_PIXEL_NOT_FOUND_IN_RANGE_LINES, re.getSpecifier());
140 }
141
142 try {
143 ruggedWith.inverseLocation(sensorName,
144 gpWithAtmosphericRefractionCorrection[0],
145 210, maxLine);
146 Assertions.fail("an exeption should have been thrown");
147 } catch (RuggedException re) {
148 Assertions.assertEquals(RuggedMessages.SENSOR_PIXEL_NOT_FOUND_IN_RANGE_LINES, re.getSpecifier());
149 }
150
151 try {
152 ruggedWith.inverseLocation(sensorName,
153 gpWithAtmosphericRefractionCorrection[0],
154 minLine, 190);
155 Assertions.fail("an exeption should have been thrown");
156 } catch (RuggedException re) {
157 Assertions.assertEquals(RuggedMessages.SENSOR_PIXEL_NOT_FOUND_IN_RANGE_LINES, re.getSpecifier());
158 }
159
160 }
161
162 private RuggedBuilder initRuggedForAtmosphericTests(final int dimension, final String sensorName) throws URISyntaxException {
163
164 String path = getClass().getClassLoader().getResource("orekit-data").toURI().getPath();
165 DataContext.getDefault().getDataProvidersManager().addProvider(new DirectoryCrawler(new File(path)));
166 final BodyShape earth = TestUtils.createEarth();
167 final Orbit orbit = TestUtils.createOrbit(Constants.EIGEN5C_EARTH_MU);
168
169
170
171
172 Vector3D position = new Vector3D(1.5, 0, -0.2);
173 TimeDependentLOS los = TestUtils.createLOSPerfectLine(new Rotation(Vector3D.PLUS_I,
174 FastMath.toRadians(5.0),
175 RotationConvention.VECTOR_OPERATOR).applyTo(Vector3D.PLUS_K),
176 Vector3D.PLUS_I,
177 FastMath.toRadians((dimension/2.) * 2.6 / 3600.0), dimension).build();
178
179
180
181
182 AbsoluteDate crossing = new AbsoluteDate("2012-01-01T12:30:00.000", TimeScalesFactory.getUTC());
183 LineDatation lineDatation = new LinearLineDatation(crossing, dimension / 2, 1.0 / 1.5e-3);
184 int firstLine = 0;
185 int lastLine = dimension;
186 LineSensor lineSensor = new LineSensor(sensorName, lineDatation, position, los);
187 AbsoluteDate minDate = lineSensor.getDate(firstLine).shiftedBy(-1.0);
188 AbsoluteDate maxDate = lineSensor.getDate(lastLine).shiftedBy(+1.0);
189
190 TileUpdater updater = new RandomLandscapeUpdater(800.0, 9000.0, 0.1, 0xf0a401650191f9f6l, FastMath.toRadians(2.0), 257);
191
192
193 RuggedBuilder builder = new RuggedBuilder().
194 setDigitalElevationModel(updater, 8).
195 setAlgorithm(AlgorithmId.DUVENHAGE).
196 setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF).
197 setTimeSpan(minDate, maxDate, 0.001, 5.0).
198 setTrajectory(InertialFrameId.EME2000,
199 TestUtils.orbitToPV(orbit, earth, minDate.shiftedBy(-1.0), maxDate.shiftedBy(+1.0), 0.25),
200 8, CartesianDerivativesFilter.USE_PV,
201 TestUtils.orbitToQ(orbit, earth, minDate.shiftedBy(-1.0), maxDate.shiftedBy(+1.0), 0.25),
202 2, AngularDerivativesFilter.USE_R).
203 setLightTimeCorrection(false).
204 setAberrationOfLightCorrection(false).
205 addLineSensor(lineSensor);
206
207 return builder;
208 }
209
210
211
212
213 @Test
214 public void testInverseLocationMargin() throws URISyntaxException {
215
216 String path = getClass().getClassLoader().getResource("orekit-data").toURI().getPath();
217 DataContext.getDefault().getDataProvidersManager().addProvider(new DirectoryCrawler(new File(path)));
218
219 RuggedBuilder builder = new RuggedBuilder();
220
221 Frame ecf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
222 builder.setEllipsoid(ReferenceEllipsoid.getWgs84(ecf));
223
224 MultiLayerModel atmosphere = new MultiLayerModel(builder.getEllipsoid());
225 builder.setRefractionCorrection(atmosphere);
226
227 builder.setLightTimeCorrection(true);
228 builder.setAberrationOfLightCorrection(true);
229 builder.setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID);
230
231 AbsoluteDate start = AbsoluteDate.ARBITRARY_EPOCH;
232 AbsoluteDate end = start.shiftedBy(10);
233 AbsoluteDate middle = start.shiftedBy(end.durationFrom(start) / 2);
234 builder.setTimeSpan(start, end, 1e-3, 1e-3);
235
236 final double h = 500e3;
237 Vector3D p = new Vector3D(6378137 + h, 0, 0);
238 Vector3D v = Vector3D.ZERO;
239 List<TimeStampedPVCoordinates> pvs = Arrays.asList(
240 new TimeStampedPVCoordinates(start, p, v),
241 new TimeStampedPVCoordinates(end, p, v));
242
243 Rotation rotation = new Rotation(Vector3D.MINUS_I, Vector3D.MINUS_K, Vector3D.PLUS_K, Vector3D.PLUS_I);
244 TimeStampedAngularCoordinates attitude =
245 new TimeStampedAngularCoordinates(
246 middle, rotation,
247 Vector3D.PLUS_I.scalarMultiply(0.1), Vector3D.ZERO);
248 List<TimeStampedAngularCoordinates> attitudes = Arrays.asList(
249 attitude.shiftedBy(start.durationFrom(attitude.getDate())),
250 attitude,
251 attitude.shiftedBy(end.durationFrom(attitude.getDate())));
252
253 builder.setTrajectory(ecf,
254 pvs, 2, CartesianDerivativesFilter.USE_P,
255 attitudes, 2, AngularDerivativesFilter.USE_R);
256
257 final double iFov = 1e-6;
258 TimeDependentLOS los = new TimeDependentLOS() {
259 @Override
260 public int getNbPixels() {
261 return 1000;
262 }
263
264 @Override
265 public Vector3D getLOS(int index, AbsoluteDate date) {
266
267 final double center = getNbPixels() / 2.0;
268 final double x = (index - center);
269 final double los = x * iFov;
270 return new Vector3D(los, 0, 1);
271 }
272
273 @Override
274 public <T extends Derivative<T>> FieldVector3D<T> getLOSDerivatives(
275 int index,
276 AbsoluteDate date,
277 DerivativeGenerator<T> generator) {
278 throw new UnsupportedOperationException("not implemented");
279 }
280
281 @Override
282 public Stream<ParameterDriver> getParametersDrivers() {
283 return Stream.empty();
284 }
285 };
286
287 LineSensor sensor = new LineSensor("sensor",
288 new LinearLineDatation(middle, 0, 1000),
289 Vector3D.ZERO,
290 los);
291 builder.addLineSensor(sensor);
292
293 Rugged ruggedWithDefaultMargin = builder.build();
294
295 GeodeticPoint point = ruggedWithDefaultMargin.directLocation(sensor.getName(), 1000)[500];
296 try {
297 final int maxLine = 4999;
298 ruggedWithDefaultMargin.inverseLocation(sensor.getName(), point, 0, maxLine);
299 Assertions.fail("An exception should have been thrown");
300
301 } catch (RuggedException re) {
302 Assertions.assertEquals(RuggedMessages.SENSOR_PIXEL_NOT_FOUND_IN_PIXELS_LINE,re.getSpecifier());
303 }
304
305
306 Assertions.assertEquals(builder.getRefractionCorrection().getComputationParameters().getDefaultInverseLocMargin(),
307 builder.getRefractionCorrection().getComputationParameters().getInverseLocMargin(),
308 1.0e-10);
309
310
311 builder.getRefractionCorrection().setInverseLocMargin(0.81);
312 Rugged ruggedWithCustomMargin = builder.build();
313
314 point = ruggedWithCustomMargin.directLocation(sensor.getName(), 1000)[500];
315 final int maxLine = 4999;
316 SensorPixel pixel = ruggedWithCustomMargin.inverseLocation(sensor.getName(), point, 0, maxLine);
317 Assertions.assertTrue(pixel != null);
318
319 }
320
321
322
323
324 @Test
325 public void testLightTimeCorrection() throws URISyntaxException {
326
327 String sensorName = "line";
328 int dimension = 4000;
329
330 RuggedBuilder builder = initRuggedForAtmosphericTests(dimension, sensorName);
331
332
333 builder.setLightTimeCorrection(true);
334 Rugged ruggedWithLightTimeWithoutRefraction = builder.build();
335
336
337 AtmosphericRefraction atmosphericRefraction = new MultiLayerModel(builder.getEllipsoid());
338 int pixelStep = 100;
339 int lineStep = 100;
340 atmosphericRefraction.setGridSteps(pixelStep, lineStep);
341
342
343 builder.setRefractionCorrection(atmosphericRefraction);
344
345
346 builder.setLightTimeCorrection(false);
347 Rugged ruggedWithoutLightTime = builder.build();
348
349
350 builder.setLightTimeCorrection(true);
351 Rugged ruggedWithLightTime = builder.build();
352
353
354
355
356
357
358
359
360 double chosenLine = 200.;
361 GeodeticPoint[] gpWithoutLightTime = ruggedWithoutLightTime.directLocation(sensorName, chosenLine);
362 GeodeticPoint[] gpWithLightTime = ruggedWithLightTime.directLocation(sensorName, chosenLine);
363
364 GeodeticPoint[] gpWithLightTimeWithoutRefraction = ruggedWithLightTimeWithoutRefraction.directLocation(sensorName, chosenLine);
365
366 double earthRadius = builder.getEllipsoid().getEquatorialRadius();
367
368
369 for (int i = 0; i < gpWithLightTime.length; i++) {
370
371 double currentRadius = earthRadius + (gpWithLightTime[i].getAltitude() + gpWithoutLightTime[i].getAltitude())/2.;
372
373 double distance = GeodeticUtilities.computeDistanceInMeter(currentRadius, gpWithLightTime[i], gpWithoutLightTime[i]);
374
375
376 Assertions.assertTrue(distance > 0.0);
377 Assertions.assertTrue(distance <= 1.2);
378
379
380 distance = GeodeticUtilities.computeDistanceInMeter(currentRadius, gpWithLightTime[i], gpWithLightTimeWithoutRefraction[i]);
381
382 Assertions.assertTrue(distance > 0.0);
383 Assertions.assertTrue(distance < 2.);
384 }
385 }
386
387 @Test
388 public void testBadConfig() {
389
390 int dimension = 400;
391
392 TimeDependentLOS los = new LOSBuilder(new ArrayList<Vector3D>(dimension)).build();
393 LineSensor lineSensor = new LineSensor("line", null, Vector3D.ZERO, los);
394
395
396 AtmosphericRefraction atmosphericRefraction = new MultiLayerModel(null);
397
398
399 atmosphericRefraction.setGridSteps(100, 100);
400 atmosphericRefraction.configureCorrectionGrid(lineSensor, 0, 300);
401 assertFalse(atmosphericRefraction.isSameContext("otherSensor", 0, 300));
402 assertFalse(atmosphericRefraction.isSameContext("line", 42, 300));
403 assertFalse(atmosphericRefraction.isSameContext("line", 0, 42));
404
405
406 try {
407 atmosphericRefraction.setGridSteps(100, 100);
408 atmosphericRefraction.configureCorrectionGrid(lineSensor, 0, 100);
409 Assertions.fail("An exception should have been thrown");
410
411 } catch (RuggedException re) {
412 Assertions.assertEquals(RuggedMessages.INVALID_RANGE_FOR_LINES,re.getSpecifier());
413 }
414
415
416 try {
417 atmosphericRefraction.setGridSteps(-5, 100);
418 atmosphericRefraction.configureCorrectionGrid(lineSensor, 0, 100);
419 Assertions.fail("An exception should have been thrown");
420
421 } catch (RuggedException re) {
422 Assertions.assertEquals(RuggedMessages.INVALID_STEP,re.getSpecifier());
423 }
424
425
426 try {
427 atmosphericRefraction.setGridSteps(10, -42);
428 atmosphericRefraction.configureCorrectionGrid(lineSensor, 0, 100);
429 Assertions.fail("An exception should have been thrown");
430
431 } catch (RuggedException re) {
432 Assertions.assertEquals(RuggedMessages.INVALID_STEP,re.getSpecifier());
433 }
434 }
435 }