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.orbits;
18  
19  import org.hipparchus.geometry.euclidean.threed.Vector3D;
20  import org.hipparchus.linear.DecompositionSolver;
21  import org.hipparchus.linear.MatrixUtils;
22  import org.hipparchus.linear.QRDecomposition;
23  import org.hipparchus.linear.RealMatrix;
24  import org.hipparchus.util.FastMath;
25  import org.hipparchus.util.MathArrays;
26  import org.orekit.errors.OrekitIllegalArgumentException;
27  import org.orekit.errors.OrekitInternalError;
28  import org.orekit.errors.OrekitMessages;
29  import org.orekit.frames.Frame;
30  import org.orekit.frames.StaticTransform;
31  import org.orekit.frames.Transform;
32  import org.orekit.time.AbsoluteDate;
33  import org.orekit.time.TimeOffset;
34  import org.orekit.time.TimeShiftable;
35  import org.orekit.time.TimeStamped;
36  import org.orekit.utils.PVCoordinates;
37  import org.orekit.utils.PVCoordinatesProvider;
38  import org.orekit.utils.TimeStampedPVCoordinates;
39  
40  /**
41   * This class handles orbital parameters.
42  
43   * <p>
44   * For user convenience, both the Cartesian and the equinoctial elements
45   * are provided by this class, regardless of the canonical representation
46   * implemented in the derived class (which may be classical Keplerian
47   * elements for example).
48   * </p>
49   * <p>
50   * The parameters are defined in a frame specified by the user. It is important
51   * to make sure this frame is consistent: it probably is inertial and centered
52   * on the central body. This information is used for example by some
53   * force models.
54   * </p>
55   * <p>
56   * Instance of this class are guaranteed to be immutable.
57   * </p>
58   * @author Luc Maisonobe
59   * @author Guylaine Prat
60   * @author Fabien Maussion
61   * @author V&eacute;ronique Pommier-Maurussane
62   */
63  public abstract class Orbit
64      implements TimeStamped, TimeShiftable<Orbit>, PVCoordinatesProvider {
65  
66      /** Absolute tolerance when checking if the rate of the position angle is Keplerian or not. */
67      protected static final double TOLERANCE_POSITION_ANGLE_RATE = 1e-15;
68  
69      /** Frame in which are defined the orbital parameters. */
70      private final Frame frame;
71  
72      /** Date of the orbital parameters. */
73      private final AbsoluteDate date;
74  
75      /** Value of mu used to compute position and velocity (m³/s²). */
76      private final double mu;
77  
78      /** Computed position.
79       * @since 12.0
80       */
81      private Vector3D position;
82  
83      /** Computed PVCoordinates. */
84      private TimeStampedPVCoordinates pvCoordinates;
85  
86      /** Jacobian of the orbital parameters with mean angle with respect to the Cartesian coordinates. */
87      private double[][] jacobianMeanWrtCartesian;
88  
89      /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with mean angle. */
90      private double[][] jacobianWrtParametersMean;
91  
92      /** Jacobian of the orbital parameters with eccentric angle with respect to the Cartesian coordinates. */
93      private double[][] jacobianEccentricWrtCartesian;
94  
95      /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with eccentric angle. */
96      private double[][] jacobianWrtParametersEccentric;
97  
98      /** Jacobian of the orbital parameters with true angle with respect to the Cartesian coordinates. */
99      private double[][] jacobianTrueWrtCartesian;
100 
101     /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with true angle. */
102     private double[][] jacobianWrtParametersTrue;
103 
104     /** Default constructor.
105      * Build a new instance with arbitrary default elements.
106      * @param frame the frame in which the parameters are defined
107      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
108      * @param date date of the orbital parameters
109      * @param mu central attraction coefficient (m^3/s^2)
110      * @exception IllegalArgumentException if frame is not a {@link
111      * Frame#isPseudoInertial pseudo-inertial frame}
112      */
113     protected Orbit(final Frame frame, final AbsoluteDate date, final double mu)
114         throws IllegalArgumentException {
115         ensurePseudoInertialFrame(frame);
116         this.date                      = date;
117         this.mu                        = mu;
118         this.pvCoordinates             = null;
119         this.frame                     = frame;
120         jacobianMeanWrtCartesian       = null;
121         jacobianWrtParametersMean      = null;
122         jacobianEccentricWrtCartesian  = null;
123         jacobianWrtParametersEccentric = null;
124         jacobianTrueWrtCartesian       = null;
125         jacobianWrtParametersTrue      = null;
126     }
127 
128     /** Set the orbit from Cartesian parameters.
129      *
130      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
131      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
132      * use {@code mu} and the position to compute the acceleration, including
133      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
134      *
135      * @param pvCoordinates the position and velocity in the inertial frame
136      * @param frame the frame in which the {@link TimeStampedPVCoordinates} are defined
137      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
138      * @param mu central attraction coefficient (m^3/s^2)
139      * @exception IllegalArgumentException if frame is not a {@link
140      * Frame#isPseudoInertial pseudo-inertial frame}
141      */
142     protected Orbit(final TimeStampedPVCoordinates pvCoordinates, final Frame frame, final double mu)
143         throws IllegalArgumentException {
144         ensurePseudoInertialFrame(frame);
145         this.date = pvCoordinates.getDate();
146         this.mu = mu;
147         if (pvCoordinates.getAcceleration().getNormSq() == 0) {
148             // the acceleration was not provided,
149             // compute it from Newtonian attraction
150             final double r2 = pvCoordinates.getPosition().getNormSq();
151             final double r3 = r2 * FastMath.sqrt(r2);
152             this.pvCoordinates = new TimeStampedPVCoordinates(pvCoordinates.getDate(),
153                                                               pvCoordinates.getPosition(),
154                                                               pvCoordinates.getVelocity(),
155                                                               new Vector3D(-mu / r3, pvCoordinates.getPosition()));
156         } else {
157             this.pvCoordinates = pvCoordinates;
158         }
159         this.frame = frame;
160     }
161 
162     /** Check if Cartesian coordinates include non-Keplerian acceleration.
163      * @param pva Cartesian coordinates
164      * @param mu central attraction coefficient
165      * @return true if Cartesian coordinates include non-Keplerian acceleration
166      */
167     protected static boolean hasNonKeplerianAcceleration(final PVCoordinates pva, final double mu) {
168 
169         final Vector3D a = pva.getAcceleration();
170         if (a == null) {
171             return false;
172         }
173 
174         final Vector3D p = pva.getPosition();
175         final double r2 = p.getNormSq();
176 
177         // Check if acceleration is relatively close to 0 compared to the Keplerian acceleration
178         final double tolerance = mu * 1e-9;
179         final Vector3D aTimesR2 = a.scalarMultiply(r2);
180         if (aTimesR2.getNorm() < tolerance) {
181             return false;
182         }
183 
184         if ((aTimesR2.add(p.normalize().scalarMultiply(mu))).getNorm() > tolerance) {
185             // we have a relevant acceleration, we can compute derivatives
186             return true;
187         } else {
188             // the provided acceleration is either too small to be reliable (probably even 0), or NaN
189             return false;
190         }
191     }
192 
193     /** Returns true if and only if the orbit is elliptical i.e. has a non-negative semi-major axis.
194      * @return true if getA() is strictly greater than 0
195      * @since 12.0
196      */
197     public boolean isElliptical() {
198         return getA() > 0.;
199     }
200 
201     /** Get the orbit type.
202      * @return orbit type
203      */
204     public abstract OrbitType getType();
205 
206     /** Ensure the defining frame is a pseudo-inertial frame.
207      * @param frame frame to check
208      * @exception IllegalArgumentException if frame is not a {@link
209      * Frame#isPseudoInertial pseudo-inertial frame}
210      */
211     private static void ensurePseudoInertialFrame(final Frame frame)
212         throws IllegalArgumentException {
213         if (!frame.isPseudoInertial()) {
214             throw new OrekitIllegalArgumentException(OrekitMessages.NON_PSEUDO_INERTIAL_FRAME,
215                                                      frame.getName());
216         }
217     }
218 
219     /** Get the frame in which the orbital parameters are defined.
220      * @return frame in which the orbital parameters are defined
221      */
222     public Frame getFrame() {
223         return frame;
224     }
225 
226     /** Get the semi-major axis.
227      * <p>Note that the semi-major axis is considered negative for hyperbolic orbits.</p>
228      * @return semi-major axis (m)
229      */
230     public abstract double getA();
231 
232     /** Get the semi-major axis derivative.
233      * <p>Note that the semi-major axis is considered negative for hyperbolic orbits.</p>
234      * <p>
235      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
236      * </p>
237      * @return semi-major axis  derivative (m/s)
238      * @since 9.0
239      */
240     public abstract double getADot();
241 
242     /** Get the first component of the equinoctial eccentricity vector.
243      * @return first component of the equinoctial eccentricity vector
244      */
245     public abstract double getEquinoctialEx();
246 
247     /** Get the first component of the equinoctial eccentricity vector derivative.
248      * <p>
249      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
250      * </p>
251      * @return first component of the equinoctial eccentricity vector derivative
252      * @since 9.0
253      */
254     public abstract double getEquinoctialExDot();
255 
256     /** Get the second component of the equinoctial eccentricity vector.
257      * @return second component of the equinoctial eccentricity vector
258      */
259     public abstract double getEquinoctialEy();
260 
261     /** Get the second component of the equinoctial eccentricity vector derivative.
262      * <p>
263      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
264      * </p>
265      * @return second component of the equinoctial eccentricity vector derivative
266      * @since 9.0
267      */
268     public abstract double getEquinoctialEyDot();
269 
270     /** Get the first component of the inclination vector.
271      * @return first component of the inclination vector
272      */
273     public abstract double getHx();
274 
275     /** Get the first component of the inclination vector derivative.
276      * <p>
277      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
278      * </p>
279      * @return first component of the inclination vector derivative
280 
281      * @since 9.0
282      */
283     public abstract double getHxDot();
284 
285     /** Get the second component of the inclination vector.
286      * @return second component of the inclination vector
287      */
288     public abstract double getHy();
289 
290     /** Get the second component of the inclination vector derivative.
291      * <p>
292      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
293      * </p>
294      * @return second component of the inclination vector derivative
295      * @since 9.0
296      */
297     public abstract double getHyDot();
298 
299     /** Get the eccentric longitude argument.
300      * @return E + ω + Ω eccentric longitude argument (rad)
301      */
302     public abstract double getLE();
303 
304     /** Get the eccentric longitude argument derivative.
305      * <p>
306      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
307      * </p>
308      * @return d(E + ω + Ω)/dt eccentric longitude argument derivative (rad/s)
309      * @since 9.0
310      */
311     public abstract double getLEDot();
312 
313     /** Get the true longitude argument.
314      * @return v + ω + Ω true longitude argument (rad)
315      */
316     public abstract double getLv();
317 
318     /** Get the true longitude argument derivative.
319      * <p>
320      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
321      * </p>
322      * @return d(v + ω + Ω)/dt true longitude argument derivative (rad/s)
323      * @since 9.0
324      */
325     public abstract double getLvDot();
326 
327     /** Get the mean longitude argument.
328      * @return M + ω + Ω mean longitude argument (rad)
329      */
330     public abstract double getLM();
331 
332     /** Get the mean longitude argument derivative.
333      * <p>
334      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
335      * </p>
336      * @return d(M + ω + Ω)/dt mean longitude argument derivative (rad/s)
337      * @since 9.0
338      */
339     public abstract double getLMDot();
340 
341     // Additional orbital elements
342 
343     /** Get the eccentricity.
344      * @return eccentricity
345      */
346     public abstract double getE();
347 
348     /** Get the eccentricity derivative.
349      * <p>
350      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
351      * </p>
352      * @return eccentricity derivative
353      * @since 9.0
354      */
355     public abstract double getEDot();
356 
357     /** Get the inclination.
358      * @return inclination (rad)
359      */
360     public abstract double getI();
361 
362     /** Get the inclination derivative.
363      * <p>
364      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
365      * </p>
366      * @return inclination derivative (rad/s)
367      * @since 9.0
368      */
369     public abstract double getIDot();
370 
371     /** Check if orbit includes non-Keplerian rates.
372      * @return true if orbit includes non-Keplerian derivatives
373      * @see #getADot()
374      * @see #getEquinoctialExDot()
375      * @see #getEquinoctialEyDot()
376      * @see #getHxDot()
377      * @see #getHyDot()
378      * @see #getLEDot()
379      * @see #getLvDot()
380      * @see #getLMDot()
381      * @see #getEDot()
382      * @see #getIDot()
383      * @since 13.0
384      */
385     public boolean hasNonKeplerianAcceleration() {
386         return hasNonKeplerianAcceleration(getPVCoordinates(), getMu());
387     }
388 
389     /** Get the central acceleration constant.
390      * @return central acceleration constant
391      */
392     public double getMu() {
393         return mu;
394     }
395 
396     /** Get the Keplerian period.
397      * <p>The Keplerian period is computed directly from semi major axis
398      * and central acceleration constant.</p>
399      * @return Keplerian period in seconds, or positive infinity for hyperbolic orbits
400      */
401     public double getKeplerianPeriod() {
402         final double a = getA();
403         return isElliptical() ? 2.0 * FastMath.PI * a * FastMath.sqrt(a / mu) : Double.POSITIVE_INFINITY;
404     }
405 
406     /** Get the Keplerian mean motion.
407      * <p>The Keplerian mean motion is computed directly from semi major axis
408      * and central acceleration constant.</p>
409      * @return Keplerian mean motion in radians per second
410      */
411     public double getKeplerianMeanMotion() {
412         final double absA = FastMath.abs(getA());
413         return FastMath.sqrt(mu / absA) / absA;
414     }
415 
416     /** Get the derivative of the mean anomaly with respect to the semi major axis.
417      * @return derivative of the mean anomaly with respect to the semi major axis
418      */
419     public double getMeanAnomalyDotWrtA() {
420         return -1.5 * getKeplerianMeanMotion() / getA();
421     }
422 
423     /** Get the date of orbital parameters.
424      * @return date of the orbital parameters
425      */
426     public AbsoluteDate getDate() {
427         return date;
428     }
429 
430     /** Get the {@link TimeStampedPVCoordinates} in a specified frame.
431      * @param outputFrame frame in which the position/velocity coordinates shall be computed
432      * @return pvCoordinates in the specified output frame
433           * @see #getPVCoordinates()
434      */
435     public TimeStampedPVCoordinates getPVCoordinates(final Frame outputFrame) {
436         if (pvCoordinates == null) {
437             pvCoordinates = initPVCoordinates();
438         }
439 
440         // If output frame requested is the same as definition frame,
441         // PV coordinates are returned directly
442         if (outputFrame == frame) {
443             return pvCoordinates;
444         }
445 
446         // Else, PV coordinates are transformed to output frame
447         final Transform t = frame.getTransformTo(outputFrame, date);
448         return t.transformPVCoordinates(pvCoordinates);
449     }
450 
451     /** {@inheritDoc} */
452     public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate otherDate, final Frame otherFrame) {
453         return shiftedBy(otherDate.durationFrom(getDate())).getPVCoordinates(otherFrame);
454     }
455 
456     /** {@inheritDoc} */
457     @Override
458     public Vector3D getPosition(final AbsoluteDate otherDate, final Frame otherFrame) {
459         return shiftedBy(otherDate.durationFrom(getDate())).getPosition(otherFrame);
460     }
461 
462     /** Get the position in a specified frame.
463      * @param outputFrame frame in which the position coordinates shall be computed
464      * @return position in the specified output frame
465      * @see #getPosition()
466      * @since 12.0
467      */
468     public Vector3D getPosition(final Frame outputFrame) {
469         if (position == null) {
470             position = initPosition();
471         }
472 
473         // If output frame requested is the same as definition frame,
474         // Position vector is returned directly
475         if (outputFrame == frame) {
476             return position;
477         }
478 
479         // Else, position vector is transformed to output frame
480         final StaticTransform t = frame.getStaticTransformTo(outputFrame, date);
481         return t.transformPosition(position);
482 
483     }
484 
485     /** Get the position in definition frame.
486      * @return position in the definition frame
487      * @see #getPVCoordinates()
488      * @since 12.0
489      */
490     public Vector3D getPosition() {
491         if (position == null) {
492             position = initPosition();
493         }
494         return position;
495     }
496 
497     /** Get the {@link TimeStampedPVCoordinates} in definition frame.
498      * @return pvCoordinates in the definition frame
499      * @see #getPVCoordinates(Frame)
500      */
501     public TimeStampedPVCoordinates getPVCoordinates() {
502         if (pvCoordinates == null) {
503             pvCoordinates = initPVCoordinates();
504             position      = pvCoordinates.getPosition();
505         }
506         return pvCoordinates;
507     }
508 
509     /** Compute the position coordinates from the canonical parameters.
510      * @return computed position coordinates
511      * @since 12.0
512      */
513     protected abstract Vector3D initPosition();
514 
515     /** Compute the position/velocity coordinates from the canonical parameters.
516      * @return computed position/velocity coordinates
517      */
518     protected abstract TimeStampedPVCoordinates initPVCoordinates();
519 
520     /**
521      * Create a new object representing the same physical orbital state, but attached to a different reference frame.
522      * If the new frame is not inertial, an exception will be thrown.
523      *
524      * @param inertialFrame reference frame of output orbit
525      * @return orbit with different frame
526      * @since 13.0
527      */
528     public abstract Orbit inFrame(Frame inertialFrame);
529 
530     /** Get a time-shifted orbit.
531      * <p>
532      * The orbit can be slightly shifted to close dates. The shifting model is a
533      * Keplerian one if no derivatives are available in the orbit, or Keplerian
534      * plus quadratic effect of the non-Keplerian acceleration if derivatives are
535      * available. Shifting is <em>not</em> intended as a replacement for proper
536      * orbit propagation but should be sufficient for small time shifts or coarse
537      * accuracy.
538      * </p>
539      * @param dt time shift in seconds
540      * @return a new orbit, shifted with respect to the instance (which is immutable)
541      */
542     @Override
543     public abstract Orbit shiftedBy(double dt);
544 
545     /** Get a time-shifted orbit.
546      * <p>
547      * The orbit can be slightly shifted to close dates. The shifting model is a
548      * Keplerian one if no derivatives are available in the orbit, or Keplerian
549      * plus quadratic effect of the non-Keplerian acceleration if derivatives are
550      * available. Shifting is <em>not</em> intended as a replacement for proper
551      * orbit propagation but should be sufficient for small time shifts or coarse
552      * accuracy.
553      * </p>
554      * @param dt time shift
555      * @return a new orbit, shifted with respect to the instance (which is immutable)
556      */
557     @Override
558     public abstract Orbit shiftedBy(TimeOffset dt);
559 
560     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
561      * <p>
562      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
563      * respect to Cartesian coordinate j. This means each row corresponds to one orbital parameter
564      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
565      * </p>
566      * @param type type of the position angle to use
567      * @param jacobian placeholder 6x6 (or larger) matrix to be filled with the Jacobian, if matrix
568      * is larger than 6x6, only the 6x6 upper left corner will be modified
569      */
570     public void getJacobianWrtCartesian(final PositionAngleType type, final double[][] jacobian) {
571 
572         final double[][] cachedJacobian;
573         synchronized (this) {
574             switch (type) {
575                 case MEAN :
576                     if (jacobianMeanWrtCartesian == null) {
577                         // first call, we need to compute the Jacobian and cache it
578                         jacobianMeanWrtCartesian = computeJacobianMeanWrtCartesian();
579                     }
580                     cachedJacobian = jacobianMeanWrtCartesian;
581                     break;
582                 case ECCENTRIC :
583                     if (jacobianEccentricWrtCartesian == null) {
584                         // first call, we need to compute the Jacobian and cache it
585                         jacobianEccentricWrtCartesian = computeJacobianEccentricWrtCartesian();
586                     }
587                     cachedJacobian = jacobianEccentricWrtCartesian;
588                     break;
589                 case TRUE :
590                     if (jacobianTrueWrtCartesian == null) {
591                         // first call, we need to compute the Jacobian and cache it
592                         jacobianTrueWrtCartesian = computeJacobianTrueWrtCartesian();
593                     }
594                     cachedJacobian = jacobianTrueWrtCartesian;
595                     break;
596                 default :
597                     throw new OrekitInternalError(null);
598             }
599         }
600 
601         // fill the user provided array
602         for (int i = 0; i < cachedJacobian.length; ++i) {
603             System.arraycopy(cachedJacobian[i], 0, jacobian[i], 0, cachedJacobian[i].length);
604         }
605 
606     }
607 
608     /** Compute the Jacobian of the Cartesian parameters with respect to the orbital parameters.
609      * <p>
610      * Element {@code jacobian[i][j]} is the derivative of Cartesian coordinate i of the orbit with
611      * respect to orbital parameter j. This means each row corresponds to one Cartesian coordinate
612      * x, y, z, xdot, ydot, zdot whereas columns 0 to 5 correspond to the orbital parameters.
613      * </p>
614      * @param type type of the position angle to use
615      * @param jacobian placeholder 6x6 (or larger) matrix to be filled with the Jacobian, if matrix
616      * is larger than 6x6, only the 6x6 upper left corner will be modified
617      */
618     public void getJacobianWrtParameters(final PositionAngleType type, final double[][] jacobian) {
619 
620         final double[][] cachedJacobian;
621         synchronized (this) {
622             switch (type) {
623                 case MEAN :
624                     if (jacobianWrtParametersMean == null) {
625                         // first call, we need to compute the Jacobian and cache it
626                         jacobianWrtParametersMean = createInverseJacobian(type);
627                     }
628                     cachedJacobian = jacobianWrtParametersMean;
629                     break;
630                 case ECCENTRIC :
631                     if (jacobianWrtParametersEccentric == null) {
632                         // first call, we need to compute the Jacobian and cache it
633                         jacobianWrtParametersEccentric = createInverseJacobian(type);
634                     }
635                     cachedJacobian = jacobianWrtParametersEccentric;
636                     break;
637                 case TRUE :
638                     if (jacobianWrtParametersTrue == null) {
639                         // first call, we need to compute the Jacobian and cache it
640                         jacobianWrtParametersTrue = createInverseJacobian(type);
641                     }
642                     cachedJacobian = jacobianWrtParametersTrue;
643                     break;
644                 default :
645                     throw new OrekitInternalError(null);
646             }
647         }
648 
649         // fill the user-provided array
650         for (int i = 0; i < cachedJacobian.length; ++i) {
651             System.arraycopy(cachedJacobian[i], 0, jacobian[i], 0, cachedJacobian[i].length);
652         }
653 
654     }
655 
656     /** Create an inverse Jacobian.
657      * @param type type of the position angle to use
658      * @return inverse Jacobian
659      */
660     private double[][] createInverseJacobian(final PositionAngleType type) {
661 
662         // get the direct Jacobian
663         final double[][] directJacobian = new double[6][6];
664         getJacobianWrtCartesian(type, directJacobian);
665 
666         // invert the direct Jacobian
667         final RealMatrix matrix = MatrixUtils.createRealMatrix(directJacobian);
668         final DecompositionSolver solver = new QRDecomposition(matrix).getSolver();
669         return solver.getInverse().getData();
670 
671     }
672 
673     /** Compute the Jacobian of the orbital parameters with mean angle with respect to the Cartesian parameters.
674      * <p>
675      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
676      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
677      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
678      * </p>
679      * <p>
680      * The array returned by this method will not be modified.
681      * </p>
682      * @return 6x6 Jacobian matrix
683      * @see #computeJacobianEccentricWrtCartesian()
684      * @see #computeJacobianTrueWrtCartesian()
685      */
686     protected abstract double[][] computeJacobianMeanWrtCartesian();
687 
688     /** Compute the Jacobian of the orbital parameters with eccentric angle with respect to the Cartesian parameters.
689      * <p>
690      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
691      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
692      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
693      * </p>
694      * <p>
695      * The array returned by this method will not be modified.
696      * </p>
697      * @return 6x6 Jacobian matrix
698      * @see #computeJacobianMeanWrtCartesian()
699      * @see #computeJacobianTrueWrtCartesian()
700      */
701     protected abstract double[][] computeJacobianEccentricWrtCartesian();
702 
703     /** Compute the Jacobian of the orbital parameters with true angle with respect to the Cartesian parameters.
704      * <p>
705      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
706      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
707      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
708      * </p>
709      * <p>
710      * The array returned by this method will not be modified.
711      * </p>
712      * @return 6x6 Jacobian matrix
713      * @see #computeJacobianMeanWrtCartesian()
714      * @see #computeJacobianEccentricWrtCartesian()
715      */
716     protected abstract double[][] computeJacobianTrueWrtCartesian();
717 
718     /** Add the contribution of the Keplerian motion to parameters derivatives
719      * <p>
720      * This method is used by integration-based propagators to evaluate the part of Keplerian
721      * motion to evolution of the orbital state.
722      * </p>
723      * @param type type of the position angle in the state
724      * @param gm attraction coefficient to use
725      * @param pDot array containing orbital state derivatives to update (the Keplerian
726      * part must be <em>added</em> to the array components, as the array may already
727      * contain some non-zero elements corresponding to non-Keplerian parts)
728      */
729     public abstract void addKeplerContribution(PositionAngleType type, double gm, double[] pDot);
730 
731         /** Fill a Jacobian half row with a single vector.
732      * @param a coefficient of the vector
733      * @param v vector
734      * @param row Jacobian matrix row
735      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
736      */
737     protected static void fillHalfRow(final double a, final Vector3D v, final double[] row, final int j) {
738         row[j]     = a * v.getX();
739         row[j + 1] = a * v.getY();
740         row[j + 2] = a * v.getZ();
741     }
742 
743     /** Fill a Jacobian half row with a linear combination of vectors.
744      * @param a1 coefficient of the first vector
745      * @param v1 first vector
746      * @param a2 coefficient of the second vector
747      * @param v2 second vector
748      * @param row Jacobian matrix row
749      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
750      */
751     protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2,
752                                       final double[] row, final int j) {
753         row[j]     = MathArrays.linearCombination(a1, v1.getX(), a2, v2.getX());
754         row[j + 1] = MathArrays.linearCombination(a1, v1.getY(), a2, v2.getY());
755         row[j + 2] = MathArrays.linearCombination(a1, v1.getZ(), a2, v2.getZ());
756     }
757 
758     /** Fill a Jacobian half row with a linear combination of vectors.
759      * @param a1 coefficient of the first vector
760      * @param v1 first vector
761      * @param a2 coefficient of the second vector
762      * @param v2 second vector
763      * @param a3 coefficient of the third vector
764      * @param v3 third vector
765      * @param row Jacobian matrix row
766      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
767      */
768     protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2,
769                                       final double a3, final Vector3D v3,
770                                       final double[] row, final int j) {
771         row[j]     = MathArrays.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX());
772         row[j + 1] = MathArrays.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY());
773         row[j + 2] = MathArrays.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ());
774     }
775 
776     /** Fill a Jacobian half row with a linear combination of vectors.
777      * @param a1 coefficient of the first vector
778      * @param v1 first vector
779      * @param a2 coefficient of the second vector
780      * @param v2 second vector
781      * @param a3 coefficient of the third vector
782      * @param v3 third vector
783      * @param a4 coefficient of the fourth vector
784      * @param v4 fourth vector
785      * @param row Jacobian matrix row
786      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
787      */
788     protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2,
789                                       final double a3, final Vector3D v3, final double a4, final Vector3D v4,
790                                       final double[] row, final int j) {
791         row[j]     = MathArrays.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX(), a4, v4.getX());
792         row[j + 1] = MathArrays.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY(), a4, v4.getY());
793         row[j + 2] = MathArrays.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ(), a4, v4.getZ());
794     }
795 
796     /** Fill a Jacobian half row with a linear combination of vectors.
797      * @param a1 coefficient of the first vector
798      * @param v1 first vector
799      * @param a2 coefficient of the second vector
800      * @param v2 second vector
801      * @param a3 coefficient of the third vector
802      * @param v3 third vector
803      * @param a4 coefficient of the fourth vector
804      * @param v4 fourth vector
805      * @param a5 coefficient of the fifth vector
806      * @param v5 fifth vector
807      * @param row Jacobian matrix row
808      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
809      */
810     protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2,
811                                       final double a3, final Vector3D v3, final double a4, final Vector3D v4,
812                                       final double a5, final Vector3D v5,
813                                       final double[] row, final int j) {
814         final double[] a = new double[] {
815             a1, a2, a3, a4, a5
816         };
817         row[j]     = MathArrays.linearCombination(a, new double[] {
818             v1.getX(), v2.getX(), v3.getX(), v4.getX(), v5.getX()
819         });
820         row[j + 1] = MathArrays.linearCombination(a, new double[] {
821             v1.getY(), v2.getY(), v3.getY(), v4.getY(), v5.getY()
822         });
823         row[j + 2] = MathArrays.linearCombination(a, new double[] {
824             v1.getZ(), v2.getZ(), v3.getZ(), v4.getZ(), v5.getZ()
825         });
826     }
827 
828     /** Fill a Jacobian half row with a linear combination of vectors.
829      * @param a1 coefficient of the first vector
830      * @param v1 first vector
831      * @param a2 coefficient of the second vector
832      * @param v2 second vector
833      * @param a3 coefficient of the third vector
834      * @param v3 third vector
835      * @param a4 coefficient of the fourth vector
836      * @param v4 fourth vector
837      * @param a5 coefficient of the fifth vector
838      * @param v5 fifth vector
839      * @param a6 coefficient of the sixth vector
840      * @param v6 sixth vector
841      * @param row Jacobian matrix row
842      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
843      */
844     protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2,
845                                       final double a3, final Vector3D v3, final double a4, final Vector3D v4,
846                                       final double a5, final Vector3D v5, final double a6, final Vector3D v6,
847                                       final double[] row, final int j) {
848         final double[] a = new double[] {
849             a1, a2, a3, a4, a5, a6
850         };
851         row[j]     = MathArrays.linearCombination(a, new double[] {
852             v1.getX(), v2.getX(), v3.getX(), v4.getX(), v5.getX(), v6.getX()
853         });
854         row[j + 1] = MathArrays.linearCombination(a, new double[] {
855             v1.getY(), v2.getY(), v3.getY(), v4.getY(), v5.getY(), v6.getY()
856         });
857         row[j + 2] = MathArrays.linearCombination(a, new double[] {
858             v1.getZ(), v2.getZ(), v3.getZ(), v4.getZ(), v5.getZ(), v6.getZ()
859         });
860     }
861 
862 }