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