FieldOrbit.java

  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;



  21. import org.hipparchus.CalculusFieldElement;
  22. import org.hipparchus.Field;
  23. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  24. import org.hipparchus.linear.FieldDecompositionSolver;
  25. import org.hipparchus.linear.FieldLUDecomposition;
  26. import org.hipparchus.linear.FieldMatrix;
  27. import org.hipparchus.linear.MatrixUtils;
  28. import org.hipparchus.util.MathArrays;
  29. import org.orekit.errors.OrekitIllegalArgumentException;
  30. import org.orekit.errors.OrekitInternalError;
  31. import org.orekit.errors.OrekitMessages;
  32. import org.orekit.frames.FieldStaticTransform;
  33. import org.orekit.frames.FieldTransform;
  34. import org.orekit.frames.Frame;
  35. import org.orekit.time.FieldAbsoluteDate;
  36. import org.orekit.time.FieldTimeShiftable;
  37. import org.orekit.time.FieldTimeStamped;
  38. import org.orekit.utils.FieldPVCoordinates;
  39. import org.orekit.utils.FieldPVCoordinatesProvider;
  40. import org.orekit.utils.TimeStampedFieldPVCoordinates;
  41. import org.orekit.utils.TimeStampedPVCoordinates;

  43. /**
  44.  * This class handles orbital parameters.

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

  71.     /** Frame in which are defined the orbital parameters. */
  72.     private final Frame frame;

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

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

  80.     /** Field-valued one.
  81.      * @since 12.0
  82.      * */
  83.     private final T one;

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

  90.     /** Field used by this class.
  91.      * @since 12.0
  92.      */
  93.     private final Field<T> field;

  95.     /** Computed position.
  96.      * @since 12.0
  97.      */
  98.     private transient FieldVector3D<T> position;

  100.     /** Computed PVCoordinates. */
  101.     private transient TimeStampedFieldPVCoordinates<T> pvCoordinates;

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

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

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

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

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

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

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

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

  185.     /** Check if Cartesian coordinates include non-Keplerian acceleration.
  186.      * @param pva Cartesian coordinates
  187.      * @param mu central attraction coefficient
  188.      * @param <T> type of the field elements
  189.      * @return true if Cartesian coordinates include non-Keplerian acceleration
  190.      */
  191.     protected static <T extends CalculusFieldElement<T>> boolean hasNonKeplerianAcceleration(final FieldPVCoordinates<T> pva, final T mu) {

  193.         final FieldVector3D<T> a = pva.getAcceleration();
  194.         if (a == null) {
  195.             return false;
  196.         }

  198.         final FieldVector3D<T> p = pva.getPosition();
  199.         final T r2 = p.getNormSq();
  200.         final T r  = r2.sqrt();
  201.         final FieldVector3D<T> keplerianAcceleration = new FieldVector3D<>(r.multiply(r2).reciprocal().multiply(mu.negate()), p);
  202.         if (a.getNorm().getReal() > 1.0e-9 * keplerianAcceleration.getNorm().getReal()) {
  203.             // we have a relevant acceleration, we can compute derivatives
  204.             return true;
  205.         } else {
  206.             // the provided acceleration is either too small to be reliable (probably even 0), or NaN
  207.             return false;
  208.         }

  210.     }

  212.     /** Returns true if and only if the orbit is elliptical i.e. has a non-negative semi-major axis.
  213.      * @return true if getA() is strictly greater than 0
  214.      * @since 12.0
  215.      */
  216.     public boolean isElliptical() {
  217.         return getA().getReal() > 0.;
  218.     }

  220.     /** Get the orbit type.
  221.      * @return orbit type
  222.      */
  223.     public abstract OrbitType getType();

  225.     /** Ensure the defining frame is a pseudo-inertial frame.
  226.      * @param frame frame to check
  227.      * @exception IllegalArgumentException if frame is not a {@link
  228.      * Frame#isPseudoInertial pseudo-inertial frame}
  229.      */
  230.     private static void ensurePseudoInertialFrame(final Frame frame)
  231.         throws IllegalArgumentException {
  232.         if (!frame.isPseudoInertial()) {
  233.             throw new OrekitIllegalArgumentException(OrekitMessages.NON_PSEUDO_INERTIAL_FRAME,
  234.                                                      frame.getName());
  235.         }
  236.     }

  238.     /** Get the frame in which the orbital parameters are defined.
  239.      * @return frame in which the orbital parameters are defined
  240.      */
  241.     public Frame getFrame() {
  242.         return frame;
  243.     }

  245.     /**Transforms the FieldOrbit instance into an Orbit instance.
  246.      * @return Orbit instance with same properties
  247.      */
  248.     public abstract Orbit toOrbit();

  250.     /** Get the semi-major axis.
  251.      * <p>Note that the semi-major axis is considered negative for hyperbolic orbits.</p>
  252.      * @return semi-major axis (m)
  253.      */
  254.     public abstract T getA();

  256.     /** Get the semi-major axis derivative.
  257.      * <p>Note that the semi-major axis is considered negative for hyperbolic orbits.</p>
  258.      * <p>
  259.      * If the orbit was created without derivatives, the value returned is null.
  260.      * </p>
  261.      * @return semi-major axis  derivative (m/s)
  262.      */
  263.     public abstract T getADot();

  265.     /** Get the first component of the equinoctial eccentricity vector.
  266.      * @return first component of the equinoctial eccentricity vector
  267.      */
  268.     public abstract T getEquinoctialEx();

  270.     /** Get the first component of the equinoctial eccentricity vector.
  271.      * <p>
  272.      * If the orbit was created without derivatives, the value returned is null.
  273.      * </p>
  274.      * @return first component of the equinoctial eccentricity vector
  275.      */
  276.     public abstract T getEquinoctialExDot();

  278.     /** Get the second component of the equinoctial eccentricity vector.
  279.      * @return second component of the equinoctial eccentricity vector
  280.      */
  281.     public abstract T getEquinoctialEy();

  283.     /** Get the second component of the equinoctial eccentricity vector.
  284.      * <p>
  285.      * If the orbit was created without derivatives, the value returned is null.
  286.      * </p>
  287.      * @return second component of the equinoctial eccentricity vector
  288.      */
  289.     public abstract T getEquinoctialEyDot();

  291.     /** Get the first component of the inclination vector.
  292.      * @return first component of the inclination vector
  293.      */
  294.     public abstract T getHx();

  296.     /** Get the first component of the inclination vector derivative.
  297.      * <p>
  298.      * If the orbit was created without derivatives, the value returned is null.
  299.      * </p>
  300.      * @return first component of the inclination vector derivative
  301.      */
  302.     public abstract T getHxDot();

  304.     /** Get the second component of the inclination vector.
  305.      * @return second component of the inclination vector
  306.      */
  307.     public abstract T getHy();

  309.     /** Get the second component of the inclination vector derivative.
  310.      * @return second component of the inclination vector derivative
  311.      */
  312.     public abstract T getHyDot();

  314.     /** Get the eccentric longitude argument.
  315.      * @return E + ω + Ω eccentric longitude argument (rad)
  316.      */
  317.     public abstract T getLE();

  319.     /** Get the eccentric longitude argument derivative.
  320.      * <p>
  321.      * If the orbit was created without derivatives, the value returned is null.
  322.      * </p>
  323.      * @return d(E + ω + Ω)/dt eccentric longitude argument derivative (rad/s)
  324.      */
  325.     public abstract T getLEDot();

  327.     /** Get the true longitude argument.
  328.      * @return v + ω + Ω true longitude argument (rad)
  329.      */
  330.     public abstract T getLv();

  332.     /** Get the true longitude argument derivative.
  333.      * <p>
  334.      * If the orbit was created without derivatives, the value returned is null.
  335.      * </p>
  336.      * @return d(v + ω + Ω)/dt true longitude argument derivative (rad/s)
  337.      */
  338.     public abstract T getLvDot();

  340.     /** Get the mean longitude argument.
  341.      * @return M + ω + Ω mean longitude argument (rad)
  342.      */
  343.     public abstract T getLM();

  345.     /** Get the mean longitude argument derivative.
  346.      * <p>
  347.      * If the orbit was created without derivatives, the value returned is null.
  348.      * </p>
  349.      * @return d(M + ω + Ω)/dt mean longitude argument derivative (rad/s)
  350.      */
  351.     public abstract T getLMDot();

  353.     // Additional orbital elements

  355.     /** Get the eccentricity.
  356.      * @return eccentricity
  357.      */
  358.     public abstract T getE();

  360.     /** Get the eccentricity derivative.
  361.      * <p>
  362.      * If the orbit was created without derivatives, the value returned is null.
  363.      * </p>
  364.      * @return eccentricity derivative
  365.      */
  366.     public abstract T getEDot();

  368.     /** Get the inclination.
  369.      * <p>
  370.      * If the orbit was created without derivatives, the value returned is null.
  371.      * </p>
  372.      * @return inclination (rad)
  373.      */
  374.     public abstract T getI();

  376.     /** Get the inclination derivative.
  377.      * @return inclination derivative (rad/s)
  378.      */
  379.     public abstract T getIDot();

  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 abstract boolean hasDerivatives();

  397.     /** Get the central attraction coefficient used for position and velocity conversions (m³/s²).
  398.      * @return central attraction coefficient used for position and velocity conversions (m³/s²)
  399.      */
  400.     public T getMu() {
  401.         return mu;
  402.     }

  404.     /** Get the Keplerian period.
  405.      * <p>The Keplerian period is computed directly from semi major axis
  406.      * and central acceleration constant.</p>
  407.      * @return Keplerian period in seconds, or positive infinity for hyperbolic orbits
  408.      */
  409.     public T getKeplerianPeriod() {
  410.         final T a = getA();
  411.         return isElliptical() ?
  412.                a.multiply(a.getPi().multiply(2.0)).multiply(a.divide(mu).sqrt()) :
  413.                zero.add(Double.POSITIVE_INFINITY);
  414.     }

  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 T getKeplerianMeanMotion() {
  422.         final T absA = getA().abs();
  423.         return absA.reciprocal().multiply(mu).sqrt().divide(absA);
  424.     }

  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 T getMeanAnomalyDotWrtA() {
  430.         return getKeplerianMeanMotion().divide(getA()).multiply(-1.5);
  431.     }

  433.     /** Get the date of orbital parameters.
  434.      * @return date of the orbital parameters
  435.      */
  436.     public FieldAbsoluteDate<T> getDate() {
  437.         return date;
  438.     }

  440.     /** Get the {@link TimeStampedPVCoordinates} in a specified frame.
  441.      * @param outputFrame frame in which the position/velocity coordinates shall be computed
  442.      * @return FieldPVCoordinates in the specified output frame
  443.           * @see #getPVCoordinates()
  444.      */
  445.     public TimeStampedFieldPVCoordinates<T> getPVCoordinates(final Frame outputFrame) {
  446.         if (pvCoordinates == null) {
  447.             pvCoordinates = initPVCoordinates();
  448.         }

  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.         }

  456.         // Else, PV coordinates are transformed to output frame
  457.         final FieldTransform<T> t = frame.getTransformTo(outputFrame, date);
  458.         return t.transformPVCoordinates(pvCoordinates);
  459.     }

  461.     /** {@inheritDoc} */
  462.     public TimeStampedFieldPVCoordinates<T> getPVCoordinates(final FieldAbsoluteDate<T> otherDate, final Frame otherFrame) {
  463.         return shiftedBy(otherDate.durationFrom(getDate())).getPVCoordinates(otherFrame);
  464.     }

  466.     /** Get the position in a specified frame.
  467.      * @param outputFrame frame in which the position coordinates shall be computed
  468.      * @return position in the specified output frame
  469.      * @see #getPosition()
  470.      * @since 12.0
  471.      */
  472.     public FieldVector3D<T> getPosition(final Frame outputFrame) {
  473.         if (position == null) {
  474.             position = initPosition();
  475.         }

  477.         // If output frame requested is the same as definition frame,
  478.         // Position vector is returned directly
  479.         if (outputFrame == frame) {
  480.             return position;
  481.         }

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

  487.     }

  489.     /** Get the position in definition frame.
  490.      * @return position in the definition frame
  491.      * @see #getPVCoordinates()
  492.      * @since 12.0
  493.      */
  494.     public FieldVector3D<T> getPosition() {
  495.         if (position == null) {
  496.             position = initPosition();
  497.         }
  498.         return position;
  499.     }

  501.     /** Get the {@link TimeStampedPVCoordinates} in definition frame.
  502.      * @return FieldPVCoordinates in the definition frame
  503.      * @see #getPVCoordinates(Frame)
  504.      */
  505.     public TimeStampedFieldPVCoordinates<T> getPVCoordinates() {
  506.         if (pvCoordinates == null) {
  507.             pvCoordinates = initPVCoordinates();
  508.             position      = pvCoordinates.getPosition();
  509.         }
  510.         return pvCoordinates;
  511.     }

  513.     /** Getter for Field-valued one.
  514.      * @return one
  515.      * */
  516.     protected T getOne() {
  517.         return one;
  518.     }

  520.     /** Getter for Field-valued zero.
  521.      * @return zero
  522.      * */
  523.     protected T getZero() {
  524.         return zero;
  525.     }

  527.     /** Getter for Field.
  528.      * @return CalculusField
  529.      * */
  530.     protected Field<T> getField() {
  531.         return field;
  532.     }

  534.     /** Compute the position coordinates from the canonical parameters.
  535.      * @return computed position coordinates
  536.      * @since 12.0
  537.      */
  538.     protected abstract FieldVector3D<T> initPosition();

  540.     /** Compute the position/velocity coordinates from the canonical parameters.
  541.      * @return computed position/velocity coordinates
  542.      */
  543.     protected abstract TimeStampedFieldPVCoordinates<T> initPVCoordinates();

  545.     /** Get a time-shifted orbit.
  546.      * <p>
  547.      * The orbit can be slightly shifted to close dates. This shift is based on
  548.      * a simple Keplerian model. It is <em>not</em> intended as a replacement
  549.      * for proper orbit and attitude propagation but should be sufficient for
  550.      * small time shifts or coarse accuracy.
  551.      * </p>
  552.      * @param dt time shift in seconds
  553.      * @return a new orbit, shifted with respect to the instance (which is immutable)
  554.      */
  555.     public abstract FieldOrbit<T> shiftedBy(T dt);

  557.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  558.      * <p>
  559.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  560.      * respect to Cartesian coordinate j. This means each row corresponds to one orbital parameter
  561.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  562.      * </p>
  563.      * @param type type of the position angle to use
  564.      * @param jacobian placeholder 6x6 (or larger) matrix to be filled with the Jacobian, if matrix
  565.      * is larger than 6x6, only the 6x6 upper left corner will be modified
  566.      */
  567.     public void getJacobianWrtCartesian(final PositionAngleType type, final T[][] jacobian) {

  569.         final T[][] cachedJacobian;
  570.         synchronized (this) {
  571.             switch (type) {
  572.                 case MEAN :
  573.                     if (jacobianMeanWrtCartesian == null) {
  574.                         // first call, we need to compute the Jacobian and cache it
  575.                         jacobianMeanWrtCartesian = computeJacobianMeanWrtCartesian();
  576.                     }
  577.                     cachedJacobian = jacobianMeanWrtCartesian;
  578.                     break;
  579.                 case ECCENTRIC :
  580.                     if (jacobianEccentricWrtCartesian == null) {
  581.                         // first call, we need to compute the Jacobian and cache it
  582.                         jacobianEccentricWrtCartesian = computeJacobianEccentricWrtCartesian();
  583.                     }
  584.                     cachedJacobian = jacobianEccentricWrtCartesian;
  585.                     break;
  586.                 case TRUE :
  587.                     if (jacobianTrueWrtCartesian == null) {
  588.                         // first call, we need to compute the Jacobian and cache it
  589.                         jacobianTrueWrtCartesian = computeJacobianTrueWrtCartesian();
  590.                     }
  591.                     cachedJacobian = jacobianTrueWrtCartesian;
  592.                     break;
  593.                 default :
  594.                     throw new OrekitInternalError(null);
  595.             }
  596.         }

  598.         // fill the user provided array
  599.         for (int i = 0; i < cachedJacobian.length; ++i) {
  600.             System.arraycopy(cachedJacobian[i], 0, jacobian[i], 0, cachedJacobian[i].length);
  601.         }

  603.     }

  605.     /** Compute the Jacobian of the Cartesian parameters with respect to the orbital parameters.
  606.      * <p>
  607.      * Element {@code jacobian[i][j]} is the derivative of Cartesian coordinate i of the orbit with
  608.      * respect to orbital parameter j. This means each row corresponds to one Cartesian coordinate
  609.      * x, y, z, xdot, ydot, zdot whereas columns 0 to 5 correspond to the orbital parameters.
  610.      * </p>
  611.      * @param type type of the position angle to use
  612.      * @param jacobian placeholder 6x6 (or larger) matrix to be filled with the Jacobian, if matrix
  613.      * is larger than 6x6, only the 6x6 upper left corner will be modified
  614.      */
  615.     public void getJacobianWrtParameters(final PositionAngleType type, final T[][] jacobian) {

  617.         final T[][] cachedJacobian;
  618.         synchronized (this) {
  619.             switch (type) {
  620.                 case MEAN :
  621.                     if (jacobianWrtParametersMean == null) {
  622.                         // first call, we need to compute the Jacobian and cache it
  623.                         jacobianWrtParametersMean = createInverseJacobian(type);
  624.                     }
  625.                     cachedJacobian = jacobianWrtParametersMean;
  626.                     break;
  627.                 case ECCENTRIC :
  628.                     if (jacobianWrtParametersEccentric == null) {
  629.                         // first call, we need to compute the Jacobian and cache it
  630.                         jacobianWrtParametersEccentric = createInverseJacobian(type);
  631.                     }
  632.                     cachedJacobian = jacobianWrtParametersEccentric;
  633.                     break;
  634.                 case TRUE :
  635.                     if (jacobianWrtParametersTrue == null) {
  636.                         // first call, we need to compute the Jacobian and cache it
  637.                         jacobianWrtParametersTrue = createInverseJacobian(type);
  638.                     }
  639.                     cachedJacobian = jacobianWrtParametersTrue;
  640.                     break;
  641.                 default :
  642.                     throw new OrekitInternalError(null);
  643.             }
  644.         }

  646.         // fill the user-provided array
  647.         for (int i = 0; i < cachedJacobian.length; ++i) {
  648.             System.arraycopy(cachedJacobian[i], 0, jacobian[i], 0, cachedJacobian[i].length);
  649.         }

  651.     }

  653.     /** Create an inverse Jacobian.
  654.      * @param type type of the position angle to use
  655.      * @return inverse Jacobian
  656.      */
  657.     private T[][] createInverseJacobian(final PositionAngleType type) {

  659.         // get the direct Jacobian
  660.         final T[][] directJacobian = MathArrays.buildArray(getA().getField(), 6, 6);
  661.         getJacobianWrtCartesian(type, directJacobian);

  663.         // invert the direct Jacobian
  664.         final FieldMatrix<T> matrix = MatrixUtils.createFieldMatrix(directJacobian);
  665.         final FieldDecompositionSolver<T> solver = new FieldLUDecomposition<>(matrix).getSolver();
  666.         return solver.getInverse().getData();

  668.     }

  670.     /** Compute the Jacobian of the orbital parameters with mean angle with respect to the Cartesian parameters.
  671.      * <p>
  672.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  673.      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
  674.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  675.      * </p>
  676.      * @return 6x6 Jacobian matrix
  677.      * @see #computeJacobianEccentricWrtCartesian()
  678.      * @see #computeJacobianTrueWrtCartesian()
  679.      */
  680.     protected abstract T[][] computeJacobianMeanWrtCartesian();

  682.     /** Compute the Jacobian of the orbital parameters with eccentric angle with respect to the Cartesian parameters.
  683.      * <p>
  684.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  685.      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
  686.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  687.      * </p>
  688.      * @return 6x6 Jacobian matrix
  689.      * @see #computeJacobianMeanWrtCartesian()
  690.      * @see #computeJacobianTrueWrtCartesian()
  691.      */
  692.     protected abstract T[][] computeJacobianEccentricWrtCartesian();

  694.     /** Compute the Jacobian of the orbital parameters with true angle with respect to the Cartesian parameters.
  695.      * <p>
  696.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  697.      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
  698.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  699.      * </p>
  700.      * @return 6x6 Jacobian matrix
  701.      * @see #computeJacobianMeanWrtCartesian()
  702.      * @see #computeJacobianEccentricWrtCartesian()
  703.      */
  704.     protected abstract T[][] computeJacobianTrueWrtCartesian();

  706.     /** Add the contribution of the Keplerian motion to parameters derivatives
  707.      * <p>
  708.      * This method is used by integration-based propagators to evaluate the part of Keplerian
  709.      * motion to evolution of the orbital state.
  710.      * </p>
  711.      * @param type type of the position angle in the state
  712.      * @param gm attraction coefficient to use
  713.      * @param pDot array containing orbital state derivatives to update (the Keplerian
  714.      * part must be <em>added</em> to the array components, as the array may already
  715.      * contain some non-zero elements corresponding to non-Keplerian parts)
  716.      */
  717.     public abstract void addKeplerContribution(PositionAngleType type, T gm, T[] pDot);

  719.         /** Fill a Jacobian half row with a single vector.
  720.      * @param a coefficient of the vector
  721.      * @param v vector
  722.      * @param row Jacobian matrix row
  723.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  724.      */

  726.     protected void fillHalfRow(final T a, final FieldVector3D<T> v, final T[] row, final int j) {
  727.         row[j]     = a.multiply(v.getX());
  728.         row[j + 1] = a.multiply(v.getY());
  729.         row[j + 2] = a.multiply(v.getZ());
  730.     }

  732.     /** Fill a Jacobian half row with a linear combination of vectors.
  733.      * @param a1 coefficient of the first vector
  734.      * @param v1 first vector
  735.      * @param a2 coefficient of the second vector
  736.      * @param v2 second vector
  737.      * @param row Jacobian matrix row
  738.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  739.      */
  740.     protected void fillHalfRow(final T a1, final FieldVector3D<T> v1, final T a2, final FieldVector3D<T> v2,
  741.                                       final T[] row, final int j) {
  742.         row[j]     = a1.linearCombination(a1, v1.getX(), a2, v2.getX());
  743.         row[j + 1] = a1.linearCombination(a1, v1.getY(), a2, v2.getY());
  744.         row[j + 2] = a1.linearCombination(a1, v1.getZ(), a2, v2.getZ());

  746.     }

  748.     /** Fill a Jacobian half row with a linear combination of vectors.
  749.      * @param a1 coefficient of the first vector
  750.      * @param v1 first vector
  751.      * @param a2 coefficient of the second vector
  752.      * @param v2 second vector
  753.      * @param a3 coefficient of the third vector
  754.      * @param v3 third vector
  755.      * @param row Jacobian matrix row
  756.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  757.      */
  758.     protected void fillHalfRow(final T a1, final FieldVector3D<T> v1,
  759.                                final T a2, final FieldVector3D<T> v2,
  760.                                final T a3, final FieldVector3D<T> v3,
  761.                                       final T[] row, final int j) {
  762.         row[j]     = a1.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX());
  763.         row[j + 1] = a1.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY());
  764.         row[j + 2] = a1.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ());

  766.     }

  768.     /** Fill a Jacobian half row with a linear combination of vectors.
  769.      * @param a1 coefficient of the first vector
  770.      * @param v1 first vector
  771.      * @param a2 coefficient of the second vector
  772.      * @param v2 second vector
  773.      * @param a3 coefficient of the third vector
  774.      * @param v3 third vector
  775.      * @param a4 coefficient of the fourth vector
  776.      * @param v4 fourth vector
  777.      * @param row Jacobian matrix row
  778.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  779.      */
  780.     protected void fillHalfRow(final T a1, final FieldVector3D<T> v1, final T a2, final FieldVector3D<T> v2,
  781.                                       final T a3, final FieldVector3D<T> v3, final T a4, final FieldVector3D<T> v4,
  782.                                       final T[] row, final int j) {
  783.         row[j]     = a1.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX(), a4, v4.getX());
  784.         row[j + 1] = a1.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY(), a4, v4.getY());
  785.         row[j + 2] = a1.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ(), a4, v4.getZ());

  787.     }

  789.     /** Fill a Jacobian half row with a linear combination of vectors.
  790.      * @param a1 coefficient of the first vector
  791.      * @param v1 first vector
  792.      * @param a2 coefficient of the second vector
  793.      * @param v2 second vector
  794.      * @param a3 coefficient of the third vector
  795.      * @param v3 third vector
  796.      * @param a4 coefficient of the fourth vector
  797.      * @param v4 fourth vector
  798.      * @param a5 coefficient of the fifth vector
  799.      * @param v5 fifth vector
  800.      * @param row Jacobian matrix row
  801.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  802.      */
  803.     protected void fillHalfRow(final T a1, final FieldVector3D<T> v1, final T a2, final FieldVector3D<T> v2,
  804.                                       final T a3, final FieldVector3D<T> v3, final T a4, final FieldVector3D<T> v4,
  805.                                       final T a5, final FieldVector3D<T> v5,
  806.                                       final T[] row, final int j) {
  807.         row[j]     = a1.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX(), a4, v4.getX()).add(a5.multiply(v5.getX()));
  808.         row[j + 1] = a1.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY(), a4, v4.getY()).add(a5.multiply(v5.getY()));
  809.         row[j + 2] = a1.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ(), a4, v4.getZ()).add(a5.multiply(v5.getZ()));

  811.     }

  813.     /** Fill a Jacobian half row with a linear combination of vectors.
  814.      * @param a1 coefficient of the first vector
  815.      * @param v1 first vector
  816.      * @param a2 coefficient of the second vector
  817.      * @param v2 second vector
  818.      * @param a3 coefficient of the third vector
  819.      * @param v3 third vector
  820.      * @param a4 coefficient of the fourth vector
  821.      * @param v4 fourth vector
  822.      * @param a5 coefficient of the fifth vector
  823.      * @param v5 fifth vector
  824.      * @param a6 coefficient of the sixth vector
  825.      * @param v6 sixth 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 a3, final FieldVector3D<T> v3, final T a4, final FieldVector3D<T> v4,
  831.                                       final T a5, final FieldVector3D<T> v5, final T a6, final FieldVector3D<T> v6,
  832.                                       final T[] row, final int j) {
  833.         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()));
  834.         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()));
  835.         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()));

  837.     }


  840. }
Created with JaCoCo 0.8.11.202310140853