1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation;
18
19 import org.hipparchus.analysis.polynomials.SmoothStepFactory;
20 import org.hipparchus.geometry.euclidean.threed.Vector3D;
21 import org.hipparchus.linear.BlockRealMatrix;
22 import org.hipparchus.linear.RealMatrix;
23 import org.hipparchus.ode.ODEIntegrator;
24 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
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.mockito.Mockito;
32 import org.orekit.Utils;
33 import org.orekit.bodies.BodyShape;
34 import org.orekit.bodies.CelestialBody;
35 import org.orekit.bodies.CelestialBodyFactory;
36 import org.orekit.bodies.OneAxisEllipsoid;
37 import org.orekit.data.DataContext;
38 import org.orekit.forces.drag.DragForce;
39 import org.orekit.forces.drag.IsotropicDrag;
40 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
41 import org.orekit.forces.gravity.OceanTides;
42 import org.orekit.forces.gravity.Relativity;
43 import org.orekit.forces.gravity.SolidTides;
44 import org.orekit.forces.gravity.ThirdBodyAttraction;
45 import org.orekit.forces.gravity.potential.AstronomicalAmplitudeReader;
46 import org.orekit.forces.gravity.potential.EGMFormatReader;
47 import org.orekit.forces.gravity.potential.FESCHatEpsilonReader;
48 import org.orekit.forces.gravity.potential.GravityFieldFactory;
49 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
50 import org.orekit.forces.gravity.potential.OceanLoadDeformationCoefficients;
51 import org.orekit.forces.radiation.IsotropicRadiationSingleCoefficient;
52 import org.orekit.forces.radiation.KnockeRediffusedForceModel;
53 import org.orekit.forces.radiation.RadiationSensitive;
54 import org.orekit.forces.radiation.SolarRadiationPressure;
55 import org.orekit.frames.Frame;
56 import org.orekit.frames.FramesFactory;
57 import org.orekit.frames.LOFType;
58 import org.orekit.models.earth.ReferenceEllipsoid;
59 import org.orekit.models.earth.atmosphere.Atmosphere;
60 import org.orekit.models.earth.atmosphere.NRLMSISE00;
61 import org.orekit.models.earth.atmosphere.data.CssiSpaceWeatherData;
62 import org.orekit.orbits.CartesianOrbit;
63 import org.orekit.orbits.Orbit;
64 import org.orekit.orbits.OrbitBlender;
65 import org.orekit.orbits.OrbitType;
66 import org.orekit.orbits.PositionAngleType;
67 import org.orekit.propagation.analytical.Ephemeris;
68 import org.orekit.propagation.analytical.KeplerianPropagator;
69 import org.orekit.propagation.numerical.NumericalPropagator;
70 import org.orekit.time.AbsoluteDate;
71 import org.orekit.time.AbstractTimeInterpolator;
72 import org.orekit.time.TimeInterpolator;
73 import org.orekit.time.TimeScalesFactory;
74 import org.orekit.time.TimeStampedPair;
75 import org.orekit.utils.CartesianDerivativesFilter;
76 import org.orekit.utils.Constants;
77 import org.orekit.utils.IERSConventions;
78 import org.orekit.utils.PVCoordinates;
79
80 import java.util.ArrayList;
81 import java.util.List;
82 import java.util.Locale;
83 import java.util.Map;
84
85 public class StateCovarianceKeplerianHermiteInterpolatorTest {
86 private static Orbit sergeiOrbit;
87 private static Frame sergeiFrame;
88 private final int DEFAULT_SERGEI_INTERPOLATION_POINTS = 2;
89 private final double DEFAULT_SERGEI_PROPAGATION_TIME = 2400;
90 private final double DEFAULT_SERGEI_TABULATED_TIMESTEP = 2400;
91
92 private static final double TOLERANCE = 1e-6;
93
94 @BeforeAll
95 public static void setUp() {
96 Utils.setDataRoot("regular-data:potential/egm-format:atmosphere:tides:regular-data/de405-ephemerides");
97 GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("EGM96-truncated-21x21", true));
98 AstronomicalAmplitudeReader aaReader =
99 new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
100 DataContext.getDefault().getDataProvidersManager().feed(aaReader.getSupportedNames(), aaReader);
101 Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
102 GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
103 0.01, FastMath.toRadians(1.0),
104 OceanLoadDeformationCoefficients.IERS_2010,
105 map));
106
107 sergeiOrbit = generateSergeiReferenceOrbit();
108 sergeiFrame = sergeiOrbit.getFrame();
109 }
110
111 public static DescriptiveStatistics[] computeStatisticsCovarianceLOFInterpolation(final double propagationDuration,
112 final double tabulatedTimeStep,
113 final LOFType lofType,
114 final TimeInterpolator<TimeStampedPair<Orbit, StateCovariance>> covarianceLOFInterpolator) {
115
116 final List<StateCovariance> referenceCovariances = new ArrayList<>();
117 final List<StateCovariance> interpolatedCovariances = new ArrayList<>();
118 final List<TimeStampedPair<Orbit, StateCovariance>> tabulatedOrbitsAndCovariances = new ArrayList<>();
119
120
121 final SpacecraftState sergeiState = generateSergeiReferenceState();
122 final Orbit sergeiOrbit = sergeiState.getOrbit();
123 final AbsoluteDate initialDate = sergeiOrbit.getDate();
124
125
126 final RealMatrix sergeiCovarianceMatrix = generateSergeiCovarianceMatrix();
127
128
129 final NumericalPropagator propagator = new NumericalPropagator(
130 generateDefaultIntegrator(sergeiOrbit, OrbitType.CARTESIAN));
131
132 propagator.setOrbitType(OrbitType.CARTESIAN);
133
134
135 final MatricesHarvester harvester =
136 propagator.setupMatricesComputation("harvester", null, null);
137
138
139 final StateCovariance sergeiCovariance =
140 new StateCovariance(sergeiCovarianceMatrix, sergeiState.getDate(), sergeiState.getFrame(),
141 OrbitType.CARTESIAN, PositionAngleType.MEAN);
142
143 final StateCovarianceMatrixProvider stateCovarianceMatrixProvider =
144 new StateCovarianceMatrixProvider("covariance", "harvester", harvester, sergeiCovariance);
145
146
147 propagator.setInitialState(sergeiState);
148
149 configurePropagatorForSergeiCase(propagator);
150
151 propagator.addAdditionalDataProvider(stateCovarianceMatrixProvider);
152
153 propagator.getMultiplexer().add(1, (currentState) -> {
154
155
156 final StateCovariance currentCovarianceFromProviderInCartesian =
157 stateCovarianceMatrixProvider.getStateCovariance(currentState);
158
159
160 final StateCovariance currentCovarianceFromProviderInLOF =
161 currentCovarianceFromProviderInCartesian.changeCovarianceFrame(currentState.getOrbit(), lofType);
162
163 referenceCovariances.add(currentCovarianceFromProviderInLOF);
164
165
166 final double durationFromStart = currentState.getDate().durationFrom(sergeiState.getDate());
167 if (durationFromStart % tabulatedTimeStep == 0) {
168 tabulatedOrbitsAndCovariances.add(new TimeStampedPair<>(currentState.getOrbit(),
169 currentCovarianceFromProviderInCartesian));
170 }
171 });
172
173
174 propagator.propagate(initialDate.shiftedBy(propagationDuration));
175
176
177 for (int dt = 0; dt < referenceCovariances.size(); dt++) {
178 final AbsoluteDate currentInterpolationDate = initialDate.shiftedBy(dt);
179 final StateCovariance currentInterpolatedCovariance =
180 covarianceLOFInterpolator.interpolate(currentInterpolationDate, tabulatedOrbitsAndCovariances)
181 .getSecond();
182
183 interpolatedCovariances.add(currentInterpolatedCovariance);
184 }
185
186
187 return computeRMSRelativeErrorsStatistics(referenceCovariances,
188 interpolatedCovariances);
189
190 }
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205 public static DescriptiveStatistics[] computeStatisticsCovarianceInterpolationOnSergeiCase(
206 final double propagationDuration,
207 final double tabulatedTimeStep,
208 final TimeInterpolator<SpacecraftState> stateInterpolator,
209 final TimeInterpolator<TimeStampedPair<Orbit, StateCovariance>> covarianceInterpolator) {
210
211
212 final List<StateCovariance> referenceCovariances = new ArrayList<>();
213 final List<StateCovariance> interpolatedCovariances = new ArrayList<>();
214 final List<SpacecraftState> tabulatedStates = new ArrayList<>();
215 final List<StateCovariance> tabulatedCovariances = new ArrayList<>();
216
217
218 final SpacecraftState sergeiState = generateSergeiReferenceState();
219 final Orbit sergeiOrbit = sergeiState.getOrbit();
220 final AbsoluteDate initialDate = sergeiOrbit.getDate();
221
222
223 final RealMatrix sergeiCovarianceMatrix = generateSergeiCovarianceMatrix();
224
225
226 final NumericalPropagator propagator = new NumericalPropagator(
227 generateDefaultIntegrator(sergeiOrbit, OrbitType.CARTESIAN));
228
229 propagator.setOrbitType(OrbitType.CARTESIAN);
230
231
232 final MatricesHarvester harvester =
233 propagator.setupMatricesComputation("harvester", null, null);
234
235
236 final StateCovariance sergeiCovariance =
237 new StateCovariance(sergeiCovarianceMatrix, sergeiState.getDate(), sergeiState.getFrame(),
238 OrbitType.CARTESIAN, PositionAngleType.MEAN);
239
240 final StateCovarianceMatrixProvider stateCovarianceMatrixProvider =
241 new StateCovarianceMatrixProvider("covariance", "harvester", harvester, sergeiCovariance);
242
243
244 propagator.setInitialState(sergeiState);
245
246 configurePropagatorForSergeiCase(propagator);
247
248 propagator.addAdditionalDataProvider(stateCovarianceMatrixProvider);
249
250 propagator.getMultiplexer().add(1, (currentState) -> {
251
252
253 final StateCovariance currentCovarianceFromProviderInCartesian =
254 stateCovarianceMatrixProvider.getStateCovariance(currentState);
255
256 referenceCovariances.add(currentCovarianceFromProviderInCartesian);
257
258
259 final double durationFromStart = currentState.getDate().durationFrom(sergeiState.getDate());
260 if (durationFromStart % tabulatedTimeStep == 0) {
261 tabulatedStates.add(currentState);
262 tabulatedCovariances.add(currentCovarianceFromProviderInCartesian);
263 }
264 });
265
266
267 propagator.propagate(initialDate.shiftedBy(propagationDuration));
268
269
270 final Ephemeris ephemeris =
271 new Ephemeris(tabulatedStates, stateInterpolator, tabulatedCovariances, covarianceInterpolator);
272
273
274 for (int dt = 0; dt < referenceCovariances.size(); dt++) {
275 final AbsoluteDate currentInterpolationDate = initialDate.shiftedBy(dt);
276
277 interpolatedCovariances.add(ephemeris.getCovariance(currentInterpolationDate).get());
278 }
279
280
281 return computeRMSRelativeErrorsStatistics(referenceCovariances,
282 interpolatedCovariances);
283
284 }
285
286
287
288
289
290
291
292
293
294 public static DescriptiveStatistics[] computeRMSRelativeErrorsStatistics(
295 final List<StateCovariance> referenceCovariances,
296 final List<StateCovariance> interpolatedCovariances) {
297 final DescriptiveStatistics[] maxRelativeRMSPosAndVelError = new DescriptiveStatistics[2];
298
299 final DescriptiveStatistics relativeRMSPosSigmaErrorStat = new DescriptiveStatistics();
300 final DescriptiveStatistics relativeRMSVelSigmaErrorStat = new DescriptiveStatistics();
301
302 for (int i = 0; i < referenceCovariances.size(); i++) {
303
304 final RealMatrix currentReferenceCovariance = referenceCovariances.get(i).getMatrix();
305 final RealMatrix currentInterpolatedCovariance = interpolatedCovariances.get(i).getMatrix();
306
307 final double[] currentReferenceSigmas = extractSigmas(currentReferenceCovariance);
308 final double[] currentInterpolatedSigmas = extractSigmas(currentInterpolatedCovariance);
309
310 final double[] currentReferencePosAndVelSigmas = computeRMSSigmaArray(currentReferenceCovariance);
311
312 final double[] deltaPositionSigmas = new double[] {
313 currentReferenceSigmas[0] - currentInterpolatedSigmas[0],
314 currentReferenceSigmas[1] - currentInterpolatedSigmas[1],
315 currentReferenceSigmas[2] - currentInterpolatedSigmas[2] };
316
317 final double[] deltaVelocitySigmas = new double[] {
318 currentReferenceSigmas[3] - currentInterpolatedSigmas[3],
319 currentReferenceSigmas[4] - currentInterpolatedSigmas[4],
320 currentReferenceSigmas[5] - currentInterpolatedSigmas[5] };
321
322
323 final double relativeRMSPosSigmaError =
324 computeRMS(deltaPositionSigmas) * 100 / currentReferencePosAndVelSigmas[0];
325 final double relativeRMSVelSigmaError =
326 computeRMS(deltaVelocitySigmas) * 100 / currentReferencePosAndVelSigmas[1];
327
328 relativeRMSPosSigmaErrorStat.addValue(relativeRMSPosSigmaError);
329 relativeRMSVelSigmaErrorStat.addValue(relativeRMSVelSigmaError);
330 }
331 maxRelativeRMSPosAndVelError[0] = relativeRMSPosSigmaErrorStat;
332 maxRelativeRMSPosAndVelError[1] = relativeRMSVelSigmaErrorStat;
333
334 return maxRelativeRMSPosAndVelError;
335 }
336
337 public static double[] computeRMSSigmaArray(final RealMatrix covarianceMatrix) {
338
339
340 final double[] rmsSigmaArray = new double[2];
341
342 final double[] positionSigmas = new double[] {
343 FastMath.sqrt(FastMath.abs(covarianceMatrix.getEntry(0, 0))),
344 FastMath.sqrt(FastMath.abs(covarianceMatrix.getEntry(1, 1))),
345 FastMath.sqrt(FastMath.abs(covarianceMatrix.getEntry(2, 2))) };
346
347 final double[] velocitySigmas = new double[] {
348 FastMath.sqrt(FastMath.abs(covarianceMatrix.getEntry(3, 3))),
349 FastMath.sqrt(FastMath.abs(covarianceMatrix.getEntry(4, 4))),
350 FastMath.sqrt(FastMath.abs(covarianceMatrix.getEntry(5, 5))) };
351
352 rmsSigmaArray[0] = computeRMS(positionSigmas);
353 rmsSigmaArray[1] = computeRMS(velocitySigmas);
354
355 return rmsSigmaArray;
356 }
357
358
359
360
361
362
363
364
365 public static double computeRMS(final double[] data) {
366 double sum = 0;
367 for (final double element : data) {
368 sum += element * element;
369 }
370 return FastMath.sqrt(sum / data.length);
371 }
372
373
374
375
376
377
378
379
380 public static double[] extractSigmas(final RealMatrix covarianceMatrix) {
381 return new double[]
382 { FastMath.sqrt(FastMath.abs(covarianceMatrix.getEntry(0, 0))),
383 FastMath.sqrt(FastMath.abs(covarianceMatrix.getEntry(1, 1))),
384 FastMath.sqrt(FastMath.abs(covarianceMatrix.getEntry(2, 2))),
385 FastMath.sqrt(FastMath.abs(covarianceMatrix.getEntry(3, 3))),
386 FastMath.sqrt(FastMath.abs(covarianceMatrix.getEntry(4, 4))),
387 FastMath.sqrt(FastMath.abs(covarianceMatrix.getEntry(5, 5))) };
388 }
389
390 public static ODEIntegrator generateDefaultIntegrator(final Orbit orbit, final OrbitType orbitType) {
391 final double dP = 1;
392 final double[][] TOLERANCEs = ToleranceProvider.getDefaultToleranceProvider(dP).getTolerances(orbit, orbitType);
393 final double minStep = 0.001;
394 final double maxStep = 300.;
395
396 return new DormandPrince853Integrator(minStep, maxStep, TOLERANCEs[0], TOLERANCEs[1]);
397 }
398
399 public static void configurePropagatorForSergeiCase(final NumericalPropagator propagator) {
400
401
402 final IERSConventions conventions = IERSConventions.IERS_2010;
403 final Frame itrf = FramesFactory.getITRF(conventions, true);
404 final BodyShape earthShape = new OneAxisEllipsoid(Constants.IERS2010_EARTH_EQUATORIAL_RADIUS,
405 Constants.IERS2010_EARTH_FLATTENING, itrf);
406
407 final CelestialBody sun = CelestialBodyFactory.getSun();
408 final CelestialBody moon = CelestialBodyFactory.getMoon();
409
410
411 final NormalizedSphericalHarmonicsProvider gravity =
412 GravityFieldFactory.getNormalizedProvider(21, 21);
413 propagator.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, gravity));
414
415
416 final double crossSection = 20;
417
418
419 final double cr = 1;
420 final RadiationSensitive radiationSensitive = new IsotropicRadiationSingleCoefficient(crossSection, cr);
421 final SolarRadiationPressure solarRadiationPressure =
422 new SolarRadiationPressure(sun,
423 ReferenceEllipsoid.getIers2010(itrf),
424 radiationSensitive);
425
426
427 solarRadiationPressure.addOccultingBody(moon, Constants.MOON_EQUATORIAL_RADIUS);
428
429 final KnockeRediffusedForceModel albedo =
430 new KnockeRediffusedForceModel(sun, radiationSensitive,
431 Constants.IERS2010_EARTH_EQUATORIAL_RADIUS,
432 FastMath.toRadians(5));
433
434 propagator.addForceModel(solarRadiationPressure);
435 propagator.addForceModel(albedo);
436
437
438 final double cd = 2.2;
439 final IsotropicDrag dragSensitive = new IsotropicDrag(crossSection, cd);
440 final Atmosphere atmosphere =
441 new NRLMSISE00(new CssiSpaceWeatherData(CssiSpaceWeatherData.DEFAULT_SUPPORTED_NAMES), sun, earthShape);
442 final DragForce dragForce = new DragForce(atmosphere, dragSensitive, false);
443 propagator.addForceModel(dragForce);
444
445
446 final ThirdBodyAttraction moonAttraction = new ThirdBodyAttraction(moon);
447 final ThirdBodyAttraction sunAttraction = new ThirdBodyAttraction(sun);
448
449 propagator.addForceModel(moonAttraction);
450 propagator.addForceModel(sunAttraction);
451
452
453 final SolidTides solidTides = new SolidTides(earthShape.getBodyFrame(), gravity.getAe(), gravity.getMu(),
454 gravity.getTideSystem(), conventions,
455 TimeScalesFactory.getUT1(conventions, true),
456 moon, sun);
457 propagator.addForceModel(solidTides);
458
459
460 final OceanTides oceanTides = new OceanTides(earthShape.getBodyFrame(), gravity.getAe(), gravity.getMu(),
461 4, 4, conventions,
462 TimeScalesFactory.getUT1(conventions, true));
463 propagator.addForceModel(oceanTides);
464
465
466 final Relativity relativity = new Relativity(Constants.IERS2010_EARTH_MU);
467 propagator.addForceModel(relativity);
468
469 }
470
471 public static SpacecraftState generateSergeiReferenceState() {
472 final Orbit orbit = generateSergeiReferenceOrbit();
473 final double sergeiMass = 0.04;
474 return new SpacecraftState(orbit).withMass(sergeiMass);
475 }
476
477 public static Orbit generateSergeiReferenceOrbit() {
478
479
480 final AbsoluteDate sergeiDate = new AbsoluteDate(2003, 11, 22, 19, 0, 0, TimeScalesFactory.getUTC());
481
482 return new CartesianOrbit(new PVCoordinates(
483 new Vector3D(-2397200, 4217850, 5317450),
484 new Vector3D(-1303.9, 5558.9, -4839.6)), FramesFactory.getGCRF(), sergeiDate, Constants.EIGEN5C_EARTH_MU);
485 }
486
487 public static RealMatrix generateSergeiCovarianceMatrix() {
488
489
490 final double sigmaPosX = 98676;
491 final double sigmaPosY = 420547;
492 final double sigmaPosZ = 366438;
493
494
495 final double sigmaVelX = 194;
496 final double sigmaVelY = 341;
497 final double sigmaVelZ = 430;
498
499
500 final double posXYCorr = -0.999985;
501 final double posXZCorr = 0.999983;
502 final double posYZCorr = -0.999997;
503
504
505 final double velXYCorr = -0.999998;
506 final double velXZCorr = -0.999997;
507 final double velYZCorr = 0.999997;
508
509
510 final double posXVelXCorr = -0.999982;
511 final double posXVelYCorr = 0.999989;
512 final double posXVelZCorr = 0.999983;
513
514 final double posYVelXCorr = 0.999997;
515 final double posYVelYCorr = -0.999999;
516 final double posYVelZCorr = -0.999995;
517
518 final double posZVelXCorr = -0.999996;
519 final double posZVelYCorr = 0.999996;
520 final double posZVelZCorr = 0.999999;
521
522 return new BlockRealMatrix(new double[][] {
523 { sigmaPosX * sigmaPosX, sigmaPosX * sigmaPosY * posXYCorr, sigmaPosX * sigmaPosZ * posXZCorr,
524 sigmaPosX * sigmaVelX * posXVelXCorr, sigmaPosX * sigmaVelY * posXVelYCorr,
525 sigmaPosX * sigmaVelZ * posXVelZCorr },
526 { sigmaPosX * sigmaPosY * posXYCorr, sigmaPosY * sigmaPosY, sigmaPosY * sigmaPosZ * posYZCorr,
527 sigmaPosY * sigmaVelX * posYVelXCorr, sigmaPosY * sigmaVelY * posYVelYCorr,
528 sigmaPosY * sigmaVelZ * posYVelZCorr },
529 { sigmaPosX * sigmaPosZ * posXZCorr, sigmaPosY * sigmaPosZ * posYZCorr, sigmaPosZ * sigmaPosZ,
530 sigmaPosZ * sigmaVelX * posZVelXCorr, sigmaPosZ * sigmaVelY * posZVelYCorr,
531 sigmaPosZ * sigmaVelZ * posZVelZCorr },
532 { sigmaPosX * sigmaVelX * posXVelXCorr, sigmaPosY * sigmaVelX * posYVelXCorr,
533 sigmaPosZ * sigmaVelX * posZVelXCorr, sigmaVelX * sigmaVelX, sigmaVelX * sigmaVelY * velXYCorr,
534 sigmaVelX * sigmaVelZ * velXZCorr },
535 { sigmaPosX * sigmaVelY * posXVelYCorr, sigmaPosY * sigmaVelY * posYVelYCorr,
536 sigmaPosZ * sigmaVelY * posZVelYCorr, sigmaVelX * sigmaVelY * velXYCorr, sigmaVelY * sigmaVelY,
537 sigmaVelY * sigmaVelZ * velYZCorr },
538 { sigmaPosX * sigmaVelZ * posXVelZCorr, sigmaPosY * sigmaVelZ * posYVelZCorr,
539 sigmaPosZ * sigmaVelZ * posZVelZCorr, sigmaVelX * sigmaVelZ * velXZCorr,
540 sigmaVelY * sigmaVelZ * velYZCorr, sigmaVelZ * sigmaVelZ } });
541 }
542
543 private void doTestInterpolation(final TimeInterpolator<SpacecraftState> stateInterpolator,
544 final TimeInterpolator<TimeStampedPair<Orbit, StateCovariance>> covarianceInterpolator,
545 final double propagationHorizon,
546 final double tabulatedTimeStep,
547 final double expectedMeanRMSPositionError,
548 final double expectedMeanRMSVelocityError,
549 final double expectedMedianRMSPositionError,
550 final double expectedMedianRMSVelocityError,
551 final double expectedMaxRMSPositionError,
552 final double expectedMaxRMSVelocityError,
553 final double TOLERANCE,
554 final boolean showResults) {
555 final DescriptiveStatistics[] relativeRMSSigmaError =
556 computeStatisticsCovarianceInterpolationOnSergeiCase(propagationHorizon, tabulatedTimeStep,
557 stateInterpolator, covarianceInterpolator);
558
559
560 if (showResults) {
561
562 final int nDigits = (int) FastMath.log10(1/TOLERANCE);
563 final String fmt = String.format(Locale.US, "%%35s = %%%02d.%02df%%n", nDigits + 4, nDigits);
564
565 System.out.format(Locale.US, fmt, "relativeRMSSigmaError[0].getMean", relativeRMSSigmaError[0].getMean());
566 System.out.format(Locale.US, fmt, "relativeRMSSigmaError[1].getMean", relativeRMSSigmaError[1].getMean());
567 System.out.format(Locale.US, fmt, "relativeRMSSigmaError[0].getMedian", relativeRMSSigmaError[0].getPercentile(50));
568 System.out.format(Locale.US, fmt, "relativeRMSSigmaError[1].getMedian", relativeRMSSigmaError[1].getPercentile(50));
569 System.out.format(Locale.US, fmt, "relativeRMSSigmaError[0].getMax", relativeRMSSigmaError[0].getMax());
570 System.out.format(Locale.US, fmt, "relativeRMSSigmaError[1].getMax", relativeRMSSigmaError[1].getMax());
571
572 }
573
574 Assertions.assertEquals(expectedMeanRMSPositionError, relativeRMSSigmaError[0].getMean(), TOLERANCE);
575 Assertions.assertEquals(expectedMeanRMSVelocityError, relativeRMSSigmaError[1].getMean(), TOLERANCE);
576 Assertions.assertEquals(expectedMedianRMSPositionError, relativeRMSSigmaError[0].getPercentile(50), TOLERANCE);
577 Assertions.assertEquals(expectedMedianRMSVelocityError, relativeRMSSigmaError[1].getPercentile(50), TOLERANCE);
578 Assertions.assertEquals(expectedMaxRMSPositionError, relativeRMSSigmaError[0].getMax(), TOLERANCE);
579 Assertions.assertEquals(expectedMaxRMSVelocityError, relativeRMSSigmaError[1].getMax(), TOLERANCE);
580 }
581
582
583
584
585
586
587
588
589
590
591
592
593
594 @Test
595 @DisplayName("test quintic Keplerian interpolation (two time derivatives used) on full force model test case from : "
596 + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
597 + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
598 void testQuinticKeplerianInterpolation() {
599
600
601 final boolean showResults = false;
602
603
604 final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
605
606 final TimeInterpolator<Orbit> orbitInterpolator =
607 new OrbitBlender(blendingFunction, new KeplerianPropagator(sergeiOrbit), sergeiFrame);
608
609 final StateCovarianceKeplerianHermiteInterpolator covarianceInterpolator =
610 new StateCovarianceKeplerianHermiteInterpolator(orbitInterpolator, sergeiFrame, OrbitType.CARTESIAN,
611 PositionAngleType.MEAN);
612
613
614 final TimeInterpolator<SpacecraftState> stateInterpolator =
615 new SpacecraftStateInterpolator(AbstractTimeInterpolator.DEFAULT_INTERPOLATION_POINTS, AbstractTimeInterpolator.DEFAULT_EXTRAPOLATION_THRESHOLD_SEC,
616 sergeiFrame, orbitInterpolator,
617 null, null, null, null);
618
619
620 doTestInterpolation(stateInterpolator, covarianceInterpolator,
621 DEFAULT_SERGEI_PROPAGATION_TIME, DEFAULT_SERGEI_TABULATED_TIMESTEP,
622 0.069908803,
623 0.207658157,
624 0.068772221,
625 0.207587317,
626 0.185452386,
627 0.609433792,
628 TOLERANCE,
629 showResults);
630
631
632
633
634
635
636
637
638
639
640 Assertions.assertEquals(CartesianDerivativesFilter.USE_PVA, covarianceInterpolator.getFilter());
641
642 }
643
644
645
646
647
648
649
650
651 @Test
652 @DisplayName("test cubic Keplerian interpolation (first time derivative used) on full force model test case from : "
653 + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
654 + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
655 void testCubicKeplerianInterpolation() {
656
657
658 final boolean showResults = false;
659
660
661 final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
662
663 final TimeInterpolator<Orbit> orbitInterpolator = new OrbitBlender(blendingFunction,
664 new KeplerianPropagator(sergeiOrbit),
665 sergeiFrame);
666
667 final StateCovarianceKeplerianHermiteInterpolator covarianceInterpolator =
668 new StateCovarianceKeplerianHermiteInterpolator(DEFAULT_SERGEI_INTERPOLATION_POINTS, orbitInterpolator,
669 CartesianDerivativesFilter.USE_PV, sergeiFrame,
670 OrbitType.CARTESIAN,
671 PositionAngleType.MEAN);
672
673
674 final TimeInterpolator<SpacecraftState> stateInterpolator =
675 new SpacecraftStateInterpolator(AbstractTimeInterpolator.DEFAULT_INTERPOLATION_POINTS, AbstractTimeInterpolator.DEFAULT_EXTRAPOLATION_THRESHOLD_SEC,
676 sergeiFrame, orbitInterpolator, null, null, null, null);
677
678
679 doTestInterpolation(stateInterpolator, covarianceInterpolator,
680 DEFAULT_SERGEI_PROPAGATION_TIME, DEFAULT_SERGEI_TABULATED_TIMESTEP,
681 0.099905557,
682 0.201661811,
683 0.107651204,
684 0.176451233,
685 0.211682516,
686 0.639107824,
687 TOLERANCE,
688 showResults);
689
690
691
692
693
694
695
696
697
698 Assertions.assertEquals(CartesianDerivativesFilter.USE_PV, covarianceInterpolator.getFilter());
699
700 }
701
702
703
704
705
706
707
708
709 @Test
710 @DisplayName("test linear Keplerian interpolation (no derivatives used) on full force model test case from : "
711 + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
712 + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
713 void testLinearKeplerianInterpolation() {
714
715
716 final boolean showResults = false;
717
718
719 final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
720
721 final TimeInterpolator<Orbit> orbitInterpolator = new OrbitBlender(blendingFunction,
722 new KeplerianPropagator(sergeiOrbit),
723 sergeiFrame);
724
725 final StateCovarianceKeplerianHermiteInterpolator covarianceInterpolator =
726 new StateCovarianceKeplerianHermiteInterpolator(DEFAULT_SERGEI_INTERPOLATION_POINTS, orbitInterpolator,
727 CartesianDerivativesFilter.USE_P, sergeiFrame,
728 OrbitType.CARTESIAN,
729 PositionAngleType.MEAN);
730
731
732 final TimeInterpolator<SpacecraftState> stateInterpolator =
733 new SpacecraftStateInterpolator(AbstractTimeInterpolator.DEFAULT_INTERPOLATION_POINTS, AbstractTimeInterpolator.DEFAULT_EXTRAPOLATION_THRESHOLD_SEC,
734 sergeiFrame, orbitInterpolator, null, null, null, null);
735
736
737 doTestInterpolation(stateInterpolator, covarianceInterpolator,
738 DEFAULT_SERGEI_PROPAGATION_TIME, DEFAULT_SERGEI_TABULATED_TIMESTEP,
739 0.315401052,
740 0.277053458,
741 0.359035707,
742 0.237517672,
743 0.493068781,
744 0.829556363,
745 TOLERANCE,
746 showResults);
747
748
749
750
751
752
753
754
755
756 Assertions.assertEquals(CartesianDerivativesFilter.USE_P, covarianceInterpolator.getFilter());
757
758 }
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774 @Test
775 @DisplayName("test quintic Keplerian interpolation (output in LOF) on full force model test case from : "
776 + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
777 + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
778 void testLOFInterpolation() {
779
780
781 final boolean showResults = false;
782
783
784 final Orbit orbit = generateSergeiReferenceOrbit();
785 final Frame frame = orbit.getFrame();
786
787
788 final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
789
790 final TimeInterpolator<Orbit> orbitInterpolator = new OrbitBlender(blendingFunction,
791 new KeplerianPropagator(orbit),
792 frame);
793
794 final LOFType DEFAULT_LOFTYPE = LOFType.TNW;
795 final TimeInterpolator<TimeStampedPair<Orbit, StateCovariance>> covarianceInterpolator =
796 new StateCovarianceKeplerianHermiteInterpolator(orbitInterpolator, DEFAULT_LOFTYPE);
797
798
799 final DescriptiveStatistics[] relativeRMSSigmaError =
800 computeStatisticsCovarianceLOFInterpolation(DEFAULT_SERGEI_PROPAGATION_TIME,
801 DEFAULT_SERGEI_TABULATED_TIMESTEP,
802 DEFAULT_LOFTYPE, covarianceInterpolator);
803
804
805 if (showResults) {
806 System.out.format(Locale.US, "%35s = %20.12f%n", "relativeRMSSigmaError[0].getMean", relativeRMSSigmaError[0].getMean());
807 System.out.format(Locale.US, "%35s = %20.12f%n", "relativeRMSSigmaError[1].getMean", relativeRMSSigmaError[1].getMean());
808 System.out.format(Locale.US, "%35s = %20.12f%n", "relativeRMSSigmaError[0].getMedian", relativeRMSSigmaError[0].getPercentile(50));
809 System.out.format(Locale.US, "%35s = %20.12f%n", "relativeRMSSigmaError[1].getMedian", relativeRMSSigmaError[1].getPercentile(50));
810 System.out.format(Locale.US, "%35s = %20.12f%n", "relativeRMSSigmaError[0].getMax", relativeRMSSigmaError[0].getMax());
811 System.out.format(Locale.US, "%35s = %20.12f%n", "relativeRMSSigmaError[1].getMax", relativeRMSSigmaError[1].getMax());
812
813 }
814 Assertions.assertEquals( 0.081480745610, relativeRMSSigmaError[0].getMean(), TOLERANCE);
815 Assertions.assertEquals(19.410354129217, relativeRMSSigmaError[1].getMean(), TOLERANCE);
816 Assertions.assertEquals( 0.084585943748, relativeRMSSigmaError[0].getPercentile(50), TOLERANCE);
817 Assertions.assertEquals(14.968739949999694, relativeRMSSigmaError[1].getPercentile(50), TOLERANCE);
818 Assertions.assertEquals( 0.165222304428, relativeRMSSigmaError[0].getMax(), TOLERANCE);
819 Assertions.assertEquals(75.48232040417717, relativeRMSSigmaError[1].getMax(), 3 * TOLERANCE);
820 }
821
822 @Test
823 @SuppressWarnings("unchecked")
824 void testIssue1521() {
825
826
827
828 final TimeInterpolator<Orbit> orbitInterpolator = Mockito.mock(TimeInterpolator.class);
829
830
831
832 final StateCovarianceKeplerianHermiteInterpolator interpolator =
833 new StateCovarianceKeplerianHermiteInterpolator(orbitInterpolator, LOFType.NTW);
834
835
836 final List<TimeInterpolator<?>> subInterpolators = interpolator.getSubInterpolators();
837
838
839 Assertions.assertEquals(1, subInterpolators.size());
840 Assertions.assertEquals(orbitInterpolator, subInterpolators.get(0));
841 }
842 }