1   /* Copyright 2002-2025 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.propagation;
18  
19  import java.util.Arrays;
20  
21  import org.hipparchus.geometry.euclidean.threed.Vector3D;
22  import org.hipparchus.linear.MatrixUtils;
23  import org.hipparchus.linear.RealMatrix;
24  import org.hipparchus.ode.ODEIntegrator;
25  import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator;
26  import org.hipparchus.util.FastMath;
27  import org.junit.jupiter.api.Assertions;
28  import org.junit.jupiter.api.BeforeEach;
29  import org.junit.jupiter.api.Test;
30  import org.junit.jupiter.params.ParameterizedTest;
31  import org.junit.jupiter.params.provider.EnumSource;
32  import org.orekit.Utils;
33  import org.orekit.forces.ForceModel;
34  import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
35  import org.orekit.forces.gravity.J2OnlyPerturbation;
36  import org.orekit.forces.gravity.potential.GravityFieldFactory;
37  import org.orekit.forces.gravity.potential.ICGEMFormatReader;
38  import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
39  import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
40  import org.orekit.frames.Frame;
41  import org.orekit.frames.FramesFactory;
42  import org.orekit.frames.LOFType;
43  import org.orekit.orbits.CartesianOrbit;
44  import org.orekit.orbits.Orbit;
45  import org.orekit.orbits.OrbitType;
46  import org.orekit.orbits.PositionAngleType;
47  import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
48  import org.orekit.propagation.analytical.KeplerianPropagator;
49  import org.orekit.propagation.numerical.NumericalPropagator;
50  import org.orekit.propagation.semianalytical.dsst.DSSTHarvester;
51  import org.orekit.propagation.semianalytical.dsst.DSSTPropagator;
52  import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
53  import org.orekit.propagation.semianalytical.dsst.forces.DSSTJ2SquaredClosedForm;
54  import org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal;
55  import org.orekit.propagation.semianalytical.dsst.forces.ZeisModel;
56  import org.orekit.time.AbsoluteDate;
57  import org.orekit.time.TimeScalesFactory;
58  import org.orekit.utils.Constants;
59  import org.orekit.utils.IERSConventions;
60  import org.orekit.utils.PVCoordinates;
61  
62  class StateCovarianceMatrixProviderTest {
63  
64      /** Initial S/C state. */
65      private SpacecraftState initialState;
66      
67      /** Initial covariance. */
68      private RealMatrix  initCov;
69      
70      /** Reference covariance after 60s prop with a J2 numerical model (computed with another solution). */
71      private RealMatrix  refCovAfter60s;
72      
73      @BeforeEach
74      void setUp() {
75          Utils.setDataRoot("orbit-determination/february-2016:potential/icgem-format");
76          GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true));
77          Orbit initialOrbit = new CartesianOrbit(
78                  new PVCoordinates(new Vector3D(7526993.581890527, -9646310.10026971, 1464110.4928112086),
79                                    new Vector3D(3033.79456099698, 1715.265069098717, -4447.658745923895)),
80                  FramesFactory.getEME2000(),
81                  new AbsoluteDate("2016-02-13T16:00:00.000", TimeScalesFactory.getUTC()),
82                  Constants.WGS84_EARTH_MU);
83          initialState = new SpacecraftState(initialOrbit);
84          initCov = MatrixUtils.createRealMatrix( new double[][] {
85              { 8.651816029e+01, 5.689987127e+01, -2.763870764e+01, -2.435617201e-02, 2.058274137e-02, -5.872883051e-03 },
86              { 5.689987127e+01, 7.070624321e+01, 1.367120909e+01, -6.112622013e-03, 7.623626008e-03, -1.239413190e-02 },
87              { -2.763870764e+01, 1.367120909e+01, 1.811858898e+02, 3.143798992e-02, -4.963106559e-02, -7.420114385e-04 },
88              { -2.435617201e-02, -6.112622013e-03, 3.143798992e-02, 4.657077389e-05, 1.469943634e-05, 3.328475593e-05 },
89              { 2.058274137e-02, 7.623626008e-03, -4.963106559e-02, 1.469943634e-05, 3.950715934e-05, 2.516044258e-05 },
90              { -5.872883051e-03, -1.239413190e-02, -7.420114385e-04, 3.328475593e-05, 2.516044258e-05, 3.547466120e-05 }
91          });
92          
93          // Reference covariance after 60s prop with a J2 numerical model (computed with another solution)
94          refCovAfter60s = MatrixUtils.createRealMatrix( new double[][] {
95              { 5.770543135e+02, 2.316979550e+02, -5.172369105e+02, -2.585893247e-01, 2.113809017e-01, -1.759509343e-01 },
96              { 2.316979550e+02, 1.182942930e+02, -1.788422178e+02, -9.570305681e-02, 7.792155309e-02, -7.435822327e-02 },
97              { -5.172369105e+02, -1.788422178e+02, 6.996248500e+02, 2.633605389e-01, -2.480144888e-01, 1.908427233e-01 },
98              { -2.585893247e-01, -9.570305681e-02, 2.633605389e-01, 1.419148897e-04, -8.715858320e-05, 1.024944399e-04 },
99              { 2.113809017e-01, 7.792155309e-02, -2.480144888e-01, -8.715858320e-05, 1.069566588e-04, -5.667563856e-05 },
100             { -1.759509343e-01, -7.435822327e-02, 1.908427233e-01, 1.024944399e-04, -5.667563856e-05, 8.178356868e-05 }
101         });
102     }
103 
104     /**
105      * Absolute comparison of two covariance matrices.
106      *
107      * @param reference reference covariance
108      * @param computed  computed covariance
109      * @param threshold threshold for comparison
110      */
111     private void compareCovariance(final RealMatrix reference, final RealMatrix computed, final double threshold) {
112         for (int row = 0; row < reference.getRowDimension(); row++) {
113             for (int column = 0; column < reference.getColumnDimension(); column++) {
114                 if (reference.getEntry(row, column) == 0) {
115                     Assertions.assertEquals(reference.getEntry(row, column), computed.getEntry(row, column),
116                                             threshold);
117                 }
118                 else {
119                     Assertions.assertEquals(reference.getEntry(row, column), computed.getEntry(row, column),
120                                             FastMath.abs(threshold * reference.getEntry(row, column)));
121                 }
122             }
123         }
124     }
125 
126     /**
127      * Unit test for covariance propagation in Cartesian elements.
128      */
129     @Test
130     void testWithNumericalPropagatorCartesian() {
131         // Integrator
132         final double        step       = 60.0;
133         final ODEIntegrator integrator = new ClassicalRungeKuttaIntegrator(step);
134 
135         // Numerical propagator
136         final OrbitType           propType   = OrbitType.CARTESIAN;
137         final PositionAngleType angleType  = PositionAngleType.MEAN;
138         final NumericalPropagator propagator = new NumericalPropagator(integrator);
139         // Add a force model
140         final NormalizedSphericalHarmonicsProvider gravity = GravityFieldFactory.getNormalizedProvider(2, 0);
141         final ForceModel holmesFeatherstone =
142                 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), gravity);
143         propagator.addForceModel(holmesFeatherstone);
144         
145         // Finalize setting
146         propagator.setOrbitType(propType);
147         propagator.setPositionAngleType(angleType);
148         propagator.setInitialState(initialState);
149 
150         // Configure covariance propagation
151         final StateCovarianceMatrixProvider provider = setUpCovariancePropagation(propagator);
152         
153 
154         // Propagate
155         final SpacecraftState propagated = propagator.propagate(initialState.getDate().shiftedBy(Constants.JULIAN_DAY));
156 
157         // Get the propagated covariance
158         final StateCovariance propagatedStateCov = provider.getStateCovariance(propagated);
159         final RealMatrix propagatedCov = propagatedStateCov.getMatrix();
160 
161         // Verify
162         compareCovariance(refCovAfter60s, propagatedCov, 4.0e-7);
163         Assertions.assertEquals(OrbitType.CARTESIAN, provider.getCovarianceOrbitType());
164         Assertions.assertEquals(OrbitType.CARTESIAN, propagatedStateCov.getOrbitType());
165         Assertions.assertNull(propagatedStateCov.getLOF());
166 
167         ///////////
168         // Test the frame transformation
169         ///////////
170 
171         // Define a new output frame
172         final Frame frameB = FramesFactory.getTEME();
173 
174         // Get the covariance in TEME frame
175         RealMatrix transformedCovA = provider.getStateCovariance(propagated, frameB).getMatrix();
176 
177         // Second transformation
178         RealMatrix transformedCovB =
179         		propagatedStateCov.changeCovarianceFrame(propagated.getOrbit(), frameB).getMatrix();
180 
181         // Verify
182         compareCovariance(transformedCovA, transformedCovB, 1.0e-15);
183 
184         ///////////
185         // Test the orbit type transformation
186         ///////////
187 
188         // Define a new output frame
189         final OrbitType     outOrbitType = OrbitType.KEPLERIAN;
190         final PositionAngleType outAngleType = PositionAngleType.MEAN;
191 
192         // Transformation using getStateJacobian() method
193         RealMatrix transformedCovC = provider.getStateCovariance(propagated, outOrbitType, outAngleType).getMatrix();
194 
195         // Second transformation
196         RealMatrix transformedCovD =
197         		propagatedStateCov.changeCovarianceType(propagated.getOrbit(), outOrbitType, outAngleType).getMatrix();
198 
199         // Verify
200         compareCovariance(transformedCovC, transformedCovD, 1.0e-15);
201 
202     }
203 
204     /**
205      * Unit test for covariance propagation in Cartesian elements.
206      * This test verifies the mechanism of AdditionalDataProvider for a RealMatrix.
207      */
208     @Test
209     void testAdditionalDataProvider() {
210         // Integrator
211         final double step = 60.0;
212         final ODEIntegrator integrator = new ClassicalRungeKuttaIntegrator(step);
213 
214         // Numerical propagator
215         final NumericalPropagator propagator = new NumericalPropagator(integrator);
216         // Add a force model
217         final NormalizedSphericalHarmonicsProvider gravity = GravityFieldFactory.getNormalizedProvider(2, 0);
218         final ForceModel j2OnlyPerturbation = new J2OnlyPerturbation(GravityFieldFactory.getUnnormalizedProvider(gravity),
219                 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
220         propagator.addForceModel(j2OnlyPerturbation);
221         propagator.setInitialState(initialState);
222 
223         // Configure covariance propagation
224         final StateCovarianceMatrixProvider provider = setUpCovariancePropagation(propagator);
225 
226         // Propagate
227         final SpacecraftState propagated = propagator.propagate(initialState.getDate().shiftedBy(Constants.JULIAN_DAY));
228 
229         // Get the propagated covariance
230         final RealMatrix propagatedCov = provider.getAdditionalData(propagated);
231 
232         // Verify
233         compareCovariance(refCovAfter60s, propagatedCov, 3.0e-5);
234         Assertions.assertEquals(OrbitType.CARTESIAN, provider.getCovarianceOrbitType());
235     }
236 
237     @ParameterizedTest
238     @EnumSource(value = LOFType.class, names = {"QSW", "NTW", "LVLH", "TNW"})
239     void testStmLof(final LOFType lofType) {
240         final NumericalPropagator propagator = new NumericalPropagator(new ClassicalRungeKuttaIntegrator(100));
241         propagator.resetInitialState(initialState);
242         propagator.setOrbitType(OrbitType.EQUINOCTIAL);
243         final MatricesHarvester harvester = propagator.setupMatricesComputation("stm", null, null);
244         final StateCovariance initialCovariance = new StateCovariance(MatrixUtils.createRealIdentityMatrix(6).scalarMultiply(1e-2),
245                 initialState.getDate(), lofType);
246         final StateCovarianceMatrixProvider provider =
247                 new StateCovarianceMatrixProvider("cov", "stm", harvester, initialCovariance);
248         propagator.setInitialState(initialState);
249         propagator.addAdditionalDataProvider(provider);
250         // WHEN
251         final SpacecraftState sameState = propagator.propagate(initialState.getDate());
252         // THEN
253         final StateCovariance actualCovariance = provider.getStateCovariance(sameState);
254         Assertions.assertEquals(initialCovariance.getDate(), actualCovariance.getDate());
255         Assertions.assertEquals(initialCovariance.getLOF(), actualCovariance.getLOF());
256         Assertions.assertEquals(initialCovariance.getFrame(), actualCovariance.getFrame());
257         Assertions.assertEquals(initialCovariance.getOrbitType(), actualCovariance.getOrbitType());
258     }
259 
260     /**
261      * Unit test for covariance propagation in Cartesian elements. The difference here is that the propagator uses its
262      * default orbit type: EQUINOCTIAL
263      */
264     @Test
265     void testWithNumericalPropagatorDefault() {
266         // Integrator
267         final double        step       = 60.0;
268         final ODEIntegrator integrator = new ClassicalRungeKuttaIntegrator(step);
269 
270         // Numerical propagator
271         final NumericalPropagator propagator = new NumericalPropagator(integrator);
272         // Add a force model
273         final NormalizedSphericalHarmonicsProvider gravity = GravityFieldFactory.getNormalizedProvider(2, 0);
274         final ForceModel holmesFeatherstone =
275                 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), gravity);
276         propagator.addForceModel(holmesFeatherstone);
277         propagator.setInitialState(initialState);
278 
279         // Configure covariance propagation
280         final StateCovarianceMatrixProvider provider = setUpCovariancePropagation(propagator);
281 
282         // Propagate
283         final SpacecraftState propagated = propagator.propagate(initialState.getDate().shiftedBy(Constants.JULIAN_DAY));
284 
285         // Get the propagated covariance
286         final StateCovariance propagatedStateCov = provider.getStateCovariance(propagated);
287         final RealMatrix propagatedCov = propagatedStateCov.getMatrix();
288 
289         // Verify
290         compareCovariance(refCovAfter60s, propagatedCov, 3.0e-5);
291         Assertions.assertEquals(OrbitType.CARTESIAN, provider.getCovarianceOrbitType());
292 
293         ///////////
294         // Test the frame transformation
295         ///////////
296 
297         // Define a new output frame
298         final Frame frameB = FramesFactory.getTEME();
299 
300         // Get the covariance in TEME frame
301         RealMatrix transformedCovA = provider.getStateCovariance(propagated, frameB).getMatrix();
302 
303         // Second transformation
304         RealMatrix transformedCovB =
305         		propagatedStateCov.changeCovarianceFrame(propagated.getOrbit(), frameB).getMatrix();
306 
307         // Verify
308         compareCovariance(transformedCovA, transformedCovB, 1.0e-15);
309 
310         ///////////
311         // Test the orbit type transformation
312         ///////////
313 
314         // Define a new output frame
315         final OrbitType     outOrbitType = OrbitType.KEPLERIAN;
316         final PositionAngleType outAngleType = PositionAngleType.MEAN;
317 
318         // Transformation using getStateJacobian() method
319         RealMatrix transformedCovC = provider.getStateCovariance(propagated, outOrbitType, outAngleType).getMatrix();
320 
321         // Second transformation
322         RealMatrix transformedCovD =
323         		propagatedStateCov.changeCovarianceType(propagated.getOrbit(), outOrbitType, outAngleType).getMatrix();
324 
325         // Verify
326         compareCovariance(transformedCovC, transformedCovD, 1.0e-15);
327 
328     }
329 
330     /**
331      * Unit test for covariance propagation in Keplerian elements. The difference here is that the propagator uses its
332      * default orbit type: EQUINOCTIAL.
333      * <p>
334      * The additional purpose of this test is to make sure that the propagated state covariance is expressed in the right
335      * orbit type.
336      */
337     @Test
338     void testWithNumericalPropagatorDefaultAndKeplerianOrbitType() {
339         // Integrator
340         final double        step       = 60.0;
341         final ODEIntegrator integrator = new ClassicalRungeKuttaIntegrator(step);
342 
343         // Numerical propagator
344         final String              stmName    = "STM";
345         final NumericalPropagator propagator = new NumericalPropagator(integrator);
346         // Add a force model
347         final NormalizedSphericalHarmonicsProvider gravity = GravityFieldFactory.getNormalizedProvider(2, 0);
348         final ForceModel holmesFeatherstone =
349                 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), gravity);
350         propagator.addForceModel(holmesFeatherstone);
351         // Finalize setting
352         final MatricesHarvester harvester = propagator.setupMatricesComputation(stmName, null, null);
353 
354         // Create additional state
355         final String     additionalName = "cartCov";
356         final StateCovariance initialStateCovariance = new StateCovariance(initCov, initialState.getDate(), initialState.getFrame(), OrbitType.CARTESIAN, PositionAngleType.MEAN);
357         final StateCovariance initialStateCovarianceInKep = initialStateCovariance.changeCovarianceType(initialState.getOrbit(), OrbitType.KEPLERIAN, PositionAngleType.MEAN);
358         final StateCovarianceMatrixProvider provider =
359                 new StateCovarianceMatrixProvider(additionalName, stmName, harvester, initialStateCovarianceInKep);
360         propagator.setInitialState(initialState);
361         propagator.addAdditionalDataProvider(provider);
362 
363         // Propagate
364         final SpacecraftState propagated = propagator.propagate(initialState.getDate().shiftedBy(Constants.JULIAN_DAY));
365 
366         // Get the propagated covariance
367         final StateCovariance propagatedStateCov = provider.getStateCovariance(propagated);
368         final StateCovariance propagatedStateCovInCart = propagatedStateCov.changeCovarianceType(propagated.getOrbit(), OrbitType.CARTESIAN, PositionAngleType.MEAN);
369         final RealMatrix propagatedCovInCart = propagatedStateCovInCart.getMatrix();
370 
371         // Verify
372         compareCovariance(refCovAfter60s, propagatedCovInCart, 3.0e-5);
373         Assertions.assertEquals(OrbitType.KEPLERIAN, provider.getCovarianceOrbitType());
374         Assertions.assertEquals(OrbitType.KEPLERIAN, propagatedStateCov.getOrbitType());
375 
376     }
377 
378     /**
379      * Unit test for covariance propagation in Cartesian elements.
380      */
381     @Test
382     void testWithAnalyticalPropagator() {
383         // Numerical propagator
384         final EcksteinHechlerPropagator propagator = new EcksteinHechlerPropagator(initialState.getOrbit(),
385                                                                                    GravityFieldFactory.getUnnormalizedProvider(6, 0));
386 
387         // Configure covariance propagation
388         final StateCovarianceMatrixProvider provider = setUpCovariancePropagation(propagator);
389 
390         // Propagate
391         final SpacecraftState propagated = propagator.propagate(initialState.getDate().shiftedBy(Constants.JULIAN_DAY));
392 
393         // Get the propagated covariance
394         final StateCovariance propagatedStateCov = provider.getStateCovariance(propagated);
395         final RealMatrix propagatedCov = propagatedStateCov.getMatrix();
396 
397         // Verify
398         compareCovariance(refCovAfter60s, propagatedCov, 5.0e-4);
399         Assertions.assertEquals(OrbitType.CARTESIAN, provider.getCovarianceOrbitType());
400 
401         ///////////
402         // Test the frame transformation
403         ///////////
404 
405         // Define a new output frame
406         final Frame frameB = FramesFactory.getTEME();
407 
408         // Get the covariance in TEME frame
409         RealMatrix transformedCovA = provider.getStateCovariance(propagated, frameB).getMatrix();
410 
411         // Second transformation
412         RealMatrix transformedCovB =
413         		propagatedStateCov.changeCovarianceFrame(propagated.getOrbit(), frameB).getMatrix();
414 
415         // Verify
416         compareCovariance(transformedCovA, transformedCovB, 1.0e-15);
417 
418         ///////////
419         // Test the orbit type transformation
420         ///////////
421 
422         // Define a new output frame
423         final OrbitType     outOrbitType = OrbitType.KEPLERIAN;
424         final PositionAngleType outAngleType = PositionAngleType.MEAN;
425 
426         // Transformation using getStateJacobian() method
427         RealMatrix transformedCovC = provider.getStateCovariance(propagated, outOrbitType, outAngleType).getMatrix();
428 
429         // Second transformation
430         RealMatrix transformedCovD =
431         		propagatedStateCov.changeCovarianceType(propagated.getOrbit(), outOrbitType, outAngleType).getMatrix();
432 
433         // Verify
434         compareCovariance(transformedCovC, transformedCovD, 1.0e-15);
435 
436     }
437 
438     /**
439      * Unit test for covariance propagation with DSST propagator.
440      */
441     @Test
442     void testWithDSSTPropagatorDefault() {
443         // Integrator
444         final double        step       = 3600.0;
445         final ODEIntegrator integrator = new ClassicalRungeKuttaIntegrator(step);
446 
447         // DSST propagator
448         final String         stmName    = "STM";
449         final DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING);
450         // Add a force model
451         final UnnormalizedSphericalHarmonicsProvider gravity = GravityFieldFactory.getUnnormalizedProvider(2, 0);
452         final DSSTForceModel zonal = new DSSTZonal(gravity);
453         propagator.addForceModel(zonal);
454         propagator.addForceModel(new DSSTJ2SquaredClosedForm(new ZeisModel(), gravity));
455         // Finalize setting
456         final DSSTHarvester harvester = (DSSTHarvester) propagator.setupMatricesComputation(stmName, null, null);
457         harvester.initializeFieldShortPeriodTerms(DSSTPropagator.computeMeanState(initialState, propagator.getAttitudeProvider(), Arrays.asList(zonal)));
458 
459         // Create additional state
460         final String     additionalName = "cartCov";
461         final StateCovariance initialStateCovariance = new StateCovariance(initCov, initialState.getDate(), initialState.getFrame(), OrbitType.CARTESIAN, PositionAngleType.MEAN);
462         final StateCovarianceMatrixProvider provider =
463                 new StateCovarianceMatrixProvider(additionalName, stmName, harvester, initialStateCovariance);
464         propagator.setInitialState(initialState);
465         propagator.addAdditionalDataProvider(provider);
466 
467         // Propagate
468         final SpacecraftState propagated = propagator.propagate(initialState.getDate().shiftedBy(Constants.JULIAN_DAY));
469 
470         // Get the propagated covariance
471         final StateCovariance propagatedStateCov = provider.getStateCovariance(propagated);
472         final RealMatrix propagatedCov = propagatedStateCov.getMatrix();
473 
474         // Verify (3% error with respect to reference)
475         compareCovariance(refCovAfter60s, propagatedCov, 0.03);
476         Assertions.assertEquals(OrbitType.CARTESIAN, provider.getCovarianceOrbitType());
477     }
478 
479     /**
480      * Unit test for shiftedBy() method.
481      * The method is compared to covariance propagation using the Keplerian propagator.
482      */
483     @Test
484     void testCovarianceShift() {
485         // Keplerian propagator
486         final KeplerianPropagator propagator = new KeplerianPropagator(initialState.getOrbit());
487         final double dt = 60.0;
488 
489         // Configure covariance propagation
490         final StateCovarianceMatrixProvider provider = setUpCovariancePropagation(propagator);
491 
492         // Propagate
493         final SpacecraftState propagated = propagator.propagate(initialState.getDate().shiftedBy(dt));
494 
495         // Get the propagated covariance
496         final StateCovariance propagatedStateCov = provider.getStateCovariance(propagated);
497         final RealMatrix propagatedCov = propagatedStateCov.getMatrix();
498 
499         // Use of shiftedBy
500         final StateCovariance initialStateCovariance = new StateCovariance(initCov, initialState.getDate(), initialState.getFrame(), OrbitType.CARTESIAN, PositionAngleType.MEAN);
501         final StateCovariance shiftedStateCov = initialStateCovariance.shiftedBy(initialState.getOrbit(), dt);
502         final RealMatrix shiftedCov = shiftedStateCov.getMatrix();
503 
504         // Verify
505         compareCovariance(propagatedCov, shiftedCov, 4.0e-12);
506         Assertions.assertEquals(propagatedStateCov.getDate(), shiftedStateCov.getDate());
507         Assertions.assertEquals(propagatedStateCov.getOrbitType(), shiftedStateCov.getOrbitType());
508         Assertions.assertEquals(propagatedStateCov.getPositionAngleType(), shiftedStateCov.getPositionAngleType());
509     }
510 
511     /**
512      * Unit test for issue 1253: incorrect covariance propagation with a bounded propagator.
513      *
514      * <p>Bug detection courtesy of Christophe Le Bris
515      *
516      * <p>The problem was that, when using ephemeris mode with an integrated propagator, the covariance propagated with the
517      * output BoundedPropagator was wrong. Actually, it was always equal to the initial covariance matrix.
518      */
519     @Test
520     void testIssue1253_IntegratedPropagator() {
521 
522         // GIVEN
523         // -----
524 
525         // dt
526         final double dt = 3600.0;
527         final AbsoluteDate propDate = initialState.getDate().shiftedBy(dt);
528 
529         // Numerical propagator
530         final OrbitType           propType   = OrbitType.CARTESIAN;
531         final PositionAngleType angleType  = PositionAngleType.MEAN;
532         NumericalPropagator propagator = new NumericalPropagator(new ClassicalRungeKuttaIntegrator(60.));
533         propagator.setOrbitType(propType);
534         propagator.setPositionAngleType(angleType);
535         propagator.setInitialState(initialState);
536         StateCovarianceMatrixProvider provider = setUpCovariancePropagation(propagator);
537 
538         // Propagate to dt/2
539         SpacecraftState propagated = propagator.propagate(propDate);
540 
541         // Get the reference propagated covariance
542         final StateCovariance refPropagatedStateCov = provider.getStateCovariance(propagated);
543         final RealMatrix refPropagatedCov = refPropagatedStateCov.getMatrix();
544         
545         // WHEN
546         // ----
547         
548         // 2. Ephemeris on [t0, t0 + dt]
549         propagator = new NumericalPropagator(new ClassicalRungeKuttaIntegrator(60.));
550         propagator.setOrbitType(propType);
551         propagator.setPositionAngleType(angleType);
552         propagator.setInitialState(initialState);
553         provider = setUpCovariancePropagation(propagator);
554 
555         // Propagate to dt and store ephemeris
556         final EphemerisGenerator generator = propagator.getEphemerisGenerator();
557         propagator.propagate(propDate);
558         final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
559         
560         // 3. Propagate ephemeris to dt
561         propagated = ephemeris.propagate(propDate);
562         final StateCovariance propagatedStateCov = provider.getStateCovariance(propagated);
563         final RealMatrix propagatedCov = provider.getStateCovariance(propagated).getMatrix();
564         
565         
566         // THEN
567         // ----
568 
569         // Verify that both covariances are equal
570         Assertions.assertEquals(refPropagatedStateCov.getDate(), propagatedStateCov.getDate());
571         Assertions.assertEquals(refPropagatedStateCov.getOrbitType(), propagatedStateCov.getOrbitType());
572         Assertions.assertEquals(refPropagatedStateCov.getPositionAngleType(), propagatedStateCov.getPositionAngleType());
573         compareCovariance(refPropagatedCov, propagatedCov, 0.);
574     }
575 
576     /**
577      * Unit test for issue 1253: incorrect covariance propagation with a bounded propagator.
578      * <p>
579      * With analytical propagators, the former code worked. Test if it still works with the new version of the code.
580      */
581     @Test
582     void testIssue1253_AnalyticalPropagator() {
583 
584         // GIVEN
585         // -----
586 
587         // dt
588         final double dt = 3600.0;
589         final AbsoluteDate propDate = initialState.getDate().shiftedBy(dt);
590 
591         // Keplerian propagator
592         Propagator propagator = new KeplerianPropagator(initialState.getOrbit());
593         StateCovarianceMatrixProvider provider = setUpCovariancePropagation(propagator);
594 
595         // Propagate to dt/2
596         SpacecraftState propagated = propagator.propagate(propDate);
597 
598         // Get the reference propagated covariance
599         final StateCovariance refPropagatedStateCov = provider.getStateCovariance(propagated);
600         final RealMatrix refPropagatedCov = refPropagatedStateCov.getMatrix();
601         
602         // WHEN
603         // ----
604         
605         // 2. Ephemeris on [t0, t0 + dt]
606         propagator = new KeplerianPropagator(initialState.getOrbit());
607         provider = setUpCovariancePropagation(propagator);
608 
609         // Propagate to dt and store ephemeris
610         final EphemerisGenerator generator = propagator.getEphemerisGenerator();
611         propagator.propagate(propDate);
612         final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
613         
614         // 3. Propagate ephemeris to dt
615         propagated = ephemeris.propagate(propDate);
616         final StateCovariance propagatedStateCov = provider.getStateCovariance(propagated);
617         final RealMatrix propagatedCov = provider.getStateCovariance(propagated).getMatrix();
618         
619         
620         // THEN
621         // ----
622 
623         // Verify that both covariances are equal
624         Assertions.assertEquals(refPropagatedStateCov.getDate(), propagatedStateCov.getDate());
625         Assertions.assertEquals(refPropagatedStateCov.getOrbitType(), propagatedStateCov.getOrbitType());
626         Assertions.assertEquals(refPropagatedStateCov.getPositionAngleType(), propagatedStateCov.getPositionAngleType());
627         compareCovariance(refPropagatedCov, propagatedCov, 1e-10);
628     }
629     
630     /** Setup a StateCovarianceMatrixProvider without changing orbit type.
631      * 
632      * @param propagator input propagator
633      * @return configured StateCovarianceMatrixProvider
634      */
635     private StateCovarianceMatrixProvider setUpCovariancePropagation(final Propagator propagator) {
636         
637         final String stmName = "STM";
638         final String     additionalName = "cartCov";
639         
640         // Finalize setting
641         final MatricesHarvester harvester = propagator.setupMatricesComputation(stmName, null, null);
642 
643         // Create additional state
644         final StateCovariance initialStateCovariance = new StateCovariance(initCov, initialState.getDate(), initialState.getFrame(), OrbitType.CARTESIAN, PositionAngleType.MEAN);
645         final StateCovarianceMatrixProvider provider =
646                 new StateCovarianceMatrixProvider(additionalName, stmName, harvester, initialStateCovariance);
647         propagator.addAdditionalDataProvider(provider);
648         
649         return provider;
650     }
651 
652     @Test
653     void testIssue1933() {
654         // Numerical propagator with EQUINOCTIAL ELEMENTS
655         NumericalPropagator propagator = new NumericalPropagator(new ClassicalRungeKuttaIntegrator(60.));
656         propagator.setOrbitType(OrbitType.EQUINOCTIAL);
657         propagator.setInitialState(initialState);
658 
659         // Create covariance matrix provider with initial covariance in CARTESIAN elements
660         final StateCovarianceMatrixProvider provider = setUpCovariancePropagation(propagator);
661 
662         // Propagate and store ephemeris on [t0 - dt, t0 + dt]
663         final AbsoluteDate initialDate = initialState.getDate();
664         final EphemerisGenerator generator = propagator.getEphemerisGenerator();
665         final double dt = 3600.0;
666         propagator.propagate(initialDate.shiftedBy(-dt), initialDate.shiftedBy(dt));
667         final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
668 
669         // Propagate ephemeris
670         SpacecraftState propagated = ephemeris.propagate(initialDate);
671         final StateCovariance propagatedStateCov = provider.getStateCovariance(propagated);
672 
673         // Verify that both covariances are equal
674         Assertions.assertEquals(initialDate, propagatedStateCov.getDate());
675         Assertions.assertEquals(provider.getCovarianceOrbitType(), propagatedStateCov.getOrbitType());
676         Assertions.assertEquals(0., initCov.subtract(propagatedStateCov.getMatrix()).getNorm1(), 1e-6);
677     }
678 }