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