1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation;
18
19 import java.util.ArrayList;
20 import java.util.List;
21 import java.util.Locale;
22 import java.util.Map;
23
24 import org.hipparchus.analysis.polynomials.SmoothStepFactory;
25 import org.hipparchus.stat.descriptive.DescriptiveStatistics;
26 import org.hipparchus.util.FastMath;
27 import org.junit.jupiter.api.Assertions;
28 import org.junit.jupiter.api.BeforeAll;
29 import org.junit.jupiter.api.DisplayName;
30 import org.junit.jupiter.api.Test;
31 import org.orekit.TestUtils;
32 import org.orekit.Utils;
33 import org.orekit.data.DataContext;
34 import org.orekit.forces.gravity.potential.AstronomicalAmplitudeReader;
35 import org.orekit.forces.gravity.potential.EGMFormatReader;
36 import org.orekit.forces.gravity.potential.FESCHatEpsilonReader;
37 import org.orekit.forces.gravity.potential.GravityFieldFactory;
38 import org.orekit.forces.gravity.potential.OceanLoadDeformationCoefficients;
39 import org.orekit.frames.Frame;
40 import org.orekit.frames.FramesFactory;
41 import org.orekit.frames.LOFType;
42 import org.orekit.orbits.Orbit;
43 import org.orekit.orbits.OrbitBlender;
44 import org.orekit.orbits.OrbitHermiteInterpolator;
45 import org.orekit.orbits.OrbitType;
46 import org.orekit.orbits.PositionAngleType;
47 import org.orekit.propagation.analytical.AbstractAnalyticalPropagator;
48 import org.orekit.propagation.analytical.BrouwerLyddanePropagator;
49 import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
50 import org.orekit.propagation.analytical.KeplerianPropagator;
51 import org.orekit.time.AbsoluteDate;
52 import org.orekit.time.TimeInterpolator;
53 import org.orekit.time.TimeStampedPair;
54 import org.orekit.utils.Constants;
55
56 class StateCovarianceBlenderTest {
57
58 private static SpacecraftState sergeiState;
59 private static Orbit sergeiOrbit;
60 private static Frame sergeiFrame;
61
62
63 private final double DEFAULT_SERGEI_PROPAGATION_TIME = 2400;
64 private final double DEFAUTL_SERGEI_TABULATED_TIMESTEP = 2400;
65
66 @BeforeAll
67 public static void setUp() {
68 Utils.setDataRoot("regular-data:potential/egm-format:atmosphere:tides:regular-data/de405-ephemerides");
69 GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("EGM96-truncated-21x21", true));
70 AstronomicalAmplitudeReader aaReader =
71 new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
72 DataContext.getDefault().getDataProvidersManager().feed(aaReader.getSupportedNames(), aaReader);
73 Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
74 GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
75 0.01, FastMath.toRadians(1.0),
76 OceanLoadDeformationCoefficients.IERS_2010,
77 map));
78
79 sergeiState = StateCovarianceKeplerianHermiteInterpolatorTest.generateSergeiReferenceState();
80 sergeiOrbit = sergeiState.getOrbit();
81 sergeiFrame = sergeiOrbit.getFrame();
82 }
83
84 private void doTestBlending(final double propagationHorizon, final double tabulatedTimeStep,
85 final SmoothStepFactory.SmoothStepFunction blendingFunction,
86 final AbstractAnalyticalPropagator analyticalPropagator,
87 final double expectedMeanRMSPositionError,
88 final double expectedMeanRMSVelocityError,
89 final double expectedMedianRMSPositionError,
90 final double expectedMedianRMSVelocityError,
91 final double expectedMaxRMSPositionError,
92 final double expectedMaxRMSVelocityError,
93 final double tolerance,
94 final boolean showResults) {
95 final TimeInterpolator<Orbit> orbitInterpolator = new OrbitBlender(blendingFunction,
96 analyticalPropagator,
97 sergeiFrame);
98
99 final TimeInterpolator<TimeStampedPair<Orbit, StateCovariance>> covarianceInterpolator =
100 new StateCovarianceBlender(blendingFunction, orbitInterpolator,
101 sergeiFrame, OrbitType.CARTESIAN,
102 PositionAngleType.MEAN);
103
104
105 final TimeInterpolator<SpacecraftState> stateInterpolator =
106 new SpacecraftStateInterpolator(2, 1.0e-3, sergeiFrame, orbitInterpolator, null, null, null, null);
107
108
109 final DescriptiveStatistics[] relativeRMSSigmaError =
110 StateCovarianceKeplerianHermiteInterpolatorTest.computeStatisticsCovarianceInterpolationOnSergeiCase(
111 propagationHorizon, tabulatedTimeStep, stateInterpolator, covarianceInterpolator);
112
113
114 if (showResults) {
115 System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMean", relativeRMSSigmaError[0].getMean());
116 System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMean", relativeRMSSigmaError[1].getMean());
117 System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMedian", relativeRMSSigmaError[0].getPercentile(50));
118 System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMedian", relativeRMSSigmaError[1].getPercentile(50));
119 System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMax", relativeRMSSigmaError[0].getMax());
120 System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMax", relativeRMSSigmaError[1].getMax());
121
122 }
123
124
125 Assertions.assertEquals(expectedMeanRMSPositionError, relativeRMSSigmaError[0].getMean(), tolerance);
126 Assertions.assertEquals(expectedMeanRMSVelocityError, relativeRMSSigmaError[1].getMean(), tolerance);
127 Assertions.assertEquals(expectedMedianRMSPositionError, relativeRMSSigmaError[0].getPercentile(50), tolerance);
128 Assertions.assertEquals(expectedMedianRMSVelocityError, relativeRMSSigmaError[1].getPercentile(50), tolerance);
129 Assertions.assertEquals(expectedMaxRMSPositionError, relativeRMSSigmaError[0].getMax(), tolerance);
130 Assertions.assertEquals(expectedMaxRMSVelocityError, relativeRMSSigmaError[1].getMax(), tolerance);
131 }
132
133
134
135
136
137
138
139
140
141
142 @Test
143 @DisplayName("test Keplerian quadratic blending interpolation on full force model test case from: "
144 + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
145 + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
146 void testKeplerianQuadraticBlending() {
147
148 final boolean showResults = false;
149 final double tolerance = 1.e-6;
150
151
152 final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
153
154
155 doTestBlending(DEFAULT_SERGEI_PROPAGATION_TIME, DEFAUTL_SERGEI_TABULATED_TIMESTEP, blendingFunction,
156 new KeplerianPropagator(sergeiOrbit),
157 0.1089229597534360,
158 0.2597137453057761,
159 0.1051914244596399,
160 0.2804213769546742,
161 0.2970138301592348,
162 0.5755370093511280,
163 tolerance,
164 showResults);
165
166
167
168
169
170
171
172
173 }
174
175
176
177
178
179
180
181
182 @Test
183 @DisplayName("test Brouwer Lyddane quadratic blending interpolation on full force model test case from: "
184 + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
185 + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
186 void testBrouwerLyddaneQuadraticBlending() {
187
188 final boolean showResults = false;
189 final double tolerance = 1.e-6;
190
191
192 final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
193
194 final AbstractAnalyticalPropagator propagator = new BrouwerLyddanePropagator(sergeiOrbit,
195 sergeiState.getMass(),
196 Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
197 Constants.EIGEN5C_EARTH_MU,
198 Constants.EIGEN5C_EARTH_C20,
199 Constants.EIGEN5C_EARTH_C30,
200 Constants.EIGEN5C_EARTH_C40,
201 Constants.EIGEN5C_EARTH_C50,
202 0);
203
204
205 doTestBlending(DEFAULT_SERGEI_PROPAGATION_TIME, DEFAUTL_SERGEI_TABULATED_TIMESTEP, blendingFunction,
206 propagator,
207 0.1385177282847412,
208 0.1559076539117118,
209 0.1518228759308863,
210 0.1537370078767798,
211 0.2189165721963918,
212 0.4694243494150299,
213 tolerance,
214 showResults);
215
216
217
218
219
220
221
222
223 }
224
225
226
227
228
229
230
231
232 @Deprecated
233 @Test
234 @DisplayName("test Brouwer Lyddane quadratic blending interpolation on full force model test case from: "
235 + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
236 + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
237 void testBrouwerLyddaneQuadraticBlendingDeprecated() {
238
239 final boolean showResults = false;
240 final double tolerance = 1.e-6;
241
242
243 final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
244
245 final AbstractAnalyticalPropagator propagator = new BrouwerLyddanePropagator(sergeiOrbit,
246 sergeiState.getMass(),
247 Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
248 Constants.EIGEN5C_EARTH_MU,
249 Constants.EIGEN5C_EARTH_C20,
250 Constants.EIGEN5C_EARTH_C30,
251 Constants.EIGEN5C_EARTH_C40,
252 Constants.EIGEN5C_EARTH_C50,
253 0);
254
255
256 doTestBlending(DEFAULT_SERGEI_PROPAGATION_TIME, DEFAUTL_SERGEI_TABULATED_TIMESTEP, blendingFunction,
257 propagator,
258 0.1385177282847412,
259 0.1559076539117118,
260 0.1518228759308863,
261 0.1537370078767798,
262 0.2189165721963918,
263 0.4694243494150299,
264 tolerance,
265 showResults);
266
267
268
269
270
271
272
273
274 }
275
276
277
278
279
280
281
282
283 @Test
284 @DisplayName("test Ekstein Hechler quadratic blending interpolation on full force model test case from: "
285 + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
286 + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
287 void testEksteinHechlerQuadraticBlending() {
288
289 final boolean showResults = false;
290 final double tolerance = 1.e-6;
291
292
293 final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
294
295 final AbstractAnalyticalPropagator propagator = new EcksteinHechlerPropagator(sergeiOrbit,
296 sergeiState.getMass(),
297 Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
298 Constants.EIGEN5C_EARTH_MU,
299 Constants.EIGEN5C_EARTH_C20,
300 Constants.EIGEN5C_EARTH_C30,
301 Constants.EIGEN5C_EARTH_C40,
302 Constants.EIGEN5C_EARTH_C50,
303 Constants.EIGEN5C_EARTH_C60);
304
305
306 doTestBlending(DEFAULT_SERGEI_PROPAGATION_TIME, DEFAUTL_SERGEI_TABULATED_TIMESTEP, blendingFunction,
307 propagator,
308 0.0992103832979174,
309 0.2033256863240140,
310 0.0996928970142715,
311 0.2152205038403896,
312 0.1819318695859633,
313 0.5458515065235870,
314 tolerance,
315 showResults);
316
317
318
319
320
321
322
323
324 }
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340 @Test
341 @DisplayName("test Keplerian quadratic blending expressed in LOF on full force model test case from : "
342 + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
343 + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
344 void testLOFKeplerianBlending() {
345
346 final boolean showResults = false;
347 final double tolerance = 1.e-6;
348
349
350 final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
351
352 final TimeInterpolator<Orbit> orbitInterpolator = new OrbitBlender(blendingFunction,
353 new KeplerianPropagator(sergeiOrbit),
354 sergeiFrame);
355
356 final LOFType DEFAULT_LOFTYPE = LOFType.TNW;
357 final AbstractStateCovarianceInterpolator covarianceInterpolator =
358 new StateCovarianceBlender(blendingFunction, orbitInterpolator, DEFAULT_LOFTYPE);
359
360
361 final DescriptiveStatistics[] relativeRMSSigmaError =
362 StateCovarianceKeplerianHermiteInterpolatorTest.computeStatisticsCovarianceLOFInterpolation(
363 DEFAULT_SERGEI_PROPAGATION_TIME, DEFAUTL_SERGEI_TABULATED_TIMESTEP,
364 DEFAULT_LOFTYPE, covarianceInterpolator);
365
366
367 if (showResults) {
368 System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMean", relativeRMSSigmaError[0].getMean());
369 System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMean", relativeRMSSigmaError[1].getMean());
370 System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMedian", relativeRMSSigmaError[0].getPercentile(50));
371 System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMedian", relativeRMSSigmaError[1].getPercentile(50));
372 System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMax", relativeRMSSigmaError[0].getMax());
373 System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMax", relativeRMSSigmaError[1].getMax());
374
375 }
376
377
378 Assertions.assertEquals( 0.1219305982260829, relativeRMSSigmaError[0].getMean(), tolerance);
379 Assertions.assertEquals(19.1744012628613100, relativeRMSSigmaError[1].getMean(), tolerance);
380 Assertions.assertEquals( 0.1268691798373903, relativeRMSSigmaError[0].getPercentile(50), tolerance);
381 Assertions.assertEquals(15.074357209468475, relativeRMSSigmaError[1].getPercentile(50), tolerance);
382 Assertions.assertEquals( 0.2395147080521065, relativeRMSSigmaError[0].getMax(), tolerance);
383 Assertions.assertEquals(75.13545251822399, relativeRMSSigmaError[1].getMax(), 3 * tolerance);
384
385
386 Assertions.assertNull(covarianceInterpolator.getOutFrame());
387 Assertions.assertEquals(DEFAULT_LOFTYPE, covarianceInterpolator.getOutLOF());
388 Assertions.assertEquals(OrbitType.CARTESIAN, covarianceInterpolator.getOutOrbitType());
389 Assertions.assertEquals(PositionAngleType.MEAN, covarianceInterpolator.getOutPositionAngleType());
390 Assertions.assertEquals(orbitInterpolator, covarianceInterpolator.getOrbitInterpolator());
391
392 }
393
394
395
396
397
398
399
400 @Test
401 void testStateCovarianceBlenderWithOrbitInterpolatorUsingMoreInterpolationPoints() {
402
403
404 final OrbitHermiteInterpolator orbitInterpolator =
405 new OrbitHermiteInterpolator(3, FramesFactory.getGCRF());
406
407
408 final StateCovarianceBlender covarianceBlender = new StateCovarianceBlender(SmoothStepFactory.getQuadratic(),
409 orbitInterpolator,
410 LOFType.TNW);
411
412
413 final AbsoluteDate interpolationDate = new AbsoluteDate();
414
415
416 final Orbit orbit1 = TestUtils.getDefaultOrbit(interpolationDate.shiftedBy(-1));
417 final Orbit orbit2 = TestUtils.getDefaultOrbit(interpolationDate);
418 final Orbit orbit3 = TestUtils.getDefaultOrbit(interpolationDate.shiftedBy(1));
419
420
421 final StateCovariance covariance1 = TestUtils.getFakeStateCovarianceInLOF(interpolationDate.shiftedBy(-1), LOFType.TNW);
422 final StateCovariance covariance2 = TestUtils.getFakeStateCovarianceInLOF(interpolationDate, LOFType.TNW);
423 final StateCovariance covariance3 = TestUtils.getFakeStateCovarianceInLOF(interpolationDate.shiftedBy(1), LOFType.TNW);
424
425
426
427 final List<TimeStampedPair<Orbit, StateCovariance>> samples = new ArrayList<>();
428 samples.add(new TimeStampedPair<>(orbit1, covariance1));
429 samples.add(new TimeStampedPair<>(orbit2, covariance2));
430 samples.add(new TimeStampedPair<>(orbit3, covariance3));
431
432
433
434 Assertions.assertDoesNotThrow(() -> covarianceBlender.interpolate(interpolationDate.shiftedBy(-1), samples));
435 Assertions.assertDoesNotThrow(() -> covarianceBlender.interpolate(interpolationDate.shiftedBy(-0.9), samples));
436 Assertions.assertDoesNotThrow(() -> covarianceBlender.interpolate(interpolationDate.shiftedBy(-0.1), samples));
437 Assertions.assertDoesNotThrow(() -> covarianceBlender.interpolate(interpolationDate, samples));
438 Assertions.assertDoesNotThrow(() -> covarianceBlender.interpolate(interpolationDate.shiftedBy(0.1), samples));
439 Assertions.assertDoesNotThrow(() -> covarianceBlender.interpolate(interpolationDate.shiftedBy(0.9), samples));
440 Assertions.assertDoesNotThrow(() -> covarianceBlender.interpolate(interpolationDate.shiftedBy(1), samples));
441
442 }
443
444
445
446
447
448
449 @Test
450 void testGetSubInterpolators() {
451
452
453 final OrbitHermiteInterpolator orbitInterpolator =
454 new OrbitHermiteInterpolator(FramesFactory.getGCRF());
455
456
457 final StateCovarianceBlender covarianceBlender = new StateCovarianceBlender(SmoothStepFactory.getQuadratic(),
458 orbitInterpolator,
459 LOFType.TNW);
460
461
462 final List<TimeInterpolator<?>> actualSubInterpolators = covarianceBlender.getSubInterpolators();
463
464
465 Assertions.assertEquals(1, actualSubInterpolators.size());
466 Assertions.assertEquals(orbitInterpolator, actualSubInterpolators.get(0));
467 }
468
469 }