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