Orbit.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;

  19. import java.io.Serializable;

  21. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  22. import org.hipparchus.linear.DecompositionSolver;
  23. import org.hipparchus.linear.MatrixUtils;
  24. import org.hipparchus.linear.QRDecomposition;
  25. import org.hipparchus.linear.RealMatrix;
  26. import org.hipparchus.util.FastMath;
  27. import org.hipparchus.util.MathArrays;
  28. import org.orekit.errors.OrekitIllegalArgumentException;
  29. import org.orekit.errors.OrekitInternalError;
  30. import org.orekit.errors.OrekitMessages;
  31. import org.orekit.frames.Frame;
  32. import org.orekit.frames.StaticTransform;
  33. import org.orekit.frames.Transform;
  34. import org.orekit.time.AbsoluteDate;
  35. import org.orekit.time.TimeShiftable;
  36. import org.orekit.time.TimeStamped;
  37. import org.orekit.utils.PVCoordinates;
  38. import org.orekit.utils.PVCoordinatesProvider;
  39. import org.orekit.utils.TimeStampedPVCoordinates;

  41. /**
  42.  * This class handles orbital parameters.

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

  67.     /** Serializable UID. */
  68.     private static final long serialVersionUID = 438733454597999578L;

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

  73.     /** Date of the orbital parameters. */
  74.     private final AbsoluteDate date;

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

  79.     /** Computed position.
  80.      * @since 12.0
  81.      */
  82.     private transient Vector3D position;

  84.     /** Computed PVCoordinates. */
  85.     private transient TimeStampedPVCoordinates pvCoordinates;

  87.     /** Jacobian of the orbital parameters with mean angle with respect to the Cartesian coordinates. */
  88.     private transient double[][] jacobianMeanWrtCartesian;

  90.     /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with mean angle. */
  91.     private transient double[][] jacobianWrtParametersMean;

  93.     /** Jacobian of the orbital parameters with eccentric angle with respect to the Cartesian coordinates. */
  94.     private transient double[][] jacobianEccentricWrtCartesian;

  96.     /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with eccentric angle. */
  97.     private transient double[][] jacobianWrtParametersEccentric;

  99.     /** Jacobian of the orbital parameters with true angle with respect to the Cartesian coordinates. */
  100.     private transient double[][] jacobianTrueWrtCartesian;

  102.     /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with true angle. */
  103.     private transient double[][] jacobianWrtParametersTrue;

  105.     /** Default constructor.
  106.      * Build a new instance with arbitrary default elements.
  107.      * @param frame the frame in which the parameters are defined
  108.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  109.      * @param date date of the orbital parameters
  110.      * @param mu central attraction coefficient (m^3/s^2)
  111.      * @exception IllegalArgumentException if frame is not a {@link
  112.      * Frame#isPseudoInertial pseudo-inertial frame}
  113.      */
  114.     protected Orbit(final Frame frame, final AbsoluteDate date, final double mu)
  115.         throws IllegalArgumentException {
  116.         ensurePseudoInertialFrame(frame);
  117.         this.date                      = date;
  118.         this.mu                        = mu;
  119.         this.pvCoordinates             = null;
  120.         this.frame                     = frame;
  121.         jacobianMeanWrtCartesian       = null;
  122.         jacobianWrtParametersMean      = null;
  123.         jacobianEccentricWrtCartesian  = null;
  124.         jacobianWrtParametersEccentric = null;
  125.         jacobianTrueWrtCartesian       = null;
  126.         jacobianWrtParametersTrue      = null;
  127.     }

  129.     /** Set the orbit from Cartesian parameters.
  130.      *
  131.      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
  132.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  133.      * use {@code mu} and the position to compute the acceleration, including
  134.      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
  135.      *
  136.      * @param pvCoordinates the position and velocity in the inertial frame
  137.      * @param frame the frame in which the {@link TimeStampedPVCoordinates} are defined
  138.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  139.      * @param mu central attraction coefficient (m^3/s^2)
  140.      * @exception IllegalArgumentException if frame is not a {@link
  141.      * Frame#isPseudoInertial pseudo-inertial frame}
  142.      */
  143.     protected Orbit(final TimeStampedPVCoordinates pvCoordinates, final Frame frame, final double mu)
  144.         throws IllegalArgumentException {
  145.         ensurePseudoInertialFrame(frame);
  146.         this.date = pvCoordinates.getDate();
  147.         this.mu = mu;
  148.         if (pvCoordinates.getAcceleration().getNormSq() == 0) {
  149.             // the acceleration was not provided,
  150.             // compute it from Newtonian attraction
  151.             final double r2 = pvCoordinates.getPosition().getNormSq();
  152.             final double r3 = r2 * FastMath.sqrt(r2);
  153.             this.pvCoordinates = new TimeStampedPVCoordinates(pvCoordinates.getDate(),
  154.                                                               pvCoordinates.getPosition(),
  155.                                                               pvCoordinates.getVelocity(),
  156.                                                               new Vector3D(-mu / r3, pvCoordinates.getPosition()));
  157.         } else {
  158.             this.pvCoordinates = pvCoordinates;
  159.         }
  160.         this.frame = frame;
  161.     }

  163.     /** Check if Cartesian coordinates include non-Keplerian acceleration.
  164.      * @param pva Cartesian coordinates
  165.      * @param mu central attraction coefficient
  166.      * @return true if Cartesian coordinates include non-Keplerian acceleration
  167.      */
  168.     protected static boolean hasNonKeplerianAcceleration(final PVCoordinates pva, final double mu) {

  170.         final Vector3D a = pva.getAcceleration();
  171.         if (a == null) {
  172.             return false;
  173.         }

  175.         final Vector3D p = pva.getPosition();
  176.         final double r2 = p.getNormSq();

  178.         // Check if acceleration is relatively close to 0 compared to the Keplerian acceleration
  179.         final double tolerance = mu * 1e-9;
  180.         final Vector3D aTimesR2 = a.scalarMultiply(r2);
  181.         if (aTimesR2.getNorm() < tolerance) {
  182.             return false;
  183.         }

  185.         if ((aTimesR2.add(p.normalize().scalarMultiply(mu))).getNorm() > tolerance) {
  186.             // we have a relevant acceleration, we can compute derivatives
  187.             return true;
  188.         } else {
  189.             // the provided acceleration is either too small to be reliable (probably even 0), or NaN
  190.             return false;
  191.         }
  192.     }

  194.     /** Returns true if and only if the orbit is elliptical i.e. has a non-negative semi-major axis.
  195.      * @return true if getA() is strictly greater than 0
  196.      * @since 12.0
  197.      */
  198.     public boolean isElliptical() {
  199.         return getA() > 0.;
  200.     }

  202.     /** Get the orbit type.
  203.      * @return orbit type
  204.      */
  205.     public abstract OrbitType getType();

  207.     /** Ensure the defining frame is a pseudo-inertial frame.
  208.      * @param frame frame to check
  209.      * @exception IllegalArgumentException if frame is not a {@link
  210.      * Frame#isPseudoInertial pseudo-inertial frame}
  211.      */
  212.     private static void ensurePseudoInertialFrame(final Frame frame)
  213.         throws IllegalArgumentException {
  214.         if (!frame.isPseudoInertial()) {
  215.             throw new OrekitIllegalArgumentException(OrekitMessages.NON_PSEUDO_INERTIAL_FRAME,
  216.                                                      frame.getName());
  217.         }
  218.     }

  220.     /** Get the frame in which the orbital parameters are defined.
  221.      * @return frame in which the orbital parameters are defined
  222.      */
  223.     public Frame getFrame() {
  224.         return frame;
  225.     }

  227.     /** Get the semi-major axis.
  228.      * <p>Note that the semi-major axis is considered negative for hyperbolic orbits.</p>
  229.      * @return semi-major axis (m)
  230.      */
  231.     public abstract double getA();

  233.     /** Get the semi-major axis derivative.
  234.      * <p>Note that the semi-major axis is considered negative for hyperbolic orbits.</p>
  235.      * <p>
  236.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  237.      * </p>
  238.      * @return semi-major axis  derivative (m/s)
  239.      * @see #hasDerivatives()
  240.      * @since 9.0
  241.      */
  242.     public abstract double getADot();

  244.     /** Get the first component of the equinoctial eccentricity vector.
  245.      * @return first component of the equinoctial eccentricity vector
  246.      */
  247.     public abstract double getEquinoctialEx();

  249.     /** Get the first component of the equinoctial eccentricity vector derivative.
  250.      * <p>
  251.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  252.      * </p>
  253.      * @return first component of the equinoctial eccentricity vector derivative
  254.      * @see #hasDerivatives()
  255.      * @since 9.0
  256.      */
  257.     public abstract double getEquinoctialExDot();

  259.     /** Get the second component of the equinoctial eccentricity vector.
  260.      * @return second component of the equinoctial eccentricity vector
  261.      */
  262.     public abstract double getEquinoctialEy();

  264.     /** Get the second component of the equinoctial eccentricity vector derivative.
  265.      * <p>
  266.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  267.      * </p>
  268.      * @return second component of the equinoctial eccentricity vector derivative
  269.      * @see #hasDerivatives()
  270.      * @since 9.0
  271.      */
  272.     public abstract double getEquinoctialEyDot();

  274.     /** Get the first component of the inclination vector.
  275.      * @return first component of the inclination vector
  276.      */
  277.     public abstract double getHx();

  279.     /** Get the first component of the inclination vector derivative.
  280.      * <p>
  281.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  282.      * </p>
  283.      * @return first component of the inclination vector derivative
  284.      * @see #hasDerivatives()
  285.      * @since 9.0
  286.      */
  287.     public abstract double getHxDot();

  289.     /** Get the second component of the inclination vector.
  290.      * @return second component of the inclination vector
  291.      */
  292.     public abstract double getHy();

  294.     /** Get the second component of the inclination vector derivative.
  295.      * <p>
  296.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  297.      * </p>
  298.      * @return second component of the inclination vector derivative
  299.      * @see #hasDerivatives()
  300.      * @since 9.0
  301.      */
  302.     public abstract double getHyDot();

  304.     /** Get the eccentric longitude argument.
  305.      * @return E + ω + Ω eccentric longitude argument (rad)
  306.      */
  307.     public abstract double getLE();

  309.     /** Get the eccentric longitude argument derivative.
  310.      * <p>
  311.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  312.      * </p>
  313.      * @return d(E + ω + Ω)/dt eccentric longitude argument derivative (rad/s)
  314.      * @see #hasDerivatives()
  315.      * @since 9.0
  316.      */
  317.     public abstract double getLEDot();

  319.     /** Get the true longitude argument.
  320.      * @return v + ω + Ω true longitude argument (rad)
  321.      */
  322.     public abstract double getLv();

  324.     /** Get the true longitude argument derivative.
  325.      * <p>
  326.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  327.      * </p>
  328.      * @return d(v + ω + Ω)/dt true longitude argument derivative (rad/s)
  329.      * @see #hasDerivatives()
  330.      * @since 9.0
  331.      */
  332.     public abstract double getLvDot();

  334.     /** Get the mean longitude argument.
  335.      * @return M + ω + Ω mean longitude argument (rad)
  336.      */
  337.     public abstract double getLM();

  339.     /** Get the mean longitude argument derivative.
  340.      * <p>
  341.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  342.      * </p>
  343.      * @return d(M + ω + Ω)/dt mean longitude argument derivative (rad/s)
  344.      * @see #hasDerivatives()
  345.      * @since 9.0
  346.      */
  347.     public abstract double getLMDot();

  349.     // Additional orbital elements

  351.     /** Get the eccentricity.
  352.      * @return eccentricity
  353.      */
  354.     public abstract double getE();

  356.     /** Get the eccentricity derivative.
  357.      * <p>
  358.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  359.      * </p>
  360.      * @return eccentricity derivative
  361.      * @see #hasDerivatives()
  362.      * @since 9.0
  363.      */
  364.     public abstract double getEDot();

  366.     /** Get the inclination.
  367.      * @return inclination (rad)
  368.      */
  369.     public abstract double getI();

  371.     /** Get the inclination derivative.
  372.      * <p>
  373.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  374.      * </p>
  375.      * @return inclination derivative (rad/s)
  376.      * @see #hasDerivatives()
  377.      * @since 9.0
  378.      */
  379.     public abstract double getIDot();

  381.     /** Check if orbit includes derivatives.
  382.      * @return true if orbit includes derivatives
  383.      * @see #getADot()
  384.      * @see #getEquinoctialExDot()
  385.      * @see #getEquinoctialEyDot()
  386.      * @see #getHxDot()
  387.      * @see #getHyDot()
  388.      * @see #getLEDot()
  389.      * @see #getLvDot()
  390.      * @see #getLMDot()
  391.      * @see #getEDot()
  392.      * @see #getIDot()
  393.      * @since 9.0
  394.      */
  395.     public boolean hasDerivatives() {
  396.         return !Double.isNaN(getADot());
  397.     }

  399.     /** Get the central acceleration constant.
  400.      * @return central acceleration constant
  401.      */
  402.     public double getMu() {
  403.         return mu;
  404.     }

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

  416.     /** Get the Keplerian mean motion.
  417.      * <p>The Keplerian mean motion is computed directly from semi major axis
  418.      * and central acceleration constant.</p>
  419.      * @return Keplerian mean motion in radians per second
  420.      */
  421.     public double getKeplerianMeanMotion() {
  422.         final double absA = FastMath.abs(getA());
  423.         return FastMath.sqrt(mu / absA) / absA;
  424.     }

  426.     /** Get the derivative of the mean anomaly with respect to the semi major axis.
  427.      * @return derivative of the mean anomaly with respect to the semi major axis
  428.      */
  429.     public double getMeanAnomalyDotWrtA() {
  430.         return -1.5 * getKeplerianMeanMotion() / getA();
  431.     }

  433.     /** Get the date of orbital parameters.
  434.      * @return date of the orbital parameters
  435.      */
  436.     public AbsoluteDate 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 pvCoordinates in the specified output frame
  443.           * @see #getPVCoordinates()
  444.      */
  445.     public TimeStampedPVCoordinates getPVCoordinates(final Frame outputFrame) {
  446.         if (pvCoordinates == null) {
  447.             pvCoordinates = initPVCoordinates();
  448.         }

  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 Transform t = frame.getTransformTo(outputFrame, date);
  458.         return t.transformPVCoordinates(pvCoordinates);
  459.     }

  461.     /** {@inheritDoc} */
  462.     public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate otherDate, final Frame otherFrame) {
  463.         return shiftedBy(otherDate.durationFrom(getDate())).getPVCoordinates(otherFrame);
  464.     }

  466.     /** {@inheritDoc} */
  467.     @Override
  468.     public Vector3D getPosition(final AbsoluteDate otherDate, final Frame otherFrame) {
  469.         return shiftedBy(otherDate.durationFrom(getDate())).getPosition(otherFrame);
  470.     }

  472.     /** Get the position in a specified frame.
  473.      * @param outputFrame frame in which the position coordinates shall be computed
  474.      * @return position in the specified output frame
  475.      * @see #getPosition()
  476.      * @since 12.0
  477.      */
  478.     public Vector3D getPosition(final Frame outputFrame) {
  479.         if (position == null) {
  480.             position = initPosition();
  481.         }

  483.         // If output frame requested is the same as definition frame,
  484.         // Position vector is returned directly
  485.         if (outputFrame == frame) {
  486.             return position;
  487.         }

  489.         // Else, position vector is transformed to output frame
  490.         final StaticTransform t = frame.getStaticTransformTo(outputFrame, date);
  491.         return t.transformPosition(position);

  493.     }

  495.     /** Get the position in definition frame.
  496.      * @return position in the definition frame
  497.      * @see #getPVCoordinates()
  498.      * @since 12.0
  499.      */
  500.     public Vector3D getPosition() {
  501.         if (position == null) {
  502.             position = initPosition();
  503.         }
  504.         return position;
  505.     }

  507.     /** Get the {@link TimeStampedPVCoordinates} in definition frame.
  508.      * @return pvCoordinates in the definition frame
  509.      * @see #getPVCoordinates(Frame)
  510.      */
  511.     public TimeStampedPVCoordinates getPVCoordinates() {
  512.         if (pvCoordinates == null) {
  513.             pvCoordinates = initPVCoordinates();
  514.             position      = pvCoordinates.getPosition();
  515.         }
  516.         return pvCoordinates;
  517.     }

  519.     /** Compute the position coordinates from the canonical parameters.
  520.      * @return computed position coordinates
  521.      * @since 12.0
  522.      */
  523.     protected abstract Vector3D initPosition();

  525.     /** Compute the position/velocity coordinates from the canonical parameters.
  526.      * @return computed position/velocity coordinates
  527.      */
  528.     protected abstract TimeStampedPVCoordinates initPVCoordinates();

  530.     /** Get a time-shifted orbit.
  531.      * <p>
  532.      * The orbit can be slightly shifted to close dates. The shifting model is a
  533.      * Keplerian one if no derivatives are available in the orbit, or Keplerian
  534.      * plus quadratic effect of the non-Keplerian acceleration if derivatives are
  535.      * available. Shifting is <em>not</em> intended as a replacement for proper
  536.      * orbit propagation but should be sufficient for small time shifts or coarse
  537.      * accuracy.
  538.      * </p>
  539.      * @param dt time shift in seconds
  540.      * @return a new orbit, shifted with respect to the instance (which is immutable)
  541.      */
  542.     public abstract Orbit shiftedBy(double dt);

  544.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  545.      * <p>
  546.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  547.      * respect to Cartesian coordinate j. This means each row corresponds to one orbital parameter
  548.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  549.      * </p>
  550.      * @param type type of the position angle to use
  551.      * @param jacobian placeholder 6x6 (or larger) matrix to be filled with the Jacobian, if matrix
  552.      * is larger than 6x6, only the 6x6 upper left corner will be modified
  553.      */
  554.     public void getJacobianWrtCartesian(final PositionAngleType type, final double[][] jacobian) {

  556.         final double[][] cachedJacobian;
  557.         synchronized (this) {
  558.             switch (type) {
  559.                 case MEAN :
  560.                     if (jacobianMeanWrtCartesian == null) {
  561.                         // first call, we need to compute the Jacobian and cache it
  562.                         jacobianMeanWrtCartesian = computeJacobianMeanWrtCartesian();
  563.                     }
  564.                     cachedJacobian = jacobianMeanWrtCartesian;
  565.                     break;
  566.                 case ECCENTRIC :
  567.                     if (jacobianEccentricWrtCartesian == null) {
  568.                         // first call, we need to compute the Jacobian and cache it
  569.                         jacobianEccentricWrtCartesian = computeJacobianEccentricWrtCartesian();
  570.                     }
  571.                     cachedJacobian = jacobianEccentricWrtCartesian;
  572.                     break;
  573.                 case TRUE :
  574.                     if (jacobianTrueWrtCartesian == null) {
  575.                         // first call, we need to compute the Jacobian and cache it
  576.                         jacobianTrueWrtCartesian = computeJacobianTrueWrtCartesian();
  577.                     }
  578.                     cachedJacobian = jacobianTrueWrtCartesian;
  579.                     break;
  580.                 default :
  581.                     throw new OrekitInternalError(null);
  582.             }
  583.         }

  585.         // fill the user provided array
  586.         for (int i = 0; i < cachedJacobian.length; ++i) {
  587.             System.arraycopy(cachedJacobian[i], 0, jacobian[i], 0, cachedJacobian[i].length);
  588.         }

  590.     }

  592.     /** Compute the Jacobian of the Cartesian parameters with respect to the orbital parameters.
  593.      * <p>
  594.      * Element {@code jacobian[i][j]} is the derivative of Cartesian coordinate i of the orbit with
  595.      * respect to orbital parameter j. This means each row corresponds to one Cartesian coordinate
  596.      * x, y, z, xdot, ydot, zdot whereas columns 0 to 5 correspond to the orbital parameters.
  597.      * </p>
  598.      * @param type type of the position angle to use
  599.      * @param jacobian placeholder 6x6 (or larger) matrix to be filled with the Jacobian, if matrix
  600.      * is larger than 6x6, only the 6x6 upper left corner will be modified
  601.      */
  602.     public void getJacobianWrtParameters(final PositionAngleType type, final double[][] jacobian) {

  604.         final double[][] cachedJacobian;
  605.         synchronized (this) {
  606.             switch (type) {
  607.                 case MEAN :
  608.                     if (jacobianWrtParametersMean == null) {
  609.                         // first call, we need to compute the Jacobian and cache it
  610.                         jacobianWrtParametersMean = createInverseJacobian(type);
  611.                     }
  612.                     cachedJacobian = jacobianWrtParametersMean;
  613.                     break;
  614.                 case ECCENTRIC :
  615.                     if (jacobianWrtParametersEccentric == null) {
  616.                         // first call, we need to compute the Jacobian and cache it
  617.                         jacobianWrtParametersEccentric = createInverseJacobian(type);
  618.                     }
  619.                     cachedJacobian = jacobianWrtParametersEccentric;
  620.                     break;
  621.                 case TRUE :
  622.                     if (jacobianWrtParametersTrue == null) {
  623.                         // first call, we need to compute the Jacobian and cache it
  624.                         jacobianWrtParametersTrue = createInverseJacobian(type);
  625.                     }
  626.                     cachedJacobian = jacobianWrtParametersTrue;
  627.                     break;
  628.                 default :
  629.                     throw new OrekitInternalError(null);
  630.             }
  631.         }

  633.         // fill the user-provided array
  634.         for (int i = 0; i < cachedJacobian.length; ++i) {
  635.             System.arraycopy(cachedJacobian[i], 0, jacobian[i], 0, cachedJacobian[i].length);
  636.         }

  638.     }

  640.     /** Create an inverse Jacobian.
  641.      * @param type type of the position angle to use
  642.      * @return inverse Jacobian
  643.      */
  644.     private double[][] createInverseJacobian(final PositionAngleType type) {

  646.         // get the direct Jacobian
  647.         final double[][] directJacobian = new double[6][6];
  648.         getJacobianWrtCartesian(type, directJacobian);

  650.         // invert the direct Jacobian
  651.         final RealMatrix matrix = MatrixUtils.createRealMatrix(directJacobian);
  652.         final DecompositionSolver solver = new QRDecomposition(matrix).getSolver();
  653.         return solver.getInverse().getData();

  655.     }

  657.     /** Compute the Jacobian of the orbital parameters with mean angle with respect to the Cartesian parameters.
  658.      * <p>
  659.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  660.      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
  661.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  662.      * </p>
  663.      * <p>
  664.      * The array returned by this method will not be modified.
  665.      * </p>
  666.      * @return 6x6 Jacobian matrix
  667.      * @see #computeJacobianEccentricWrtCartesian()
  668.      * @see #computeJacobianTrueWrtCartesian()
  669.      */
  670.     protected abstract double[][] computeJacobianMeanWrtCartesian();

  672.     /** Compute the Jacobian of the orbital parameters with eccentric angle with respect to the Cartesian parameters.
  673.      * <p>
  674.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  675.      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
  676.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  677.      * </p>
  678.      * <p>
  679.      * The array returned by this method will not be modified.
  680.      * </p>
  681.      * @return 6x6 Jacobian matrix
  682.      * @see #computeJacobianMeanWrtCartesian()
  683.      * @see #computeJacobianTrueWrtCartesian()
  684.      */
  685.     protected abstract double[][] computeJacobianEccentricWrtCartesian();

  687.     /** Compute the Jacobian of the orbital parameters with true angle with respect to the Cartesian parameters.
  688.      * <p>
  689.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  690.      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
  691.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  692.      * </p>
  693.      * <p>
  694.      * The array returned by this method will not be modified.
  695.      * </p>
  696.      * @return 6x6 Jacobian matrix
  697.      * @see #computeJacobianMeanWrtCartesian()
  698.      * @see #computeJacobianEccentricWrtCartesian()
  699.      */
  700.     protected abstract double[][] computeJacobianTrueWrtCartesian();

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

  715.         /** Fill a Jacobian half row with a single vector.
  716.      * @param a coefficient of the vector
  717.      * @param v vector
  718.      * @param row Jacobian matrix row
  719.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  720.      */
  721.     protected static void fillHalfRow(final double a, final Vector3D v, final double[] row, final int j) {
  722.         row[j]     = a * v.getX();
  723.         row[j + 1] = a * v.getY();
  724.         row[j + 2] = a * v.getZ();
  725.     }

  727.     /** Fill a Jacobian half row with a linear combination of vectors.
  728.      * @param a1 coefficient of the first vector
  729.      * @param v1 first vector
  730.      * @param a2 coefficient of the second vector
  731.      * @param v2 second vector
  732.      * @param row Jacobian matrix row
  733.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  734.      */
  735.     protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2,
  736.                                       final double[] row, final int j) {
  737.         row[j]     = MathArrays.linearCombination(a1, v1.getX(), a2, v2.getX());
  738.         row[j + 1] = MathArrays.linearCombination(a1, v1.getY(), a2, v2.getY());
  739.         row[j + 2] = MathArrays.linearCombination(a1, v1.getZ(), a2, v2.getZ());
  740.     }

  742.     /** Fill a Jacobian half row with a linear combination of vectors.
  743.      * @param a1 coefficient of the first vector
  744.      * @param v1 first vector
  745.      * @param a2 coefficient of the second vector
  746.      * @param v2 second vector
  747.      * @param a3 coefficient of the third vector
  748.      * @param v3 third vector
  749.      * @param row Jacobian matrix row
  750.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  751.      */
  752.     protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2,
  753.                                       final double a3, final Vector3D v3,
  754.                                       final double[] row, final int j) {
  755.         row[j]     = MathArrays.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX());
  756.         row[j + 1] = MathArrays.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY());
  757.         row[j + 2] = MathArrays.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ());
  758.     }

  760.     /** Fill a Jacobian half row with a linear combination of vectors.
  761.      * @param a1 coefficient of the first vector
  762.      * @param v1 first vector
  763.      * @param a2 coefficient of the second vector
  764.      * @param v2 second vector
  765.      * @param a3 coefficient of the third vector
  766.      * @param v3 third vector
  767.      * @param a4 coefficient of the fourth vector
  768.      * @param v4 fourth vector
  769.      * @param row Jacobian matrix row
  770.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  771.      */
  772.     protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2,
  773.                                       final double a3, final Vector3D v3, final double a4, final Vector3D v4,
  774.                                       final double[] row, final int j) {
  775.         row[j]     = MathArrays.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX(), a4, v4.getX());
  776.         row[j + 1] = MathArrays.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY(), a4, v4.getY());
  777.         row[j + 2] = MathArrays.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ(), a4, v4.getZ());
  778.     }

  780.     /** Fill a Jacobian half row with a linear combination of vectors.
  781.      * @param a1 coefficient of the first vector
  782.      * @param v1 first vector
  783.      * @param a2 coefficient of the second vector
  784.      * @param v2 second vector
  785.      * @param a3 coefficient of the third vector
  786.      * @param v3 third vector
  787.      * @param a4 coefficient of the fourth vector
  788.      * @param v4 fourth vector
  789.      * @param a5 coefficient of the fifth vector
  790.      * @param v5 fifth vector
  791.      * @param row Jacobian matrix row
  792.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  793.      */
  794.     protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2,
  795.                                       final double a3, final Vector3D v3, final double a4, final Vector3D v4,
  796.                                       final double a5, final Vector3D v5,
  797.                                       final double[] row, final int j) {
  798.         final double[] a = new double[] {
  799.             a1, a2, a3, a4, a5
  800.         };
  801.         row[j]     = MathArrays.linearCombination(a, new double[] {
  802.             v1.getX(), v2.getX(), v3.getX(), v4.getX(), v5.getX()
  803.         });
  804.         row[j + 1] = MathArrays.linearCombination(a, new double[] {
  805.             v1.getY(), v2.getY(), v3.getY(), v4.getY(), v5.getY()
  806.         });
  807.         row[j + 2] = MathArrays.linearCombination(a, new double[] {
  808.             v1.getZ(), v2.getZ(), v3.getZ(), v4.getZ(), v5.getZ()
  809.         });
  810.     }

  812.     /** Fill a Jacobian half row with a linear combination of vectors.
  813.      * @param a1 coefficient of the first vector
  814.      * @param v1 first vector
  815.      * @param a2 coefficient of the second vector
  816.      * @param v2 second vector
  817.      * @param a3 coefficient of the third vector
  818.      * @param v3 third vector
  819.      * @param a4 coefficient of the fourth vector
  820.      * @param v4 fourth vector
  821.      * @param a5 coefficient of the fifth vector
  822.      * @param v5 fifth vector
  823.      * @param a6 coefficient of the sixth vector
  824.      * @param v6 sixth vector
  825.      * @param row Jacobian matrix row
  826.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  827.      */
  828.     protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2,
  829.                                       final double a3, final Vector3D v3, final double a4, final Vector3D v4,
  830.                                       final double a5, final Vector3D v5, final double a6, final Vector3D v6,
  831.                                       final double[] row, final int j) {
  832.         final double[] a = new double[] {
  833.             a1, a2, a3, a4, a5, a6
  834.         };
  835.         row[j]     = MathArrays.linearCombination(a, new double[] {
  836.             v1.getX(), v2.getX(), v3.getX(), v4.getX(), v5.getX(), v6.getX()
  837.         });
  838.         row[j + 1] = MathArrays.linearCombination(a, new double[] {
  839.             v1.getY(), v2.getY(), v3.getY(), v4.getY(), v5.getY(), v6.getY()
  840.         });
  841.         row[j + 2] = MathArrays.linearCombination(a, new double[] {
  842.             v1.getZ(), v2.getZ(), v3.getZ(), v4.getZ(), v5.getZ(), v6.getZ()
  843.         });
  844.     }

  846. }
Created with JaCoCo 0.8.12.202403310830