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