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 org.hipparchus.geometry.euclidean.threed.Vector3D;
20  import org.hipparchus.linear.BlockRealMatrix;
21  import org.hipparchus.linear.MatrixUtils;
22  import org.hipparchus.linear.RealMatrix;
23  import org.hipparchus.ode.ODEIntegrator;
24  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
25  import org.hipparchus.util.FastMath;
26  import org.junit.jupiter.api.Assertions;
27  import org.junit.jupiter.api.DisplayName;
28  import org.junit.jupiter.api.Test;
29  import org.junit.jupiter.params.ParameterizedTest;
30  import org.junit.jupiter.params.provider.EnumSource;
31  import org.orekit.Utils;
32  import org.orekit.errors.OrekitException;
33  import org.orekit.forces.gravity.potential.GravityFieldFactory;
34  import org.orekit.forces.gravity.potential.ICGEMFormatReader;
35  import org.orekit.frames.Frame;
36  import org.orekit.frames.FramesFactory;
37  import org.orekit.frames.LOF;
38  import org.orekit.frames.LOFType;
39  import org.orekit.orbits.CartesianOrbit;
40  import org.orekit.orbits.EquinoctialOrbit;
41  import org.orekit.orbits.KeplerianOrbit;
42  import org.orekit.orbits.Orbit;
43  import org.orekit.orbits.OrbitType;
44  import org.orekit.orbits.PositionAngleType;
45  import org.orekit.propagation.numerical.NumericalPropagator;
46  import org.orekit.time.AbsoluteDate;
47  import org.orekit.time.TimeScalesFactory;
48  import org.orekit.utils.Constants;
49  import org.orekit.utils.IERSConventions;
50  import org.orekit.utils.PVCoordinates;
51  
52  import static org.junit.jupiter.api.Assertions.assertDoesNotThrow;
53  import static org.junit.jupiter.api.Assertions.assertEquals;
54  
55  public class StateCovarianceTest {
56  
57      final private static double          DEFAULT_VALLADO_THRESHOLD = 1e-6;
58      private       SpacecraftState initialState;
59      private       double[][]      initCov;
60  
61      /**
62       * Unit test for the covariance frame transformation.
63       */
64      @Test
65      void testFrameConversion() {
66  
67          // Initialization
68          setUp();
69  
70          // Define frames
71          final Frame frameA = FramesFactory.getEME2000();
72          final Frame frameB = FramesFactory.getTEME();
73  
74          // State covariance
75          final StateCovariance reference = new StateCovariance(MatrixUtils.createRealMatrix(initCov), initialState.getDate(), frameA, OrbitType.CARTESIAN, PositionAngleType.MEAN);
76  
77          // First transformation
78          StateCovariance transformed = reference.changeCovarianceFrame(initialState.getOrbit(), frameB);
79  
80          // Second transformation
81          transformed = transformed.changeCovarianceFrame(initialState.getOrbit(), frameA);
82  
83          // Verify
84          compareCovariance(reference.getMatrix(), transformed.getMatrix(), 5.2e-15);
85  
86      }
87  
88      @Test
89      void testConstructor() {
90          final AbsoluteDate initialDate          = new AbsoluteDate();
91          final Frame        initialInertialFrame = FramesFactory.getGCRF();
92      	final StateCovariance stateCovariance = new StateCovariance(getValladoInitialCovarianceMatrix(), initialDate, initialInertialFrame, OrbitType.CARTESIAN, PositionAngleType.MEAN);
93      	assertEquals(OrbitType.CARTESIAN, stateCovariance.getOrbitType());
94      	assertEquals(PositionAngleType.MEAN, stateCovariance.getPositionAngleType());
95      	assertEquals(initialInertialFrame, stateCovariance.getFrame());
96      	Assertions.assertNull(stateCovariance.getLOF());
97      	assertEquals(initialDate, stateCovariance.getDate());
98      }
99  
100     void setUp() {
101         Utils.setDataRoot("orbit-determination/february-2016:potential/icgem-format");
102         GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true));
103         Orbit initialOrbit = new CartesianOrbit(
104                 new PVCoordinates(new Vector3D(7526993.581890527, -9646310.10026971, 1464110.4928112086),
105                                   new Vector3D(3033.79456099698, 1715.265069098717, -4447.658745923895)),
106                 FramesFactory.getEME2000(),
107                 new AbsoluteDate("2016-02-13T16:00:00.000", TimeScalesFactory.getUTC()),
108                 Constants.WGS84_EARTH_MU);
109         initialState = new SpacecraftState(initialOrbit);
110         initCov = new double[][] {
111                 { 8.651816029e+01, 5.689987127e+01, -2.763870764e+01, -2.435617201e-02, 2.058274137e-02, -5.872883051e-03 },
112                 { 5.689987127e+01, 7.070624321e+01, 1.367120909e+01, -6.112622013e-03, 7.623626008e-03, -1.239413190e-02 },
113                 { -2.763870764e+01, 1.367120909e+01, 1.811858898e+02, 3.143798992e-02, -4.963106559e-02, -7.420114385e-04 },
114                 { -2.435617201e-02, -6.112622013e-03, 3.143798992e-02, 4.657077389e-05, 1.469943634e-05, 3.328475593e-05 },
115                 { 2.058274137e-02, 7.623626008e-03, -4.963106559e-02, 1.469943634e-05, 3.950715934e-05, 2.516044258e-05 },
116                 { -5.872883051e-03, -1.239413190e-02, -7.420114385e-04, 3.328475593e-05, 2.516044258e-05, 3.547466120e-05 }
117         };
118     }
119 
120     /**
121      * Compare two covariance matrices
122      *
123      * @param reference reference covariance
124      * @param computed  computed covariance
125      * @param threshold threshold for comparison
126      */
127     public static void compareCovariance(final RealMatrix reference, final RealMatrix computed, final double threshold) {
128         for (int row = 0; row < reference.getRowDimension(); row++) {
129             for (int column = 0; column < reference.getColumnDimension(); column++) {
130                 if (reference.getEntry(row, column) == 0) {
131                     assertEquals(reference.getEntry(row, column), computed.getEntry(row, column),
132                                             threshold);
133                 }
134                 else {
135                     assertEquals(reference.getEntry(row, column), computed.getEntry(row, column),
136                                             FastMath.abs(threshold * reference.getEntry(row, column)));
137                 }
138             }
139         }
140     }
141 
142     @Test
143     @DisplayName("Test conversion from inertial frame to RTN local orbital frame")
144     void should_return_same_covariance_matrix() {
145 
146         // Given
147         final AbsoluteDate initialDate          = new AbsoluteDate();
148         final Frame        initialInertialFrame = FramesFactory.getGCRF();
149         final double       mu                   = 398600e9;
150 
151         final PVCoordinates initialPV = new PVCoordinates(
152                 new Vector3D(6778000, 0, 0),
153                 new Vector3D(0, 7668.63, 0));
154 
155         final Orbit initialOrbit = new CartesianOrbit(initialPV, initialInertialFrame, initialDate, mu);
156 
157         final RealMatrix initialCovarianceInInertialFrame = new BlockRealMatrix(new double[][] {
158                 { 1, 0, 0, 0, 1e-5, 1e-4 },
159                 { 0, 1, 0, 0, 0, 1e-5 },
160                 { 0, 0, 1, 0, 0, 0 },
161                 { 0, 0, 0, 1e-3, 0, 0 },
162                 { 1e-5, 0, 0, 0, 1e-3, 0 },
163                 { 1e-4, 1e-5, 0, 0, 0, 1e-3 }
164         });
165 
166         final StateCovariance stateCovariance =
167                 new StateCovariance(initialCovarianceInInertialFrame, initialDate, initialInertialFrame,
168                         OrbitType.CARTESIAN, PositionAngleType.MEAN);
169         // When
170         final RealMatrix covarianceMatrixInRTN =
171                 stateCovariance.changeCovarianceFrame(initialOrbit, LOFType.QSW_INERTIAL).getMatrix();
172 
173         // Then
174         compareCovariance(initialCovarianceInInertialFrame, covarianceMatrixInRTN, DEFAULT_VALLADO_THRESHOLD);
175 
176     }
177 
178     /**
179      * Unit test for the covariance type transformation.
180      */
181     @Test
182     void testTypeConversion() {
183 
184         // Initialization
185         setUp();
186 
187         // Define orbit types
188         final OrbitType cart = OrbitType.CARTESIAN;
189         final OrbitType kep  = OrbitType.KEPLERIAN;
190 
191         // State covariance
192         final StateCovariance reference =
193                 new StateCovariance(MatrixUtils.createRealMatrix(initCov), initialState.getDate(), initialState.getFrame(),
194                         cart, PositionAngleType.MEAN);
195 
196         // First transformation
197         StateCovariance transformed = reference.changeCovarianceType(initialState.getOrbit(), kep, PositionAngleType.MEAN);
198 
199         // Second transformation
200         transformed = transformed.changeCovarianceType(initialState.getOrbit(), cart, PositionAngleType.MEAN);
201 
202         // Verify
203         compareCovariance(reference.getMatrix(), transformed.getMatrix(), 3.5e-12);
204 
205     }
206 
207     @Test
208     @DisplayName("Test covariance conversion from inertial frame to RTN local orbital frame")
209     void should_rotate_covariance_matrix_by_ninety_degrees() {
210 
211         // Given
212         final AbsoluteDate initialDate          = new AbsoluteDate();
213         final Frame        initialInertialFrame = FramesFactory.getGCRF();
214         final double       mu                   = 398600e9;
215 
216         final PVCoordinates initialPV = new PVCoordinates(
217                 new Vector3D(0, 6778000, 0),
218                 new Vector3D(-7668.63, 0, 0));
219 
220         final Orbit initialOrbit = new CartesianOrbit(initialPV, initialInertialFrame, initialDate, mu);
221 
222         final RealMatrix initialCovarianceInInertialFrame = new BlockRealMatrix(new double[][] {
223                 { 1, 0, 0, 0, 0, 1e-5 },
224                 { 0, 1, 0, 0, 0, 0 },
225                 { 0, 0, 1, 0, 0, 0 },
226                 { 0, 0, 0, 1e-3, 0, 0 },
227                 { 0, 0, 0, 0, 1e-3, 0 },
228                 { 1e-5, 0, 0, 0, 0, 1e-3 }
229         });
230 
231         final StateCovariance stateCovariance =
232                 new StateCovariance(initialCovarianceInInertialFrame, initialDate, initialInertialFrame,
233                         OrbitType.CARTESIAN, PositionAngleType.MEAN);
234 
235         // When
236         final RealMatrix convertedCovarianceMatrixInRTN =
237                 stateCovariance.changeCovarianceFrame(initialOrbit, LOFType.QSW_INERTIAL).getMatrix();
238 
239         // Then
240         // Expected covariance matrix obtained by rotation initial covariance matrix by 90 degrees
241         final RealMatrix expectedMatrixInRTN = new BlockRealMatrix(new double[][] {
242         	{ 1.000000e+00,  0.000000e+00, 0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00 },
243         	{ 0.000000e+00,  1.000000e+00, 0.000000e+00,  0.000000e+00,  0.000000e+00, -1.000000e-05 },
244         	{ 0.000000e+00,  0.000000e+00, 1.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00 },
245         	{ 0.000000e+00,  0.000000e+00, 0.000000e+00,  1.000000e-03,  0.000000e+00,  0.000000e+00 },
246         	{ 0.000000e+00,  0.000000e+00, 0.000000e+00,  0.000000e+00,  1.000000e-03,  0.000000e+00 },
247         	{ 0.000000e+00, -1.000000e-05, 0.000000e+00,  0.000000e+00,  0.000000e+00,  1.000000e-03 },
248         });
249 
250         compareCovariance(expectedMatrixInRTN, convertedCovarianceMatrixInRTN, DEFAULT_VALLADO_THRESHOLD);
251     }
252 
253     @Test
254     @DisplayName("Test covariance conversion from RTN local orbital frame to inertial frame")
255     void should_rotate_covariance_matrix_by_minus_ninety_degrees() {
256 
257         // Given
258         final AbsoluteDate initialDate          = new AbsoluteDate();
259         final Frame        initialInertialFrame = FramesFactory.getGCRF();
260         final double       mu                   = 398600e9;
261 
262         final PVCoordinates initialPV = new PVCoordinates(
263                 new Vector3D(0, 6778000, 0),
264                 new Vector3D(-7668.63, 0, 0));
265 
266         final Orbit initialOrbit = new CartesianOrbit(initialPV, initialInertialFrame, initialDate, mu);
267 
268         final RealMatrix initialCovarianceInRTN = new BlockRealMatrix(new double[][] {
269                 { 1, 0, 0, 0, 0, 0 },
270                 { 0, 1, 0, 0, 0, -1e-5 },
271                 { 0, 0, 1, 0, 0, 0 },
272                 { 0, 0, 0, 1e-3, 0, 0 },
273                 { 0, 0, 0, 0, 1e-3, 0 },
274                 { 0, -1e-5, 0, 0, 0, 1e-3 }
275         });
276 
277         final StateCovariance stateCovariance = new StateCovariance(initialCovarianceInRTN, initialDate, LOFType.QSW_INERTIAL);
278 
279         // When
280         final RealMatrix convertedCovarianceMatrixInInertialFrame = stateCovariance.changeCovarianceFrame(initialOrbit, initialInertialFrame).getMatrix();
281 
282         // Then
283 
284         // Expected covariance matrix obtained by rotation initial covariance matrix by -90 degrees
285         final RealMatrix expectedMatrixInInertialFrame = new BlockRealMatrix(new double[][] {
286         	{1.000000e+00,  0.000000e+00, 0.000000e+00,  0.000000e+00, 0.000000e+00, 1.000000e-05 },
287         	{0.000000e+00,  1.000000e+00, 0.000000e+00,  0.000000e+00, 0.000000e+00, 0.000000e+00 },
288         	{0.000000e+00,  0.000000e+00, 1.000000e+00,  0.000000e+00, 0.000000e+00, 0.000000e+00 },
289         	{0.000000e+00,  0.000000e+00, 0.000000e+00,  1.000000e-03, 0.000000e+00, 0.000000e+00 },
290         	{0.000000e+00,  0.000000e+00, 0.000000e+00,  0.000000e+00, 1.000000e-03, 0.000000e+00 },
291         	{1.000000e-05,  0.000000e+00, 0.000000e+00,  0.000000e+00, 0.000000e+00, 1.000000e-03 },
292         });
293 
294         compareCovariance(expectedMatrixInInertialFrame, convertedCovarianceMatrixInInertialFrame, DEFAULT_VALLADO_THRESHOLD);
295 
296     }
297 
298     /**
299      * This test is based on the following paper : Covariance Transformations for Satellite Flight Dynamics Operations
300      * from David A. Vallado.
301      * <p>
302      * More specifically, we're using the initial covariance matrix from p.14 and compare the computed result with the
303      * cartesian covariance in RSW from p.19.
304      * <p>
305      * In this case, the same transformation as the one in Vallado's paper is applied (only rotation) so the same values
306      * are expected.
307      * <p>
308      * Note that the followings local orbital frame are equivalent RSW=RTN=QSW.
309      */
310     @Test
311     @DisplayName("Test covariance conversion Vallado test case : ECI cartesian to RTN")
312     void should_return_Vallado_RSW_covariance_matrix_from_ECI() {
313 
314         // Initialize Orekit
315         Utils.setDataRoot("regular-data");
316 
317         // Given
318         final AbsoluteDate  initialDate          = getValladoInitialDate();
319         final PVCoordinates initialPV            = getValladoInitialPV();
320         final Frame         initialInertialFrame = FramesFactory.getGCRF();
321         final Orbit initialOrbit =
322                 new CartesianOrbit(initialPV, initialInertialFrame, initialDate, getValladoMu());
323 
324         final RealMatrix initialCovarianceMatrix = getValladoInitialCovarianceMatrix();
325 
326         final StateCovariance stateCovariance =
327                 new StateCovariance(initialCovarianceMatrix, initialDate, initialInertialFrame, OrbitType.CARTESIAN,
328                                     PositionAngleType.MEAN);
329         // When
330         final RealMatrix convertedCovarianceMatrixInRTN =
331                 stateCovariance.changeCovarianceFrame(initialOrbit, LOFType.QSW_INERTIAL).getMatrix();
332 
333         // Then
334         final RealMatrix expectedCovarianceMatrixInRTN = new BlockRealMatrix(new double[][] {
335                 { 9.918921e-001, 6.700644e-003, -2.878187e-003, 1.892086e-005, 6.700644e-005, -2.878187e-005 },
336                 { 6.700644e-003, 1.013730e+000, -1.019283e-002, 6.700644e-005, 2.372970e-004, -1.019283e-004 },
337                 { -2.878187e-003, -1.019283e-002, 9.943782e-001, -2.878187e-005, -1.019283e-004, 4.378217e-005 },
338                 { 1.892086e-005, 6.700644e-005, -2.878187e-005, 1.892086e-007, 6.700644e-007, -2.878187e-007 },
339                 { 6.700644e-005, 2.372970e-004, -1.019283e-004, 6.700644e-007, 2.372970e-006, -1.019283e-006 },
340                 { -2.878187e-005, -1.019283e-004, 4.378217e-005, -2.878187e-007, -1.019283e-006, 4.378217e-007 } });
341 
342         compareCovariance(expectedCovarianceMatrixInRTN, convertedCovarianceMatrixInRTN, DEFAULT_VALLADO_THRESHOLD);
343 
344     }
345 
346     /**
347      * This test is based on the following paper : Covariance Transformations for Satellite Flight Dynamics Operations
348      * from David A. Vallado.
349      * <p>
350      * More specifically, we're using the initial covariance matrix from p.14 and compare the computed result with the
351      * cartesian covariance in RSW from p.19.
352      * <p>
353      * In this case, Orekit applies the full frame transformation while Vallado's paper only take into account the
354      * rotation part. Therefore, some values are different with respect to the reference ones in the paper.
355      * <p>
356      * Note that the followings local orbital frame are equivalent RSW=RTN=QSW.
357      */
358     @Test
359     @DisplayName("Test covariance conversion Vallado test case : ECI cartesian to RTN")
360     void should_return_Vallado_RSW_non_inertial_covariance_matrix_from_ECI() {
361 
362         // Initialize Orekit
363         Utils.setDataRoot("regular-data");
364 
365         // Given
366         final AbsoluteDate  initialDate          = getValladoInitialDate();
367         final PVCoordinates initialPV            = getValladoInitialPV();
368         final Frame         initialInertialFrame = FramesFactory.getGCRF();
369         final Orbit initialOrbit =
370                 new CartesianOrbit(initialPV, initialInertialFrame, initialDate, getValladoMu());
371 
372         final RealMatrix initialCovarianceMatrix = getValladoInitialCovarianceMatrix();
373 
374         final StateCovariance stateCovariance =
375                 new StateCovariance(initialCovarianceMatrix, initialDate, initialInertialFrame, OrbitType.CARTESIAN,
376                         PositionAngleType.MEAN);
377         // When
378         final RealMatrix convertedCovarianceMatrixInRTN =
379                 stateCovariance.changeCovarianceFrame(initialOrbit, LOFType.QSW).getMatrix();
380 
381         // Then
382         final RealMatrix expectedCovarianceMatrixInRTN = new BlockRealMatrix(new double[][] {
383                 { 9.918921e-01, 6.700644e-03, -2.878187e-03, 2.637186e-05, -1.035961e-03, -2.878187e-05 },
384                 { 6.700644e-03, 1.013730e+00, -1.019283e-02, 1.194257e-03, 2.298460e-04, -1.019283e-04 },
385                 { -2.878187e-03, -1.019283e-02, 9.943782e-01, -4.011613e-05, -9.872780e-05, 4.378217e-05 },
386                 { 2.637186e-05, 1.194257e-03, -4.011613e-05, 1.591713e-06, 9.046096e-07, -4.011613e-07 },
387                 { -1.035961e-03, 2.298460e-04, -9.872780e-05, 9.046096e-07, 3.450431e-06, -9.872780e-07 },
388                 { -2.878187e-05, -1.019283e-04, 4.378217e-05, -4.011613e-07, -9.872780e-07, 4.378217e-07 }
389         });
390 
391         compareCovariance(expectedCovarianceMatrixInRTN, convertedCovarianceMatrixInRTN, DEFAULT_VALLADO_THRESHOLD);
392 
393     }
394 
395     private AbsoluteDate getValladoInitialDate() {
396         return new AbsoluteDate(2000, 12, 15, 16, 58, 50.208, TimeScalesFactory.getUTC());
397     }
398 
399     private PVCoordinates getValladoInitialPV() {
400         return new PVCoordinates(
401                 new Vector3D(-605792.21660, -5870229.51108, 3493053.19896),
402                 new Vector3D(-1568.25429, -3702.34891, -6479.48395));
403     }
404 
405     private double getValladoMu() {
406         return Constants.IERS2010_EARTH_MU;
407     }
408 
409     private RealMatrix getValladoInitialCovarianceMatrix() {
410         return new BlockRealMatrix(new double[][] {
411                 { 1, 1e-2, 1e-2, 1e-4, 1e-4, 1e-4 },
412                 { 1e-2, 1, 1e-2, 1e-4, 1e-4, 1e-4 },
413                 { 1e-2, 1e-2, 1, 1e-4, 1e-4, 1e-4 },
414                 { 1e-4, 1e-4, 1e-4, 1e-6, 1e-6, 1e-6 },
415                 { 1e-4, 1e-4, 1e-4, 1e-6, 1e-6, 1e-6 },
416                 { 1e-4, 1e-4, 1e-4, 1e-6, 1e-6, 1e-6 }
417         });
418     }
419 
420     /**
421      * This test is based on the following paper : Covariance Transformations for Satellite Flight Dynamics Operations
422      * from David A. Vallado.
423      * <p>
424      * More specifically, we're using the initial covariance matrix from p.14 and compare the computed result with the
425      * cartesian covariance in NTW from p.19.
426      * <p>
427      * In this case, the same transformation as the one in Vallado's paper is applied (only rotation) so the same values
428      * are expected.
429      * <p>
430      */
431     @Test
432     @DisplayName("Test covariance conversion Vallado test case : ECI cartesian to NTW ( considered inertial)")
433     void should_return_Vallado_NTW_covariance_matrix_from_ECI() {
434 
435         // Initialize orekit
436         Utils.setDataRoot("regular-data");
437 
438         // Given
439         final AbsoluteDate  initialDate          = getValladoInitialDate();
440         final PVCoordinates initialPV            = getValladoInitialPV();
441         final Frame         initialInertialFrame = FramesFactory.getGCRF();
442         final Orbit initialOrbit =
443                 new CartesianOrbit(initialPV, initialInertialFrame, initialDate, getValladoMu());
444 
445         final RealMatrix initialCovarianceMatrix = getValladoInitialCovarianceMatrix();
446 
447         final StateCovariance stateCovariance =
448                 new StateCovariance(initialCovarianceMatrix, initialDate, initialInertialFrame,
449                         OrbitType.CARTESIAN, PositionAngleType.MEAN);
450         // When
451         final RealMatrix convertedCovarianceMatrixInNTW =
452         		stateCovariance.changeCovarianceFrame(initialOrbit, LOFType.NTW_INERTIAL).getMatrix();
453 
454         // Then
455         final RealMatrix expectedCovarianceMatrixInNTW = new BlockRealMatrix(new double[][] {
456                 { 9.918792e-001, 6.679546e-003, -2.868345e-003, 1.879167e-005, 6.679546e-005, -2.868345e-005 },
457                 { 6.679546e-003, 1.013743e+000, -1.019560e-002, 6.679546e-005, 2.374262e-004, -1.019560e-004 },
458                 { -2.868345e-003, -1.019560e-002, 9.943782e-001, -2.868345e-005, -1.019560e-004, 4.378217e-005 },
459                 { 1.879167e-005, 6.679546e-005, -2.868345e-005, 1.879167e-007, 6.679546e-007, -2.868345e-007 },
460                 { 6.679546e-005, 2.374262e-004, -1.019560e-004, 6.679546e-007, 2.374262e-006, -1.019560e-006 },
461                 { -2.868345e-005, -1.019560e-004, 4.378217e-005, -2.868345e-007, -1.019560e-006, 4.378217e-007 } });
462 
463         compareCovariance(expectedCovarianceMatrixInNTW, convertedCovarianceMatrixInNTW, DEFAULT_VALLADO_THRESHOLD);
464 
465     }
466 
467     /**
468      * This test is based on the following paper : Covariance Transformations for Satellite Flight Dynamics Operations
469      * from David A. Vallado.
470      * <p>
471      * More specifically, we're using the initial covariance matrix from p.14 and compare the computed result with the
472      * cartesian covariance in NTW from p.19.
473      * <p>
474      * In this case, Orekit applies the full frame transformation while Vallado's paper only take into account the
475      * rotation part. Therefore, some values are different with respect to the reference ones in the paper.
476      * <p>
477      */
478     @Test
479     @DisplayName("Test covariance conversion Vallado test case : ECI cartesian to NTW (non inertial)")
480     void should_return_Vallado_NTW_non_inertial_covariance_matrix_from_ECI() {
481 
482         // Initialize orekit
483         Utils.setDataRoot("regular-data");
484 
485         // Given
486         final AbsoluteDate  initialDate          = getValladoInitialDate();
487         final PVCoordinates initialPV            = getValladoInitialPV();
488         final Frame         initialInertialFrame = FramesFactory.getGCRF();
489         final Orbit initialOrbit =
490                 new CartesianOrbit(initialPV, initialInertialFrame, initialDate, getValladoMu());
491 
492         final RealMatrix initialCovarianceMatrix = getValladoInitialCovarianceMatrix();
493 
494         final StateCovariance stateCovariance =
495                 new StateCovariance(initialCovarianceMatrix, initialDate, initialInertialFrame,
496                         OrbitType.CARTESIAN, PositionAngleType.MEAN);
497 
498         // When
499         final RealMatrix convertedCovarianceMatrixInNTW =
500                 stateCovariance.changeCovarianceFrame(initialOrbit, LOFType.NTW).getMatrix();
501 
502         // Then
503         final RealMatrix expectedCovarianceMatrixInNTW = new BlockRealMatrix(new double[][] {
504                 { 9.918792e-01,  6.679546e-03, -2.868345e-03, 2.6215894e-05, -1.035665e-03, -2.868345e-05 },
505                 { 6.679546e-03,  1.013743e+00, -1.019560e-02,  1.193556e-03,  2.300019e-04, -1.019560e-04 },
506                 {-2.868345e-03, -1.019560e-02,  9.943782e-01, -4.0015724e-05, -9.8767909e-05,  4.378217e-05 },
507                 { 2.6215894e-05,  1.193556e-03, -4.0015724e-05,  1.58878e-06,  9.0271200e-07, -4.0015724e-07 },
508                 {-1.035665e-03,  2.300019e-04, -9.8767909e-05,  9.0271200e-07,  3.4511471e-06, -9.8767909e-07 },
509                 {-2.868345e-05, -1.019560e-04,  4.378217e-05, -4.0015724e-07, -9.8767909e-07,  4.378217e-07 },
510         });
511 
512         compareCovariance(expectedCovarianceMatrixInNTW, convertedCovarianceMatrixInNTW, DEFAULT_VALLADO_THRESHOLD);
513 
514     }
515 
516     /**
517      * This test is based on the following paper : Covariance Transformations for Satellite Flight Dynamics Operations
518      * from David A. Vallado.
519      * <p>
520      * More specifically, we're using the initial covariance matrix from p.14 as a reference to test multiple
521      * conversions.
522      * <p>
523      * This test aims to verify the numerical precision after various conversions and serves as a non regression test
524      * for future updates.
525      * <p>
526      * Also, note that the conversion from the RTN to TEME tests the fact that the orbit is initially expressed in GCRF
527      * while we want the covariance expressed in TEME. Hence, it tests that the rotation from RTN to TEME needs to be
528      * obtained by expressing the orbit PVCoordinates in the TEME frame (hence the use of orbit.gtPVCoordinates(frameOut)
529      * ,see relevant changeCovarianceFrame method).
530      */
531     @Test
532     @DisplayName("Test custom covariance conversion Vallado test case : GCRF -> TEME -> IRTF -> NTW -> RTN -> ITRF -> GCRF")
533     void should_return_initial_covariance_after_multiple_conversion() {
534 
535         // Initialize orekit
536         Utils.setDataRoot("regular-data");
537 
538         // Given
539         final AbsoluteDate  initialDate          = getValladoInitialDate();
540         final PVCoordinates initialPV            = getValladoInitialPV();
541         final Frame         initialInertialFrame = FramesFactory.getGCRF();
542         final Orbit initialOrbit =
543                 new CartesianOrbit(initialPV, initialInertialFrame, initialDate, getValladoMu());
544 
545         final RealMatrix initialCovarianceMatrixInGCRF = getValladoInitialCovarianceMatrix();
546 
547         final Frame teme = FramesFactory.getTEME();
548 
549         final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, false);
550 
551         final StateCovariance stateCovariance =
552                 new StateCovariance(initialCovarianceMatrixInGCRF, initialDate, initialInertialFrame,
553                         OrbitType.CARTESIAN, PositionAngleType.MEAN);
554 
555         // When
556         // GCRF -> TEME
557         final StateCovariance convertedCovarianceMatrixInTEME =
558         		stateCovariance.changeCovarianceFrame(initialOrbit, teme);
559 
560         // TEME -> ITRF
561         final StateCovariance convertedCovarianceMatrixInITRF =
562         		convertedCovarianceMatrixInTEME.changeCovarianceFrame(initialOrbit, itrf);
563 
564         // ITRF -> NTW
565         final StateCovariance convertedCovarianceMatrixInNTW =
566         		convertedCovarianceMatrixInITRF.changeCovarianceFrame(initialOrbit, LOFType.NTW);
567 
568         // NTW -> RTN
569         final StateCovariance convertedCovarianceMatrixInRTN =
570         		convertedCovarianceMatrixInNTW.changeCovarianceFrame(initialOrbit, LOFType.QSW);
571 
572         // RTN -> ITRF
573         final StateCovariance convertedCovarianceMatrixBackInITRF =
574         		convertedCovarianceMatrixInRTN.changeCovarianceFrame(initialOrbit, itrf);
575 
576         // ITRF -> TEME
577         final StateCovariance convertedCovarianceMatrixBackInTEME =
578         		convertedCovarianceMatrixBackInITRF.changeCovarianceFrame(initialOrbit, teme);
579 
580         // TEME -> GCRF
581         final StateCovariance convertedCovarianceMatrixInGCRF =
582         		convertedCovarianceMatrixBackInTEME.changeCovarianceFrame(initialOrbit, initialInertialFrame);
583 
584         // Then
585 
586         compareCovariance(initialCovarianceMatrixInGCRF, convertedCovarianceMatrixInGCRF.getMatrix(), 1e-12);
587 
588     }
589 
590     /**
591      * This test is based on the following paper : Covariance Transformations for Satellite Flight Dynamics Operations
592      * from David A. Vallado.
593      * <p>
594      * More specifically, we're using the initial covariance matrix from p.14 and compare the computed result with the
595      * cartesian covariance in ECEF from p.18.
596      * </p>
597      * <p>
598      * <b>BEWARE: It has been found that the earth rotation in this Vallado's case is given 1 million times slower than
599      * the expected value. This has been corrected and the expected covariance matrix is now the covariance matrix
600      * computed by Orekit given the similarities with Vallado's results. In addition, the small differences potentially
601      * come from the custom EOP that Vallado uses. Hence, this test can be considered as a <u>non regression
602      * test</u>.</b>
603      * </p>
604      */
605     @Test
606     @DisplayName("Test covariance conversion Vallado test case : ECI cartesian to PEF")
607     void should_return_Vallado_PEF_covariance_matrix_from_ECI() {
608 
609         // Initialize orekit
610         Utils.setDataRoot("regular-data");
611 
612         // Given
613         // Initialize orbit
614         final AbsoluteDate  initialDate          = getValladoInitialDate();
615         final PVCoordinates initialPV            = getValladoInitialPV();
616         final Frame         initialInertialFrame = FramesFactory.getGCRF();
617 
618         final Orbit initialOrbit = new CartesianOrbit(initialPV, initialInertialFrame, initialDate, getValladoMu());
619 
620         // Initialize input covariance matrix
621         final RealMatrix initialCovarianceMatrix = getValladoInitialCovarianceMatrix();
622 
623         // Initialize output frame
624         final Frame outputFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
625 
626         final StateCovariance stateCovariance = new StateCovariance(initialCovarianceMatrix, initialDate, initialInertialFrame, OrbitType.CARTESIAN, PositionAngleType.MEAN);
627         // When
628         final RealMatrix convertedCovarianceMatrixInITRF =
629         		stateCovariance.changeCovarianceFrame(initialOrbit, outputFrame).getMatrix();
630 
631         // Then
632         final RealMatrix expectedCovarianceMatrixInITRF = new BlockRealMatrix(new double[][] {
633                 { 9.9340005761276870e-01, 7.5124999798868530e-03, 5.8312675007359050e-03, 3.4548396261054936e-05, 2.6851237046859200e-06, 5.8312677693153940e-05 },
634                 { 7.5124999798868025e-03, 1.0065990293034541e+00, 1.2884310200351924e-02, 1.4852736004690684e-04, 1.6544247282904867e-04, 1.2884310644320954e-04 },
635                 { 5.8312675007359040e-03, 1.2884310200351924e-02, 1.0000009130837746e+00, 5.9252211072590390e-05, 1.2841787487219444e-04, 1.0000913090989617e-04 },
636                 { 3.4548396261054936e-05, 1.4852736004690686e-04, 5.9252211072590403e-05, 3.5631474857130520e-07, 7.6083489184819870e-07, 5.9252213790760030e-07 },
637                 { 2.6851237046859150e-06, 1.6544247282904864e-04, 1.2841787487219447e-04, 7.6083489184819880e-07, 1.6542289254142709e-06, 1.2841787929229964e-06 },
638                 { 5.8312677693153934e-05, 1.2884310644320950e-04, 1.0000913090989616e-04, 5.9252213790760020e-07, 1.2841787929229960e-06, 1.0000913098203875e-06 }
639         });
640 
641         compareCovariance(expectedCovarianceMatrixInITRF, convertedCovarianceMatrixInITRF, 1e-20);
642 
643     }
644 
645     /**
646      * This test is based on the following paper : Covariance Transformations for Satellite Flight Dynamics Operations
647      * from David A. Vallado.
648      * <p>
649      * More specifically, we're using the initial covariance matrix from p.14 and compare the computed result with the
650      * cartesian covariance in PEF from p.18.
651      */
652     @Test
653     @DisplayName("Test covariance conversion Vallado test case : PEF cartesian to ECI")
654     void should_return_Vallado_ECI_covariance_matrix_from_PEF() {
655 
656         // Initialize orekit
657         Utils.setDataRoot("regular-data");
658 
659         // Given
660         final AbsoluteDate  initialDate          = getValladoInitialDate();
661         final PVCoordinates initialPV            = getValladoInitialPV();
662         final Frame         initialInertialFrame = FramesFactory.getGCRF();
663         final Orbit initialOrbit =
664                 new CartesianOrbit(initialPV, initialInertialFrame, initialDate, getValladoMu());
665 
666         final RealMatrix initialCovarianceMatrixInPEF = new BlockRealMatrix(new double[][] {
667                 { 9.9340005761276870e-01, 7.5124999798868530e-03, 5.8312675007359050e-03, 3.4548396261054936e-05, 2.6851237046859200e-06, 5.8312677693153940e-05 },
668                 { 7.5124999798868025e-03, 1.0065990293034541e+00, 1.2884310200351924e-02, 1.4852736004690684e-04, 1.6544247282904867e-04, 1.2884310644320954e-04 },
669                 { 5.8312675007359040e-03, 1.2884310200351924e-02, 1.0000009130837746e+00, 5.9252211072590390e-05, 1.2841787487219444e-04, 1.0000913090989617e-04 },
670                 { 3.4548396261054936e-05, 1.4852736004690686e-04, 5.9252211072590403e-05, 3.5631474857130520e-07, 7.6083489184819870e-07, 5.9252213790760030e-07 },
671                 { 2.6851237046859150e-06, 1.6544247282904864e-04, 1.2841787487219447e-04, 7.6083489184819880e-07, 1.6542289254142709e-06, 1.2841787929229964e-06 },
672                 { 5.8312677693153934e-05, 1.2884310644320950e-04, 1.0000913090989616e-04, 5.9252213790760020e-07, 1.2841787929229960e-06, 1.0000913098203875e-06 }
673         });
674 
675         final Frame inputFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, false);
676 
677         // State covariance
678         final StateCovariance stateCovariance = new StateCovariance(initialCovarianceMatrixInPEF, initialDate, inputFrame, OrbitType.CARTESIAN, PositionAngleType.MEAN);
679 
680         // When
681         final RealMatrix convertedCovarianceMatrixInECI = stateCovariance.changeCovarianceFrame(initialOrbit, initialInertialFrame).getMatrix();
682 
683         // Then
684         final RealMatrix expectedCovarianceMatrixInECI = getValladoInitialCovarianceMatrix();
685 
686         compareCovariance(expectedCovarianceMatrixInECI, convertedCovarianceMatrixInECI, 1e-7);
687 
688     }
689 
690     /**
691      * This test is based on the following paper : Covariance Transformations for Satellite Flight Dynamics Operations
692      * from David A. Vallado.
693      * <p>
694      * More specifically, we're using the initial covariance matrix from p.14 and compare the computed result with the
695      * cartesian covariance in MOD from p.17.
696      * </p>
697      */
698     @Test
699     @DisplayName("Test covariance conversion Vallado test case : ECI cartesian to MOD")
700     void should_return_Vallado_MOD_covariance_matrix_from_ECI() {
701 
702         // Initialize orekit
703         Utils.setDataRoot("regular-data");
704 
705         // Given
706         final AbsoluteDate  initialDate          = getValladoInitialDate();
707         final PVCoordinates initialPV            = getValladoInitialPV();
708         final Frame         initialInertialFrame = FramesFactory.getGCRF();
709         final Orbit initialOrbit =
710                 new CartesianOrbit(initialPV, initialInertialFrame, initialDate, getValladoMu());
711 
712         final RealMatrix initialCovarianceMatrix = getValladoInitialCovarianceMatrix();
713 
714         final Frame outputFrame = FramesFactory.getMOD(IERSConventions.IERS_2010);
715 
716         final StateCovariance stateCovariance =
717                 new StateCovariance(initialCovarianceMatrix, initialDate, initialInertialFrame,
718                         OrbitType.CARTESIAN, PositionAngleType.MEAN);
719         // When
720         final RealMatrix convertedCovarianceMatrixInMOD =
721                 stateCovariance.changeCovarianceFrame(initialOrbit, outputFrame).getMatrix();
722 
723         // Then
724         final RealMatrix expectedCovarianceMatrixInMOD = new BlockRealMatrix(new double[][] {
725                 { 9.999939e-001, 9.999070e-003, 9.997861e-003, 9.993866e-005, 9.999070e-005, 9.997861e-005 },
726                 { 9.999070e-003, 1.000004e+000, 1.000307e-002, 9.999070e-005, 1.000428e-004, 1.000307e-004 },
727                 { 9.997861e-003, 1.000307e-002, 1.000002e+000, 9.997861e-005, 1.000307e-004, 1.000186e-004 },
728                 { 9.993866e-005, 9.999070e-005, 9.997861e-005, 9.993866e-007, 9.999070e-007, 9.997861e-007 },
729                 { 9.999070e-005, 1.000428e-004, 1.000307e-004, 9.999070e-007, 1.000428e-006, 1.000307e-006 },
730                 { 9.997861e-005, 1.000307e-004, 1.000186e-004, 9.997861e-007, 1.000307e-006, 1.000186e-006 },
731         });
732 
733         compareCovariance(expectedCovarianceMatrixInMOD, convertedCovarianceMatrixInMOD, DEFAULT_VALLADO_THRESHOLD);
734 
735     }
736 
737     /**
738      * This test is based on the following paper : Covariance Transformations for Satellite Flight Dynamics Operations
739      * from David A. Vallado.
740      * <p>
741      * More specifically, we're using the initial NTW covariance matrix from p.19 and compare the computed result with
742      * the cartesian covariance in RSW from the same page.
743      * </p>
744      */
745     @Test
746     @DisplayName("Test covariance conversion from Vallado test case NTW to RSW")
747     void should_convert_Vallado_NTW_to_RSW() {
748 
749         // Initialize orekit
750         Utils.setDataRoot("regular-data");
751 
752         // Given
753         final AbsoluteDate  initialDate          = getValladoInitialDate();
754         final PVCoordinates initialPV            = getValladoInitialPV();
755         final Frame         initialInertialFrame = FramesFactory.getGCRF();
756         final Orbit initialOrbit =
757                 new CartesianOrbit(initialPV, initialInertialFrame, initialDate, getValladoMu());
758 
759         final RealMatrix initialCovarianceMatrixInNTW = new BlockRealMatrix(new double[][] {
760                 { 9.918792e-001, 6.679546e-003, -2.868345e-003, 1.879167e-005, 6.679546e-005, -2.868345e-005 },
761                 { 6.679546e-003, 1.013743e+000, -1.019560e-002, 6.679546e-005, 2.374262e-004, -1.019560e-004 },
762                 { -2.868345e-003, -1.019560e-002, 9.943782e-001, -2.868345e-005, -1.019560e-004, 4.378217e-005 },
763                 { 1.879167e-005, 6.679546e-005, -2.868345e-005, 1.879167e-007, 6.679546e-007, -2.868345e-007 },
764                 { 6.679546e-005, 2.374262e-004, -1.019560e-004, 6.679546e-007, 2.374262e-006, -1.019560e-006 },
765                 { -2.868345e-005, -1.019560e-004, 4.378217e-005, -2.868345e-007, -1.019560e-006, 4.378217e-007 }
766         });
767 
768         final StateCovariance stateCovariance =
769                 new StateCovariance(initialCovarianceMatrixInNTW, initialDate, LOFType.NTW);
770 
771         // When
772         final RealMatrix convertedCovarianceMatrixInRTN =
773         		stateCovariance.changeCovarianceFrame(initialOrbit, LOFType.QSW).getMatrix();
774 
775         // Then
776         final RealMatrix expectedCovarianceMatrixInRTN = new BlockRealMatrix(new double[][] {
777                 { 9.918921e-001, 6.700644e-003, -2.878187e-003, 1.8924189e-005, 6.651362e-005, -2.878187e-005 },
778                 { 6.700644e-003, 1.013730e+000, -1.019283e-002, 6.7510094e-005, 2.3729368e-004, -1.01928257e-004 },
779                 { -2.878187e-003, -1.019283e-002, 9.943782e-001, -2.8786942e-005, -1.01926827e-004, 4.378217e-005 },
780                 { 1.8924189e-005, 6.7510094e-005, -2.8786942e-005, 1.89275434e-007, 6.7017283e-007, -2.8786942e-007 },
781                 { 6.651362e-005, 2.3729368e-004, -1.01926827e-004, 6.7017283e-007, 2.372903e-006, -1.0192682e-006 },
782                 { -2.878187e-005, -1.01928257e-004, 4.378217e-005, -2.87869424e-007, -1.0192682e-006, 4.378217e-007 }
783         });
784 
785         compareCovariance(expectedCovarianceMatrixInRTN, convertedCovarianceMatrixInRTN, DEFAULT_VALLADO_THRESHOLD);
786 
787     }
788 
789     @Test
790     @Deprecated
791     void testGetStm() {
792         // Given
793         final AbsoluteDate  initialDate      = new AbsoluteDate();
794         final Frame         inertialFrame    = FramesFactory.getGCRF();
795         final PVCoordinates pv               = getValladoInitialPV();
796         final double        mu               = Constants.IERS2010_EARTH_MU;
797         final Orbit initialOrbit = new CartesianOrbit(pv, inertialFrame, initialDate, mu);
798         final double dt = 100;
799 
800         // When
801         final RealMatrix stm = StateCovariance.getStm(initialOrbit, dt);
802 
803         // Then
804         final RealMatrix expectedStm          = MatrixUtils.createRealIdentityMatrix(6);
805         final double     sma          = initialOrbit.getA();
806         final double     contribution = -1.5 * dt * FastMath.sqrt(mu / FastMath.pow(sma, 5));
807         stm.setEntry(5, 0, contribution);
808         Assertions.assertEquals(0., expectedStm.subtract(stm).getNorm1(), 1e-7);
809     }
810 
811     @ParameterizedTest
812     @EnumSource(OrbitType.class)
813     void testGetKeplerianStm(final OrbitType orbitType) {
814         // GIVEN
815         final AbsoluteDate  initialDate      = new AbsoluteDate();
816         final Frame         inertialFrame    = FramesFactory.getGCRF();
817         final PVCoordinates pv               = getValladoInitialPV();
818         final Orbit initialOrbit = new CartesianOrbit(pv, inertialFrame, initialDate, Constants.EGM96_EARTH_MU);
819         final StateCovariance cartesianCovariance = new StateCovariance(getValladoInitialCovarianceMatrix(), initialDate,
820                 inertialFrame, OrbitType.CARTESIAN, PositionAngleType.MEAN);
821         final StateCovariance covariance = cartesianCovariance.changeCovarianceType(initialOrbit, orbitType,
822                 cartesianCovariance.getPositionAngleType());
823         final double dt = 1e3;
824 
825         // WHEN
826         final RealMatrix actualStm = covariance.getKeplerianStm(initialOrbit, dt);
827 
828         // THEN
829         final CartesianOrbit orbitWithNonKeplerianAcceleration = new CartesianOrbit(new PVCoordinates(pv.getPosition(),
830                 pv.getVelocity(), new Vector3D(1, 2, 3)), inertialFrame, initialDate, initialOrbit.getMu());
831         final RealMatrix expectedStm = covariance.getKeplerianStm(orbitWithNonKeplerianAcceleration, dt);
832         Assertions.assertEquals(0., expectedStm.subtract(actualStm).getNorm1());
833 
834     }
835 
836     /**
837      * The goal of this test is to check the shiftedBy method of {@link StateCovariance} by creating one state
838      * covariance expressed in 3 different ways (inertial Equinoctial, LOF cartesian and non inertial cartesian) -> shift
839      * them -> reconvert them back to the same initial frame & type -> compare them with expected, manually computed
840      * covariance matrix.
841      */
842     @Test
843     @DisplayName("Test shiftedBy method of StateCovariance")
844     void should_return_expected_shifted_state_covariance() {
845 
846         // Initialize orekit
847         Utils.setDataRoot("regular-data");
848 
849         // Given
850         final AbsoluteDate  initialDate      = new AbsoluteDate();
851         final Frame         inertialFrame    = FramesFactory.getGCRF();
852         final Frame         nonInertialFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
853         final PVCoordinates pv               = getValladoInitialPV();
854         final double        mu               = Constants.IERS2010_EARTH_MU;
855         final Orbit initialOrbit = new CartesianOrbit(pv, inertialFrame, initialDate, mu);
856         final EquinoctialOrbit equinoctialOrbit = ((EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(initialOrbit)).withCachedPositionAngleType(PositionAngleType.MEAN);
857         final KeplerianOrbit keplerianOrbit = ((KeplerianOrbit) OrbitType.KEPLERIAN.convertType(initialOrbit)).withCachedPositionAngleType(PositionAngleType.MEAN);
858 
859         final double timeShift = 300; // In s
860 
861         // Initializing initial covariance matrix common to all
862         final StateCovariance initialCovarianceInCartesian =
863                 new StateCovariance(getValladoInitialCovarianceMatrix(), initialDate, inertialFrame,
864                                     OrbitType.CARTESIAN, PositionAngleType.MEAN);
865 
866         final StateCovariance covarianceInEquinoctial =
867                 initialCovarianceInCartesian.changeCovarianceType(equinoctialOrbit, OrbitType.EQUINOCTIAL, PositionAngleType.MEAN);
868 
869         final StateCovariance covarianceInCartesianInLOF =
870                 initialCovarianceInCartesian.changeCovarianceFrame(initialOrbit, LOFType.QSW);
871 
872         final StateCovariance covarianceInCartesianInNonInertial =
873                 initialCovarianceInCartesian.changeCovarianceFrame(initialOrbit, nonInertialFrame);
874 
875         // When
876         final StateCovariance shiftedCovarianceInEquinoctial =
877                 covarianceInEquinoctial.shiftedBy(equinoctialOrbit, timeShift);
878         final RealMatrix shiftedCovarianceInEquinoctialBackToInitial =
879                 shiftedCovarianceInEquinoctial.changeCovarianceType(equinoctialOrbit.shiftedBy(timeShift),
880                                                                   OrbitType.CARTESIAN, PositionAngleType.MEAN)
881                         .getMatrix();
882 
883         final StateCovariance shiftedCovarianceInCartesianInLOF =
884                 covarianceInCartesianInLOF.shiftedBy(initialOrbit, timeShift);
885         final RealMatrix shiftedCovarianceInCartesianInLOFBackToInitial =
886                 shiftedCovarianceInCartesianInLOF.changeCovarianceFrame(initialOrbit.shiftedBy(timeShift),
887                                                                         inertialFrame)
888                         .getMatrix();
889 
890         final StateCovariance shiftedCovarianceInCartesianInNonInertial =
891                 covarianceInCartesianInNonInertial.shiftedBy(initialOrbit, timeShift);
892         final RealMatrix shiftedCovarianceInCartesianInNonInertialBackToInitial =
893                 shiftedCovarianceInCartesianInNonInertial.changeCovarianceFrame(initialOrbit.shiftedBy(timeShift),
894                                                                                 inertialFrame)
895                         .getMatrix();
896 
897         // Then
898 
899         final OrbitType keplerianType = OrbitType.KEPLERIAN;
900         final StateCovariance initialCovarianceInKeplerian =
901                 initialCovarianceInCartesian.changeCovarianceType(keplerianOrbit, keplerianType, PositionAngleType.MEAN);
902 
903         final RealMatrix stm          = initialCovarianceInKeplerian.getKeplerianStm(keplerianOrbit, timeShift);
904         final RealMatrix referenceCovarianceMatrixInKeplerian =
905                 stm.multiply(initialCovarianceInKeplerian.getMatrix().multiplyTransposed(stm));
906 
907         final RealMatrix referenceCovarianceMatrixInCartesian =
908                 new StateCovariance(referenceCovarianceMatrixInKeplerian, initialDate.shiftedBy(timeShift),
909                                     inertialFrame, keplerianType, PositionAngleType.MEAN).changeCovarianceType(
910                         keplerianOrbit.shiftedBy(timeShift), OrbitType.CARTESIAN, PositionAngleType.MEAN).getMatrix();
911 
912         // Compare with results
913         compareCovariance(referenceCovarianceMatrixInCartesian, shiftedCovarianceInEquinoctialBackToInitial, 1e-6);
914         compareCovariance(referenceCovarianceMatrixInCartesian, shiftedCovarianceInCartesianInLOFBackToInitial, 1e-7);
915         compareCovariance(referenceCovarianceMatrixInCartesian, shiftedCovarianceInCartesianInNonInertialBackToInitial, 1e-7);
916 
917     }
918 
919     /**
920      * This test is based on the following paper : Covariance Transformations for Satellite Flight Dynamics Operations
921      * from David A. Vallado.
922      * <p>
923      * More specifically, we're using the initial NTW covariance matrix from p.19 and compare the computed result with
924      * the cartesian covariance in RSW from the same page.
925      * </p>
926      */
927     @Test
928     @DisplayName("Test thrown error if input frame is not pseudo-inertial and "
929             + "the covariance matrix is not expressed in cartesian elements")
930     void should_return_orekit_exception() {
931 
932         // Initialize orekit
933         Utils.setDataRoot("regular-data");
934 
935         // Given
936         final AbsoluteDate  initialDate          = getValladoInitialDate();
937         final PVCoordinates initialPV            = getValladoInitialPV();
938         final Frame         initialInertialFrame = FramesFactory.getGCRF();
939         final Orbit initialOrbit =
940                 new CartesianOrbit(initialPV, initialInertialFrame, initialDate, getValladoMu());
941 
942         final RealMatrix randomCovarianceMatrix = new BlockRealMatrix(new double[][] {
943                 { 9.918792e-001, 6.679546e-003, -2.868345e-003, 1.879167e-005, 6.679546e-005, -2.868345e-005 },
944                 { 6.679546e-003, 1.013743e+000, -1.019560e-002, 6.679546e-005, 2.374262e-004, -1.019560e-004 },
945                 { -2.868345e-003, -1.019560e-002, 9.943782e-001, -2.868345e-005, -1.019560e-004, 4.378217e-005 },
946                 { 1.879167e-005, 6.679546e-005, -2.868345e-005, 1.879167e-007, 6.679546e-007, -2.868345e-007 },
947                 { 6.679546e-005, 2.374262e-004, -1.019560e-004, 6.679546e-007, 2.374262e-006, -1.019560e-006 },
948                 { -2.868345e-005, -1.019560e-004, 4.378217e-005, -2.868345e-007, -1.019560e-006, 4.378217e-007 } });
949 
950         final Frame nonInertialFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
951 
952         final Frame inertialFrame = FramesFactory.getGCRF();
953 
954         // When & Then
955         Assertions.assertThrows(OrekitException.class,
956                 () -> new StateCovariance(randomCovarianceMatrix, initialDate, nonInertialFrame, OrbitType.CIRCULAR,
957                         PositionAngleType.MEAN).changeCovarianceFrame(initialOrbit, inertialFrame));
958 
959         Assertions.assertThrows(OrekitException.class,
960                 () -> new StateCovariance(randomCovarianceMatrix, initialDate, nonInertialFrame, OrbitType.EQUINOCTIAL,
961                         PositionAngleType.MEAN).changeCovarianceFrame(initialOrbit, LOFType.QSW));
962 
963         Assertions.assertThrows(OrekitException.class,
964                 () -> new StateCovariance(randomCovarianceMatrix, initialDate, nonInertialFrame, OrbitType.EQUINOCTIAL,
965                         PositionAngleType.MEAN).changeCovarianceType(initialOrbit, OrbitType.KEPLERIAN,
966                         PositionAngleType.MEAN));
967 
968         Assertions.assertThrows(OrekitException.class,
969                 () -> new StateCovariance(randomCovarianceMatrix, initialDate, LOFType.QSW).changeCovarianceType(
970                         initialOrbit, OrbitType.KEPLERIAN, PositionAngleType.MEAN));
971 
972         Assertions.assertThrows(OrekitException.class,
973                 () -> new StateCovariance(randomCovarianceMatrix, initialDate, nonInertialFrame, OrbitType.CARTESIAN,
974                         PositionAngleType.MEAN).changeCovarianceType(initialOrbit, OrbitType.KEPLERIAN,
975                         PositionAngleType.MEAN));
976 
977     }
978 
979     @Test
980     @DisplayName("test issue 1052 : Error when propagating covariance defined in non-inertial frame")
981     void testIssue1052() {
982 
983         // Initialize orekit
984         Utils.setDataRoot("regular-data");
985 
986         // Initialize state and covariance
987         setUp();
988 
989         // Given
990         final AbsoluteDate  initialDate      = initialState.getDate();
991         final OrbitType     stmOrbitType     = OrbitType.CARTESIAN;
992         final PositionAngleType stmPositionAngleType = PositionAngleType.MEAN;
993 
994         final NumericalPropagator propagator = buildDefaultPropagator(initialState, stmOrbitType, stmPositionAngleType);
995 
996         // Initialize harvester
997         final String            stmAdditionalName = "stm";
998         final MatricesHarvester harvester         = propagator.setupMatricesComputation(stmAdditionalName, null, null);
999 
1000         // Initialize covariance
1001         final StateCovariance initialCovariance = new StateCovariance(new BlockRealMatrix(initCov),
1002                                                                       initialDate, LOFType.QSW);
1003 
1004         // Initialize covariance provider
1005         final StateCovarianceMatrixProvider provider =
1006                 new StateCovarianceMatrixProvider("covariance", stmAdditionalName, harvester, initialCovariance);
1007 
1008         propagator.addAdditionalDataProvider(provider);
1009 
1010         // When
1011         final SpacecraftState propagatedState      = propagator.propagate(initialDate.shiftedBy(1));
1012         final StateCovariance propagatedCovariance = provider.getStateCovariance(propagatedState);
1013 
1014         // Assert that the error message is not thrown anymore (cannot change covariance type if defined in LOF)
1015         Assertions.assertDoesNotThrow(() -> propagator.propagate(initialDate.shiftedBy(1)));
1016 
1017         // Assert that propagated covariance is in the same LOF as the initial covariance
1018         assertEquals(LOFType.QSW, propagatedCovariance.getLOF());
1019 
1020     }
1021 
1022     @Test
1023     void testIssue1485CartesianOrbitTypeAndNullAngleType() {
1024         // GIVEN
1025         final AbsoluteDate      date           = AbsoluteDate.ARBITRARY_EPOCH;
1026         final RealMatrix        expectedMatrix = MatrixUtils.createRealIdentityMatrix(6);
1027         final PositionAngleType anomalyType    = PositionAngleType.MEAN;
1028         final Frame             inFrame        = FramesFactory.getEME2000();
1029         final Frame             outFrame       = FramesFactory.getGCRF();
1030         final double            mu             = Constants.EGM96_EARTH_MU;
1031 
1032         final KeplerianOrbit orbit = new KeplerianOrbit(1e7, 0.01, 1., 2., 3., 4., anomalyType, inFrame, date, mu);
1033 
1034         final StateCovariance originalCovariance =
1035                 new StateCovariance(expectedMatrix, date, inFrame, OrbitType.CARTESIAN, null);
1036 
1037         // WHEN & THEN
1038         assertDoesNotThrow(()-> originalCovariance.changeCovarianceFrame(orbit, outFrame));
1039 
1040         // THEN
1041         final StateCovariance convertedCovariance = originalCovariance.changeCovarianceFrame(orbit, outFrame);
1042         assertEquals(originalCovariance.getMatrix().getNorm1(), convertedCovariance.getMatrix().getNorm1(), 1e-15);
1043     }
1044 
1045     @Test
1046     void testIssue1485SameLOF() {
1047         // GIVEN
1048         final LOF lof = LOFType.TNW;
1049 
1050         // WHEN & THEN
1051         doTestIssue1485SameFrameDefinition(lof);
1052     }
1053 
1054     @Test
1055     void testIssue1485SameFrame() {
1056         // GIVEN
1057         final Frame eme2000 = FramesFactory.getEME2000();
1058 
1059         // WHEN & THEN
1060         doTestIssue1485SameFrameDefinition(eme2000);
1061     }
1062 
1063     <T> void doTestIssue1485SameFrameDefinition(T frameOrLOF) {
1064         // GIVEN
1065         final AbsoluteDate      date           = AbsoluteDate.ARBITRARY_EPOCH;
1066         final RealMatrix        expectedMatrix = MatrixUtils.createRealIdentityMatrix(6);
1067         final PositionAngleType anomalyType    = PositionAngleType.MEAN;
1068         final Frame             frame          = FramesFactory.getEME2000();
1069         final double            mu             = Constants.EGM96_EARTH_MU;
1070 
1071         final KeplerianOrbit orbit = new KeplerianOrbit(1e7, 0.01, 1., 2., 3., 4., anomalyType, frame, date, mu);
1072 
1073         // WHEN
1074         final StateCovariance originalCovariance;
1075         final StateCovariance covariance;
1076         if (frameOrLOF instanceof LOF) {
1077             originalCovariance = new StateCovariance(expectedMatrix, date, (LOF) frameOrLOF);
1078             covariance         = originalCovariance.changeCovarianceFrame(orbit, (LOF) frameOrLOF);
1079         } else {
1080             originalCovariance = new StateCovariance(expectedMatrix, date, (Frame) frameOrLOF, orbit.getType(), anomalyType);
1081             covariance         = originalCovariance.changeCovarianceFrame(orbit, (Frame) frameOrLOF);
1082         }
1083 
1084         // THEN
1085         final RealMatrix actualMatrix = covariance.getMatrix();
1086         assertEquals(0., actualMatrix.subtract(expectedMatrix).getNorm1());
1087     }
1088 
1089     private NumericalPropagator buildDefaultPropagator(final SpacecraftState state,
1090                                                        final OrbitType orbitType,
1091                                                        final PositionAngleType positionAngleType) {
1092 
1093         // Build default ODEIntegrator
1094         final ODEIntegrator integrator = buildDefaultODEIntegrator(state.getOrbit(), orbitType);
1095 
1096         // Build and initialize numerical propagator
1097         final NumericalPropagator propagator = new NumericalPropagator(integrator);
1098         propagator.setInitialState(state);
1099         propagator.setOrbitType(orbitType);
1100         propagator.setPositionAngleType(positionAngleType);
1101 
1102         return propagator;
1103     }
1104 
1105     private ODEIntegrator buildDefaultODEIntegrator(final Orbit orbit, final OrbitType orbitType) {
1106         final double     dP         = 1;
1107         final double     minStep    = 0.001;
1108         final double     maxStep    = 60;
1109         final double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(dP).getTolerances(orbit, orbitType);
1110 
1111         return new DormandPrince853Integrator(minStep, maxStep, tolerances[0], tolerances[1]);
1112     }
1113 
1114 }
1115