1   /* Copyright 2002-2025 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.propagation;
18  
19  import org.hipparchus.linear.Array2DRowRealMatrix;
20  import org.hipparchus.linear.MatrixUtils;
21  import org.hipparchus.linear.RealMatrix;
22  import org.orekit.errors.OrekitException;
23  import org.orekit.errors.OrekitMessages;
24  import org.orekit.frames.Frame;
25  import org.orekit.frames.KinematicTransform;
26  import org.orekit.frames.LOF;
27  import org.orekit.orbits.Orbit;
28  import org.orekit.orbits.OrbitType;
29  import org.orekit.orbits.PositionAngleType;
30  import org.orekit.time.AbsoluteDate;
31  import org.orekit.time.TimeStamped;
32  import org.orekit.utils.CartesianCovarianceUtils;
33  
34  /** This class is the representation of a covariance matrix at a given date.
35   * <p>
36   * Currently, the covariance only represents the orbital elements.
37   * <p>
38   * It is possible to change the covariance frame by using the
39   * {@link #changeCovarianceFrame(Orbit, Frame)} or {@link #changeCovarianceFrame(Orbit, LOF)} method.
40   * These methods are based on Equations (18) and (20) of <i>Covariance Transformations for Satellite
41   * Flight Dynamics Operations</i> by David A. SVallado.
42   * <p>
43   * Finally, covariance orbit type can be changed using the
44   * {@link #changeCovarianceType(Orbit, OrbitType, PositionAngleType)} method.
45   *
46   * @author Bryan Cazabonne
47   * @author Vincent Cucchietti
48   * @since 11.3
49   */
50  public class StateCovariance implements TimeStamped {
51  
52      /** State dimension. */
53      public static final int STATE_DIMENSION = 6;
54  
55      /** Default position angle for covariance expressed in Cartesian elements. */
56      private static final PositionAngleType DEFAULT_POSITION_ANGLE = PositionAngleType.TRUE;
57  
58      /** Orbital covariance [6x6]. */
59      private final RealMatrix orbitalCovariance;
60  
61      /** Covariance frame (can be null if LOF is defined). */
62      private final Frame frame;
63  
64      /** Covariance LOF type (can be null if frame is defined). */
65      private final LOF lof;
66  
67      /** Covariance epoch. */
68      private final AbsoluteDate epoch;
69  
70      /** Covariance orbit type. */
71      private final OrbitType orbitType;
72  
73      /** Covariance position angle type (not used if orbitType is CARTESIAN). */
74      private final PositionAngleType angleType;
75  
76      /**
77       * Constructor.
78       * @param orbitalCovariance 6x6 orbital parameters covariance
79       * @param epoch epoch of the covariance
80       * @param lof covariance LOF type
81       */
82      public StateCovariance(final RealMatrix orbitalCovariance, final AbsoluteDate epoch, final LOF lof) {
83          this(orbitalCovariance, epoch, null, lof, OrbitType.CARTESIAN, DEFAULT_POSITION_ANGLE);
84      }
85  
86      /**
87       * Constructor.
88       * @param orbitalCovariance 6x6 orbital parameters covariance
89       * @param epoch epoch of the covariance
90       * @param covarianceFrame covariance frame (inertial or Earth fixed)
91       * @param orbitType orbit type of the covariance (CARTESIAN if covarianceFrame is not pseudo-inertial)
92       * @param angleType position angle type of the covariance (not used if orbitType is CARTESIAN)
93       */
94      public StateCovariance(final RealMatrix orbitalCovariance, final AbsoluteDate epoch,
95                             final Frame covarianceFrame,
96                             final OrbitType orbitType, final PositionAngleType angleType) {
97          this(orbitalCovariance, epoch, covarianceFrame, null, orbitType, angleType);
98      }
99  
100     /**
101      * Private constructor.
102      * @param orbitalCovariance 6x6 orbital parameters covariance
103      * @param epoch epoch of the covariance
104      * @param covarianceFrame covariance frame (inertial or Earth fixed)
105      * @param lof covariance LOF type
106      * @param orbitType orbit type of the covariance
107      * @param angleType position angle type of the covariance (not used if orbitType is CARTESIAN)
108      */
109     private StateCovariance(final RealMatrix orbitalCovariance, final AbsoluteDate epoch,
110                             final Frame covarianceFrame, final LOF lof,
111                             final OrbitType orbitType, final PositionAngleType angleType) {
112 
113         checkFrameAndOrbitTypeConsistency(covarianceFrame, orbitType);
114 
115         this.orbitalCovariance = orbitalCovariance;
116         this.epoch = epoch;
117         this.frame     = covarianceFrame;
118         this.lof       = lof;
119         this.orbitType = orbitType;
120         this.angleType = angleType;
121 
122     }
123 
124     /**
125      * Check constructor's inputs consistency.
126      *
127      * @param covarianceFrame covariance frame (inertial or Earth fixed)
128      * @param inputType orbit type of the covariance
129      *
130      * @throws OrekitException if input frame is not pseudo-inertial AND the orbit type is not Cartesian
131      */
132     public static void checkFrameAndOrbitTypeConsistency(final Frame covarianceFrame, final OrbitType inputType) {
133 
134         // State covariance expressed in a celestial body frame
135         if (covarianceFrame != null) {
136 
137             // Input frame is not pseudo-inertial
138             if (!covarianceFrame.isPseudoInertial() && inputType != OrbitType.CARTESIAN) {
139                 throw new OrekitException(OrekitMessages.WRONG_ORBIT_PARAMETERS_TYPE,
140                                           inputType.name(),
141                                           OrbitType.CARTESIAN.name());
142             }
143         }
144     }
145 
146     /**
147      * Checks if input/output orbit and angle types are identical.
148      *
149      * @param inOrbitType input orbit type
150      * @param inAngleType input angle type
151      * @param outOrbitType output orbit type
152      * @param outAngleType output angle type
153      * @return flag defining if input/output orbit and angle types are identical
154      */
155     public static boolean inputAndOutputAreIdentical(final OrbitType inOrbitType, final PositionAngleType inAngleType,
156                                                      final OrbitType outOrbitType, final PositionAngleType outAngleType) {
157         return inOrbitType == outOrbitType && inAngleType == outAngleType;
158     }
159 
160     /**
161      * Checks if input and output orbit types are both {@code OrbitType.CARTESIAN}.
162      *
163      * @param inOrbitType input orbit type
164      * @param outOrbitType output orbit type
165      * @return flag defining if input and output orbit types are both {@code OrbitType.CARTESIAN}
166      */
167     public static boolean inputAndOutputOrbitTypesAreCartesian(final OrbitType inOrbitType, final OrbitType outOrbitType) {
168         return inOrbitType == OrbitType.CARTESIAN && outOrbitType == OrbitType.CARTESIAN;
169     }
170 
171     /** {@inheritDoc}. */
172     @Override
173     public AbsoluteDate getDate() {
174         return epoch;
175     }
176 
177     /**
178      * Get the covariance matrix.
179      * @return the covariance matrix
180      */
181     public RealMatrix getMatrix() {
182         return orbitalCovariance;
183     }
184 
185     /**
186      * Get the covariance orbit type.
187      * @return the covariance orbit type
188      */
189     public OrbitType getOrbitType() {
190         return orbitType;
191     }
192 
193     /**
194      * Get the covariance angle type.
195      * @return the covariance angle type
196      */
197     public PositionAngleType getPositionAngleType() {
198         return angleType;
199     }
200 
201     /**
202      * Get the covariance frame.
203      * @return the covariance frame (can be null)
204      * @see #getLOF()
205      */
206     public Frame getFrame() {
207         return frame;
208     }
209 
210     /**
211      * Get the covariance LOF type.
212      * @return the covariance LOF type (can be null)
213      * @see #getFrame()
214      */
215     public LOF getLOF() {
216         return lof;
217     }
218 
219     /**
220      * Get the covariance matrix in another orbit type.
221      * <p>
222      * The covariance orbit type <b>cannot</b> be changed if the covariance
223      * matrix is expressed in a {@link LOF local orbital frame} or a
224      * non-pseudo inertial frame.
225      * <p>
226      * As this type change uses the jacobian matrix of the transformation, it introduces a linear approximation.
227      * Hence, the current covariance matrix <b>will not exactly match</b> the new linearized case and the
228      * distribution will not follow a generalized Gaussian distribution anymore.
229      * <p>
230      * This is based on equation (1) to (6) from "Vallado, D. A. (2004). Covariance transformations for satellite flight
231      * dynamics operations."
232      *
233      * @param orbit orbit to which the covariance matrix should correspond
234      * @param outOrbitType target orbit type of the state covariance matrix
235      * @param outAngleType target position angle type of the state covariance matrix
236      * @return a new covariance state, expressed in the target orbit type with the target position angle
237      * @see #changeCovarianceFrame(Orbit, Frame)
238      */
239     public StateCovariance changeCovarianceType(final Orbit orbit, final OrbitType outOrbitType,
240                                                 final PositionAngleType outAngleType) {
241 
242         // Handle case where the covariance is already expressed in the output type
243         if (outOrbitType == orbitType && (outOrbitType == OrbitType.CARTESIAN || outAngleType == angleType)) {
244             if (lof == null) {
245                 return new StateCovariance(orbitalCovariance, epoch, frame, orbitType, angleType);
246             }
247             else {
248                 return new StateCovariance(orbitalCovariance, epoch, lof);
249             }
250         }
251 
252         // Check if the covariance is expressed in a celestial body frame
253         if (frame != null) {
254 
255             // Check if the covariance is defined in an inertial frame
256             if (frame.isPseudoInertial()) {
257                 return changeTypeAndCreate(orbit, epoch, frame, orbitType, angleType, outOrbitType, outAngleType,
258                                            orbitalCovariance);
259             }
260 
261             // The covariance is not defined in an inertial frame. The orbit type cannot be changed
262             throw new OrekitException(OrekitMessages.CANNOT_CHANGE_COVARIANCE_TYPE_IF_DEFINED_IN_NON_INERTIAL_FRAME);
263 
264         }
265 
266         // The covariance is not expressed in a celestial body frame. The orbit type cannot be changed
267         throw new OrekitException(OrekitMessages.CANNOT_CHANGE_COVARIANCE_TYPE_IF_DEFINED_IN_LOF);
268 
269     }
270 
271     /**
272      * Get the covariance in a given local orbital frame.
273      * <p>
274      * Changing the covariance frame is a linear process, this method does not introduce approximation unless a change
275      * in covariance orbit type is required.
276      * <p>
277      * This is based on equation (18) to (20) "from Vallado, D. A. (2004). Covariance transformations for satellite
278      * flight dynamics operations."
279      *
280      * @param orbit orbit to which the covariance matrix should correspond
281      * @param lofOut output local orbital frame
282      * @return a new covariance state, expressed in the output local orbital frame
283      */
284     public StateCovariance changeCovarianceFrame(final Orbit orbit, final LOF lofOut) {
285 
286         // Verify current covariance frame
287         if (lof != null) {
288 
289             // Check specific case where output covariance will be the same
290             if (lofOut == lof) {
291                 return new StateCovariance(orbitalCovariance, epoch, lof);
292             }
293 
294             // Change the covariance local orbital frame
295             return changeFrameAndCreate(orbit, epoch, lof, lofOut, orbitalCovariance);
296 
297         } else {
298 
299             // Covariance is expressed in celestial body frame
300             return changeFrameAndCreate(orbit, epoch, frame, lofOut, orbitalCovariance, orbitType, angleType);
301 
302         }
303 
304     }
305 
306     /**
307      * Get the covariance in the output frame.
308      * <p>
309      * Changing the covariance frame is a linear process, this method does not introduce approximation unless a change
310      * in covariance orbit type is required.
311      * <p>
312      * This is based on equation (18) to (20) "from Vallado, D. A. (2004). Covariance transformations for satellite
313      * flight dynamics operations."
314      *
315      * @param orbit orbit to which the covariance matrix should correspond
316      * @param frameOut output frame
317      * @return a new covariance state, expressed in the output frame
318      */
319     public StateCovariance changeCovarianceFrame(final Orbit orbit, final Frame frameOut) {
320 
321         // Verify current covariance frame
322         if (lof != null) {
323 
324             // Covariance is expressed in local orbital frame
325             return changeFrameAndCreate(orbit, epoch, lof, frameOut, orbitalCovariance);
326 
327         } else {
328 
329             // Check specific case where output covariance will be the same
330             if (frame == frameOut) {
331                 return new StateCovariance(orbitalCovariance, epoch, frame, orbitType, angleType);
332             }
333 
334             // Change covariance frame
335             return changeFrameAndCreate(orbit, epoch, frame, frameOut, orbitalCovariance, orbitType, angleType);
336 
337         }
338 
339     }
340 
341     /**
342      * Get a time-shifted covariance matrix.
343      * <p>
344      * The shifting model is a linearized, Keplerian one. In other words, it is based on a state transition matrix that
345      * is computed assuming Keplerian motion.
346      * <p>
347      * Shifting is <em>not</em> intended as a replacement for proper covariance propagation, but should be sufficient
348      * for small time shifts or coarse accuracy.
349      *
350      * @param orbit orbit to which the covariance matrix should correspond
351      * @param dt time shift in seconds
352      * @return a new covariance state, shifted with respect to the instance
353      */
354     public StateCovariance shiftedBy(final Orbit orbit, final double dt) {
355 
356         // Shifted orbit
357         final Orbit shifted = orbit.shiftedBy(dt);
358 
359         // State covariance expressed in celestial body frame
360         if (frame != null) {
361 
362             // State covariance expressed in a pseudo-inertial frame
363             if (frame.isPseudoInertial()) {
364 
365                 // Compute STM
366                 final RealMatrix stm = getStm(orbit, dt);
367 
368                 // Convert covariance in STM type (i.e., Equinoctial elements)
369                 final StateCovariance inStmType = changeTypeAndCreate(orbit, epoch, frame, orbitType, angleType,
370                                                                       OrbitType.EQUINOCTIAL, PositionAngleType.MEAN,
371                                                                       orbitalCovariance);
372 
373                 // Shift covariance by applying the STM
374                 final RealMatrix shiftedCov = stm.multiply(inStmType.getMatrix().multiplyTransposed(stm));
375 
376                 // Restore the initial covariance type
377                 return changeTypeAndCreate(shifted, shifted.getDate(), frame,
378                                            OrbitType.EQUINOCTIAL, PositionAngleType.MEAN,
379                                            orbitType, angleType, shiftedCov);
380             }
381 
382             // State covariance expressed in a non pseudo-inertial frame
383             else {
384 
385                 // Convert state covariance to orbit pseudo-inertial frame
386                 final StateCovariance inOrbitFrame = changeCovarianceFrame(orbit, orbit.getFrame());
387 
388                 // Shift the state covariance
389                 final StateCovariance shiftedCovariance = inOrbitFrame.shiftedBy(orbit, dt);
390 
391                 // Restore the initial covariance frame
392                 return shiftedCovariance.changeCovarianceFrame(shifted, frame);
393             }
394         }
395 
396         // State covariance expressed in a commonly used local orbital frame (LOF)
397         else {
398 
399             // Convert state covariance to orbit pseudo-inertial frame
400             final StateCovariance inOrbitFrame = changeCovarianceFrame(orbit, orbit.getFrame());
401 
402             // Shift the state covariance
403             final StateCovariance shiftedCovariance = inOrbitFrame.shiftedBy(orbit, dt);
404 
405             // Restore the initial covariance frame
406             return shiftedCovariance.changeCovarianceFrame(shifted, lof);
407         }
408 
409     }
410 
411     /**
412      * Create a covariance matrix in another orbit type.
413      * <p>
414      * As this type change uses the jacobian matrix of the transformation, it introduces a linear approximation.
415      * Hence, the input covariance matrix <b>will not exactly match</b> the new linearized case and the
416      * distribution will not follow a generalized Gaussian distribution anymore.
417      * <p>
418      * This is based on equation (1) to (6) from "Vallado, D. A. (2004). Covariance transformations for satellite flight
419      * dynamics operations."
420      *
421      * @param orbit orbit to which the covariance matrix should correspond
422      * @param date covariance epoch
423      * @param covFrame covariance frame
424      * @param inOrbitType initial orbit type of the state covariance matrix
425      * @param inAngleType initial position angle type of the state covariance matrix
426      * @param outOrbitType target orbit type of the state covariance matrix
427      * @param outAngleType target position angle type of the state covariance matrix
428      * @param inputCov input covariance
429      * @return the covariance expressed in the target orbit type with the target position angle
430      */
431     private static StateCovariance changeTypeAndCreate(final Orbit orbit,
432                                                        final AbsoluteDate date,
433                                                        final Frame covFrame,
434                                                        final OrbitType inOrbitType, final PositionAngleType inAngleType,
435                                                        final OrbitType outOrbitType, final PositionAngleType outAngleType,
436                                                        final RealMatrix inputCov) {
437 
438         // Check if type change is really necessary, if not then return input covariance
439         if (inputAndOutputOrbitTypesAreCartesian(inOrbitType, outOrbitType) ||
440             inputAndOutputAreIdentical(inOrbitType, inAngleType, outOrbitType, outAngleType)) {
441             return new StateCovariance(inputCov, date, covFrame, inOrbitType, inAngleType);
442         }
443 
444         // Notations:
445         // I: Input orbit type
446         // O: Output orbit type
447         // C: Cartesian parameters
448 
449         // Compute dOutputdCartesian
450         final double[][] aOC               = new double[STATE_DIMENSION][STATE_DIMENSION];
451         final Orbit      orbitInOutputType = outOrbitType.convertType(orbit);
452         orbitInOutputType.getJacobianWrtCartesian(outAngleType, aOC);
453         final RealMatrix dOdC = new Array2DRowRealMatrix(aOC, false);
454 
455         // Compute dCartesiandInput
456         final double[][] aCI              = new double[STATE_DIMENSION][STATE_DIMENSION];
457         final Orbit      orbitInInputType = inOrbitType.convertType(orbit);
458         orbitInInputType.getJacobianWrtParameters(inAngleType, aCI);
459         final RealMatrix dCdI = new Array2DRowRealMatrix(aCI, false);
460 
461         // Compute dOutputdInput
462         final RealMatrix dOdI = dOdC.multiply(dCdI);
463 
464         // Conversion of the state covariance matrix in target orbit type
465         final RealMatrix outputCov = dOdI.multiply(inputCov.multiplyTransposed(dOdI));
466 
467         // Return the converted covariance
468         return new StateCovariance(outputCov, date, covFrame, outOrbitType, outAngleType);
469 
470     }
471 
472     /**
473      * Create a covariance matrix from a {@link LOF local orbital frame} to another
474      * {@link LOF local orbital frame}.
475      * <p>
476      * Changing the covariance frame is a linear process, this method does not introduce approximation.
477      * <p>
478      * The transformation is based on Equation (18) to (20) of "Covariance Transformations for Satellite Flight Dynamics
479      * Operations" by David A. Vallado.
480      *
481      * @param orbit orbit to which the covariance matrix should correspond
482      * @param date covariance epoch
483      * @param lofIn the local orbital frame in which the input covariance matrix is expressed
484      * @param lofOut the target local orbital frame
485      * @param inputCartesianCov input covariance {@code CARTESIAN})
486      * @return the covariance matrix expressed in the target commonly used local orbital frame in Cartesian elements
487      */
488     private static StateCovariance changeFrameAndCreate(final Orbit orbit, final AbsoluteDate date,
489                                                         final LOF lofIn, final LOF lofOut,
490                                                         final RealMatrix inputCartesianCov) {
491 
492         // Builds the matrix to perform covariance transformation
493         final RealMatrix jacobianFromLofInToLofOut =
494                 getJacobian(LOF.transformFromLOFInToLOFOut(lofIn, lofOut, date, orbit.getPVCoordinates()));
495 
496         // Get the Cartesian covariance matrix converted to frameOut
497         final RealMatrix cartesianCovarianceOut =
498                 jacobianFromLofInToLofOut.multiply(inputCartesianCov.multiplyTransposed(jacobianFromLofInToLofOut));
499 
500         // Output converted covariance
501         return new StateCovariance(cartesianCovarianceOut, date, lofOut);
502 
503     }
504 
505     /**
506      * Convert the covariance matrix from a {@link Frame frame} to a {@link LOF local orbital frame}.
507      * <p>
508      * Changing the covariance frame is a linear process, this method does not introduce approximation unless a change
509      * in covariance orbit type is required.
510      * <p>
511      * The transformation is based on Equation (18) to (20) of "Covariance Transformations for Satellite Flight Dynamics
512      * Operations" by David A. Vallado.
513      * <p>
514      * As the frame transformation must be performed with the covariance expressed in Cartesian elements, both the orbit
515      * and position angle types of the input covariance must be provided.
516      * <p>
517      * <b>The output covariance matrix will provided in Cartesian orbit type and not converted back to
518      * its original expression (if input different from Cartesian elements).</b>
519      *
520      * @param orbit orbit to which the covariance matrix should correspond
521      * @param date covariance epoch
522      * @param frameIn the frame in which the input covariance matrix is expressed. In case the frame is <b>not</b>
523      * pseudo-inertial, the input covariance matrix is expected to be expressed in <b>Cartesian elements</b>.
524      * @param lofOut the target local orbital frame
525      * @param inputCov input covariance
526      * @param covOrbitType orbit type of the covariance matrix (used if frameIn is pseudo-inertial)
527      * @param covAngleType position angle type of the covariance matrix (used if frameIn is pseudo-inertial) (not used
528      * if covOrbitType equals {@code CARTESIAN})
529      * @return the covariance matrix expressed in the target local orbital frame in Cartesian elements
530      * @throws OrekitException if input frame is <b>not</b> pseudo-inertial <b>and</b> the input covariance is
531      * <b>not</b> expressed in Cartesian elements.
532      */
533     private static StateCovariance changeFrameAndCreate(final Orbit orbit, final AbsoluteDate date,
534                                                         final Frame frameIn, final LOF lofOut,
535                                                         final RealMatrix inputCov,
536                                                         final OrbitType covOrbitType,
537                                                         final PositionAngleType covAngleType) {
538 
539         // Input frame is inertial
540         if (frameIn.isPseudoInertial()) {
541 
542             // Convert input matrix to Cartesian parameters in input frame
543             final RealMatrix cartesianCovarianceIn =
544                     changeTypeAndCreate(orbit, date, frameIn, covOrbitType, covAngleType,
545                                         OrbitType.CARTESIAN, PositionAngleType.MEAN,
546                                         inputCov).getMatrix();
547 
548             // Builds the matrix to perform covariance transformation
549             final RealMatrix jacobianFromFrameInToLofOut =
550                             getJacobian(lofOut.transformFromInertial(date, orbit.getPVCoordinates(frameIn)));
551 
552             // Get the Cartesian covariance matrix converted to frameOut
553             final RealMatrix cartesianCovarianceOut =
554                     jacobianFromFrameInToLofOut.multiply(cartesianCovarianceIn.multiplyTransposed(jacobianFromFrameInToLofOut));
555 
556             // Return converted covariance matrix expressed in Cartesian elements
557             return new StateCovariance(cartesianCovarianceOut, date, lofOut);
558 
559         }
560 
561         // Input frame is not inertial so the covariance matrix is expected to be in Cartesian elements
562         else {
563 
564             // Method checkInputConsistency already verify that input covariance is defined in CARTESIAN type
565             final Frame orbitInertialFrame = orbit.getFrame();
566 
567             // Compute rotation matrix from frameIn to orbit inertial frame
568             final RealMatrix cartesianCovarianceInOrbitFrame =
569                    changeFrameAndCreate(orbit, date, frameIn, orbitInertialFrame, inputCov,
570                                          OrbitType.CARTESIAN, PositionAngleType.MEAN).getMatrix();
571 
572             // Convert from orbit inertial frame to lofOut
573             return changeFrameAndCreate(orbit, date, orbitInertialFrame, lofOut, cartesianCovarianceInOrbitFrame,
574                                         OrbitType.CARTESIAN, PositionAngleType.MEAN);
575 
576         }
577 
578     }
579 
580     /**
581      * Convert the covariance matrix from a {@link LOF  local orbital frame} to a {@link Frame frame}.
582      * <p>
583      * Changing the covariance frame is a linear process, this method does not introduce approximation.
584      * <p>
585      * The transformation is based on Equation (18) to (20) of "Covariance Transformations for Satellite Flight Dynamics
586      * Operations" by David A. Vallado.
587      * <p>
588      * The <u>input</u> covariance matrix is necessarily provided in <b>Cartesian orbit type</b>.
589      * <p>
590      * The <u>output</u> covariance matrix will be expressed in <b>Cartesian elements</b>.
591      *
592      * @param orbit orbit to which the covariance matrix should correspond
593      * @param date covariance epoch
594      * @param lofIn the local orbital frame in which the input covariance matrix is expressed
595      * @param frameOut the target frame
596      * @param inputCartesianCov input covariance ({@code CARTESIAN})
597      * @return the covariance matrix expressed in the target frame in Cartesian elements
598      */
599     private static StateCovariance changeFrameAndCreate(final Orbit orbit, final AbsoluteDate date,
600                                                         final LOF lofIn, final Frame frameOut,
601                                                         final RealMatrix inputCartesianCov) {
602 
603         // Output frame is pseudo-inertial
604         if (frameOut.isPseudoInertial()) {
605 
606             // Builds the matrix to perform covariance transformation
607             final RealMatrix jacobianFromLofInToFrameOut =
608                     getJacobian(lofIn.transformFromInertial(date, orbit.getPVCoordinates(frameOut)).getInverse());
609 
610             // Transform covariance
611             final RealMatrix transformedCovariance =
612                     jacobianFromLofInToFrameOut.multiply(inputCartesianCov.multiplyTransposed(jacobianFromLofInToFrameOut));
613 
614             // Get the Cartesian covariance matrix converted to frameOut
615             return new StateCovariance(transformedCovariance, date, frameOut, OrbitType.CARTESIAN,
616                                        DEFAULT_POSITION_ANGLE);
617 
618         }
619 
620         // Output frame is not pseudo-inertial
621         else {
622 
623             // Builds the matrix to perform covariance transformation
624             final RealMatrix jacobianFromLofInToOrbitFrame =
625                     getJacobian(lofIn.transformFromInertial(date, orbit.getPVCoordinates()).getInverse());
626 
627             // Get the Cartesian covariance matrix converted to orbit inertial frame
628             final RealMatrix cartesianCovarianceInOrbitFrame = jacobianFromLofInToOrbitFrame.multiply(
629                     inputCartesianCov.multiplyTransposed(jacobianFromLofInToOrbitFrame));
630 
631             // Get the Cartesian covariance matrix converted to frameOut
632             return changeFrameAndCreate(orbit, date, orbit.getFrame(), frameOut, cartesianCovarianceInOrbitFrame,
633                                         OrbitType.CARTESIAN, PositionAngleType.MEAN);
634         }
635 
636     }
637 
638     /**
639      * Get the covariance matrix in another frame.
640      * <p>
641      * The transformation is based on Equation (18) of "Covariance Transformations for Satellite Flight Dynamics
642      * Operations" by David A. Vallado.
643      * <p>
644      * Changing the covariance frame is a linear process, this method does not introduce approximation unless a change
645      * in covariance orbit type is required.
646      * <p>
647      * As the frame transformation must be performed with the covariance expressed in Cartesian elements, both the orbit
648      * and position angle types of the input covariance must be provided.
649      * <p>
650      * In case the <u>input</u> frame is <b>not</b> pseudo-inertial, the <u>input</u> covariance matrix is necessarily
651      * expressed in <b>Cartesian elements</b>.
652      * <p>
653      * In case the <u>output</u> frame is <b>not</b> pseudo-inertial, the <u>output</u> covariance matrix will be
654      * expressed in <b>Cartesian elements</b>.
655      *
656      * @param orbit orbit to which the covariance matrix should correspond
657      * @param date covariance epoch
658      * @param frameIn the frame in which the input covariance matrix is expressed
659      * @param frameOut the target frame
660      * @param inputCov input covariance
661      * @param covOrbitType orbit type of the covariance matrix (used if frameIn is pseudo-inertial)
662      * @param covAngleType position angle type of the covariance matrix (used if frameIn is pseudo-inertial) (<b>not</b>
663      * used if covOrbitType equals {@code CARTESIAN})
664      * @return the covariance matrix expressed in the target frame
665      * @throws OrekitException if input frame is <b>not</b> pseudo-inertial <b>and</b> the input covariance is
666      * <b>not</b> expressed in Cartesian elements.
667      */
668     private static StateCovariance changeFrameAndCreate(final Orbit orbit, final AbsoluteDate date,
669                                                         final Frame frameIn, final Frame frameOut,
670                                                         final RealMatrix inputCov,
671                                                         final OrbitType covOrbitType,
672                                                         final PositionAngleType covAngleType) {
673 
674         // Input frame pseudo-inertial
675         if (frameIn.isPseudoInertial()) {
676 
677             // Convert input matrix to Cartesian parameters in input frame
678             final RealMatrix cartesianCovarianceIn =
679                     changeTypeAndCreate(orbit, date, frameIn, covOrbitType, covAngleType,
680                                         OrbitType.CARTESIAN, PositionAngleType.MEAN,
681                                         inputCov).getMatrix();
682 
683             // Get the Cartesian covariance matrix converted to frameOut
684             final RealMatrix cartesianCovarianceOut = CartesianCovarianceUtils.changeReferenceFrame(frameIn,
685                     cartesianCovarianceIn, date, frameOut);
686 
687             // Output frame is pseudo-inertial
688             if (frameOut.isPseudoInertial()) {
689 
690                 // Convert output Cartesian matrix to initial orbit type and position angle
691                 return changeTypeAndCreate(orbit, date, frameOut, OrbitType.CARTESIAN, PositionAngleType.MEAN,
692                                            covOrbitType, covAngleType, cartesianCovarianceOut);
693 
694             }
695 
696             // Output frame is not pseudo-inertial
697             else {
698 
699                 // Output Cartesian matrix
700                 return new StateCovariance(cartesianCovarianceOut, date, frameOut, OrbitType.CARTESIAN,
701                                            DEFAULT_POSITION_ANGLE);
702 
703             }
704         }
705 
706         // Input frame is not pseudo-inertial
707         else {
708 
709             // Method checkInputConsistency already verify that input covariance is defined in CARTESIAN type
710 
711             // Convert covariance matrix to frameOut
712             final RealMatrix covInFrameOut = CartesianCovarianceUtils.changeReferenceFrame(frameIn,
713                     inputCov, date, frameOut);
714 
715             // Output the Cartesian covariance matrix converted to frameOut
716             return new StateCovariance(covInFrameOut, date, frameOut, OrbitType.CARTESIAN, DEFAULT_POSITION_ANGLE);
717 
718         }
719 
720     }
721 
722     /**
723      * Builds the matrix to perform covariance frame transformation.
724      * @param transform input transformation
725      * @return the matrix to perform the covariance frame transformation
726      */
727     private static RealMatrix getJacobian(final KinematicTransform transform) {
728         return MatrixUtils.createRealMatrix(transform.getPVJacobian());
729     }
730 
731     /**
732      * Get the state transition matrix considering Keplerian contribution only.
733      *
734      * @param initialOrbit orbit to which the initial covariance matrix should correspond
735      * @param dt time difference between the two orbits
736      * @return the state transition matrix used to shift the covariance matrix
737      */
738     public static RealMatrix getStm(final Orbit initialOrbit, final double dt) {
739 
740         // initialize the STM
741         final RealMatrix stm = MatrixUtils.createRealIdentityMatrix(STATE_DIMENSION);
742 
743         // State transition matrix using Keplerian contribution only
744         final double contribution = initialOrbit.getMeanAnomalyDotWrtA() * dt;
745         stm.setEntry(5, 0, contribution);
746 
747         // Return
748         return stm;
749 
750     }
751 
752 }