FieldOrbit.java

  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;



  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.FieldStaticTransform;
  35. import org.orekit.frames.FieldTransform;
  36. import org.orekit.frames.Frame;
  37. import org.orekit.time.FieldAbsoluteDate;
  38. import org.orekit.time.FieldTimeShiftable;
  39. import org.orekit.time.FieldTimeStamped;
  40. import org.orekit.utils.FieldPVCoordinates;
  41. import org.orekit.utils.FieldPVCoordinatesProvider;
  42. import org.orekit.utils.TimeStampedFieldPVCoordinates;
  43. import org.orekit.utils.TimeStampedPVCoordinates;

  45. /**
  46.  * This class handles orbital parameters.

  48.  * <p>
  49.  * For user convenience, both the Cartesian and the equinoctial elements
  50.  * are provided by this class, regardless of the canonical representation
  51.  * implemented in the derived class (which may be classical Keplerian
  52.  * elements for example).
  53.  * </p>
  54.  * <p>
  55.  * The parameters are defined in a frame specified by the user. It is important
  56.  * to make sure this frame is consistent: it probably is inertial and centered
  57.  * on the central body. This information is used for example by some
  58.  * force models.
  59.  * </p>
  60.  * <p>
  61.  * Instance of this class are guaranteed to be immutable.
  62.  * </p>
  63.  * @author Luc Maisonobe
  64.  * @author Guylaine Prat
  65.  * @author Fabien Maussion
  66.  * @author V&eacute;ronique Pommier-Maurussane
  67.  * @since 9.0
  68.  * @see Orbit
  69.  * @param <T> type of the field elements
  70.  */
  71. public abstract class FieldOrbit<T extends CalculusFieldElement<T>>
  72.     implements FieldPVCoordinatesProvider<T>, FieldTimeStamped<T>, FieldTimeShiftable<FieldOrbit<T>, T> {

  74.     /** Absolute tolerance when checking if the rate of the position angle is Keplerian or not. */
  75.     protected static final double TOLERANCE_POSITION_ANGLE_RATE = 1e-15;

  77.     /** Frame in which are defined the orbital parameters. */
  78.     private final Frame frame;

  80.     /** Date of the orbital parameters. */
  81.     private final FieldAbsoluteDate<T> date;

  83.     /** Value of mu used to compute position and velocity (m³/s²). */
  84.     private final T mu;

  86.     /** Field-valued one.
  87.      * @since 12.0
  88.      * */
  89.     private final T one;

  91.     /** Field-valued zero.
  92.      * @since 12.0
  93.      * */
  94.     private final T zero;

  96.     /** Field used by this class.
  97.      * @since 12.0
  98.      */
  99.     private final Field<T> field;

  101.     /** Computed position.
  102.      * @since 12.0
  103.      */
  104.     private FieldVector3D<T> position;

  106.     /** Computed PVCoordinates. */
  107.     private TimeStampedFieldPVCoordinates<T> pvCoordinates;

  109.     /** Jacobian of the orbital parameters with mean angle with respect to the Cartesian coordinates. */
  110.     private T[][] jacobianMeanWrtCartesian;

  112.     /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with mean angle. */
  113.     private T[][] jacobianWrtParametersMean;

  115.     /** Jacobian of the orbital parameters with eccentric angle with respect to the Cartesian coordinates. */
  116.     private T[][] jacobianEccentricWrtCartesian;

  118.     /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with eccentric angle. */
  119.     private T[][] jacobianWrtParametersEccentric;

  121.     /** Jacobian of the orbital parameters with true angle with respect to the Cartesian coordinates. */
  122.     private T[][] jacobianTrueWrtCartesian;

  124.     /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with true angle. */
  125.     private T[][] jacobianWrtParametersTrue;

  127.     /** Default constructor.
  128.      * Build a new instance with arbitrary default elements.
  129.      * @param frame the frame in which the parameters are defined
  130.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  131.      * @param date date of the orbital parameters
  132.      * @param mu central attraction coefficient (m^3/s^2)
  133.      * @exception IllegalArgumentException if frame is not a {@link
  134.      * Frame#isPseudoInertial pseudo-inertial frame}
  135.      */
  136.     protected FieldOrbit(final Frame frame, final FieldAbsoluteDate<T> date, final T mu)
  137.         throws IllegalArgumentException {
  138.         ensurePseudoInertialFrame(frame);
  139.         this.field = mu.getField();
  140.         this.zero = this.field.getZero();
  141.         this.one = this.field.getOne();
  142.         this.date                      = date;
  143.         this.mu                        = mu;
  144.         this.pvCoordinates             = null;
  145.         this.frame                     = frame;
  146.         jacobianMeanWrtCartesian       = null;
  147.         jacobianWrtParametersMean      = null;
  148.         jacobianEccentricWrtCartesian  = null;
  149.         jacobianWrtParametersEccentric = null;
  150.         jacobianTrueWrtCartesian       = null;
  151.         jacobianWrtParametersTrue      = null;
  152.     }

  154.     /** Set the orbit from Cartesian parameters.
  155.      *
  156.      * <p> The acceleration provided in {@code FieldPVCoordinates} is accessible using
  157.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  158.      * use {@code mu} and the position to compute the acceleration, including
  159.      * {@link #shiftedBy(CalculusFieldElement)} and {@link #getPVCoordinates(FieldAbsoluteDate, Frame)}.
  160.      *
  161.      * @param fieldPVCoordinates the position and velocity in the inertial frame
  162.      * @param frame the frame in which the {@link TimeStampedPVCoordinates} are defined
  163.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  164.      * @param mu central attraction coefficient (m^3/s^2)
  165.      * @exception IllegalArgumentException if frame is not a {@link
  166.      * Frame#isPseudoInertial pseudo-inertial frame}
  167.      */
  168.     protected FieldOrbit(final TimeStampedFieldPVCoordinates<T> fieldPVCoordinates, final Frame frame, final T mu)
  169.         throws IllegalArgumentException {
  170.         ensurePseudoInertialFrame(frame);
  171.         this.field = mu.getField();
  172.         this.zero = this.field.getZero();
  173.         this.one = this.field.getOne();
  174.         this.date = fieldPVCoordinates.getDate();
  175.         this.mu = mu;
  176.         if (fieldPVCoordinates.getAcceleration().getNormSq().getReal() == 0.0) {
  177.             // the acceleration was not provided,
  178.             // compute it from Newtonian attraction
  179.             final T r2 = fieldPVCoordinates.getPosition().getNormSq();
  180.             final T r3 = r2.multiply(r2.sqrt());
  181.             this.pvCoordinates = new TimeStampedFieldPVCoordinates<>(fieldPVCoordinates.getDate(),
  182.                                                                      fieldPVCoordinates.getPosition(),
  183.                                                                      fieldPVCoordinates.getVelocity(),
  184.                                                                      new FieldVector3D<>(r3.reciprocal().multiply(mu.negate()), fieldPVCoordinates.getPosition()));
  185.         } else {
  186.             this.pvCoordinates = fieldPVCoordinates;
  187.         }
  188.         this.frame = frame;
  189.     }

  191.     /** Compute non-Keplerian part of the acceleration from first time derivatives.
  192.      * @return non-Keplerian part of the acceleration
  193.      * @since 13.1
  194.      */
  195.     protected FieldVector3D<T> nonKeplerianAcceleration() {

  197.         final T[][] dPdC = MathArrays.buildArray(getField(), 6, 6);
  198.         final PositionAngleType positionAngleType = PositionAngleType.MEAN;
  199.         getJacobianWrtCartesian(positionAngleType, dPdC);
  200.         final FieldMatrix<T> subMatrix = MatrixUtils.createFieldMatrix(dPdC);

  202.         final FieldDecompositionSolver<T> solver = getDecompositionSolver(subMatrix);

  204.         final T[] derivatives = dPdC[0].clone();
  205.         getType().mapOrbitToArray(this, positionAngleType, derivatives.clone(), derivatives);
  206.         derivatives[5] = derivatives[5].subtract(getKeplerianMeanMotion());

  208.         final FieldVector<T> solution = solver.solve(MatrixUtils.createFieldVector(derivatives));
  209.         // solution is vector-acceleration vector
  210.         return new FieldVector3D<>(solution.getEntry(3), solution.getEntry(4), solution.getEntry(5));

  212.     }

  214.     /** Check if Cartesian coordinates include non-Keplerian acceleration.
  215.      * @param pva Cartesian coordinates
  216.      * @param mu central attraction coefficient
  217.      * @param <T> type of the field elements
  218.      * @return true if Cartesian coordinates include non-Keplerian acceleration
  219.      */
  220.     protected static <T extends CalculusFieldElement<T>> boolean hasNonKeplerianAcceleration(final FieldPVCoordinates<T> pva, final T mu) {

  222.         if (mu.getField().getZero() instanceof Derivative<?>) {
  223.             return Orbit.hasNonKeplerianAcceleration(pva.toPVCoordinates(), mu.getReal()); // for performance

  225.         } else {
  226.             final FieldVector3D<T> a = pva.getAcceleration();
  227.             if (a == null) {
  228.                 return false;
  229.             }

  231.             final FieldVector3D<T> p = pva.getPosition();
  232.             final T r2 = p.getNormSq();

  234.             // Check if acceleration is relatively close to 0 compared to the Keplerian acceleration
  235.             final double tolerance = mu.getReal() * 1e-9;
  236.             final FieldVector3D<T> aTimesR2 = a.scalarMultiply(r2);
  237.             if (aTimesR2.getNorm().getReal() < tolerance) {
  238.                 return false;
  239.             }

  241.             if ((aTimesR2.add(p.normalize().scalarMultiply(mu))).getNorm().getReal() > tolerance) {
  242.                 // we have a relevant acceleration, we can compute derivatives
  243.                 return true;
  244.             } else {
  245.                 // the provided acceleration is either too small to be reliable (probably even 0), or NaN
  246.                 return false;
  247.             }
  248.         }

  250.     }

  252.     /** Returns true if and only if the orbit is elliptical i.e. has a non-negative semi-major axis.
  253.      * @return true if getA() is strictly greater than 0
  254.      * @since 12.0
  255.      */
  256.     public boolean isElliptical() {
  257.         return getA().getReal() > 0.;
  258.     }

  260.     /** Get the orbit type.
  261.      * @return orbit type
  262.      */
  263.     public abstract OrbitType getType();

  265.     /** Ensure the defining frame is a pseudo-inertial frame.
  266.      * @param frame frame to check
  267.      * @exception IllegalArgumentException if frame is not a {@link
  268.      * Frame#isPseudoInertial pseudo-inertial frame}
  269.      */
  270.     private static void ensurePseudoInertialFrame(final Frame frame)
  271.         throws IllegalArgumentException {
  272.         if (!frame.isPseudoInertial()) {
  273.             throw new OrekitIllegalArgumentException(OrekitMessages.NON_PSEUDO_INERTIAL_FRAME,
  274.                                                      frame.getName());
  275.         }
  276.     }

  278.     /** Get the frame in which the orbital parameters are defined.
  279.      * @return frame in which the orbital parameters are defined
  280.      */
  281.     public Frame getFrame() {
  282.         return frame;
  283.     }

  285.     /**Transforms the FieldOrbit instance into an Orbit instance.
  286.      * @return Orbit instance with same properties
  287.      */
  288.     public abstract Orbit toOrbit();

  290.     /** Get the semi-major axis.
  291.      * <p>Note that the semi-major axis is considered negative for hyperbolic orbits.</p>
  292.      * @return semi-major axis (m)
  293.      */
  294.     public abstract T getA();

  296.     /** Get the semi-major axis derivative.
  297.      * <p>Note that the semi-major axis is considered negative for hyperbolic orbits.</p>
  298.      * <p>
  299.      * If the orbit was created without derivatives, the value returned is null.
  300.      * </p>
  301.      * @return semi-major axis  derivative (m/s)
  302.      */
  303.     public abstract T getADot();

  305.     /** Get the first component of the equinoctial eccentricity vector.
  306.      * @return first component of the equinoctial eccentricity vector
  307.      */
  308.     public abstract T getEquinoctialEx();

  310.     /** Get the first component of the equinoctial eccentricity vector derivative.
  311.      * <p>
  312.      * If the orbit was created without derivatives, the value returned is null.
  313.      * </p>
  314.      * @return first component of the equinoctial eccentricity vector derivative
  315.      */
  316.     public abstract T getEquinoctialExDot();

  318.     /** Get the second component of the equinoctial eccentricity vector.
  319.      * @return second component of the equinoctial eccentricity vector
  320.      */
  321.     public abstract T getEquinoctialEy();

  323.     /** Get the second component of the equinoctial eccentricity vector derivative.
  324.      * <p>
  325.      * If the orbit was created without derivatives, the value returned is null.
  326.      * </p>
  327.      * @return second component of the equinoctial eccentricity vector derivative
  328.      */
  329.     public abstract T getEquinoctialEyDot();

  331.     /** Get the first component of the inclination vector.
  332.      * @return first component of the inclination vector
  333.      */
  334.     public abstract T getHx();

  336.     /** Get the first component of the inclination vector derivative.
  337.      * <p>
  338.      * If the orbit was created without derivatives, the value returned is null.
  339.      * </p>
  340.      * @return first component of the inclination vector derivative
  341.      */
  342.     public abstract T getHxDot();

  344.     /** Get the second component of the inclination vector.
  345.      * @return second component of the inclination vector
  346.      */
  347.     public abstract T getHy();

  349.     /** Get the second component of the inclination vector derivative.
  350.      * @return second component of the inclination vector derivative
  351.      */
  352.     public abstract T getHyDot();

  354.     /** Get the eccentric longitude argument.
  355.      * @return E + ω + Ω eccentric longitude argument (rad)
  356.      */
  357.     public abstract T getLE();

  359.     /** Get the eccentric longitude argument derivative.
  360.      * <p>
  361.      * If the orbit was created without derivatives, the value returned is null.
  362.      * </p>
  363.      * @return d(E + ω + Ω)/dt eccentric longitude argument derivative (rad/s)
  364.      */
  365.     public abstract T getLEDot();

  367.     /** Get the true longitude argument.
  368.      * @return v + ω + Ω true longitude argument (rad)
  369.      */
  370.     public abstract T getLv();

  372.     /** Get the true longitude argument derivative.
  373.      * <p>
  374.      * If the orbit was created without derivatives, the value returned is null.
  375.      * </p>
  376.      * @return d(v + ω + Ω)/dt true longitude argument derivative (rad/s)
  377.      */
  378.     public abstract T getLvDot();

  380.     /** Get the mean longitude argument.
  381.      * @return M + ω + Ω mean longitude argument (rad)
  382.      */
  383.     public abstract T getLM();

  385.     /** Get the mean longitude argument derivative.
  386.      * <p>
  387.      * If the orbit was created without derivatives, the value returned is null.
  388.      * </p>
  389.      * @return d(M + ω + Ω)/dt mean longitude argument derivative (rad/s)
  390.      */
  391.     public abstract T getLMDot();

  393.     // Additional orbital elements

  395.     /** Get the eccentricity.
  396.      * @return eccentricity
  397.      */
  398.     public abstract T getE();

  400.     /** Get the eccentricity derivative.
  401.      * <p>
  402.      * If the orbit was created without derivatives, the value returned is null.
  403.      * </p>
  404.      * @return eccentricity derivative
  405.      */
  406.     public abstract T getEDot();

  408.     /** Get the inclination.
  409.      * <p>
  410.      * If the orbit was created without derivatives, the value returned is null.
  411.      * </p>
  412.      * @return inclination (rad)
  413.      */
  414.     public abstract T getI();

  416.     /** Get the inclination derivative.
  417.      * @return inclination derivative (rad/s)
  418.      */
  419.     public abstract T getIDot();

  421.     /** Check if orbit includes non-Keplerian rates.
  422.      * @return true if orbit includes non-Keplerian derivatives
  423.      * @see #getADot()
  424.      * @see #getEquinoctialExDot()
  425.      * @see #getEquinoctialEyDot()
  426.      * @see #getHxDot()
  427.      * @see #getHyDot()
  428.      * @see #getLEDot()
  429.      * @see #getLvDot()
  430.      * @see #getLMDot()
  431.      * @see #getEDot()
  432.      * @see #getIDot()
  433.      * @since 13.0
  434.      */
  435.     public boolean hasNonKeplerianAcceleration() {
  436.         return hasNonKeplerianAcceleration(getPVCoordinates(), getMu());
  437.     }

  439.     /** Get the central attraction coefficient used for position and velocity conversions (m³/s²).
  440.      * @return central attraction coefficient used for position and velocity conversions (m³/s²)
  441.      */
  442.     public T getMu() {
  443.         return mu;
  444.     }

  446.     /** Get the Keplerian period.
  447.      * <p>The Keplerian period is computed directly from semi major axis
  448.      * and central acceleration constant.</p>
  449.      * @return Keplerian period in seconds, or positive infinity for hyperbolic orbits
  450.      */
  451.     public T getKeplerianPeriod() {
  452.         final T a = getA();
  453.         return isElliptical() ?
  454.                a.multiply(a.getPi().multiply(2.0)).multiply(a.divide(mu).sqrt()) :
  455.                zero.add(Double.POSITIVE_INFINITY);
  456.     }

  458.     /** Get the Keplerian mean motion.
  459.      * <p>The Keplerian mean motion is computed directly from semi major axis
  460.      * and central acceleration constant.</p>
  461.      * @return Keplerian mean motion in radians per second
  462.      */
  463.     public T getKeplerianMeanMotion() {
  464.         final T absA = getA().abs();
  465.         return absA.reciprocal().multiply(mu).sqrt().divide(absA);
  466.     }

  468.     /** Get the derivative of the mean anomaly with respect to the semi major axis.
  469.      * @return derivative of the mean anomaly with respect to the semi major axis
  470.      */
  471.     public T getMeanAnomalyDotWrtA() {
  472.         return getKeplerianMeanMotion().divide(getA()).multiply(-1.5);
  473.     }

  475.     /** Get the date of orbital parameters.
  476.      * @return date of the orbital parameters
  477.      */
  478.     public FieldAbsoluteDate<T> getDate() {
  479.         return date;
  480.     }

  482.     /** Get the {@link TimeStampedPVCoordinates} in a specified frame.
  483.      * @param outputFrame frame in which the position/velocity coordinates shall be computed
  484.      * @return FieldPVCoordinates in the specified output frame
  485.           * @see #getPVCoordinates()
  486.      */
  487.     public TimeStampedFieldPVCoordinates<T> getPVCoordinates(final Frame outputFrame) {
  488.         if (pvCoordinates == null) {
  489.             pvCoordinates = initPVCoordinates();
  490.         }

  492.         // If output frame requested is the same as definition frame,
  493.         // PV coordinates are returned directly
  494.         if (outputFrame == frame) {
  495.             return pvCoordinates;
  496.         }

  498.         // Else, PV coordinates are transformed to output frame
  499.         final FieldTransform<T> t = frame.getTransformTo(outputFrame, date);
  500.         return t.transformPVCoordinates(pvCoordinates);
  501.     }

  503.     /** {@inheritDoc} */
  504.     public TimeStampedFieldPVCoordinates<T> getPVCoordinates(final FieldAbsoluteDate<T> otherDate, final Frame otherFrame) {
  505.         return shiftedBy(otherDate.durationFrom(getDate())).getPVCoordinates(otherFrame);
  506.     }

  508.     /** {@inheritDoc} */
  509.     @Override
  510.     public FieldVector3D<T> getPosition(final FieldAbsoluteDate<T> otherDate, final Frame otherFrame) {
  511.         return shiftedBy(otherDate.durationFrom(getDate())).getPosition(otherFrame);
  512.     }

  514.     /** Get the position in a specified frame.
  515.      * @param outputFrame frame in which the position coordinates shall be computed
  516.      * @return position in the specified output frame
  517.      * @see #getPosition()
  518.      * @since 12.0
  519.      */
  520.     public FieldVector3D<T> getPosition(final Frame outputFrame) {
  521.         if (position == null) {
  522.             position = initPosition();
  523.         }

  525.         // If output frame requested is the same as definition frame,
  526.         // Position vector is returned directly
  527.         if (outputFrame == frame) {
  528.             return position;
  529.         }

  531.         // Else, position vector is transformed to output frame
  532.         final FieldStaticTransform<T> t = frame.getStaticTransformTo(outputFrame, date);
  533.         return t.transformPosition(position);

  535.     }

  537.     /** Get the position in definition frame.
  538.      * @return position in the definition frame
  539.      * @see #getPVCoordinates()
  540.      * @since 12.0
  541.      */
  542.     public FieldVector3D<T> getPosition() {
  543.         if (position == null) {
  544.             position = initPosition();
  545.         }
  546.         return position;
  547.     }

  549.     /** Get the {@link TimeStampedPVCoordinates} in definition frame.
  550.      * @return FieldPVCoordinates in the definition frame
  551.      * @see #getPVCoordinates(Frame)
  552.      */
  553.     public TimeStampedFieldPVCoordinates<T> getPVCoordinates() {
  554.         if (pvCoordinates == null) {
  555.             pvCoordinates = initPVCoordinates();
  556.             position      = pvCoordinates.getPosition();
  557.         }
  558.         return pvCoordinates;
  559.     }

  561.     /** Getter for Field-valued one.
  562.      * @return one
  563.      * */
  564.     protected T getOne() {
  565.         return one;
  566.     }

  568.     /** Getter for Field-valued zero.
  569.      * @return zero
  570.      * */
  571.     protected T getZero() {
  572.         return zero;
  573.     }

  575.     /** Getter for Field.
  576.      * @return CalculusField
  577.      * */
  578.     protected Field<T> getField() {
  579.         return field;
  580.     }

  582.     /** Compute the position coordinates from the canonical parameters.
  583.      * @return computed position coordinates
  584.      * @since 12.0
  585.      */
  586.     protected abstract FieldVector3D<T> initPosition();

  588.     /** Compute the position/velocity coordinates from the canonical parameters.
  589.      * @return computed position/velocity coordinates
  590.      */
  591.     protected abstract TimeStampedFieldPVCoordinates<T> initPVCoordinates();

  593.     /**
  594.      * Create a new object representing the same physical orbital state, but attached to a different reference frame.
  595.      * If the new frame is not inertial, an exception will be thrown.
  596.      *
  597.      * @param inertialFrame reference frame of output orbit
  598.      * @return orbit with different frame
  599.      * @since 13.0
  600.      */
  601.     public abstract FieldOrbit<T> inFrame(Frame inertialFrame);

  603.     /** Get a time-shifted orbit.
  604.      * <p>
  605.      * The orbit can be slightly shifted to close dates. This shift is based on
  606.      * a simple Keplerian model. It is <em>not</em> intended as a replacement
  607.      * for proper orbit and attitude propagation but should be sufficient for
  608.      * small time shifts or coarse accuracy.
  609.      * </p>
  610.      * @param dt time shift in seconds
  611.      * @return a new orbit, shifted with respect to the instance (which is immutable)
  612.      */
  613.     public abstract FieldOrbit<T> shiftedBy(T dt);

  615.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  616.      * <p>
  617.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  618.      * respect to Cartesian coordinate j. This means each row corresponds to one orbital parameter
  619.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  620.      * </p>
  621.      * @param type type of the position angle to use
  622.      * @param jacobian placeholder 6x6 (or larger) matrix to be filled with the Jacobian, if matrix
  623.      * is larger than 6x6, only the 6x6 upper left corner will be modified
  624.      */
  625.     public void getJacobianWrtCartesian(final PositionAngleType type, final T[][] jacobian) {

  627.         final T[][] cachedJacobian;
  628.         synchronized (this) {
  629.             switch (type) {
  630.                 case MEAN :
  631.                     if (jacobianMeanWrtCartesian == null) {
  632.                         // first call, we need to compute the Jacobian and cache it
  633.                         jacobianMeanWrtCartesian = computeJacobianMeanWrtCartesian();
  634.                     }
  635.                     cachedJacobian = jacobianMeanWrtCartesian;
  636.                     break;
  637.                 case ECCENTRIC :
  638.                     if (jacobianEccentricWrtCartesian == null) {
  639.                         // first call, we need to compute the Jacobian and cache it
  640.                         jacobianEccentricWrtCartesian = computeJacobianEccentricWrtCartesian();
  641.                     }
  642.                     cachedJacobian = jacobianEccentricWrtCartesian;
  643.                     break;
  644.                 case TRUE :
  645.                     if (jacobianTrueWrtCartesian == null) {
  646.                         // first call, we need to compute the Jacobian and cache it
  647.                         jacobianTrueWrtCartesian = computeJacobianTrueWrtCartesian();
  648.                     }
  649.                     cachedJacobian = jacobianTrueWrtCartesian;
  650.                     break;
  651.                 default :
  652.                     throw new OrekitInternalError(null);
  653.             }
  654.         }

  656.         // fill the user provided array
  657.         for (int i = 0; i < cachedJacobian.length; ++i) {
  658.             System.arraycopy(cachedJacobian[i], 0, jacobian[i], 0, cachedJacobian[i].length);
  659.         }

  661.     }

  663.     /** Compute the Jacobian of the Cartesian parameters with respect to the orbital parameters.
  664.      * <p>
  665.      * Element {@code jacobian[i][j]} is the derivative of Cartesian coordinate i of the orbit with
  666.      * respect to orbital parameter j. This means each row corresponds to one Cartesian coordinate
  667.      * x, y, z, xdot, ydot, zdot whereas columns 0 to 5 correspond to the orbital parameters.
  668.      * </p>
  669.      * @param type type of the position angle to use
  670.      * @param jacobian placeholder 6x6 (or larger) matrix to be filled with the Jacobian, if matrix
  671.      * is larger than 6x6, only the 6x6 upper left corner will be modified
  672.      */
  673.     public void getJacobianWrtParameters(final PositionAngleType type, final T[][] jacobian) {

  675.         final T[][] cachedJacobian;
  676.         synchronized (this) {
  677.             switch (type) {
  678.                 case MEAN :
  679.                     if (jacobianWrtParametersMean == null) {
  680.                         // first call, we need to compute the Jacobian and cache it
  681.                         jacobianWrtParametersMean = createInverseJacobian(type);
  682.                     }
  683.                     cachedJacobian = jacobianWrtParametersMean;
  684.                     break;
  685.                 case ECCENTRIC :
  686.                     if (jacobianWrtParametersEccentric == null) {
  687.                         // first call, we need to compute the Jacobian and cache it
  688.                         jacobianWrtParametersEccentric = createInverseJacobian(type);
  689.                     }
  690.                     cachedJacobian = jacobianWrtParametersEccentric;
  691.                     break;
  692.                 case TRUE :
  693.                     if (jacobianWrtParametersTrue == null) {
  694.                         // first call, we need to compute the Jacobian and cache it
  695.                         jacobianWrtParametersTrue = createInverseJacobian(type);
  696.                     }
  697.                     cachedJacobian = jacobianWrtParametersTrue;
  698.                     break;
  699.                 default :
  700.                     throw new OrekitInternalError(null);
  701.             }
  702.         }

  704.         // fill the user-provided array
  705.         for (int i = 0; i < cachedJacobian.length; ++i) {
  706.             System.arraycopy(cachedJacobian[i], 0, jacobian[i], 0, cachedJacobian[i].length);
  707.         }

  709.     }

  711.     /** Create an inverse Jacobian.
  712.      * @param type type of the position angle to use
  713.      * @return inverse Jacobian
  714.      */
  715.     private T[][] createInverseJacobian(final PositionAngleType type) {

  717.         // get the direct Jacobian
  718.         final T[][] directJacobian = MathArrays.buildArray(getA().getField(), 6, 6);
  719.         getJacobianWrtCartesian(type, directJacobian);

  721.         // invert the direct Jacobian
  722.         final FieldMatrix<T> matrix = MatrixUtils.createFieldMatrix(directJacobian);
  723.         final FieldDecompositionSolver<T> solver = getDecompositionSolver(matrix);
  724.         return solver.getInverse().getData();

  726.     }

  728.     /**
  729.      * Method to build a matrix decomposition solver.
  730.      * @param matrix matrix
  731.      * @return solver
  732.      * @since 13.1
  733.      */
  734.     protected FieldDecompositionSolver<T> getDecompositionSolver(final FieldMatrix<T> matrix) {
  735.         return new FieldLUDecomposition<>(matrix).getSolver();
  736.     }

  738.     /** Compute the Jacobian of the orbital parameters with mean angle with respect to the Cartesian parameters.
  739.      * <p>
  740.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  741.      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
  742.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  743.      * </p>
  744.      * @return 6x6 Jacobian matrix
  745.      * @see #computeJacobianEccentricWrtCartesian()
  746.      * @see #computeJacobianTrueWrtCartesian()
  747.      */
  748.     protected abstract T[][] computeJacobianMeanWrtCartesian();

  750.     /** Compute the Jacobian of the orbital parameters with eccentric angle with respect to the Cartesian parameters.
  751.      * <p>
  752.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  753.      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
  754.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  755.      * </p>
  756.      * @return 6x6 Jacobian matrix
  757.      * @see #computeJacobianMeanWrtCartesian()
  758.      * @see #computeJacobianTrueWrtCartesian()
  759.      */
  760.     protected abstract T[][] computeJacobianEccentricWrtCartesian();

  762.     /** Compute the Jacobian of the orbital parameters with true angle with respect to the Cartesian parameters.
  763.      * <p>
  764.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  765.      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
  766.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  767.      * </p>
  768.      * @return 6x6 Jacobian matrix
  769.      * @see #computeJacobianMeanWrtCartesian()
  770.      * @see #computeJacobianEccentricWrtCartesian()
  771.      */
  772.     protected abstract T[][] computeJacobianTrueWrtCartesian();

  774.     /** Add the contribution of the Keplerian motion to parameters derivatives
  775.      * <p>
  776.      * This method is used by integration-based propagators to evaluate the part of Keplerian
  777.      * motion to evolution of the orbital state.
  778.      * </p>
  779.      * @param type type of the position angle in the state
  780.      * @param gm attraction coefficient to use
  781.      * @param pDot array containing orbital state derivatives to update (the Keplerian
  782.      * part must be <em>added</em> to the array components, as the array may already
  783.      * contain some non-zero elements corresponding to non-Keplerian parts)
  784.      */
  785.     public abstract void addKeplerContribution(PositionAngleType type, T gm, T[] pDot);

  787.         /** Fill a Jacobian half row with a single vector.
  788.      * @param a coefficient of the vector
  789.      * @param v vector
  790.      * @param row Jacobian matrix row
  791.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  792.      */

  794.     protected void fillHalfRow(final T a, final FieldVector3D<T> v, final T[] row, final int j) {
  795.         row[j]     = a.multiply(v.getX());
  796.         row[j + 1] = a.multiply(v.getY());
  797.         row[j + 2] = a.multiply(v.getZ());
  798.     }

  800.     /** Fill a Jacobian half row with a linear combination of vectors.
  801.      * @param a1 coefficient of the first vector
  802.      * @param v1 first vector
  803.      * @param a2 coefficient of the second vector
  804.      * @param v2 second vector
  805.      * @param row Jacobian matrix row
  806.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  807.      */
  808.     protected void fillHalfRow(final T a1, final FieldVector3D<T> v1, final T a2, final FieldVector3D<T> v2,
  809.                                       final T[] row, final int j) {
  810.         row[j]     = a1.linearCombination(a1, v1.getX(), a2, v2.getX());
  811.         row[j + 1] = a1.linearCombination(a1, v1.getY(), a2, v2.getY());
  812.         row[j + 2] = a1.linearCombination(a1, v1.getZ(), a2, v2.getZ());

  814.     }

  816.     /** Fill a Jacobian half row with a linear combination of vectors.
  817.      * @param a1 coefficient of the first vector
  818.      * @param v1 first vector
  819.      * @param a2 coefficient of the second vector
  820.      * @param v2 second vector
  821.      * @param a3 coefficient of the third vector
  822.      * @param v3 third vector
  823.      * @param row Jacobian matrix row
  824.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  825.      */
  826.     protected void fillHalfRow(final T a1, final FieldVector3D<T> v1,
  827.                                final T a2, final FieldVector3D<T> v2,
  828.                                final T a3, final FieldVector3D<T> v3,
  829.                                       final T[] row, final int j) {
  830.         row[j]     = a1.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX());
  831.         row[j + 1] = a1.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY());
  832.         row[j + 2] = a1.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ());

  834.     }

  836.     /** Fill a Jacobian half row with a linear combination of vectors.
  837.      * @param a1 coefficient of the first vector
  838.      * @param v1 first vector
  839.      * @param a2 coefficient of the second vector
  840.      * @param v2 second vector
  841.      * @param a3 coefficient of the third vector
  842.      * @param v3 third vector
  843.      * @param a4 coefficient of the fourth vector
  844.      * @param v4 fourth vector
  845.      * @param row Jacobian matrix row
  846.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  847.      */
  848.     protected void fillHalfRow(final T a1, final FieldVector3D<T> v1, final T a2, final FieldVector3D<T> v2,
  849.                                       final T a3, final FieldVector3D<T> v3, final T a4, final FieldVector3D<T> v4,
  850.                                       final T[] row, final int j) {
  851.         row[j]     = a1.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX(), a4, v4.getX());
  852.         row[j + 1] = a1.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY(), a4, v4.getY());
  853.         row[j + 2] = a1.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ(), a4, v4.getZ());

  855.     }

  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 a5 coefficient of the fifth vector
  867.      * @param v5 fifth vector
  868.      * @param row Jacobian matrix row
  869.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  870.      */
  871.     protected void fillHalfRow(final T a1, final FieldVector3D<T> v1, final T a2, final FieldVector3D<T> v2,
  872.                                       final T a3, final FieldVector3D<T> v3, final T a4, final FieldVector3D<T> v4,
  873.                                       final T a5, final FieldVector3D<T> v5,
  874.                                       final T[] row, final int j) {
  875.         row[j]     = a1.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX(), a4, v4.getX()).add(a5.multiply(v5.getX()));
  876.         row[j + 1] = a1.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY(), a4, v4.getY()).add(a5.multiply(v5.getY()));
  877.         row[j + 2] = a1.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ(), a4, v4.getZ()).add(a5.multiply(v5.getZ()));

  879.     }

  881.     /** Fill a Jacobian half row with a linear combination of vectors.
  882.      * @param a1 coefficient of the first vector
  883.      * @param v1 first vector
  884.      * @param a2 coefficient of the second vector
  885.      * @param v2 second vector
  886.      * @param a3 coefficient of the third vector
  887.      * @param v3 third vector
  888.      * @param a4 coefficient of the fourth vector
  889.      * @param v4 fourth vector
  890.      * @param a5 coefficient of the fifth vector
  891.      * @param v5 fifth vector
  892.      * @param a6 coefficient of the sixth vector
  893.      * @param v6 sixth vector
  894.      * @param row Jacobian matrix row
  895.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  896.      */
  897.     protected void fillHalfRow(final T a1, final FieldVector3D<T> v1, final T a2, final FieldVector3D<T> v2,
  898.                                       final T a3, final FieldVector3D<T> v3, final T a4, final FieldVector3D<T> v4,
  899.                                       final T a5, final FieldVector3D<T> v5, final T a6, final FieldVector3D<T> v6,
  900.                                       final T[] row, final int j) {
  901.         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()));
  902.         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()));
  903.         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()));

  905.     }


  908. }
Created with JaCoCo 0.8.12.202403310830