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

  19. import org.hipparchus.analysis.differentiation.UnivariateDerivative1;
  20. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  21. import org.hipparchus.util.FastMath;
  22. import org.hipparchus.util.SinCos;
  23. import org.orekit.errors.OrekitException;
  24. import org.orekit.errors.OrekitIllegalArgumentException;
  25. import org.orekit.errors.OrekitInternalError;
  26. import org.orekit.errors.OrekitMessages;
  27. import org.orekit.frames.Frame;
  28. import org.orekit.frames.KinematicTransform;
  29. import org.orekit.time.AbsoluteDate;
  30. import org.orekit.time.TimeOffset;
  31. import org.orekit.utils.PVCoordinates;
  32. import org.orekit.utils.TimeStampedPVCoordinates;


  35. /**
  36.  * This class handles traditional Keplerian orbital parameters.

  38.  * <p>
  39.  * The parameters used internally are the classical Keplerian elements:
  40.  *   <pre>
  41.  *     a
  42.  *     e
  43.  *     i
  44.  *     ω
  45.  *     Ω
  46.  *     v
  47.  *   </pre>
  48.  * where ω stands for the Perigee Argument, Ω stands for the
  49.  * Right Ascension of the Ascending Node and v stands for the true anomaly.
  50.  *
  51.  * <p>
  52.  * This class supports hyperbolic orbits, using the convention that semi major
  53.  * axis is negative for such orbits (and of course eccentricity is greater than 1).
  54.  * </p>
  55.  * <p>
  56.  * When orbit is either equatorial or circular, some Keplerian elements
  57.  * (more precisely ω and Ω) become ambiguous so this class should not
  58.  * be used for such orbits. For this reason, {@link EquinoctialOrbit equinoctial
  59.  * orbits} is the recommended way to represent orbits.
  60.  * </p>

  62.  * <p>
  63.  * The instance <code>KeplerianOrbit</code> is guaranteed to be immutable.
  64.  * </p>
  65.  * @see     Orbit
  66.  * @see    CircularOrbit
  67.  * @see    CartesianOrbit
  68.  * @see    EquinoctialOrbit
  69.  * @author Luc Maisonobe
  70.  * @author Guylaine Prat
  71.  * @author Fabien Maussion
  72.  * @author V&eacute;ronique Pommier-Maurussane
  73.  */
  74. public class KeplerianOrbit extends Orbit implements PositionAngleBased<KeplerianOrbit> {

  76.     /** Name of the eccentricity parameter. */
  77.     private static final String ECCENTRICITY = "eccentricity";

  79.     /** Semi-major axis (m). */
  80.     private final double a;

  82.     /** Eccentricity. */
  83.     private final double e;

  85.     /** Inclination (rad). */
  86.     private final double i;

  88.     /** Perigee Argument (rad). */
  89.     private final double pa;

  91.     /** Right Ascension of Ascending Node (rad). */
  92.     private final double raan;

  94.     /** Cached anomaly (rad). */
  95.     private final double cachedAnomaly;

  97.     /** Semi-major axis derivative (m/s). */
  98.     private final double aDot;

  100.     /** Eccentricity derivative. */
  101.     private final double eDot;

  103.     /** Inclination derivative (rad/s). */
  104.     private final double iDot;

  106.     /** Perigee Argument derivative (rad/s). */
  107.     private final double paDot;

  109.     /** Right Ascension of Ascending Node derivative (rad/s). */
  110.     private final double raanDot;

  112.     /** Derivative of cached anomaly (rad/s). */
  113.     private final double cachedAnomalyDot;

  115.     /** Cached type of position angle. */
  116.     private final PositionAngleType cachedPositionAngleType;

  118.     /** Partial Cartesian coordinates (position and velocity are valid, acceleration may be missing). */
  119.     private transient PVCoordinates partialPV;

  121.     /** Creates a new instance.
  122.      * @param a  semi-major axis (m), negative for hyperbolic orbits
  123.      * @param e eccentricity (positive or equal to 0)
  124.      * @param i inclination (rad)
  125.      * @param pa perigee argument (ω, rad)
  126.      * @param raan right ascension of ascending node (Ω, rad)
  127.      * @param anomaly mean, eccentric or true anomaly (rad)
  128.      * @param type type of anomaly
  129.      * @param cachedPositionAngleType type of cached anomaly
  130.      * @param frame the frame in which the parameters are defined
  131.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  132.      * @param date date of the orbital parameters
  133.      * @param mu central attraction coefficient (m³/s²)
  134.      * @exception IllegalArgumentException if frame is not a {@link
  135.      * Frame#isPseudoInertial pseudo-inertial frame} or a and e don't match for hyperbolic orbits,
  136.      * or v is out of range for hyperbolic orbits
  137.      * @since 12.1
  138.      */
  139.     public KeplerianOrbit(final double a, final double e, final double i,
  140.                           final double pa, final double raan, final double anomaly,
  141.                           final PositionAngleType type, final PositionAngleType cachedPositionAngleType,
  142.                           final Frame frame, final AbsoluteDate date, final double mu)
  143.             throws IllegalArgumentException {
  144.         this(a, e, i, pa, raan, anomaly,
  145.                 0., 0., 0., 0., 0., computeKeplerianAnomalyDot(type, a, e, mu, anomaly, type),
  146.                 type, cachedPositionAngleType, frame, date, mu);
  147.     }

  149.     /** Creates a new instance without derivatives and with cached position angle same as value inputted.
  150.      * @param a  semi-major axis (m), negative for hyperbolic orbits
  151.      * @param e eccentricity (positive or equal to 0)
  152.      * @param i inclination (rad)
  153.      * @param pa perigee argument (ω, rad)
  154.      * @param raan right ascension of ascending node (Ω, rad)
  155.      * @param anomaly mean, eccentric or true anomaly (rad)
  156.      * @param type type of anomaly
  157.      * @param frame the frame in which the parameters are defined
  158.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  159.      * @param date date of the orbital parameters
  160.      * @param mu central attraction coefficient (m³/s²)
  161.      * @exception IllegalArgumentException if frame is not a {@link
  162.      * Frame#isPseudoInertial pseudo-inertial frame} or a and e don't match for hyperbolic orbits,
  163.      * or v is out of range for hyperbolic orbits
  164.      */
  165.     public KeplerianOrbit(final double a, final double e, final double i,
  166.                           final double pa, final double raan, final double anomaly,
  167.                           final PositionAngleType type,
  168.                           final Frame frame, final AbsoluteDate date, final double mu)
  169.             throws IllegalArgumentException {
  170.         this(a, e, i, pa, raan, anomaly, type, type, frame, date, mu);
  171.     }

  173.     /** Creates a new instance.
  174.      * @param a  semi-major axis (m), negative for hyperbolic orbits
  175.      * @param e eccentricity (positive or equal to 0)
  176.      * @param i inclination (rad)
  177.      * @param pa perigee argument (ω, rad)
  178.      * @param raan right ascension of ascending node (Ω, rad)
  179.      * @param anomaly mean, eccentric or true anomaly (rad)
  180.      * @param aDot  semi-major axis derivative (m/s)
  181.      * @param eDot eccentricity derivative
  182.      * @param iDot inclination derivative (rad/s)
  183.      * @param paDot perigee argument derivative (rad/s)
  184.      * @param raanDot right ascension of ascending node derivative (rad/s)
  185.      * @param anomalyDot mean, eccentric or true anomaly derivative (rad/s)
  186.      * @param type type of input anomaly
  187.      * @param cachedPositionAngleType type of cached anomaly
  188.      * @param frame the frame in which the parameters are defined
  189.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  190.      * @param date date of the orbital parameters
  191.      * @param mu central attraction coefficient (m³/s²)
  192.      * @exception IllegalArgumentException if frame is not a {@link
  193.      * Frame#isPseudoInertial pseudo-inertial frame} or a and e don't match for hyperbolic orbits,
  194.      * or v is out of range for hyperbolic orbits
  195.      * @since 12.1
  196.      */
  197.     public KeplerianOrbit(final double a, final double e, final double i,
  198.                           final double pa, final double raan, final double anomaly,
  199.                           final double aDot, final double eDot, final double iDot,
  200.                           final double paDot, final double raanDot, final double anomalyDot,
  201.                           final PositionAngleType type, final PositionAngleType cachedPositionAngleType,
  202.                           final Frame frame, final AbsoluteDate date, final double mu)
  203.             throws IllegalArgumentException {
  204.         super(frame, date, mu);
  205.         this.cachedPositionAngleType = cachedPositionAngleType;

  207.         if (a * (1 - e) < 0) {
  208.             throw new OrekitIllegalArgumentException(OrekitMessages.ORBIT_A_E_MISMATCH_WITH_CONIC_TYPE, a, e);
  209.         }

  211.         // Checking eccentricity range
  212.         checkParameterRangeInclusive(ECCENTRICITY, e, 0.0, Double.POSITIVE_INFINITY);

  214.         this.a       = a;
  215.         this.aDot    = aDot;
  216.         this.e       = e;
  217.         this.eDot    = eDot;
  218.         this.i       = i;
  219.         this.iDot    = iDot;
  220.         this.pa      = pa;
  221.         this.paDot   = paDot;
  222.         this.raan    = raan;
  223.         this.raanDot = raanDot;

  225.         final UnivariateDerivative1 cachedAnomalyUD = initializeCachedAnomaly(anomaly, anomalyDot, type);
  226.         this.cachedAnomaly = cachedAnomalyUD.getValue();
  227.         this.cachedAnomalyDot = cachedAnomalyUD.getFirstDerivative();

  229.         // check true anomaly range
  230.         if (!isElliptical()) {
  231.             final double trueAnomaly = getTrueAnomaly();
  232.             if (1 + e * FastMath.cos(trueAnomaly) <= 0) {
  233.                 final double vMax = FastMath.acos(-1 / e);
  234.                 throw new OrekitIllegalArgumentException(OrekitMessages.ORBIT_ANOMALY_OUT_OF_HYPERBOLIC_RANGE,
  235.                         trueAnomaly, e, -vMax, vMax);
  236.             }
  237.         }

  239.         this.partialPV = null;

  241.     }

  243.     /** Creates a new instance with cached position angle same as value inputted.
  244.      * @param a  semi-major axis (m), negative for hyperbolic orbits
  245.      * @param e eccentricity (positive or equal to 0)
  246.      * @param i inclination (rad)
  247.      * @param pa perigee argument (ω, rad)
  248.      * @param raan right ascension of ascending node (Ω, rad)
  249.      * @param anomaly mean, eccentric or true anomaly (rad)
  250.      * @param aDot  semi-major axis derivative (m/s)
  251.      * @param eDot eccentricity derivative
  252.      * @param iDot inclination derivative (rad/s)
  253.      * @param paDot perigee argument derivative (rad/s)
  254.      * @param raanDot right ascension of ascending node derivative (rad/s)
  255.      * @param anomalyDot mean, eccentric or true anomaly derivative (rad/s)
  256.      * @param type type of anomaly
  257.      * @param frame the frame in which the parameters are defined
  258.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  259.      * @param date date of the orbital parameters
  260.      * @param mu central attraction coefficient (m³/s²)
  261.      * @exception IllegalArgumentException if frame is not a {@link
  262.      * Frame#isPseudoInertial pseudo-inertial frame} or a and e don't match for hyperbolic orbits,
  263.      * or v is out of range for hyperbolic orbits
  264.      * @since 9.0
  265.      */
  266.     public KeplerianOrbit(final double a, final double e, final double i,
  267.                           final double pa, final double raan, final double anomaly,
  268.                           final double aDot, final double eDot, final double iDot,
  269.                           final double paDot, final double raanDot, final double anomalyDot,
  270.                           final PositionAngleType type,
  271.                           final Frame frame, final AbsoluteDate date, final double mu)
  272.             throws IllegalArgumentException {
  273.         this(a, e, i, pa, raan, anomaly, aDot, eDot, iDot, paDot, raanDot, anomalyDot, type, type,
  274.                 frame, date, mu);
  275.     }

  277.     /** Constructor from Cartesian parameters.
  278.      *
  279.      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
  280.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  281.      * use {@code mu} and the position to compute the acceleration, including
  282.      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
  283.      *
  284.      * @param pvCoordinates the PVCoordinates of the satellite
  285.      * @param frame the frame in which are defined the {@link PVCoordinates}
  286.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  287.      * @param mu central attraction coefficient (m³/s²)
  288.      * @exception IllegalArgumentException if frame is not a {@link
  289.      * Frame#isPseudoInertial pseudo-inertial frame}
  290.      */
  291.     public KeplerianOrbit(final TimeStampedPVCoordinates pvCoordinates,
  292.                           final Frame frame, final double mu)
  293.             throws IllegalArgumentException {
  294.         this(pvCoordinates, frame, mu, hasNonKeplerianAcceleration(pvCoordinates, mu));
  295.     }

  297.     /** Constructor from Cartesian parameters.
  298.      *
  299.      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
  300.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  301.      * use {@code mu} and the position to compute the acceleration, including
  302.      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
  303.      *
  304.      * @param pvCoordinates the PVCoordinates of the satellite
  305.      * @param frame the frame in which are defined the {@link PVCoordinates}
  306.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  307.      * @param mu central attraction coefficient (m³/s²)
  308.      * @param reliableAcceleration if true, the acceleration is considered to be reliable
  309.      * @exception IllegalArgumentException if frame is not a {@link
  310.      * Frame#isPseudoInertial pseudo-inertial frame}
  311.      */
  312.     private KeplerianOrbit(final TimeStampedPVCoordinates pvCoordinates,
  313.                            final Frame frame, final double mu,
  314.                            final boolean reliableAcceleration)
  315.             throws IllegalArgumentException {
  316.         super(pvCoordinates, frame, mu);

  318.         // compute inclination
  319.         final Vector3D momentum = pvCoordinates.getMomentum();
  320.         final double m2 = momentum.getNormSq();
  321.         i = Vector3D.angle(momentum, Vector3D.PLUS_K);

  323.         // compute right ascension of ascending node
  324.         raan = Vector3D.crossProduct(Vector3D.PLUS_K, momentum).getAlpha();

  326.         // preliminary computations for parameters depending on orbit shape (elliptic or hyperbolic)
  327.         final Vector3D pvP     = pvCoordinates.getPosition();
  328.         final Vector3D pvV     = pvCoordinates.getVelocity();
  329.         final Vector3D pvA     = pvCoordinates.getAcceleration();
  330.         final double   r2      = pvP.getNormSq();
  331.         final double   r       = FastMath.sqrt(r2);
  332.         final double   V2      = pvV.getNormSq();
  333.         final double   rV2OnMu = r * V2 / mu;

  335.         // compute semi-major axis (will be negative for hyperbolic orbits)
  336.         a = r / (2 - rV2OnMu);
  337.         final double muA = mu * a;

  339.         // compute cached anomaly
  340.         if (isElliptical()) {
  341.             // elliptic or circular orbit
  342.             final double eSE = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(muA);
  343.             final double eCE = rV2OnMu - 1;
  344.             e = FastMath.sqrt(eSE * eSE + eCE * eCE);
  345.             this.cachedPositionAngleType = PositionAngleType.ECCENTRIC;
  346.             cachedAnomaly = FastMath.atan2(eSE, eCE);
  347.         } else {
  348.             // hyperbolic orbit
  349.             final double eSH = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(-muA);
  350.             final double eCH = rV2OnMu - 1;
  351.             e = FastMath.sqrt(1 - m2 / muA);
  352.             this.cachedPositionAngleType = PositionAngleType.TRUE;
  353.             cachedAnomaly = KeplerianAnomalyUtility.hyperbolicEccentricToTrue(e, FastMath.log((eCH + eSH) / (eCH - eSH)) / 2);
  354.         }

  356.         // Checking eccentricity range
  357.         checkParameterRangeInclusive(ECCENTRICITY, e, 0.0, Double.POSITIVE_INFINITY);

  359.         // compute perigee argument
  360.         final Vector3D node = new Vector3D(raan, 0.0);
  361.         final double px = Vector3D.dotProduct(pvP, node);
  362.         final double py = Vector3D.dotProduct(pvP, Vector3D.crossProduct(momentum, node)) / FastMath.sqrt(m2);
  363.         pa = FastMath.atan2(py, px) - getTrueAnomaly();

  365.         partialPV = pvCoordinates;

  367.         if (reliableAcceleration) {
  368.             // we have a relevant acceleration, we can compute derivatives

  370.             final double[][] jacobian = new double[6][6];
  371.             getJacobianWrtCartesian(PositionAngleType.MEAN, jacobian);

  373.             final Vector3D keplerianAcceleration    = new Vector3D(-mu / (r * r2), pvP);
  374.             final Vector3D nonKeplerianAcceleration = pvA.subtract(keplerianAcceleration);
  375.             final double   aX                       = nonKeplerianAcceleration.getX();
  376.             final double   aY                       = nonKeplerianAcceleration.getY();
  377.             final double   aZ                       = nonKeplerianAcceleration.getZ();
  378.             aDot    = jacobian[0][3] * aX + jacobian[0][4] * aY + jacobian[0][5] * aZ;
  379.             eDot    = jacobian[1][3] * aX + jacobian[1][4] * aY + jacobian[1][5] * aZ;
  380.             iDot    = jacobian[2][3] * aX + jacobian[2][4] * aY + jacobian[2][5] * aZ;
  381.             paDot   = jacobian[3][3] * aX + jacobian[3][4] * aY + jacobian[3][5] * aZ;
  382.             raanDot = jacobian[4][3] * aX + jacobian[4][4] * aY + jacobian[4][5] * aZ;

  384.             // in order to compute cached anomaly derivative, we must compute
  385.             // mean anomaly derivative including Keplerian motion and convert to required anomaly
  386.             final double MDot = getKeplerianMeanMotion() +
  387.                     jacobian[5][3] * aX + jacobian[5][4] * aY + jacobian[5][5] * aZ;
  388.             final UnivariateDerivative1 eUD = new UnivariateDerivative1(e, eDot);
  389.             final UnivariateDerivative1 MUD = new UnivariateDerivative1(getMeanAnomaly(), MDot);
  390.             if (cachedPositionAngleType == PositionAngleType.ECCENTRIC) {
  391.                 final UnivariateDerivative1 EUD = (a < 0) ?
  392.                         FieldKeplerianAnomalyUtility.hyperbolicMeanToEccentric(eUD, MUD) :
  393.                         FieldKeplerianAnomalyUtility.ellipticMeanToEccentric(eUD, MUD);
  394.                 cachedAnomalyDot = EUD.getFirstDerivative();
  395.             } else { // TRUE
  396.                 final UnivariateDerivative1 vUD = (a < 0) ?
  397.                         FieldKeplerianAnomalyUtility.hyperbolicMeanToTrue(eUD, MUD) :
  398.                         FieldKeplerianAnomalyUtility.ellipticMeanToTrue(eUD, MUD);
  399.                 cachedAnomalyDot = vUD.getFirstDerivative();
  400.             }

  402.         } else {
  403.             // acceleration is either almost zero or NaN,
  404.             // we assume acceleration was not known
  405.             aDot    = 0.;
  406.             eDot    = 0.;
  407.             iDot    = 0.;
  408.             paDot   = 0.;
  409.             raanDot = 0.;
  410.             cachedAnomalyDot = computeKeplerianAnomalyDot(cachedPositionAngleType, a, e, mu, cachedAnomaly, cachedPositionAngleType);
  411.         }

  413.     }

  415.     /** Constructor from Cartesian parameters.
  416.      *
  417.      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
  418.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  419.      * use {@code mu} and the position to compute the acceleration, including
  420.      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
  421.      *
  422.      * @param pvCoordinates the PVCoordinates of the satellite
  423.      * @param frame the frame in which are defined the {@link PVCoordinates}
  424.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  425.      * @param date date of the orbital parameters
  426.      * @param mu central attraction coefficient (m³/s²)
  427.      * @exception IllegalArgumentException if frame is not a {@link
  428.      * Frame#isPseudoInertial pseudo-inertial frame}
  429.      */
  430.     public KeplerianOrbit(final PVCoordinates pvCoordinates,
  431.                           final Frame frame, final AbsoluteDate date, final double mu)
  432.             throws IllegalArgumentException {
  433.         this(new TimeStampedPVCoordinates(date, pvCoordinates), frame, mu);
  434.     }

  436.     /** Constructor from any kind of orbital parameters.
  437.      * @param op orbital parameters to copy
  438.      */
  439.     public KeplerianOrbit(final Orbit op) {
  440.         this(op.getPVCoordinates(), op.getFrame(), op.getMu(), op.hasNonKeplerianAcceleration());
  441.     }

  443.     /** {@inheritDoc} */
  444.     @Override
  445.     public boolean hasNonKeplerianAcceleration() {
  446.         return aDot != 0. || eDot != 0. || paDot != 0. || iDot != 0. || raanDot != 0. ||
  447.                 FastMath.abs(cachedAnomalyDot - computeKeplerianAnomalyDot(cachedPositionAngleType, a, e, getMu(), cachedAnomaly, cachedPositionAngleType)) > TOLERANCE_POSITION_ANGLE_RATE;
  448.     }

  450.     /** {@inheritDoc} */
  451.     @Override
  452.     public OrbitType getType() {
  453.         return OrbitType.KEPLERIAN;
  454.     }

  456.     /** {@inheritDoc} */
  457.     @Override
  458.     public double getA() {
  459.         return a;
  460.     }

  462.     /** {@inheritDoc} */
  463.     @Override
  464.     public double getADot() {
  465.         return aDot;
  466.     }

  468.     /** {@inheritDoc} */
  469.     @Override
  470.     public double getE() {
  471.         return e;
  472.     }

  474.     /** {@inheritDoc} */
  475.     @Override
  476.     public double getEDot() {
  477.         return eDot;
  478.     }

  480.     /** {@inheritDoc} */
  481.     @Override
  482.     public double getI() {
  483.         return i;
  484.     }

  486.     /** {@inheritDoc} */
  487.     @Override
  488.     public double getIDot() {
  489.         return iDot;
  490.     }

  492.     /** Get the perigee argument.
  493.      * @return perigee argument (rad)
  494.      */
  495.     public double getPerigeeArgument() {
  496.         return pa;
  497.     }

  499.     /** Get the perigee argument derivative.
  500.      * <p>
  501.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  502.      * </p>
  503.      * @return perigee argument derivative (rad/s)
  504.      * @since 9.0
  505.      */
  506.     public double getPerigeeArgumentDot() {
  507.         return paDot;
  508.     }

  510.     /** Get the right ascension of the ascending node.
  511.      * @return right ascension of the ascending node (rad)
  512.      */
  513.     public double getRightAscensionOfAscendingNode() {
  514.         return raan;
  515.     }

  517.     /** Get the right ascension of the ascending node derivative.
  518.      * <p>
  519.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  520.      * </p>
  521.      * @return right ascension of the ascending node derivative (rad/s)
  522.      * @since 9.0
  523.      */
  524.     public double getRightAscensionOfAscendingNodeDot() {
  525.         return raanDot;
  526.     }

  528.     /** Get the true anomaly.
  529.      * @return true anomaly (rad)
  530.      */
  531.     public double getTrueAnomaly() {
  532.         switch (cachedPositionAngleType) {
  533.             case MEAN: return (a < 0.) ? KeplerianAnomalyUtility.hyperbolicMeanToTrue(e, cachedAnomaly) :
  534.                     KeplerianAnomalyUtility.ellipticMeanToTrue(e, cachedAnomaly);

  536.             case TRUE: return cachedAnomaly;

  538.             case ECCENTRIC: return (a < 0.) ? KeplerianAnomalyUtility.hyperbolicEccentricToTrue(e, cachedAnomaly) :
  539.                     KeplerianAnomalyUtility.ellipticEccentricToTrue(e, cachedAnomaly);

  541.             default: throw new OrekitInternalError(null);
  542.         }
  543.     }

  545.     /** Get the true anomaly derivative.
  546.      * @return true anomaly derivative (rad/s)
  547.      */
  548.     public double getTrueAnomalyDot() {
  549.         switch (cachedPositionAngleType) {
  550.             case MEAN:
  551.                 final UnivariateDerivative1 eUD = new UnivariateDerivative1(e, eDot);
  552.                 final UnivariateDerivative1 MUD = new UnivariateDerivative1(cachedAnomaly, cachedAnomalyDot);
  553.                 final UnivariateDerivative1 vUD = (a < 0) ?
  554.                         FieldKeplerianAnomalyUtility.hyperbolicMeanToTrue(eUD, MUD) :
  555.                         FieldKeplerianAnomalyUtility.ellipticMeanToTrue(eUD, MUD);
  556.                 return vUD.getFirstDerivative();

  558.             case TRUE:
  559.                 return cachedAnomalyDot;

  561.             case ECCENTRIC:
  562.                 final UnivariateDerivative1 eUD2 = new UnivariateDerivative1(e, eDot);
  563.                 final UnivariateDerivative1 EUD = new UnivariateDerivative1(cachedAnomaly, cachedAnomalyDot);
  564.                 final UnivariateDerivative1 vUD2 = (a < 0) ?
  565.                         FieldKeplerianAnomalyUtility.hyperbolicEccentricToTrue(eUD2, EUD) :
  566.                         FieldKeplerianAnomalyUtility.ellipticEccentricToTrue(eUD2, EUD);
  567.                 return vUD2.getFirstDerivative();

  569.             default:
  570.                 throw new OrekitInternalError(null);
  571.         }
  572.     }

  574.     /** Get the eccentric anomaly.
  575.      * @return eccentric anomaly (rad)
  576.      */
  577.     public double getEccentricAnomaly() {
  578.         switch (cachedPositionAngleType) {
  579.             case MEAN:
  580.                 return (a < 0.) ? KeplerianAnomalyUtility.hyperbolicMeanToEccentric(e, cachedAnomaly) :
  581.                     KeplerianAnomalyUtility.ellipticMeanToEccentric(e, cachedAnomaly);

  583.             case ECCENTRIC:
  584.                 return cachedAnomaly;

  586.             case TRUE:
  587.                 return (a < 0.) ? KeplerianAnomalyUtility.hyperbolicTrueToEccentric(e, cachedAnomaly) :
  588.                     KeplerianAnomalyUtility.ellipticTrueToEccentric(e, cachedAnomaly);

  590.             default:
  591.                 throw new OrekitInternalError(null);
  592.         }
  593.     }

  595.     /** Get the eccentric anomaly derivative.
  596.      * @return eccentric anomaly derivative (rad/s)
  597.      * @since 9.0
  598.      */
  599.     public double getEccentricAnomalyDot() {
  600.         switch (cachedPositionAngleType) {
  601.             case ECCENTRIC:
  602.                 return cachedAnomalyDot;

  604.             case TRUE:
  605.                 final UnivariateDerivative1 eUD = new UnivariateDerivative1(e, eDot);
  606.                 final UnivariateDerivative1 vUD = new UnivariateDerivative1(cachedAnomaly, cachedAnomalyDot);
  607.                 final UnivariateDerivative1 EUD = (a < 0) ?
  608.                         FieldKeplerianAnomalyUtility.hyperbolicTrueToEccentric(eUD, vUD) :
  609.                         FieldKeplerianAnomalyUtility.ellipticTrueToEccentric(eUD, vUD);
  610.                 return EUD.getFirstDerivative();

  612.             case MEAN:
  613.                 final UnivariateDerivative1 eUD2 = new UnivariateDerivative1(e, eDot);
  614.                 final UnivariateDerivative1 MUD = new UnivariateDerivative1(cachedAnomaly, cachedAnomalyDot);
  615.                 final UnivariateDerivative1 EUD2 = (a < 0) ?
  616.                         FieldKeplerianAnomalyUtility.hyperbolicMeanToEccentric(eUD2, MUD) :
  617.                         FieldKeplerianAnomalyUtility.ellipticMeanToEccentric(eUD2, MUD);
  618.                 return EUD2.getFirstDerivative();

  620.             default:
  621.                 throw new OrekitInternalError(null);
  622.         }
  623.     }

  625.     /** Get the mean anomaly.
  626.      * @return mean anomaly (rad)
  627.      */
  628.     public double getMeanAnomaly() {
  629.         switch (cachedPositionAngleType) {
  630.             case ECCENTRIC: return (a < 0.) ? KeplerianAnomalyUtility.hyperbolicEccentricToMean(e, cachedAnomaly) :
  631.                     KeplerianAnomalyUtility.ellipticEccentricToMean(e, cachedAnomaly);

  633.             case MEAN: return cachedAnomaly;

  635.             case TRUE: return (a < 0.) ? KeplerianAnomalyUtility.hyperbolicTrueToMean(e, cachedAnomaly) :
  636.                     KeplerianAnomalyUtility.ellipticTrueToMean(e, cachedAnomaly);

  638.             default: throw new OrekitInternalError(null);
  639.         }
  640.     }

  642.     /** Get the mean anomaly derivative.
  643.      * @return mean anomaly derivative (rad/s)
  644.      * @since 9.0
  645.      */
  646.     public double getMeanAnomalyDot() {
  647.         switch (cachedPositionAngleType) {
  648.             case MEAN:
  649.                 return cachedAnomalyDot;

  651.             case ECCENTRIC:
  652.                 final UnivariateDerivative1 eUD = new UnivariateDerivative1(e, eDot);
  653.                 final UnivariateDerivative1 EUD = new UnivariateDerivative1(cachedAnomaly, cachedAnomalyDot);
  654.                 final UnivariateDerivative1 MUD = (a < 0) ?
  655.                         FieldKeplerianAnomalyUtility.hyperbolicEccentricToMean(eUD, EUD) :
  656.                         FieldKeplerianAnomalyUtility.ellipticEccentricToMean(eUD, EUD);
  657.                 return MUD.getFirstDerivative();

  659.             case TRUE:
  660.                 final UnivariateDerivative1 eUD2 = new UnivariateDerivative1(e, eDot);
  661.                 final UnivariateDerivative1 vUD = new UnivariateDerivative1(cachedAnomaly, cachedAnomalyDot);
  662.                 final UnivariateDerivative1 MUD2 = (a < 0) ?
  663.                         FieldKeplerianAnomalyUtility.hyperbolicTrueToMean(eUD2, vUD) :
  664.                         FieldKeplerianAnomalyUtility.ellipticTrueToMean(eUD2, vUD);
  665.                 return MUD2.getFirstDerivative();

  667.             default:
  668.                 throw new OrekitInternalError(null);
  669.         }
  670.     }

  672.     /** Get the anomaly.
  673.      * @param type type of the angle
  674.      * @return anomaly (rad)
  675.      */
  676.     public double getAnomaly(final PositionAngleType type) {
  677.         return (type == PositionAngleType.MEAN) ? getMeanAnomaly() :
  678.                 ((type == PositionAngleType.ECCENTRIC) ? getEccentricAnomaly() :
  679.                         getTrueAnomaly());
  680.     }

  682.     /** Get the anomaly derivative.
  683.      * @param type type of the angle
  684.      * @return anomaly derivative (rad/s)
  685.      * @since 9.0
  686.      */
  687.     public double getAnomalyDot(final PositionAngleType type) {
  688.         return (type == PositionAngleType.MEAN) ? getMeanAnomalyDot() :
  689.                 ((type == PositionAngleType.ECCENTRIC) ? getEccentricAnomalyDot() :
  690.                         getTrueAnomalyDot());
  691.     }

  693.     /** {@inheritDoc} */
  694.     @Override
  695.     public double getEquinoctialEx() {
  696.         return e * FastMath.cos(pa + raan);
  697.     }

  699.     /** {@inheritDoc} */
  700.     @Override
  701.     public double getEquinoctialExDot() {
  702.         if (!hasNonKeplerianAcceleration()) {
  703.             return 0.;
  704.         }
  705.         final double paPraan = pa + raan;
  706.         final SinCos sc      = FastMath.sinCos(paPraan);
  707.         return eDot * sc.cos() - e * sc.sin() * (paDot + raanDot);
  708.     }

  710.     /** {@inheritDoc} */
  711.     @Override
  712.     public double getEquinoctialEy() {
  713.         return e * FastMath.sin(pa + raan);
  714.     }

  716.     /** {@inheritDoc} */
  717.     @Override
  718.     public double getEquinoctialEyDot() {
  719.         if (!hasNonKeplerianAcceleration()) {
  720.             return 0.;
  721.         }
  722.         final double paPraan = pa + raan;
  723.         final SinCos sc      = FastMath.sinCos(paPraan);
  724.         return eDot * sc.sin() + e * sc.cos() * (paDot + raanDot);
  725.     }

  727.     /** {@inheritDoc} */
  728.     @Override
  729.     public double getHx() {
  730.         // Check for equatorial retrograde orbit
  731.         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
  732.             return Double.NaN;
  733.         }
  734.         return FastMath.cos(raan) * FastMath.tan(0.5 * i);
  735.     }

  737.     /** {@inheritDoc} */
  738.     @Override
  739.     public double getHxDot() {
  740.         // Check for equatorial retrograde orbit
  741.         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
  742.             return Double.NaN;
  743.         }
  744.         if (!hasNonKeplerianAcceleration()) {
  745.             return 0.;
  746.         }
  747.         final SinCos sc      = FastMath.sinCos(raan);
  748.         final double tan     = FastMath.tan(0.5 * i);
  749.         return 0.5 * (1 + tan * tan) * sc.cos() * iDot - tan * sc.sin() * raanDot;
  750.     }

  752.     /** {@inheritDoc} */
  753.     @Override
  754.     public double getHy() {
  755.         // Check for equatorial retrograde orbit
  756.         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
  757.             return Double.NaN;
  758.         }
  759.         return FastMath.sin(raan) * FastMath.tan(0.5 * i);
  760.     }

  762.     /** {@inheritDoc} */
  763.     @Override
  764.     public double getHyDot() {
  765.         // Check for equatorial retrograde orbit
  766.         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
  767.             return Double.NaN;
  768.         }
  769.         if (!hasNonKeplerianAcceleration()) {
  770.             return 0.;
  771.         }
  772.         final SinCos sc      = FastMath.sinCos(raan);
  773.         final double tan     = FastMath.tan(0.5 * i);
  774.         return 0.5 * (1 + tan * tan) * sc.sin() * iDot + tan * sc.cos() * raanDot;
  775.     }

  777.     /** {@inheritDoc} */
  778.     @Override
  779.     public double getLv() {
  780.         return pa + raan + getTrueAnomaly();
  781.     }

  783.     /** {@inheritDoc} */
  784.     @Override
  785.     public double getLvDot() {
  786.         return paDot + raanDot + getTrueAnomalyDot();
  787.     }

  789.     /** {@inheritDoc} */
  790.     @Override
  791.     public double getLE() {
  792.         return pa + raan + getEccentricAnomaly();
  793.     }

  795.     /** {@inheritDoc} */
  796.     @Override
  797.     public double getLEDot() {
  798.         return paDot + raanDot + getEccentricAnomalyDot();
  799.     }

  801.     /** {@inheritDoc} */
  802.     @Override
  803.     public double getLM() {
  804.         return pa + raan + getMeanAnomaly();
  805.     }

  807.     /** {@inheritDoc} */
  808.     @Override
  809.     public double getLMDot() {
  810.         return paDot + raanDot + getMeanAnomalyDot();
  811.     }

  813.     /** Initialize cached anomaly with rate.
  814.      * @param anomaly input anomaly
  815.      * @param anomalyDot rate of input anomaly
  816.      * @param inputType position angle type passed as input
  817.      * @return anomaly to cache with rate
  818.      * @since 12.1
  819.      */
  820.     private UnivariateDerivative1 initializeCachedAnomaly(final double anomaly, final double anomalyDot,
  821.                                                           final PositionAngleType inputType) {
  822.         if (cachedPositionAngleType == inputType) {
  823.             return new UnivariateDerivative1(anomaly, anomalyDot);

  825.         } else {
  826.             final UnivariateDerivative1 eUD = new UnivariateDerivative1(e, eDot);
  827.             final UnivariateDerivative1 anomalyUD = new UnivariateDerivative1(anomaly, anomalyDot);

  829.             if (a < 0) {
  830.                 switch (cachedPositionAngleType) {
  831.                     case MEAN:
  832.                         if (inputType == PositionAngleType.ECCENTRIC) {
  833.                             return FieldKeplerianAnomalyUtility.hyperbolicEccentricToMean(eUD, anomalyUD);
  834.                         } else {
  835.                             return FieldKeplerianAnomalyUtility.hyperbolicTrueToMean(eUD, anomalyUD);
  836.                         }

  838.                     case ECCENTRIC:
  839.                         if (inputType == PositionAngleType.MEAN) {
  840.                             return FieldKeplerianAnomalyUtility.hyperbolicMeanToEccentric(eUD, anomalyUD);
  841.                         } else {
  842.                             return FieldKeplerianAnomalyUtility.hyperbolicTrueToEccentric(eUD, anomalyUD);
  843.                         }

  845.                     case TRUE:
  846.                         if (inputType == PositionAngleType.MEAN) {
  847.                             return FieldKeplerianAnomalyUtility.hyperbolicMeanToTrue(eUD, anomalyUD);
  848.                         } else {
  849.                             return FieldKeplerianAnomalyUtility.hyperbolicEccentricToTrue(eUD, anomalyUD);
  850.                         }

  852.                     default:
  853.                         break;
  854.                 }

  856.             } else {
  857.                 switch (cachedPositionAngleType) {
  858.                     case MEAN:
  859.                         if (inputType == PositionAngleType.ECCENTRIC) {
  860.                             return FieldKeplerianAnomalyUtility.ellipticEccentricToMean(eUD, anomalyUD);
  861.                         } else {
  862.                             return FieldKeplerianAnomalyUtility.ellipticTrueToMean(eUD, anomalyUD);
  863.                         }

  865.                     case ECCENTRIC:
  866.                         if (inputType == PositionAngleType.MEAN) {
  867.                             return FieldKeplerianAnomalyUtility.ellipticMeanToEccentric(eUD, anomalyUD);
  868.                         } else {
  869.                             return FieldKeplerianAnomalyUtility.ellipticTrueToEccentric(eUD, anomalyUD);
  870.                         }

  872.                     case TRUE:
  873.                         if (inputType == PositionAngleType.MEAN) {
  874.                             return FieldKeplerianAnomalyUtility.ellipticMeanToTrue(eUD, anomalyUD);
  875.                         } else {
  876.                             return FieldKeplerianAnomalyUtility.ellipticEccentricToTrue(eUD, anomalyUD);
  877.                         }

  879.                     default:
  880.                         break;
  881.                 }

  883.             }
  884.             throw new OrekitInternalError(null);
  885.         }

  887.     }

  889.     /** Compute reference axes.
  890.      * @return reference axes
  891.      * @since 12.0
  892.      */
  893.     private Vector3D[] referenceAxes() {
  894.         // preliminary variables
  895.         final SinCos scRaan  = FastMath.sinCos(raan);
  896.         final SinCos scPa    = FastMath.sinCos(pa);
  897.         final SinCos scI     = FastMath.sinCos(i);
  898.         final double cosRaan = scRaan.cos();
  899.         final double sinRaan = scRaan.sin();
  900.         final double cosPa   = scPa.cos();
  901.         final double sinPa   = scPa.sin();
  902.         final double cosI    = scI.cos();
  903.         final double sinI    = scI.sin();

  905.         final double crcp    = cosRaan * cosPa;
  906.         final double crsp    = cosRaan * sinPa;
  907.         final double srcp    = sinRaan * cosPa;
  908.         final double srsp    = sinRaan * sinPa;

  910.         // reference axes defining the orbital plane
  911.         return new Vector3D[] {
  912.             new Vector3D( crcp - cosI * srsp,  srcp + cosI * crsp, sinI * sinPa),
  913.             new Vector3D(-crsp - cosI * srcp, -srsp + cosI * crcp, sinI * cosPa)
  914.         };

  916.     }

  918.     /** Compute position and velocity but not acceleration.
  919.      */
  920.     private void computePVWithoutA() {

  922.         if (partialPV != null) {
  923.             // already computed
  924.             return;
  925.         }

  927.         final Vector3D[] axes = referenceAxes();

  929.         if (isElliptical()) {

  931.             // elliptical case

  933.             // elliptic eccentric anomaly
  934.             final double uME2   = (1 - e) * (1 + e);
  935.             final double s1Me2  = FastMath.sqrt(uME2);
  936.             final SinCos scE    = FastMath.sinCos(getEccentricAnomaly());
  937.             final double cosE   = scE.cos();
  938.             final double sinE   = scE.sin();

  940.             // coordinates of position and velocity in the orbital plane
  941.             final double x      = a * (cosE - e);
  942.             final double y      = a * sinE * s1Me2;
  943.             final double factor = FastMath.sqrt(getMu() / a) / (1 - e * cosE);
  944.             final double xDot   = -sinE * factor;
  945.             final double yDot   =  cosE * s1Me2 * factor;

  947.             final Vector3D position = new Vector3D(x, axes[0], y, axes[1]);
  948.             final Vector3D velocity = new Vector3D(xDot, axes[0], yDot, axes[1]);
  949.             partialPV = new PVCoordinates(position, velocity);

  951.         } else {

  953.             // hyperbolic case

  955.             // compute position and velocity factors
  956.             final SinCos scV       = FastMath.sinCos(getTrueAnomaly());
  957.             final double sinV      = scV.sin();
  958.             final double cosV      = scV.cos();
  959.             final double f         = a * (1 - e * e);
  960.             final double posFactor = f / (1 + e * cosV);
  961.             final double velFactor = FastMath.sqrt(getMu() / f);

  963.             final double   x            =  posFactor * cosV;
  964.             final double   y            =  posFactor * sinV;
  965.             final double   xDot         = -velFactor * sinV;
  966.             final double   yDot         =  velFactor * (e + cosV);

  968.             final Vector3D position = new Vector3D(x, axes[0], y, axes[1]);
  969.             final Vector3D velocity = new Vector3D(xDot, axes[0], yDot, axes[1]);
  970.             partialPV = new PVCoordinates(position, velocity);

  972.         }

  974.     }

  976.     /** Compute non-Keplerian part of the acceleration from first time derivatives.
  977.      * @return non-Keplerian part of the acceleration
  978.      */
  979.     private Vector3D nonKeplerianAcceleration() {

  981.         final double[][] dCdP = new double[6][6];
  982.         getJacobianWrtParameters(PositionAngleType.MEAN, dCdP);

  984.         final double nonKeplerianMeanMotion = getMeanAnomalyDot() - getKeplerianMeanMotion();
  985.         final double nonKeplerianAx = dCdP[3][0] * aDot    + dCdP[3][1] * eDot    + dCdP[3][2] * iDot    +
  986.                 dCdP[3][3] * paDot   + dCdP[3][4] * raanDot + dCdP[3][5] * nonKeplerianMeanMotion;
  987.         final double nonKeplerianAy = dCdP[4][0] * aDot    + dCdP[4][1] * eDot    + dCdP[4][2] * iDot    +
  988.                 dCdP[4][3] * paDot   + dCdP[4][4] * raanDot + dCdP[4][5] * nonKeplerianMeanMotion;
  989.         final double nonKeplerianAz = dCdP[5][0] * aDot    + dCdP[5][1] * eDot    + dCdP[5][2] * iDot    +
  990.                 dCdP[5][3] * paDot   + dCdP[5][4] * raanDot + dCdP[5][5] * nonKeplerianMeanMotion;

  992.         return new Vector3D(nonKeplerianAx, nonKeplerianAy, nonKeplerianAz);

  994.     }

  996.     /** {@inheritDoc} */
  997.     @Override
  998.     protected Vector3D initPosition() {

  1000.         final Vector3D[] axes = referenceAxes();

  1002.         if (isElliptical()) {

  1004.             // elliptical case

  1006.             // elliptic eccentric anomaly
  1007.             final double uME2   = (1 - e) * (1 + e);
  1008.             final double s1Me2  = FastMath.sqrt(uME2);
  1009.             final SinCos scE    = FastMath.sinCos(getEccentricAnomaly());
  1010.             final double cosE   = scE.cos();
  1011.             final double sinE   = scE.sin();

  1013.             return new Vector3D(a * (cosE - e), axes[0], a * sinE * s1Me2, axes[1]);

  1015.         } else {

  1017.             // hyperbolic case

  1019.             // compute position and velocity factors
  1020.             final SinCos scV       = FastMath.sinCos(getTrueAnomaly());
  1021.             final double sinV      = scV.sin();
  1022.             final double cosV      = scV.cos();
  1023.             final double f         = a * (1 - e * e);
  1024.             final double posFactor = f / (1 + e * cosV);

  1026.             return new Vector3D(posFactor * cosV, axes[0], posFactor * sinV, axes[1]);

  1028.         }

  1030.     }

  1032.     /** {@inheritDoc} */
  1033.     @Override
  1034.     protected TimeStampedPVCoordinates initPVCoordinates() {

  1036.         // position and velocity
  1037.         computePVWithoutA();

  1039.         // acceleration
  1040.         final double r2 = partialPV.getPosition().getNormSq();
  1041.         final Vector3D keplerianAcceleration = new Vector3D(-getMu() / (r2 * FastMath.sqrt(r2)), partialPV.getPosition());
  1042.         final Vector3D acceleration = hasNonKeplerianAcceleration() ?
  1043.                 keplerianAcceleration.add(nonKeplerianAcceleration()) :
  1044.                 keplerianAcceleration;

  1046.         return new TimeStampedPVCoordinates(getDate(), partialPV.getPosition(), partialPV.getVelocity(), acceleration);

  1048.     }

  1050.     /** {@inheritDoc} */
  1051.     @Override
  1052.     public KeplerianOrbit inFrame(final Frame inertialFrame) {
  1053.         final PVCoordinates pvCoordinates;
  1054.         if (hasNonKeplerianAcceleration()) {
  1055.             pvCoordinates = getPVCoordinates(inertialFrame);
  1056.         } else {
  1057.             final KinematicTransform transform = getFrame().getKinematicTransformTo(inertialFrame, getDate());
  1058.             pvCoordinates = transform.transformOnlyPV(getPVCoordinates());
  1059.         }
  1060.         final KeplerianOrbit keplerianOrbit = new KeplerianOrbit(pvCoordinates, inertialFrame, getDate(), getMu());
  1061.         if (keplerianOrbit.getCachedPositionAngleType() == getCachedPositionAngleType()) {
  1062.             return keplerianOrbit;
  1063.         } else {
  1064.             return keplerianOrbit.withCachedPositionAngleType(getCachedPositionAngleType());
  1065.         }
  1066.     }

  1068.     /** {@inheritDoc} */
  1069.     @Override
  1070.     public KeplerianOrbit withCachedPositionAngleType(final PositionAngleType positionAngleType) {
  1071.         return new KeplerianOrbit(a, e, i, pa, raan, getAnomaly(positionAngleType), aDot, eDot, iDot, paDot, raanDot,
  1072.                 getAnomalyDot(positionAngleType), positionAngleType, getFrame(), getDate(), getMu());
  1073.     }

  1075.     /** {@inheritDoc} */
  1076.     @Override
  1077.     public KeplerianOrbit shiftedBy(final double dt) {
  1078.         return shiftedBy(new TimeOffset(dt));
  1079.     }

  1081.     /** {@inheritDoc} */
  1082.     @Override
  1083.     public KeplerianOrbit shiftedBy(final TimeOffset dt) {

  1085.         final double dtS = dt.toDouble();

  1087.         // use Keplerian-only motion
  1088.         final KeplerianOrbit keplerianShifted = new KeplerianOrbit(a, e, i, pa, raan,
  1089.                 getMeanAnomaly() + getKeplerianMeanMotion() * dtS, PositionAngleType.MEAN,
  1090.                 cachedPositionAngleType, getFrame(), getDate().shiftedBy(dt), getMu());

  1092.         if (hasNonKeplerianAcceleration()) {

  1094.             // extract non-Keplerian acceleration from first time derivatives
  1095.             final Vector3D nonKeplerianAcceleration = nonKeplerianAcceleration();

  1097.             // add quadratic effect of non-Keplerian acceleration to Keplerian-only shift
  1098.             keplerianShifted.computePVWithoutA();
  1099.             final Vector3D fixedP   = new Vector3D(1, keplerianShifted.partialPV.getPosition(),
  1100.                     0.5 * dtS * dtS, nonKeplerianAcceleration);
  1101.             final double   fixedR2 = fixedP.getNormSq();
  1102.             final double   fixedR  = FastMath.sqrt(fixedR2);
  1103.             final Vector3D fixedV  = new Vector3D(1, keplerianShifted.partialPV.getVelocity(),
  1104.                     dtS, nonKeplerianAcceleration);
  1105.             final Vector3D fixedA  = new Vector3D(-getMu() / (fixedR2 * fixedR), keplerianShifted.partialPV.getPosition(),
  1106.                     1, nonKeplerianAcceleration);

  1108.             // build a new orbit, taking non-Keplerian acceleration into account
  1109.             return new KeplerianOrbit(new TimeStampedPVCoordinates(keplerianShifted.getDate(),
  1110.                     fixedP, fixedV, fixedA),
  1111.                     keplerianShifted.getFrame(), keplerianShifted.getMu());

  1113.         } else {
  1114.             // Keplerian-only motion is all we can do
  1115.             return keplerianShifted;
  1116.         }

  1118.     }

  1120.     /** {@inheritDoc} */
  1121.     @Override
  1122.     protected double[][] computeJacobianMeanWrtCartesian() {
  1123.         if (isElliptical()) {
  1124.             return computeJacobianMeanWrtCartesianElliptical();
  1125.         } else {
  1126.             return computeJacobianMeanWrtCartesianHyperbolic();
  1127.         }
  1128.     }

  1130.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  1131.      * <p>
  1132.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  1133.      * respect to Cartesian coordinate j (x for j=0, y for j=1, z for j=2, xDot for j=3,
  1134.      * yDot for j=4, zDot for j=5).
  1135.      * </p>
  1136.      * @return 6x6 Jacobian matrix
  1137.      */
  1138.     private double[][] computeJacobianMeanWrtCartesianElliptical() {

  1140.         final double[][] jacobian = new double[6][6];

  1142.         // compute various intermediate parameters
  1143.         computePVWithoutA();
  1144.         final Vector3D position = partialPV.getPosition();
  1145.         final Vector3D velocity = partialPV.getVelocity();
  1146.         final Vector3D momentum = partialPV.getMomentum();
  1147.         final double v2         = velocity.getNormSq();
  1148.         final double r2         = position.getNormSq();
  1149.         final double r          = FastMath.sqrt(r2);
  1150.         final double r3         = r * r2;

  1152.         final double px         = position.getX();
  1153.         final double py         = position.getY();
  1154.         final double pz         = position.getZ();
  1155.         final double vx         = velocity.getX();
  1156.         final double vy         = velocity.getY();
  1157.         final double vz         = velocity.getZ();
  1158.         final double mx         = momentum.getX();
  1159.         final double my         = momentum.getY();
  1160.         final double mz         = momentum.getZ();

  1162.         final double mu         = getMu();
  1163.         final double sqrtMuA    = FastMath.sqrt(a * mu);
  1164.         final double sqrtAoMu   = FastMath.sqrt(a / mu);
  1165.         final double a2         = a * a;
  1166.         final double twoA       = 2 * a;
  1167.         final double rOnA       = r / a;

  1169.         final double oMe2       = 1 - e * e;
  1170.         final double epsilon    = FastMath.sqrt(oMe2);
  1171.         final double sqrtRec    = 1 / epsilon;

  1173.         final SinCos scI        = FastMath.sinCos(i);
  1174.         final SinCos scPA       = FastMath.sinCos(pa);
  1175.         final double cosI       = scI.cos();
  1176.         final double sinI       = scI.sin();
  1177.         final double cosPA      = scPA.cos();
  1178.         final double sinPA      = scPA.sin();

  1180.         final double pv         = Vector3D.dotProduct(position, velocity);
  1181.         final double cosE       = (a - r) / (a * e);
  1182.         final double sinE       = pv / (e * sqrtMuA);

  1184.         // da
  1185.         final Vector3D vectorAR = new Vector3D(2 * a2 / r3, position);
  1186.         final Vector3D vectorARDot = velocity.scalarMultiply(2 * a2 / mu);
  1187.         fillHalfRow(1, vectorAR,    jacobian[0], 0);
  1188.         fillHalfRow(1, vectorARDot, jacobian[0], 3);

  1190.         // de
  1191.         final double factorER3 = pv / twoA;
  1192.         final Vector3D vectorER   = new Vector3D(cosE * v2 / (r * mu), position,
  1193.                 sinE / sqrtMuA, velocity,
  1194.                 -factorER3 * sinE / sqrtMuA, vectorAR);
  1195.         final Vector3D vectorERDot = new Vector3D(sinE / sqrtMuA, position,
  1196.                 cosE * 2 * r / mu, velocity,
  1197.                 -factorER3 * sinE / sqrtMuA, vectorARDot);
  1198.         fillHalfRow(1, vectorER,    jacobian[1], 0);
  1199.         fillHalfRow(1, vectorERDot, jacobian[1], 3);

  1201.         // dE / dr (Eccentric anomaly)
  1202.         final double coefE = cosE / (e * sqrtMuA);
  1203.         final Vector3D  vectorEAnR =
  1204.                 new Vector3D(-sinE * v2 / (e * r * mu), position, coefE, velocity,
  1205.                         -factorER3 * coefE, vectorAR);

  1207.         // dE / drDot
  1208.         final Vector3D  vectorEAnRDot =
  1209.                 new Vector3D(-sinE * 2 * r / (e * mu), velocity, coefE, position,
  1210.                         -factorER3 * coefE, vectorARDot);

  1212.         // precomputing some more factors
  1213.         final double s1 = -sinE * pz / r - cosE * vz * sqrtAoMu;
  1214.         final double s2 = -cosE * pz / r3;
  1215.         final double s3 = -sinE * vz / (2 * sqrtMuA);
  1216.         final double t1 = sqrtRec * (cosE * pz / r - sinE * vz * sqrtAoMu);
  1217.         final double t2 = sqrtRec * (-sinE * pz / r3);
  1218.         final double t3 = sqrtRec * (cosE - e) * vz / (2 * sqrtMuA);
  1219.         final double t4 = sqrtRec * (e * sinI * cosPA * sqrtRec - vz * sqrtAoMu);
  1220.         final Vector3D s = new Vector3D(cosE / r, Vector3D.PLUS_K,
  1221.                 s1,       vectorEAnR,
  1222.                 s2,       position,
  1223.                 s3,       vectorAR);
  1224.         final Vector3D sDot = new Vector3D(-sinE * sqrtAoMu, Vector3D.PLUS_K,
  1225.                 s1,               vectorEAnRDot,
  1226.                 s3,               vectorARDot);
  1227.         final Vector3D t =
  1228.                 new Vector3D(sqrtRec * sinE / r, Vector3D.PLUS_K).add(new Vector3D(t1, vectorEAnR,
  1229.                         t2, position,
  1230.                         t3, vectorAR,
  1231.                         t4, vectorER));
  1232.         final Vector3D tDot = new Vector3D(sqrtRec * (cosE - e) * sqrtAoMu, Vector3D.PLUS_K,
  1233.                 t1,                              vectorEAnRDot,
  1234.                 t3,                              vectorARDot,
  1235.                 t4,                              vectorERDot);

  1237.         // di
  1238.         final double factorI1 = -sinI * sqrtRec / sqrtMuA;
  1239.         final double i1 =  factorI1;
  1240.         final double i2 = -factorI1 * mz / twoA;
  1241.         final double i3 =  factorI1 * mz * e / oMe2;
  1242.         final double i4 = cosI * sinPA;
  1243.         final double i5 = cosI * cosPA;
  1244.         fillHalfRow(i1, new Vector3D(vy, -vx, 0), i2, vectorAR, i3, vectorER, i4, s, i5, t,
  1245.                 jacobian[2], 0);
  1246.         fillHalfRow(i1, new Vector3D(-py, px, 0), i2, vectorARDot, i3, vectorERDot, i4, sDot, i5, tDot,
  1247.                 jacobian[2], 3);

  1249.         // dpa
  1250.         fillHalfRow(cosPA / sinI, s,    -sinPA / sinI, t,    jacobian[3], 0);
  1251.         fillHalfRow(cosPA / sinI, sDot, -sinPA / sinI, tDot, jacobian[3], 3);

  1253.         // dRaan
  1254.         final double factorRaanR = 1 / (mu * a * oMe2 * sinI * sinI);
  1255.         fillHalfRow(-factorRaanR * my, new Vector3D(  0, vz, -vy),
  1256.                 factorRaanR * mx, new Vector3D(-vz,  0,  vx),
  1257.                 jacobian[4], 0);
  1258.         fillHalfRow(-factorRaanR * my, new Vector3D( 0, -pz,  py),
  1259.                 factorRaanR * mx, new Vector3D(pz,   0, -px),
  1260.                 jacobian[4], 3);

  1262.         // dM
  1263.         fillHalfRow(rOnA, vectorEAnR,    -sinE, vectorER,    jacobian[5], 0);
  1264.         fillHalfRow(rOnA, vectorEAnRDot, -sinE, vectorERDot, jacobian[5], 3);

  1266.         return jacobian;

  1268.     }

  1270.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  1271.      * <p>
  1272.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  1273.      * respect to Cartesian coordinate j (x for j=0, y for j=1, z for j=2, xDot for j=3,
  1274.      * yDot for j=4, zDot for j=5).
  1275.      * </p>
  1276.      * @return 6x6 Jacobian matrix
  1277.      */
  1278.     private double[][] computeJacobianMeanWrtCartesianHyperbolic() {

  1280.         final double[][] jacobian = new double[6][6];

  1282.         // compute various intermediate parameters
  1283.         computePVWithoutA();
  1284.         final Vector3D position = partialPV.getPosition();
  1285.         final Vector3D velocity = partialPV.getVelocity();
  1286.         final Vector3D momentum = partialPV.getMomentum();
  1287.         final double r2         = position.getNormSq();
  1288.         final double r          = FastMath.sqrt(r2);
  1289.         final double r3         = r * r2;

  1291.         final double x          = position.getX();
  1292.         final double y          = position.getY();
  1293.         final double z          = position.getZ();
  1294.         final double vx         = velocity.getX();
  1295.         final double vy         = velocity.getY();
  1296.         final double vz         = velocity.getZ();
  1297.         final double mx         = momentum.getX();
  1298.         final double my         = momentum.getY();
  1299.         final double mz         = momentum.getZ();

  1301.         final double mu         = getMu();
  1302.         final double absA       = -a;
  1303.         final double sqrtMuA    = FastMath.sqrt(absA * mu);
  1304.         final double a2         = a * a;
  1305.         final double rOa        = r / absA;

  1307.         final SinCos scI        = FastMath.sinCos(i);
  1308.         final double cosI       = scI.cos();
  1309.         final double sinI       = scI.sin();

  1311.         final double pv         = Vector3D.dotProduct(position, velocity);

  1313.         // da
  1314.         final Vector3D vectorAR = new Vector3D(-2 * a2 / r3, position);
  1315.         final Vector3D vectorARDot = velocity.scalarMultiply(-2 * a2 / mu);
  1316.         fillHalfRow(-1, vectorAR,    jacobian[0], 0);
  1317.         fillHalfRow(-1, vectorARDot, jacobian[0], 3);

  1319.         // differentials of the momentum
  1320.         final double m      = momentum.getNorm();
  1321.         final double oOm    = 1 / m;
  1322.         final Vector3D dcXP = new Vector3D(  0,  vz, -vy);
  1323.         final Vector3D dcYP = new Vector3D(-vz,   0,  vx);
  1324.         final Vector3D dcZP = new Vector3D( vy, -vx,   0);
  1325.         final Vector3D dcXV = new Vector3D(  0,  -z,   y);
  1326.         final Vector3D dcYV = new Vector3D(  z,   0,  -x);
  1327.         final Vector3D dcZV = new Vector3D( -y,   x,   0);
  1328.         final Vector3D dCP  = new Vector3D(mx * oOm, dcXP, my * oOm, dcYP, mz * oOm, dcZP);
  1329.         final Vector3D dCV  = new Vector3D(mx * oOm, dcXV, my * oOm, dcYV, mz * oOm, dcZV);

  1331.         // dp
  1332.         final double mOMu   = m / mu;
  1333.         final Vector3D dpP  = new Vector3D(2 * mOMu, dCP);
  1334.         final Vector3D dpV  = new Vector3D(2 * mOMu, dCV);

  1336.         // de
  1337.         final double p      = m * mOMu;
  1338.         final double moO2ae = 1 / (2 * absA * e);
  1339.         final double m2OaMu = -p / absA;
  1340.         fillHalfRow(moO2ae, dpP, m2OaMu * moO2ae, vectorAR,    jacobian[1], 0);
  1341.         fillHalfRow(moO2ae, dpV, m2OaMu * moO2ae, vectorARDot, jacobian[1], 3);

  1343.         // di
  1344.         final double cI1 = 1 / (m * sinI);
  1345.         final double cI2 = cosI * cI1;
  1346.         fillHalfRow(cI2, dCP, -cI1, dcZP, jacobian[2], 0);
  1347.         fillHalfRow(cI2, dCV, -cI1, dcZV, jacobian[2], 3);

  1349.         // dPA
  1350.         final double cP1     =  y * oOm;
  1351.         final double cP2     = -x * oOm;
  1352.         final double cP3     = -(mx * cP1 + my * cP2);
  1353.         final double cP4     = cP3 * oOm;
  1354.         final double cP5     = -1 / (r2 * sinI * sinI);
  1355.         final double cP6     = z  * cP5;
  1356.         final double cP7     = cP3 * cP5;
  1357.         final Vector3D dacP  = new Vector3D(cP1, dcXP, cP2, dcYP, cP4, dCP, oOm, new Vector3D(-my, mx, 0));
  1358.         final Vector3D dacV  = new Vector3D(cP1, dcXV, cP2, dcYV, cP4, dCV);
  1359.         final Vector3D dpoP  = new Vector3D(cP6, dacP, cP7, Vector3D.PLUS_K);
  1360.         final Vector3D dpoV  = new Vector3D(cP6, dacV);

  1362.         final double re2     = r2 * e * e;
  1363.         final double recOre2 = (p - r) / re2;
  1364.         final double resOre2 = (pv * mOMu) / re2;
  1365.         final Vector3D dreP  = new Vector3D(mOMu, velocity, pv / mu, dCP);
  1366.         final Vector3D dreV  = new Vector3D(mOMu, position, pv / mu, dCV);
  1367.         final Vector3D davP  = new Vector3D(-resOre2, dpP, recOre2, dreP, resOre2 / r, position);
  1368.         final Vector3D davV  = new Vector3D(-resOre2, dpV, recOre2, dreV);
  1369.         fillHalfRow(1, dpoP, -1, davP, jacobian[3], 0);
  1370.         fillHalfRow(1, dpoV, -1, davV, jacobian[3], 3);

  1372.         // dRAAN
  1373.         final double cO0 = cI1 * cI1;
  1374.         final double cO1 =  mx * cO0;
  1375.         final double cO2 = -my * cO0;
  1376.         fillHalfRow(cO1, dcYP, cO2, dcXP, jacobian[4], 0);
  1377.         fillHalfRow(cO1, dcYV, cO2, dcXV, jacobian[4], 3);

  1379.         // dM
  1380.         final double s2a    = pv / (2 * absA);
  1381.         final double oObux  = 1 / FastMath.sqrt(m * m + mu * absA);
  1382.         final double scasbu = pv * oObux;
  1383.         final Vector3D dauP = new Vector3D(1 / sqrtMuA, velocity, -s2a / sqrtMuA, vectorAR);
  1384.         final Vector3D dauV = new Vector3D(1 / sqrtMuA, position, -s2a / sqrtMuA, vectorARDot);
  1385.         final Vector3D dbuP = new Vector3D(oObux * mu / 2, vectorAR,    m * oObux, dCP);
  1386.         final Vector3D dbuV = new Vector3D(oObux * mu / 2, vectorARDot, m * oObux, dCV);
  1387.         final Vector3D dcuP = new Vector3D(oObux, velocity, -scasbu * oObux, dbuP);
  1388.         final Vector3D dcuV = new Vector3D(oObux, position, -scasbu * oObux, dbuV);
  1389.         fillHalfRow(1, dauP, -e / (1 + rOa), dcuP, jacobian[5], 0);
  1390.         fillHalfRow(1, dauV, -e / (1 + rOa), dcuV, jacobian[5], 3);

  1392.         return jacobian;

  1394.     }

  1396.     /** {@inheritDoc} */
  1397.     @Override
  1398.     protected double[][] computeJacobianEccentricWrtCartesian() {
  1399.         if (isElliptical()) {
  1400.             return computeJacobianEccentricWrtCartesianElliptical();
  1401.         } else {
  1402.             return computeJacobianEccentricWrtCartesianHyperbolic();
  1403.         }
  1404.     }

  1406.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  1407.      * <p>
  1408.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  1409.      * respect to Cartesian coordinate j (x for j=0, y for j=1, z for j=2, xDot for j=3,
  1410.      * yDot for j=4, zDot for j=5).
  1411.      * </p>
  1412.      * @return 6x6 Jacobian matrix
  1413.      */
  1414.     private double[][] computeJacobianEccentricWrtCartesianElliptical() {

  1416.         // start by computing the Jacobian with mean angle
  1417.         final double[][] jacobian = computeJacobianMeanWrtCartesianElliptical();

  1419.         // Differentiating the Kepler equation M = E - e sin E leads to:
  1420.         // dM = (1 - e cos E) dE - sin E de
  1421.         // which is inverted and rewritten as:
  1422.         // dE = a/r dM + sin E a/r de
  1423.         final SinCos scE              = FastMath.sinCos(getEccentricAnomaly());
  1424.         final double aOr              = 1 / (1 - e * scE.cos());

  1426.         // update anomaly row
  1427.         final double[] eRow           = jacobian[1];
  1428.         final double[] anomalyRow     = jacobian[5];
  1429.         for (int j = 0; j < anomalyRow.length; ++j) {
  1430.             anomalyRow[j] = aOr * (anomalyRow[j] + scE.sin() * eRow[j]);
  1431.         }

  1433.         return jacobian;

  1435.     }

  1437.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  1438.      * <p>
  1439.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  1440.      * respect to Cartesian coordinate j (x for j=0, y for j=1, z for j=2, xDot for j=3,
  1441.      * yDot for j=4, zDot for j=5).
  1442.      * </p>
  1443.      * @return 6x6 Jacobian matrix
  1444.      */
  1445.     private double[][] computeJacobianEccentricWrtCartesianHyperbolic() {

  1447.         // start by computing the Jacobian with mean angle
  1448.         final double[][] jacobian = computeJacobianMeanWrtCartesianHyperbolic();

  1450.         // Differentiating the Kepler equation M = e sinh H - H leads to:
  1451.         // dM = (e cosh H - 1) dH + sinh H de
  1452.         // which is inverted and rewritten as:
  1453.         // dH = 1 / (e cosh H - 1) dM - sinh H / (e cosh H - 1) de
  1454.         final double H      = getEccentricAnomaly();
  1455.         final double coshH  = FastMath.cosh(H);
  1456.         final double sinhH  = FastMath.sinh(H);
  1457.         final double absaOr = 1 / (e * coshH - 1);

  1459.         // update anomaly row
  1460.         final double[] eRow       = jacobian[1];
  1461.         final double[] anomalyRow = jacobian[5];
  1462.         for (int j = 0; j < anomalyRow.length; ++j) {
  1463.             anomalyRow[j] = absaOr * (anomalyRow[j] - sinhH * eRow[j]);
  1464.         }

  1466.         return jacobian;

  1468.     }

  1470.     /** {@inheritDoc} */
  1471.     @Override
  1472.     protected double[][] computeJacobianTrueWrtCartesian() {
  1473.         if (isElliptical()) {
  1474.             return computeJacobianTrueWrtCartesianElliptical();
  1475.         } else {
  1476.             return computeJacobianTrueWrtCartesianHyperbolic();
  1477.         }
  1478.     }

  1480.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  1481.      * <p>
  1482.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  1483.      * respect to Cartesian coordinate j (x for j=0, y for j=1, z for j=2, xDot for j=3,
  1484.      * yDot for j=4, zDot for j=5).
  1485.      * </p>
  1486.      * @return 6x6 Jacobian matrix
  1487.      */
  1488.     private double[][] computeJacobianTrueWrtCartesianElliptical() {

  1490.         // start by computing the Jacobian with eccentric angle
  1491.         final double[][] jacobian = computeJacobianEccentricWrtCartesianElliptical();

  1493.         // Differentiating the eccentric anomaly equation sin E = sqrt(1-e^2) sin v / (1 + e cos v)
  1494.         // and using cos E = (e + cos v) / (1 + e cos v) to get rid of cos E leads to:
  1495.         // dE = [sqrt (1 - e^2) / (1 + e cos v)] dv - [sin E / (1 - e^2)] de
  1496.         // which is inverted and rewritten as:
  1497.         // dv = sqrt (1 - e^2) a/r dE + [sin E / sqrt (1 - e^2)] a/r de
  1498.         final double e2           = e * e;
  1499.         final double oMe2         = 1 - e2;
  1500.         final double epsilon      = FastMath.sqrt(oMe2);
  1501.         final SinCos scE          = FastMath.sinCos(getEccentricAnomaly());
  1502.         final double aOr          = 1 / (1 - e * scE.cos());
  1503.         final double aFactor      = epsilon * aOr;
  1504.         final double eFactor      = scE.sin() * aOr / epsilon;

  1506.         // update anomaly row
  1507.         final double[] eRow       = jacobian[1];
  1508.         final double[] anomalyRow = jacobian[5];
  1509.         for (int j = 0; j < anomalyRow.length; ++j) {
  1510.             anomalyRow[j] = aFactor * anomalyRow[j] + eFactor * eRow[j];
  1511.         }

  1513.         return jacobian;

  1515.     }

  1517.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  1518.      * <p>
  1519.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  1520.      * respect to Cartesian coordinate j (x for j=0, y for j=1, z for j=2, xDot for j=3,
  1521.      * yDot for j=4, zDot for j=5).
  1522.      * </p>
  1523.      * @return 6x6 Jacobian matrix
  1524.      */
  1525.     private double[][] computeJacobianTrueWrtCartesianHyperbolic() {

  1527.         // start by computing the Jacobian with eccentric angle
  1528.         final double[][] jacobian = computeJacobianEccentricWrtCartesianHyperbolic();

  1530.         // Differentiating the eccentric anomaly equation sinh H = sqrt(e^2-1) sin v / (1 + e cos v)
  1531.         // and using cosh H = (e + cos v) / (1 + e cos v) to get rid of cosh H leads to:
  1532.         // dH = [sqrt (e^2 - 1) / (1 + e cos v)] dv + [sinh H / (e^2 - 1)] de
  1533.         // which is inverted and rewritten as:
  1534.         // dv = sqrt (1 - e^2) a/r dH - [sinh H / sqrt (e^2 - 1)] a/r de
  1535.         final double e2       = e * e;
  1536.         final double e2Mo     = e2 - 1;
  1537.         final double epsilon  = FastMath.sqrt(e2Mo);
  1538.         final double H        = getEccentricAnomaly();
  1539.         final double coshH    = FastMath.cosh(H);
  1540.         final double sinhH    = FastMath.sinh(H);
  1541.         final double aOr      = 1 / (e * coshH - 1);
  1542.         final double aFactor  = epsilon * aOr;
  1543.         final double eFactor  = sinhH * aOr / epsilon;

  1545.         // update anomaly row
  1546.         final double[] eRow           = jacobian[1];
  1547.         final double[] anomalyRow     = jacobian[5];
  1548.         for (int j = 0; j < anomalyRow.length; ++j) {
  1549.             anomalyRow[j] = aFactor * anomalyRow[j] - eFactor * eRow[j];
  1550.         }

  1552.         return jacobian;

  1554.     }

  1556.     /** {@inheritDoc} */
  1557.     @Override
  1558.     public void addKeplerContribution(final PositionAngleType type, final double gm,
  1559.                                       final double[] pDot) {
  1560.         pDot[5] += computeKeplerianAnomalyDot(type, a, e, gm, cachedAnomaly, cachedPositionAngleType);
  1561.     }

  1563.     /**
  1564.      * Compute rate of argument of latitude.
  1565.      * @param type position angle type of rate
  1566.      * @param a semi major axis
  1567.      * @param e eccentricity
  1568.      * @param mu mu
  1569.      * @param anomaly anomaly
  1570.      * @param cachedType position angle type of passed anomaly
  1571.      * @return first-order time derivative for anomaly
  1572.      * @since 12.2
  1573.      */
  1574.     private static double computeKeplerianAnomalyDot(final PositionAngleType type, final double a, final double e,
  1575.                                                      final double mu, final double anomaly, final PositionAngleType cachedType) {
  1576.         final double absA = FastMath.abs(a);
  1577.         final double n    = FastMath.sqrt(mu / absA) / absA;
  1578.         if (type == PositionAngleType.MEAN) {
  1579.             return n;
  1580.         }
  1581.         final double oMe2 = FastMath.abs(1 - e * e);
  1582.         final double ksi = 1 + e * FastMath.cos(KeplerianAnomalyUtility.convertAnomaly(cachedType, anomaly, e, PositionAngleType.TRUE));
  1583.         if (type == PositionAngleType.ECCENTRIC) {
  1584.             return n * ksi / oMe2;
  1585.         } else { // TRUE
  1586.             return n * ksi * ksi / (oMe2 * FastMath.sqrt(oMe2));
  1587.         }
  1588.     }

  1590.     /**  Returns a string representation of this Keplerian parameters object.
  1591.      * @return a string representation of this object
  1592.      */
  1593.     public String toString() {
  1594.         return new StringBuilder().append("Keplerian parameters: ").append('{').
  1595.                 append("a: ").append(a).
  1596.                 append("; e: ").append(e).
  1597.                 append("; i: ").append(FastMath.toDegrees(i)).
  1598.                 append("; pa: ").append(FastMath.toDegrees(pa)).
  1599.                 append("; raan: ").append(FastMath.toDegrees(raan)).
  1600.                 append("; v: ").append(FastMath.toDegrees(getTrueAnomaly())).
  1601.                 append(";}").toString();
  1602.     }

  1604.     /** {@inheritDoc} */
  1605.     @Override
  1606.     public PositionAngleType getCachedPositionAngleType() {
  1607.         return cachedPositionAngleType;
  1608.     }

  1610.     /** {@inheritDoc} */
  1611.     @Override
  1612.     public boolean hasNonKeplerianRates() {
  1613.         return hasNonKeplerianAcceleration();
  1614.     }

  1616.     /** {@inheritDoc} */
  1617.     @Override
  1618.     public KeplerianOrbit withKeplerianRates() {
  1619.         final PositionAngleType positionAngleType = getCachedPositionAngleType();
  1620.         return new KeplerianOrbit(a, e, i, pa, raan, cachedAnomaly, positionAngleType, positionAngleType,
  1621.                 getFrame(), getDate(), getMu());
  1622.     }

  1624.     /** Check if the given parameter is within an acceptable range.
  1625.      * The bounds are inclusive: an exception is raised when either of those conditions are met:
  1626.      * <ul>
  1627.      *     <li>The parameter is strictly greater than upperBound</li>
  1628.      *     <li>The parameter is strictly lower than lowerBound</li>
  1629.      * </ul>
  1630.      * <p>
  1631.      * In either of these cases, an OrekitException is raised.
  1632.      * </p>
  1633.      * @param parameterName name of the parameter
  1634.      * @param parameter value of the parameter
  1635.      * @param lowerBound lower bound of the acceptable range (inclusive)
  1636.      * @param upperBound upper bound of the acceptable range (inclusive)
  1637.      */
  1638.     private void checkParameterRangeInclusive(final String parameterName, final double parameter,
  1639.                                               final double lowerBound, final double upperBound) {
  1640.         if (parameter < lowerBound || parameter > upperBound) {
  1641.             throw new OrekitException(OrekitMessages.INVALID_PARAMETER_RANGE, parameterName,
  1642.                     parameter, lowerBound, upperBound);
  1643.         }
  1644.     }

  1646. }
Created with JaCoCo 0.8.12.202403310830