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