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