1   /* Copyright 2002-2023 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(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-16;
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.1133302045524815,
152                        0.2351882330196802,
153                        0.1111607965207670,
154                        0.2621621038900783,
155                        0.2687823732965123,
156                        0.4022217682474207,
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-16;
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.0823442943514750,
202                        0.1976467714037895,
203                        0.0888590153696339,
204                        0.2020905044160413,
205                        0.1446159909426887,
206                        0.3894610461790040,
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 Eckstein Hechler propagation. This is a
224      * non-regression test.
225      */
226     @Test
227     @DisplayName("test Ekstein Hechler quadratic blending interpolation on full force model test case from: "
228                     + "TANYGIN, Sergei. Efficient covariance interpolation using blending of approximate covariance propagations. "
229                     + "The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.")
230     void testEksteinHechlerQuadraticBlending() {
231         // Given
232         final boolean showResults = false; // Show results?
233         final double tolerance = 1.e-16;
234 
235         // Create state covariance interpolator
236         final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
237 
238         final AbstractAnalyticalPropagator propagator = new EcksteinHechlerPropagator(sergeiOrbit,
239                                                                                       sergeiState.getMass(),
240                                                                                       Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
241                                                                                       Constants.EIGEN5C_EARTH_MU,
242                                                                                       Constants.EIGEN5C_EARTH_C20,
243                                                                                       Constants.EIGEN5C_EARTH_C30,
244                                                                                       Constants.EIGEN5C_EARTH_C40,
245                                                                                       Constants.EIGEN5C_EARTH_C50,
246                                                                                       Constants.EIGEN5C_EARTH_C60);
247 
248         // When & Then
249         doTestBlending(DEFAULT_SERGEI_PROPAGATION_TIME, DEFAUTL_SERGEI_TABULATED_TIMESTEP, blendingFunction,
250                        propagator,
251                        0.0920227696921130,
252                        0.1753289773020222,
253                        0.0857543808467324,
254                        0.1933199872752203,
255                        0.1693483908858563,
256                        0.3473020830532024,
257                        tolerance,
258                        showResults);
259 
260         // Results obtained when using Sergei reference date
261         /*        Assertions.assertEquals(0.08688580997030507, relativeRMSSigmaError[0].getMean(), 1e-17);
262         Assertions.assertEquals(0.15630296926592135, relativeRMSSigmaError[1].getMean(), 1e-17);
263         Assertions.assertEquals(0.07771211628338093, relativeRMSSigmaError[0].getPercentile(50), 1e-17);
264         Assertions.assertEquals(0.17056867629454775, relativeRMSSigmaError[1].getPercentile(50), 1e-17);
265         Assertions.assertEquals(0.16747846953637413, relativeRMSSigmaError[0].getMax(), 1e-17);
266         Assertions.assertEquals(0.26557409638784585, relativeRMSSigmaError[1].getMax(), 1e-16);*/
267     }
268 
269     /**
270      * Test based on the full force model test case from TANYGIN, Sergei. Efficient covariance interpolation using blending
271      * of approximate covariance propagations. The Journal of the Astronautical Sciences, 2014, vol. 61, no 1, p. 107-132.
272      * <p>
273      * This instance of the test is a non regression test aiming to test the results obtained when interpolating in a local
274      * orbital frame.
275      * <p>
276      * Bad results are seen regarding velocity interpolation. This has been checked to be independent of the method used to
277      * interpolate the state covariance. Moreover, the same results are achieved if we interpolate in an inertial frame
278      * (which has been seen to give very good results as in {@link #testKeplerianQuadraticBlending()}),and then express it in
279      * a non-inertial local orbital frame. It has also been verified that it was not (mainly) due to errors in orbit
280      * interpolation as even using exact orbit for frame conversion would only slightly improve the results (<0.1%
281      * improvements). Hence, the only explanation found is a sensitivity issue linked to non-inertial local orbital frame.
282      */
283     @Test
284     @DisplayName("test Keplerian quadratic blending expressed in LOF 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 testLOFKeplerianBlending() {
288         // Given
289         final boolean showResults = false; // Show results?
290         final double tolerance = 1.e-16;
291 
292         // Create state covariance interpolator
293         final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
294 
295         final TimeInterpolator<Orbit> orbitInterpolator = new OrbitBlender(blendingFunction,
296                                                                            new KeplerianPropagator(sergeiOrbit),
297                                                                            sergeiFrame);
298 
299         final LOFType DEFAULT_LOFTYPE = LOFType.TNW;
300         final AbstractStateCovarianceInterpolator covarianceInterpolator =
301                         new StateCovarianceBlender(blendingFunction, orbitInterpolator, DEFAULT_LOFTYPE);
302 
303         // When
304         final DescriptiveStatistics[] relativeRMSSigmaError =
305                         StateCovarianceKeplerianHermiteInterpolatorTest.computeStatisticsCovarianceLOFInterpolation(
306                                                                                                                     DEFAULT_SERGEI_PROPAGATION_TIME, DEFAUTL_SERGEI_TABULATED_TIMESTEP,
307                                                                                                                     DEFAULT_LOFTYPE, covarianceInterpolator);
308 
309         // Then
310         if (showResults) {
311             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMean", relativeRMSSigmaError[0].getMean());
312             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMean", relativeRMSSigmaError[1].getMean());
313             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMedian", relativeRMSSigmaError[0].getPercentile(50));
314             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMedian", relativeRMSSigmaError[1].getPercentile(50));
315             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[0].getMax", relativeRMSSigmaError[0].getMax());
316             System.out.format(Locale.US, "%35s = %20.16f%n", "relativeRMSSigmaError[1].getMax", relativeRMSSigmaError[1].getMax());
317 
318         }
319 
320         // Results obtained when using modified orbit date to use truncated JPL test resource file
321         Assertions.assertEquals( 0.1190324153130111, relativeRMSSigmaError[0].getMean(), tolerance);
322         Assertions.assertEquals( 7.2800694213600920, relativeRMSSigmaError[1].getMean(), tolerance);
323         Assertions.assertEquals( 0.1221432756730858, relativeRMSSigmaError[0].getPercentile(50), tolerance);
324         Assertions.assertEquals( 7.4882845444744826, relativeRMSSigmaError[1].getPercentile(50), tolerance);
325         Assertions.assertEquals( 0.2282143889000270, relativeRMSSigmaError[0].getMax(), tolerance);
326         Assertions.assertEquals(16.0468464408560370, relativeRMSSigmaError[1].getMax(), tolerance);
327 
328         // Assert getters as well
329         Assertions.assertNull(covarianceInterpolator.getOutFrame());
330         Assertions.assertEquals(DEFAULT_LOFTYPE, covarianceInterpolator.getOutLOF());
331         Assertions.assertEquals(OrbitType.CARTESIAN, covarianceInterpolator.getOutOrbitType());
332         Assertions.assertEquals(PositionAngleType.MEAN, covarianceInterpolator.getOutPositionAngleType());
333         Assertions.assertEquals(orbitInterpolator, covarianceInterpolator.getOrbitInterpolator());
334 
335     }
336 
337 }