1   /* Copyright 2002-2025 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.propagation;
18  
19  import java.util.Locale;
20  import java.util.Map;
21  
22  import org.hipparchus.analysis.polynomials.SmoothStepFactory;
23  import org.hipparchus.stat.descriptive.DescriptiveStatistics;
24  import org.hipparchus.util.FastMath;
25  import org.junit.jupiter.api.Assertions;
26  import org.junit.jupiter.api.BeforeAll;
27  import org.junit.jupiter.api.DisplayName;
28  import org.junit.jupiter.api.Test;
29  import org.orekit.Utils;
30  import org.orekit.data.DataContext;
31  import org.orekit.forces.gravity.potential.AstronomicalAmplitudeReader;
32  import org.orekit.forces.gravity.potential.EGMFormatReader;
33  import org.orekit.forces.gravity.potential.FESCHatEpsilonReader;
34  import org.orekit.forces.gravity.potential.GravityFieldFactory;
35  import org.orekit.forces.gravity.potential.OceanLoadDeformationCoefficients;
36  import org.orekit.frames.Frame;
37  import org.orekit.frames.LOFType;
38  import org.orekit.orbits.Orbit;
39  import org.orekit.orbits.OrbitBlender;
40  import org.orekit.orbits.OrbitType;
41  import org.orekit.orbits.PositionAngleType;
42  import org.orekit.propagation.analytical.AbstractAnalyticalPropagator;
43  import org.orekit.propagation.analytical.BrouwerLyddanePropagator;
44  import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
45  import org.orekit.propagation.analytical.KeplerianPropagator;
46  import org.orekit.time.TimeInterpolator;
47  import org.orekit.time.TimeStampedPair;
48  import org.orekit.utils.Constants;
49  
50  class StateCovarianceBlenderTest {
51  
52      private static SpacecraftState sergeiState;
53      private static Orbit           sergeiOrbit;
54      private static Frame           sergeiFrame;
55  
56      // Constants
57      private final double DEFAULT_SERGEI_PROPAGATION_TIME   = 2400;
58      private final double DEFAUTL_SERGEI_TABULATED_TIMESTEP = 2400;
59  
60      @BeforeAll
61      public static void setUp() {
62          Utils.setDataRoot("regular-data:potential/egm-format:atmosphere:tides:regular-data/de405-ephemerides");
63          GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("EGM96-truncated-21x21", true));
64          AstronomicalAmplitudeReader aaReader =
65                          new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
66          DataContext.getDefault().getDataProvidersManager().feed(aaReader.getSupportedNames(), aaReader);
67          Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
68          GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
69                                                                           0.01, FastMath.toRadians(1.0),
70                                                                           OceanLoadDeformationCoefficients.IERS_2010,
71                                                                           map));
72          // Default orbit case
73          sergeiState = StateCovarianceKeplerianHermiteInterpolatorTest.generateSergeiReferenceState();
74          sergeiOrbit = sergeiState.getOrbit();
75          sergeiFrame = sergeiOrbit.getFrame();
76      }
77  
78      private void doTestBlending(final double propagationHorizon, final double tabulatedTimeStep,
79                                  final SmoothStepFactory.SmoothStepFunction blendingFunction,
80                                  final AbstractAnalyticalPropagator analyticalPropagator,
81                                  final double expectedMeanRMSPositionError,
82                                  final double expectedMeanRMSVelocityError,
83                                  final double expectedMedianRMSPositionError,
84                                  final double expectedMedianRMSVelocityError,
85                                  final double expectedMaxRMSPositionError,
86                                  final double expectedMaxRMSVelocityError,
87                                  final double tolerance,
88                                  final boolean showResults) {
89          final TimeInterpolator<Orbit> orbitInterpolator = new OrbitBlender(blendingFunction,
90                                                                             analyticalPropagator,
91                                                                             sergeiFrame);
92  
93          final TimeInterpolator<TimeStampedPair<Orbit, StateCovariance>> covarianceInterpolator =
94                          new StateCovarianceBlender(blendingFunction, orbitInterpolator,
95                                                     sergeiFrame, OrbitType.CARTESIAN,
96                                                     PositionAngleType.MEAN);
97  
98          // Create state interpolator
99          final TimeInterpolator<SpacecraftState> stateInterpolator =
100                         new SpacecraftStateInterpolator(2, 1.0e-3, sergeiFrame, orbitInterpolator, null, null, null, null);
101 
102         // When
103         final DescriptiveStatistics[] relativeRMSSigmaError =
104                         StateCovarianceKeplerianHermiteInterpolatorTest.computeStatisticsCovarianceInterpolationOnSergeiCase(
105                                                                                                                              propagationHorizon, tabulatedTimeStep, stateInterpolator, covarianceInterpolator);
106 
107         // Then
108         if (showResults) {
109             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMean", relativeRMSSigmaError[0].getMean());
110             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMean", relativeRMSSigmaError[1].getMean());
111             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMedian", relativeRMSSigmaError[0].getPercentile(50));
112             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMedian", relativeRMSSigmaError[1].getPercentile(50));
113             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMax", relativeRMSSigmaError[0].getMax());
114             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMax", relativeRMSSigmaError[1].getMax());
115 
116         }
117 
118         // Results obtained when using modified orbit date to use truncated JPL test resource file
119         Assertions.assertEquals(expectedMeanRMSPositionError, relativeRMSSigmaError[0].getMean(), tolerance);
120         Assertions.assertEquals(expectedMeanRMSVelocityError, relativeRMSSigmaError[1].getMean(), tolerance);
121         Assertions.assertEquals(expectedMedianRMSPositionError, relativeRMSSigmaError[0].getPercentile(50), tolerance);
122         Assertions.assertEquals(expectedMedianRMSVelocityError, relativeRMSSigmaError[1].getPercentile(50), tolerance);
123         Assertions.assertEquals(expectedMaxRMSPositionError, relativeRMSSigmaError[0].getMax(), tolerance);
124         Assertions.assertEquals(expectedMaxRMSVelocityError, relativeRMSSigmaError[1].getMax(), tolerance);
125     }
126 
127     /**
128      * Test based on the full force model test case from TANYGIN, Sergei. Efficient covariance interpolation using blending
129      * of approximate covariance propagations. The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.
130      * <p>
131      * However, note that the exact result is not known and only an approximated value can be deduced from the available
132      * graph in the aforementioned paper. That is why it is both a validation test and a non regression test.
133      * <p>
134      * This instance of the test aims to test the Quadratic Blending method with standard Keplerian dynamics.
135      */
136     @Test
137     @DisplayName("test Keplerian quadratic blending interpolation on full force model test case from: "
138                     + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
139                     + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
140     void testKeplerianQuadraticBlending() {
141         // Given
142         final boolean showResults = false; // Show results?
143         final double tolerance = 1.e-11;
144 
145         // Create state covariance interpolator
146         final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
147 
148         // When & Then
149         doTestBlending(DEFAULT_SERGEI_PROPAGATION_TIME, DEFAUTL_SERGEI_TABULATED_TIMESTEP, blendingFunction,
150                        new KeplerianPropagator(sergeiOrbit),
151                        0.11333019740,
152                        0.23518823987,
153                        0.11116079511,
154                        0.26216208074,
155                        0.268782359940,
156                        0.402221757858,
157                        tolerance,
158                        showResults);
159 
160         // Results obtained when using Sergei reference date
161         /*        Assertions.assertEquals(0.12086964077199346, relativeRMSSigmaError[0].getMean(), 1e-17);
162         Assertions.assertEquals(0.21871712884727984, relativeRMSSigmaError[1].getMean(), 1e-17);
163         Assertions.assertEquals(0.12172900928265865, relativeRMSSigmaError[0].getPercentile(50), 1e-17);
164         Assertions.assertEquals(0.24456923477457276, relativeRMSSigmaError[1].getPercentile(50), 1e-17);
165         Assertions.assertEquals(0.24004509397020612, relativeRMSSigmaError[0].getMax(), 1e-17);
166         Assertions.assertEquals(0.3375639988469634, relativeRMSSigmaError[1].getMax(), 1e-16);*/
167     }
168 
169     /**
170      * Test based on the full force model test case from TANYGIN, Sergei. Efficient covariance interpolation using blending
171      * of approximate covariance propagations. The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.
172      * <p>
173      * This instance of the test aims to test the Quadratic Blending method using Brouwer Lyddane propagation. This is a
174      * non-regression test.
175      */
176     @Test
177     @DisplayName("test Brouwer Lyddane quadratic blending interpolation on full force model test case from: "
178                     + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
179                     + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
180     void testBrouwerLyddaneQuadraticBlending() {
181         // Given
182         final boolean showResults = false; // Show results?
183         final double tolerance = 1.e-11;
184 
185         // Create state covariance interpolator
186         final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
187 
188         final AbstractAnalyticalPropagator propagator = new BrouwerLyddanePropagator(sergeiOrbit,
189                                                                                      sergeiState.getMass(),
190                                                                                      Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
191                                                                                      Constants.EIGEN5C_EARTH_MU,
192                                                                                      Constants.EIGEN5C_EARTH_C20,
193                                                                                      Constants.EIGEN5C_EARTH_C30,
194                                                                                      Constants.EIGEN5C_EARTH_C40,
195                                                                                      Constants.EIGEN5C_EARTH_C50,
196                                                                                      0);
197 
198         // When & Then
199         doTestBlending(DEFAULT_SERGEI_PROPAGATION_TIME, DEFAUTL_SERGEI_TABULATED_TIMESTEP, blendingFunction,
200                        propagator,
201                        0.13366703460,
202                        0.13160856461,
203                        0.14205498786,
204                        0.13656689649,
205                        0.21941543409,
206                        0.25091794341,
207                        tolerance,
208                        showResults);
209 
210         // Results obtained when using Sergei reference date
211         /*        Assertions.assertEquals(0.07645785479359624, relativeRMSSigmaError[0].getMean(), 1e-17);
212         Assertions.assertEquals(0.17941792898602038, relativeRMSSigmaError[1].getMean(), 1e-17);
213         Assertions.assertEquals(0.08259069655149026, relativeRMSSigmaError[0].getPercentile(50), 1e-17);
214         Assertions.assertEquals(0.18352413417267063, relativeRMSSigmaError[1].getPercentile(50), 1e-17);
215         Assertions.assertEquals(0.13164670404592496, relativeRMSSigmaError[0].getMax(), 1e-17);
216         Assertions.assertEquals(0.32564919981018114, relativeRMSSigmaError[1].getMax(), 1e-16);*/
217     }
218 
219     /**
220      * Test based on the full force model test case from TANYGIN, Sergei. Efficient covariance interpolation using blending
221      * of approximate covariance propagations. The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.
222      * <p>
223      * This instance of the test aims to test the Quadratic Blending method using Brouwer Lyddane propagation. This is a
224      * non-regression test.
225      */
226     @Deprecated
227     @Test
228     @DisplayName("test Brouwer Lyddane quadratic blending interpolation on full force model test case from: "
229                     + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
230                     + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
231     void testBrouwerLyddaneQuadraticBlendingDeprecated() {
232         // Given
233         final boolean showResults = false; // Show results?
234         final double tolerance = 1.e-11;
235 
236         // Create state covariance interpolator
237         final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
238 
239         final AbstractAnalyticalPropagator propagator = new BrouwerLyddanePropagator(sergeiOrbit,
240                                                                                      sergeiState.getMass(),
241                                                                                      Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
242                                                                                      Constants.EIGEN5C_EARTH_MU,
243                                                                                      Constants.EIGEN5C_EARTH_C20,
244                                                                                      Constants.EIGEN5C_EARTH_C30,
245                                                                                      Constants.EIGEN5C_EARTH_C40,
246                                                                                      Constants.EIGEN5C_EARTH_C50,
247                                                                                      0);
248 
249         // When & Then
250         doTestBlending(DEFAULT_SERGEI_PROPAGATION_TIME, DEFAUTL_SERGEI_TABULATED_TIMESTEP, blendingFunction,
251                        propagator,
252                        0.13366703460,
253                        0.13160856461,
254                        0.14205498786,
255                        0.13656689649,
256                        0.21941543409,
257                        0.25091794341,
258                        tolerance,
259                        showResults);
260 
261         // Results obtained when using Sergei reference date
262         /*        Assertions.assertEquals(0.07645785479359624, relativeRMSSigmaError[0].getMean(), 1e-17);
263         Assertions.assertEquals(0.17941792898602038, relativeRMSSigmaError[1].getMean(), 1e-17);
264         Assertions.assertEquals(0.08259069655149026, relativeRMSSigmaError[0].getPercentile(50), 1e-17);
265         Assertions.assertEquals(0.18352413417267063, relativeRMSSigmaError[1].getPercentile(50), 1e-17);
266         Assertions.assertEquals(0.13164670404592496, relativeRMSSigmaError[0].getMax(), 1e-17);
267         Assertions.assertEquals(0.32564919981018114, relativeRMSSigmaError[1].getMax(), 1e-16);*/
268     }
269 
270     /**
271      * Test based on the full force model test case from TANYGIN, Sergei. Efficient covariance interpolation using blending
272      * of approximate covariance propagations. The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.
273      * <p>
274      * This instance of the test aims to test the Quadratic Blending method using Eckstein Hechler propagation. This is a
275      * non-regression test.
276      */
277     @Test
278     @DisplayName("test Ekstein Hechler quadratic blending interpolation on full force model test case from: "
279                     + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
280                     + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
281     void testEksteinHechlerQuadraticBlending() {
282         // Given
283         final boolean showResults = false; // Show results?
284         final double tolerance = 1.e-11;
285 
286         // Create state covariance interpolator
287         final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
288 
289         final AbstractAnalyticalPropagator propagator = new EcksteinHechlerPropagator(sergeiOrbit,
290                                                                                       sergeiState.getMass(),
291                                                                                       Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
292                                                                                       Constants.EIGEN5C_EARTH_MU,
293                                                                                       Constants.EIGEN5C_EARTH_C20,
294                                                                                       Constants.EIGEN5C_EARTH_C30,
295                                                                                       Constants.EIGEN5C_EARTH_C40,
296                                                                                       Constants.EIGEN5C_EARTH_C50,
297                                                                                       Constants.EIGEN5C_EARTH_C60);
298 
299         // When & Then
300         doTestBlending(DEFAULT_SERGEI_PROPAGATION_TIME, DEFAUTL_SERGEI_TABULATED_TIMESTEP, blendingFunction,
301                        propagator,
302                        0.092022772744,
303                        0.175328981237,
304                        0.085754375985,
305                        0.193319970362,
306                        0.169348422233,
307                        0.347302066023,
308                        tolerance,
309                        showResults);
310 
311         // Results obtained when using Sergei reference date
312         /*        Assertions.assertEquals(0.08688580997030507, relativeRMSSigmaError[0].getMean(), 1e-17);
313         Assertions.assertEquals(0.15630296926592135, relativeRMSSigmaError[1].getMean(), 1e-17);
314         Assertions.assertEquals(0.07771211628338093, relativeRMSSigmaError[0].getPercentile(50), 1e-17);
315         Assertions.assertEquals(0.17056867629454775, relativeRMSSigmaError[1].getPercentile(50), 1e-17);
316         Assertions.assertEquals(0.16747846953637413, relativeRMSSigmaError[0].getMax(), 1e-17);
317         Assertions.assertEquals(0.26557409638784585, relativeRMSSigmaError[1].getMax(), 1e-16);*/
318     }
319 
320     /**
321      * Test based on the full force model test case from TANYGIN, Sergei. Efficient covariance interpolation using blending
322      * of approximate covariance propagations. The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.
323      * <p>
324      * This instance of the test is a non regression test aiming to test the results obtained when interpolating in a local
325      * orbital frame.
326      * <p>
327      * Bad results are seen regarding velocity interpolation. This has been checked to be independent of the method used to
328      * interpolate the state covariance. Moreover, the same results are achieved if we interpolate in an inertial frame
329      * (which has been seen to give very good results as in {@link #testKeplerianQuadraticBlending()}),and then express it in
330      * a non-inertial local orbital frame. It has also been verified that it was not (mainly) due to errors in orbit
331      * interpolation as even using exact orbit for frame conversion would only slightly improve the results (<0.1%
332      * improvements). Hence, the only explanation found is a sensitivity issue linked to non-inertial local orbital frame.
333      */
334     @Test
335     @DisplayName("test Keplerian quadratic blending expressed in LOF on full force model test case from : "
336                     + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
337                     + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
338     void testLOFKeplerianBlending() {
339         // Given
340         final boolean showResults = false; // Show results?
341         final double tolerance = 1.e-9;
342 
343         // Create state covariance interpolator
344         final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
345 
346         final TimeInterpolator<Orbit> orbitInterpolator = new OrbitBlender(blendingFunction,
347                                                                            new KeplerianPropagator(sergeiOrbit),
348                                                                            sergeiFrame);
349 
350         final LOFType DEFAULT_LOFTYPE = LOFType.TNW;
351         final AbstractStateCovarianceInterpolator covarianceInterpolator =
352                         new StateCovarianceBlender(blendingFunction, orbitInterpolator, DEFAULT_LOFTYPE);
353 
354         // When
355         final DescriptiveStatistics[] relativeRMSSigmaError =
356                         StateCovarianceKeplerianHermiteInterpolatorTest.computeStatisticsCovarianceLOFInterpolation(
357                                                                                                                     DEFAULT_SERGEI_PROPAGATION_TIME, DEFAUTL_SERGEI_TABULATED_TIMESTEP,
358                                                                                                                     DEFAULT_LOFTYPE, covarianceInterpolator);
359 
360         // Then
361         if (showResults) {
362             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMean", relativeRMSSigmaError[0].getMean());
363             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMean", relativeRMSSigmaError[1].getMean());
364             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMedian", relativeRMSSigmaError[0].getPercentile(50));
365             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMedian", relativeRMSSigmaError[1].getPercentile(50));
366             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMax", relativeRMSSigmaError[0].getMax());
367             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMax", relativeRMSSigmaError[1].getMax());
368 
369         }
370 
371         // Results obtained when using modified orbit date to use truncated JPL test resource file
372         Assertions.assertEquals( 0.1190324127, relativeRMSSigmaError[0].getMean(), tolerance);
373         Assertions.assertEquals( 19.9401863789, relativeRMSSigmaError[1].getMean(), tolerance);
374         Assertions.assertEquals( 0.1221432731, relativeRMSSigmaError[0].getPercentile(50), tolerance);
375         Assertions.assertEquals( 14.0023883813, relativeRMSSigmaError[1].getPercentile(50), tolerance);
376         Assertions.assertEquals( 0.2282143786, relativeRMSSigmaError[0].getMax(), tolerance);
377         Assertions.assertEquals(99.776271400, relativeRMSSigmaError[1].getMax(), tolerance);
378 
379         // Assert getters as well
380         Assertions.assertNull(covarianceInterpolator.getOutFrame());
381         Assertions.assertEquals(DEFAULT_LOFTYPE, covarianceInterpolator.getOutLOF());
382         Assertions.assertEquals(OrbitType.CARTESIAN, covarianceInterpolator.getOutOrbitType());
383         Assertions.assertEquals(PositionAngleType.MEAN, covarianceInterpolator.getOutPositionAngleType());
384         Assertions.assertEquals(orbitInterpolator, covarianceInterpolator.getOrbitInterpolator());
385 
386     }
387 
388 }