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


  20. import java.lang.reflect.Array;

  22. import org.hipparchus.CalculusFieldElement;
  23. import org.hipparchus.Field;
  24. import org.hipparchus.analysis.differentiation.FieldUnivariateDerivative1;
  25. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  26. import org.hipparchus.util.FastMath;
  27. import org.hipparchus.util.FieldSinCos;
  28. import org.hipparchus.util.MathArrays;
  29. import org.orekit.errors.OrekitException;
  30. import org.orekit.errors.OrekitIllegalArgumentException;
  31. import org.orekit.errors.OrekitInternalError;
  32. import org.orekit.errors.OrekitMessages;
  33. import org.orekit.frames.Frame;
  34. import org.orekit.time.FieldAbsoluteDate;
  35. import org.orekit.utils.FieldPVCoordinates;
  36. import org.orekit.utils.TimeStampedFieldPVCoordinates;


  39. /**
  40.  * This class handles traditional Keplerian orbital parameters.

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

  65.  * <p>
  66.  * The instance <code>KeplerianOrbit</code> is guaranteed to be immutable.
  67.  * </p>
  68.  * @see     Orbit
  69.  * @see    CircularOrbit
  70.  * @see    CartesianOrbit
  71.  * @see    EquinoctialOrbit
  72.  * @author Luc Maisonobe
  73.  * @author Guylaine Prat
  74.  * @author Fabien Maussion
  75.  * @author V&eacute;ronique Pommier-Maurussane
  76.  * @author Andrea Antolino
  77.  * @since 9.0
  78.  * @param <T> type of the field elements
  79.  */
  80. public class FieldKeplerianOrbit<T extends CalculusFieldElement<T>> extends FieldOrbit<T>
  81.         implements PositionAngleBased {

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

  86.     /** Semi-major axis (m). */
  87.     private final T a;

  89.     /** Eccentricity. */
  90.     private final T e;

  92.     /** Inclination (rad). */
  93.     private final T i;

  95.     /** Perigee Argument (rad). */
  96.     private final T pa;

  98.     /** Right Ascension of Ascending Node (rad). */
  99.     private final T raan;

  101.     /** True anomaly (rad). */
  102.     private final T v;

  104.     /** Semi-major axis derivative (m/s). */
  105.     private final T aDot;

  107.     /** Eccentricity derivative. */
  108.     private final T eDot;

  110.     /** Inclination derivative (rad/s). */
  111.     private final T iDot;

  113.     /** Perigee Argument derivative (rad/s). */
  114.     private final T paDot;

  116.     /** Right Ascension of Ascending Node derivative (rad/s). */
  117.     private final T raanDot;

  119.     /** True anomaly derivative (rad/s). */
  120.     private final T vDot;

  122.     /** Partial Cartesian coordinates (position and velocity are valid, acceleration may be missing). */
  123.     private FieldPVCoordinates<T> partialPV;

  125.     /** PThird Canonical Vector. */
  126.     private final FieldVector3D<T> PLUS_K;

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

  154.     /** Creates a new instance.
  155.      * @param a  semi-major axis (m), negative for hyperbolic orbits
  156.      * @param e eccentricity (positive or equal to 0)
  157.      * @param i inclination (rad)
  158.      * @param pa perigee argument (ω, rad)
  159.      * @param raan right ascension of ascending node (Ω, rad)
  160.      * @param anomaly mean, eccentric or true anomaly (rad)
  161.      * @param aDot  semi-major axis derivative, null if unknown (m/s)
  162.      * @param eDot eccentricity derivative, null if unknown
  163.      * @param iDot inclination derivative, null if unknown (rad/s)
  164.      * @param paDot perigee argument derivative, null if unknown (rad/s)
  165.      * @param raanDot right ascension of ascending node derivative, null if unknown (rad/s)
  166.      * @param anomalyDot mean, eccentric or true anomaly derivative, null if unknown (rad/s)
  167.      * @param type type of anomaly
  168.      * @param frame the frame in which the parameters are defined
  169.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  170.      * @param date date of the orbital parameters
  171.      * @param mu central attraction coefficient (m³/s²)
  172.      * @exception IllegalArgumentException if frame is not a {@link
  173.      * Frame#isPseudoInertial pseudo-inertial frame} or a and e don't match for hyperbolic orbits,
  174.      * or v is out of range for hyperbolic orbits
  175.      */
  176.     public FieldKeplerianOrbit(final T a, final T e, final T i,
  177.                                final T pa, final T raan, final T anomaly,
  178.                                final T aDot, final T eDot, final T iDot,
  179.                                final T paDot, final T raanDot, final T anomalyDot,
  180.                                final PositionAngleType type,
  181.                                final Frame frame, final FieldAbsoluteDate<T> date, final T mu)
  182.         throws IllegalArgumentException {
  183.         super(frame, date, mu);
  184.         if (a.multiply(e.negate().add(1)).getReal() < 0) {
  185.             throw new OrekitIllegalArgumentException(OrekitMessages.ORBIT_A_E_MISMATCH_WITH_CONIC_TYPE, a.getReal(), e.getReal());
  186.         }

  188.         // Checking eccentricity range
  189.         checkParameterRangeInclusive(ECCENTRICITY, e.getReal(), 0.0, Double.POSITIVE_INFINITY);

  191.         this.a       =    a;
  192.         this.aDot    =    aDot;
  193.         this.e       =    e;
  194.         this.eDot    =    eDot;
  195.         this.i       =    i;
  196.         this.iDot    =    iDot;
  197.         this.pa      =    pa;
  198.         this.paDot   =    paDot;
  199.         this.raan    =    raan;
  200.         this.raanDot =    raanDot;

  202.         this.PLUS_K = FieldVector3D.getPlusK(a.getField());  // third canonical vector

  204.         if (hasDerivatives()) {
  205.             final FieldUnivariateDerivative1<T> eUD        = new FieldUnivariateDerivative1<>(e, eDot);
  206.             final FieldUnivariateDerivative1<T> anomalyUD  = new FieldUnivariateDerivative1<>(anomaly,  anomalyDot);
  207.             final FieldUnivariateDerivative1<T> vUD;
  208.             switch (type) {
  209.                 case MEAN :
  210.                     vUD = (a.getReal() < 0) ?
  211.                           FieldKeplerianAnomalyUtility.hyperbolicMeanToTrue(eUD, anomalyUD) :
  212.                           FieldKeplerianAnomalyUtility.ellipticMeanToTrue(eUD, anomalyUD);
  213.                     break;
  214.                 case ECCENTRIC :
  215.                     vUD = (a.getReal() < 0) ?
  216.                           FieldKeplerianAnomalyUtility.hyperbolicEccentricToTrue(eUD, anomalyUD) :
  217.                           FieldKeplerianAnomalyUtility.ellipticEccentricToTrue(eUD, anomalyUD);
  218.                     break;
  219.                 case TRUE :
  220.                     vUD = anomalyUD;
  221.                     break;
  222.                 default : // this should never happen
  223.                     throw new OrekitInternalError(null);
  224.             }
  225.             this.v    = vUD.getValue();
  226.             this.vDot = vUD.getDerivative(1);
  227.         } else {
  228.             switch (type) {
  229.                 case MEAN :

  231.                     this.v = (a.getReal() < 0) ?
  232.                              FieldKeplerianAnomalyUtility.hyperbolicMeanToTrue(e, anomaly) :
  233.                              FieldKeplerianAnomalyUtility.ellipticMeanToTrue(e, anomaly);

  235.                     break;
  236.                 case ECCENTRIC :
  237.                     this.v = (a.getReal() < 0) ?
  238.                              FieldKeplerianAnomalyUtility.hyperbolicEccentricToTrue(e, anomaly) :
  239.                              FieldKeplerianAnomalyUtility.ellipticEccentricToTrue(e, anomaly);

  241.                     break;
  242.                 case TRUE :
  243.                     this.v = anomaly;
  244.                     break;
  245.                 default : // this should never happen
  246.                     throw new OrekitInternalError(null);
  247.             }
  248.             this.vDot = null;
  249.         }

  251.         // check true anomaly range
  252.         if (e.multiply(v.cos()).add(1).getReal() <= 0) {
  253.             final double vMax = e.reciprocal().negate().acos().getReal();
  254.             throw new OrekitIllegalArgumentException(OrekitMessages.ORBIT_ANOMALY_OUT_OF_HYPERBOLIC_RANGE,
  255.                                                      v.getReal(), e.getReal(), -vMax, vMax);
  256.         }

  258.         this.partialPV = null;

  260.     }

  262.     /** Constructor from Cartesian parameters.
  263.      *
  264.      * <p> The acceleration provided in {@code FieldPVCoordinates} is accessible using
  265.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  266.      * use {@code mu} and the position to compute the acceleration, including
  267.      * {@link #shiftedBy(CalculusFieldElement)} and {@link #getPVCoordinates(FieldAbsoluteDate, Frame)}.
  268.      *
  269.      * @param pvCoordinates the PVCoordinates of the satellite
  270.      * @param frame the frame in which are defined the {@link FieldPVCoordinates}
  271.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  272.      * @param mu central attraction coefficient (m³/s²)
  273.      * @exception IllegalArgumentException if frame is not a {@link
  274.      * Frame#isPseudoInertial pseudo-inertial frame}
  275.      */
  276.     public FieldKeplerianOrbit(final TimeStampedFieldPVCoordinates<T> pvCoordinates,
  277.                                final Frame frame, final T mu)
  278.         throws IllegalArgumentException {
  279.         this(pvCoordinates, frame, mu, hasNonKeplerianAcceleration(pvCoordinates, mu));
  280.     }

  282.     /** Constructor from Cartesian parameters.
  283.      *
  284.      * <p> The acceleration provided in {@code FieldPVCoordinates} is accessible using
  285.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  286.      * use {@code mu} and the position to compute the acceleration, including
  287.      * {@link #shiftedBy(CalculusFieldElement)} and {@link #getPVCoordinates(FieldAbsoluteDate, Frame)}.
  288.      *
  289.      * @param pvCoordinates the PVCoordinates of the satellite
  290.      * @param frame the frame in which are defined the {@link FieldPVCoordinates}
  291.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  292.      * @param mu central attraction coefficient (m³/s²)
  293.      * @param reliableAcceleration if true, the acceleration is considered to be reliable
  294.      * @exception IllegalArgumentException if frame is not a {@link
  295.      * Frame#isPseudoInertial pseudo-inertial frame}
  296.      */
  297.     private FieldKeplerianOrbit(final TimeStampedFieldPVCoordinates<T> pvCoordinates,
  298.                                 final Frame frame, final T mu,
  299.                                 final boolean reliableAcceleration)
  300.         throws IllegalArgumentException {

  302.         super(pvCoordinates, frame, mu);

  304.         // third canonical vector
  305.         this.PLUS_K = FieldVector3D.getPlusK(getOne().getField());

  307.         // compute inclination
  308.         final FieldVector3D<T> momentum = pvCoordinates.getMomentum();
  309.         final T m2 = momentum.getNormSq();

  311.         i = FieldVector3D.angle(momentum, PLUS_K);
  312.         // compute right ascension of ascending node
  313.         raan = FieldVector3D.crossProduct(PLUS_K, momentum).getAlpha();
  314.         // preliminary computations for parameters depending on orbit shape (elliptic or hyperbolic)
  315.         final FieldVector3D<T> pvP     = pvCoordinates.getPosition();
  316.         final FieldVector3D<T> pvV     = pvCoordinates.getVelocity();
  317.         final FieldVector3D<T> pvA     = pvCoordinates.getAcceleration();

  319.         final T   r2      = pvP.getNormSq();
  320.         final T   r       = r2.sqrt();
  321.         final T   V2      = pvV.getNormSq();
  322.         final T   rV2OnMu = r.multiply(V2).divide(mu);

  324.         // compute semi-major axis (will be negative for hyperbolic orbits)
  325.         a = r.divide(rV2OnMu.negate().add(2.0));
  326.         final T muA = a.multiply(mu);

  328.         // compute true anomaly
  329.         if (isElliptical()) {
  330.             // elliptic or circular orbit
  331.             final T eSE = FieldVector3D.dotProduct(pvP, pvV).divide(muA.sqrt());
  332.             final T eCE = rV2OnMu.subtract(1);
  333.             e = (eSE.multiply(eSE).add(eCE.multiply(eCE))).sqrt();
  334.             v = FieldKeplerianAnomalyUtility.ellipticEccentricToTrue(e, eSE.atan2(eCE));
  335.         } else {
  336.             // hyperbolic orbit
  337.             final T eSH = FieldVector3D.dotProduct(pvP, pvV).divide(muA.negate().sqrt());
  338.             final T eCH = rV2OnMu.subtract(1);
  339.             e = (m2.negate().divide(muA).add(1)).sqrt();
  340.             v = FieldKeplerianAnomalyUtility.hyperbolicEccentricToTrue(e, (eCH.add(eSH)).divide(eCH.subtract(eSH)).log().divide(2));
  341.         }

  343.         // Checking eccentricity range
  344.         checkParameterRangeInclusive(ECCENTRICITY, e.getReal(), 0.0, Double.POSITIVE_INFINITY);

  346.         // compute perigee argument
  347.         final FieldVector3D<T> node = new FieldVector3D<>(raan, getZero());
  348.         final T px = FieldVector3D.dotProduct(pvP, node);
  349.         final T py = FieldVector3D.dotProduct(pvP, FieldVector3D.crossProduct(momentum, node)).divide(m2.sqrt());
  350.         pa = py.atan2(px).subtract(v);

  352.         partialPV = pvCoordinates;

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

  357.             final T[][] jacobian = MathArrays.buildArray(a.getField(), 6, 6);
  358.             getJacobianWrtCartesian(PositionAngleType.MEAN, jacobian);

  360.             final FieldVector3D<T> keplerianAcceleration    = new FieldVector3D<>(r.multiply(r2).reciprocal().multiply(mu.negate()), pvP);
  361.             final FieldVector3D<T> nonKeplerianAcceleration = pvA.subtract(keplerianAcceleration);
  362.             final T   aX                       = nonKeplerianAcceleration.getX();
  363.             final T   aY                       = nonKeplerianAcceleration.getY();
  364.             final T   aZ                       = nonKeplerianAcceleration.getZ();
  365.             aDot    = jacobian[0][3].multiply(aX).add(jacobian[0][4].multiply(aY)).add(jacobian[0][5].multiply(aZ));
  366.             eDot    = jacobian[1][3].multiply(aX).add(jacobian[1][4].multiply(aY)).add(jacobian[1][5].multiply(aZ));
  367.             iDot    = jacobian[2][3].multiply(aX).add(jacobian[2][4].multiply(aY)).add(jacobian[2][5].multiply(aZ));
  368.             paDot   = jacobian[3][3].multiply(aX).add(jacobian[3][4].multiply(aY)).add(jacobian[3][5].multiply(aZ));
  369.             raanDot = jacobian[4][3].multiply(aX).add(jacobian[4][4].multiply(aY)).add(jacobian[4][5].multiply(aZ));

  371.             // in order to compute true anomaly derivative, we must compute
  372.             // mean anomaly derivative including Keplerian motion and convert to true anomaly
  373.             final T MDot = getKeplerianMeanMotion().
  374.                            add(jacobian[5][3].multiply(aX)).add(jacobian[5][4].multiply(aY)).add(jacobian[5][5].multiply(aZ));
  375.             final FieldUnivariateDerivative1<T> eUD = new FieldUnivariateDerivative1<>(e, eDot);
  376.             final FieldUnivariateDerivative1<T> MUD = new FieldUnivariateDerivative1<>(getMeanAnomaly(), MDot);
  377.             final FieldUnivariateDerivative1<T> vUD = (a.getReal() < 0) ?
  378.                                             FieldKeplerianAnomalyUtility.hyperbolicMeanToTrue(eUD, MUD) :
  379.                                             FieldKeplerianAnomalyUtility.ellipticMeanToTrue(eUD, MUD);
  380.             vDot = vUD.getDerivative(1);

  382.         } else {
  383.             // acceleration is either almost zero or NaN,
  384.             // we assume acceleration was not known
  385.             // we don't set up derivatives
  386.             aDot    = null;
  387.             eDot    = null;
  388.             iDot    = null;
  389.             paDot   = null;
  390.             raanDot = null;
  391.             vDot    = null;
  392.         }

  394.     }

  396.     /** Constructor from Cartesian parameters.
  397.      *
  398.      * <p> The acceleration provided in {@code FieldPVCoordinates} is accessible using
  399.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  400.      * use {@code mu} and the position to compute the acceleration, including
  401.      * {@link #shiftedBy(CalculusFieldElement)} and {@link #getPVCoordinates(FieldAbsoluteDate, Frame)}.
  402.      *
  403.      * @param FieldPVCoordinates the PVCoordinates of the satellite
  404.      * @param frame the frame in which are defined the {@link FieldPVCoordinates}
  405.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  406.      * @param date date of the orbital parameters
  407.      * @param mu central attraction coefficient (m³/s²)
  408.      * @exception IllegalArgumentException if frame is not a {@link
  409.      * Frame#isPseudoInertial pseudo-inertial frame}
  410.      */
  411.     public FieldKeplerianOrbit(final FieldPVCoordinates<T> FieldPVCoordinates,
  412.                                final Frame frame, final FieldAbsoluteDate<T> date, final T mu)
  413.         throws IllegalArgumentException {
  414.         this(new TimeStampedFieldPVCoordinates<>(date, FieldPVCoordinates), frame, mu);
  415.     }

  417.     /** Constructor from any kind of orbital parameters.
  418.      * @param op orbital parameters to copy
  419.      */
  420.     public FieldKeplerianOrbit(final FieldOrbit<T> op) {
  421.         this(op.getPVCoordinates(), op.getFrame(), op.getMu(), op.hasDerivatives());
  422.     }

  424.     /** Constructor from Field and KeplerianOrbit.
  425.      * <p>Build a FieldKeplerianOrbit from non-Field KeplerianOrbit.</p>
  426.      * @param field CalculusField to base object on
  427.      * @param op non-field orbit with only "constant" terms
  428.      * @since 12.0
  429.      */
  430.     public FieldKeplerianOrbit(final Field<T> field, final KeplerianOrbit op) {
  431.         this(field.getZero().add(op.getA()), field.getZero().add(op.getE()), field.getZero().add(op.getI()),
  432.                 field.getZero().add(op.getPerigeeArgument()), field.getZero().add(op.getRightAscensionOfAscendingNode()),
  433.                 field.getZero().add(op.getTrueAnomaly()),
  434.                 (op.hasDerivatives()) ? field.getZero().add(op.getADot()) : null,
  435.                 (op.hasDerivatives()) ? field.getZero().add(op.getEDot()) : null,
  436.                 (op.hasDerivatives()) ? field.getZero().add(op.getIDot()) : null,
  437.                 (op.hasDerivatives()) ? field.getZero().add(op.getPerigeeArgumentDot()) : null,
  438.                 (op.hasDerivatives()) ? field.getZero().add(op.getRightAscensionOfAscendingNodeDot()) : null,
  439.                 (op.hasDerivatives()) ? field.getZero().add(op.getTrueAnomalyDot()) : null, PositionAngleType.TRUE,
  440.                 op.getFrame(), new FieldAbsoluteDate<>(field, op.getDate()), field.getZero().add(op.getMu()));
  441.     }

  443.     /** Constructor from Field and Orbit.
  444.      * <p>Build a FieldKeplerianOrbit from any non-Field Orbit.</p>
  445.      * @param field CalculusField to base object on
  446.      * @param op non-field orbit with only "constant" terms
  447.      * @since 12.0
  448.      */
  449.     public FieldKeplerianOrbit(final Field<T> field, final Orbit op) {
  450.         this(field, new KeplerianOrbit(op));
  451.     }

  453.     /** {@inheritDoc} */
  454.     public OrbitType getType() {
  455.         return OrbitType.KEPLERIAN;
  456.     }

  458.     /** {@inheritDoc} */
  459.     public T getA() {
  460.         return a;
  461.     }

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

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

  473.     /** {@inheritDoc} */
  474.     public T getEDot() {
  475.         return eDot;
  476.     }

  478.     /** {@inheritDoc} */
  479.     public T getI() {
  480.         return i;
  481.     }

  483.     /** {@inheritDoc} */
  484.     public T getIDot() {
  485.         return iDot;
  486.     }

  488.     /** Get the perigee argument.
  489.      * @return perigee argument (rad)
  490.      */
  491.     public T getPerigeeArgument() {
  492.         return pa;
  493.     }

  495.     /** Get the perigee argument derivative.
  496.      * <p>
  497.      * If the orbit was created without derivatives, the value returned is null.
  498.      * </p>
  499.      * @return perigee argument derivative (rad/s)
  500.      */
  501.     public T getPerigeeArgumentDot() {
  502.         return paDot;
  503.     }

  505.     /** Get the right ascension of the ascending node.
  506.      * @return right ascension of the ascending node (rad)
  507.      */
  508.     public T getRightAscensionOfAscendingNode() {
  509.         return raan;
  510.     }

  512.     /** Get the right ascension of the ascending node derivative.
  513.      * <p>
  514.      * If the orbit was created without derivatives, the value returned is null.
  515.      * </p>
  516.      * @return right ascension of the ascending node derivative (rad/s)
  517.      */
  518.     public T getRightAscensionOfAscendingNodeDot() {
  519.         return raanDot;
  520.     }

  522.     /** Get the true anomaly.
  523.      * @return true anomaly (rad)
  524.      */
  525.     public T getTrueAnomaly() {
  526.         return v;
  527.     }

  529.     /** Get the true anomaly derivative.
  530.      * <p>
  531.      * If the orbit was created without derivatives, the value returned is null.
  532.      * </p>
  533.      * @return true anomaly derivative (rad/s)
  534.      */
  535.     public T getTrueAnomalyDot() {
  536.         return vDot;
  537.     }

  539.     /** Get the eccentric anomaly.
  540.      * @return eccentric anomaly (rad)
  541.      */
  542.     public T getEccentricAnomaly() {
  543.         return (a.getReal() < 0) ?
  544.                FieldKeplerianAnomalyUtility.hyperbolicTrueToEccentric(e, v) :
  545.                FieldKeplerianAnomalyUtility.ellipticTrueToEccentric(e, v);
  546.     }

  548.     /** Get the eccentric anomaly derivative.
  549.      * <p>
  550.      * If the orbit was created without derivatives, the value returned is null.
  551.      * </p>
  552.      * @return eccentric anomaly derivative (rad/s)
  553.      */
  554.     public T getEccentricAnomalyDot() {

  556.         if (!hasDerivatives()) {
  557.             return null;
  558.         }

  560.         final FieldUnivariateDerivative1<T> eUD = new FieldUnivariateDerivative1<>(e, eDot);
  561.         final FieldUnivariateDerivative1<T> vUD = new FieldUnivariateDerivative1<>(v, vDot);
  562.         final FieldUnivariateDerivative1<T> EUD = (a.getReal() < 0) ?
  563.                                                 FieldKeplerianAnomalyUtility.hyperbolicTrueToEccentric(eUD, vUD) :
  564.                                                 FieldKeplerianAnomalyUtility.ellipticTrueToEccentric(eUD, vUD);
  565.         return EUD.getDerivative(1);

  567.     }

  569.     /** Get the mean anomaly.
  570.      * @return mean anomaly (rad)
  571.      */
  572.     public T getMeanAnomaly() {
  573.         return (a.getReal() < 0) ?
  574.                FieldKeplerianAnomalyUtility.hyperbolicTrueToMean(e, v) :
  575.                FieldKeplerianAnomalyUtility.ellipticTrueToMean(e, v);
  576.     }

  578.     /** Get the mean anomaly derivative.
  579.      * <p>
  580.      * If the orbit was created without derivatives, the value returned is null.
  581.      * </p>
  582.      * @return mean anomaly derivative (rad/s)
  583.      */
  584.     public T getMeanAnomalyDot() {

  586.         if (!hasDerivatives()) {
  587.             return null;
  588.         }

  590.         final FieldUnivariateDerivative1<T> eUD = new FieldUnivariateDerivative1<>(e, eDot);
  591.         final FieldUnivariateDerivative1<T> vUD = new FieldUnivariateDerivative1<>(v, vDot);
  592.         final FieldUnivariateDerivative1<T> MUD = (a.getReal() < 0) ?
  593.                                                 FieldKeplerianAnomalyUtility.hyperbolicTrueToMean(eUD, vUD) :
  594.                                                 FieldKeplerianAnomalyUtility.ellipticTrueToMean(eUD, vUD);
  595.         return MUD.getDerivative(1);

  597.     }

  599.     /** Get the anomaly.
  600.      * @param type type of the angle
  601.      * @return anomaly (rad)
  602.      */
  603.     public T getAnomaly(final PositionAngleType type) {
  604.         return (type == PositionAngleType.MEAN) ? getMeanAnomaly() :
  605.                                               ((type == PositionAngleType.ECCENTRIC) ? getEccentricAnomaly() :
  606.                                                                                    getTrueAnomaly());
  607.     }

  609.     /** Get the anomaly derivative.
  610.      * <p>
  611.      * If the orbit was created without derivatives, the value returned is null.
  612.      * </p>
  613.      * @param type type of the angle
  614.      * @return anomaly derivative (rad/s)
  615.      */
  616.     public T getAnomalyDot(final PositionAngleType type) {
  617.         return (type == PositionAngleType.MEAN) ? getMeanAnomalyDot() :
  618.                                               ((type == PositionAngleType.ECCENTRIC) ? getEccentricAnomalyDot() :
  619.                                                                                    getTrueAnomalyDot());
  620.     }

  622.     /** {@inheritDoc} */
  623.     @Override
  624.     public boolean hasDerivatives() {
  625.         return aDot != null;
  626.     }

  628.     /** {@inheritDoc} */
  629.     public T getEquinoctialEx() {
  630.         return e.multiply(pa.add(raan).cos());
  631.     }

  633.     /** {@inheritDoc} */
  634.     public T getEquinoctialExDot() {

  636.         if (!hasDerivatives()) {
  637.             return null;
  638.         }

  640.         final FieldUnivariateDerivative1<T> eUD    = new FieldUnivariateDerivative1<>(e,    eDot);
  641.         final FieldUnivariateDerivative1<T> paUD   = new FieldUnivariateDerivative1<>(pa,   paDot);
  642.         final FieldUnivariateDerivative1<T> raanUD = new FieldUnivariateDerivative1<>(raan, raanDot);
  643.         return eUD.multiply(paUD.add(raanUD).cos()).getDerivative(1);

  645.     }

  647.     /** {@inheritDoc} */
  648.     public T getEquinoctialEy() {
  649.         return  e.multiply((pa.add(raan)).sin());
  650.     }

  652.     /** {@inheritDoc} */
  653.     public T getEquinoctialEyDot() {

  655.         if (!hasDerivatives()) {
  656.             return null;
  657.         }

  659.         final FieldUnivariateDerivative1<T> eUD    = new FieldUnivariateDerivative1<>(e,    eDot);
  660.         final FieldUnivariateDerivative1<T> paUD   = new FieldUnivariateDerivative1<>(pa,   paDot);
  661.         final FieldUnivariateDerivative1<T> raanUD = new FieldUnivariateDerivative1<>(raan, raanDot);
  662.         return eUD.multiply(paUD.add(raanUD).sin()).getDerivative(1);

  664.     }

  666.     /** {@inheritDoc} */
  667.     public T getHx() {
  668.         // Check for equatorial retrograde orbit
  669.         if (FastMath.abs(i.subtract(i.getPi()).getReal()) < 1.0e-10) {
  670.             return getZero().add(Double.NaN);
  671.         }
  672.         return  raan.cos().multiply(i.divide(2).tan());
  673.     }

  675.     /** {@inheritDoc} */
  676.     public T getHxDot() {

  678.         if (!hasDerivatives()) {
  679.             return null;
  680.         }

  682.         // Check for equatorial retrograde orbit
  683.         if (FastMath.abs(i.subtract(i.getPi()).getReal()) < 1.0e-10) {
  684.             return getZero().add(Double.NaN);
  685.         }

  687.         final FieldUnivariateDerivative1<T> iUD    = new FieldUnivariateDerivative1<>(i,    iDot);
  688.         final FieldUnivariateDerivative1<T> raanUD = new FieldUnivariateDerivative1<>(raan, raanDot);
  689.         return raanUD.cos().multiply(iUD.multiply(0.5).tan()).getDerivative(1);

  691.     }

  693.     /** {@inheritDoc} */
  694.     public T getHy() {
  695.         // Check for equatorial retrograde orbit
  696.         if (FastMath.abs(i.subtract(i.getPi()).getReal()) < 1.0e-10) {
  697.             return getZero().add(Double.NaN);
  698.         }
  699.         return  raan.sin().multiply(i.divide(2).tan());
  700.     }

  702.     /** {@inheritDoc} */
  703.     public T getHyDot() {

  705.         if (!hasDerivatives()) {
  706.             return null;
  707.         }

  709.         // Check for equatorial retrograde orbit
  710.         if (FastMath.abs(i.subtract(i.getPi()).getReal()) < 1.0e-10) {
  711.             return getZero().add(Double.NaN);
  712.         }

  714.         final FieldUnivariateDerivative1<T> iUD    = new FieldUnivariateDerivative1<>(i,    iDot);
  715.         final FieldUnivariateDerivative1<T> raanUD = new FieldUnivariateDerivative1<>(raan, raanDot);
  716.         return raanUD.sin().multiply(iUD.multiply(0.5).tan()).getDerivative(1);

  718.     }

  720.     /** {@inheritDoc} */
  721.     public T getLv() {
  722.         return pa.add(raan).add(v);
  723.     }

  725.     /** {@inheritDoc} */
  726.     public T getLvDot() {
  727.         return hasDerivatives() ?
  728.                paDot.add(raanDot).add(vDot) :
  729.                null;
  730.     }

  732.     /** {@inheritDoc} */
  733.     public T getLE() {
  734.         return pa.add(raan).add(getEccentricAnomaly());
  735.     }

  737.     /** {@inheritDoc} */
  738.     public T getLEDot() {
  739.         return hasDerivatives() ?
  740.                paDot.add(raanDot).add(getEccentricAnomalyDot()) :
  741.                null;
  742.     }

  744.     /** {@inheritDoc} */
  745.     public T getLM() {
  746.         return pa.add(raan).add(getMeanAnomaly());
  747.     }

  749.     /** {@inheritDoc} */
  750.     public T getLMDot() {
  751.         return hasDerivatives() ?
  752.                paDot.add(raanDot).add(getMeanAnomalyDot()) :
  753.                null;
  754.     }

  756.     /** Compute reference axes.
  757.      * @return referecne axes
  758.      * @since 12.0
  759.      */
  760.     private FieldVector3D<T>[] referenceAxes() {

  762.         // preliminary variables
  763.         final FieldSinCos<T> scRaan = FastMath.sinCos(raan);
  764.         final FieldSinCos<T> scPa   = FastMath.sinCos(pa);
  765.         final FieldSinCos<T> scI    = FastMath.sinCos(i);
  766.         final T cosRaan = scRaan.cos();
  767.         final T sinRaan = scRaan.sin();
  768.         final T cosPa   = scPa.cos();
  769.         final T sinPa   = scPa.sin();
  770.         final T cosI    = scI.cos();
  771.         final T sinI    = scI.sin();
  772.         final T crcp    = cosRaan.multiply(cosPa);
  773.         final T crsp    = cosRaan.multiply(sinPa);
  774.         final T srcp    = sinRaan.multiply(cosPa);
  775.         final T srsp    = sinRaan.multiply(sinPa);

  777.         // reference axes defining the orbital plane
  778.         @SuppressWarnings("unchecked")
  779.         final FieldVector3D<T>[] axes = (FieldVector3D<T>[]) Array.newInstance(FieldVector3D.class, 2);
  780.         axes[0] = new FieldVector3D<>(crcp.subtract(cosI.multiply(srsp)),  srcp.add(cosI.multiply(crsp)), sinI.multiply(sinPa));
  781.         axes[1] = new FieldVector3D<>(crsp.add(cosI.multiply(srcp)).negate(), cosI.multiply(crcp).subtract(srsp), sinI.multiply(cosPa));

  783.         return axes;

  785.     }

  787.     /** Compute position and velocity but not acceleration.
  788.      */
  789.     private void computePVWithoutA() {

  791.         if (partialPV != null) {
  792.             // already computed
  793.             return;
  794.         }

  796.         final FieldVector3D<T>[] axes = referenceAxes();

  798.         if (isElliptical()) {

  800.             // elliptical case

  802.             // elliptic eccentric anomaly
  803.             final T uME2             = e.negate().add(1).multiply(e.add(1));
  804.             final T s1Me2            = uME2.sqrt();
  805.             final FieldSinCos<T> scE = FastMath.sinCos(getEccentricAnomaly());
  806.             final T cosE             = scE.cos();
  807.             final T sinE             = scE.sin();

  809.             // coordinates of position and velocity in the orbital plane
  810.             final T x      = a.multiply(cosE.subtract(e));
  811.             final T y      = a.multiply(sinE).multiply(s1Me2);
  812.             final T factor = FastMath.sqrt(getMu().divide(a)).divide(e.negate().multiply(cosE).add(1));
  813.             final T xDot   = sinE.negate().multiply(factor);
  814.             final T yDot   = cosE.multiply(s1Me2).multiply(factor);

  816.             final FieldVector3D<T> position = new FieldVector3D<>(x, axes[0], y, axes[1]);
  817.             final FieldVector3D<T> velocity = new FieldVector3D<>(xDot, axes[0], yDot, axes[1]);
  818.             partialPV = new FieldPVCoordinates<>(position, velocity);

  820.         } else {

  822.             // hyperbolic case

  824.             // compute position and velocity factors
  825.             final FieldSinCos<T> scV = FastMath.sinCos(v);
  826.             final T sinV             = scV.sin();
  827.             final T cosV             = scV.cos();
  828.             final T f                = a.multiply(e.multiply(e).negate().add(1));
  829.             final T posFactor        = f.divide(e.multiply(cosV).add(1));
  830.             final T velFactor        = FastMath.sqrt(getMu().divide(f));

  832.             final FieldVector3D<T> position     = new FieldVector3D<>(posFactor.multiply(cosV), axes[0], posFactor.multiply(sinV), axes[1]);
  833.             final FieldVector3D<T> velocity     = new FieldVector3D<>(velFactor.multiply(sinV).negate(), axes[0], velFactor.multiply(e.add(cosV)), axes[1]);
  834.             partialPV = new FieldPVCoordinates<>(position, velocity);

  836.         }

  838.     }

  840.     /** Compute non-Keplerian part of the acceleration from first time derivatives.
  841.      * <p>
  842.      * This method should be called only when {@link #hasDerivatives()} returns true.
  843.      * </p>
  844.      * @return non-Keplerian part of the acceleration
  845.      */
  846.     private FieldVector3D<T> nonKeplerianAcceleration() {

  848.         final T[][] dCdP = MathArrays.buildArray(a.getField(), 6, 6);
  849.         getJacobianWrtParameters(PositionAngleType.MEAN, dCdP);

  851.         final T nonKeplerianMeanMotion = getMeanAnomalyDot().subtract(getKeplerianMeanMotion());
  852.         final T nonKeplerianAx =     dCdP[3][0].multiply(aDot).
  853.                                  add(dCdP[3][1].multiply(eDot)).
  854.                                  add(dCdP[3][2].multiply(iDot)).
  855.                                  add(dCdP[3][3].multiply(paDot)).
  856.                                  add(dCdP[3][4].multiply(raanDot)).
  857.                                  add(dCdP[3][5].multiply(nonKeplerianMeanMotion));
  858.         final T nonKeplerianAy =     dCdP[4][0].multiply(aDot).
  859.                                  add(dCdP[4][1].multiply(eDot)).
  860.                                  add(dCdP[4][2].multiply(iDot)).
  861.                                  add(dCdP[4][3].multiply(paDot)).
  862.                                  add(dCdP[4][4].multiply(raanDot)).
  863.                                  add(dCdP[4][5].multiply(nonKeplerianMeanMotion));
  864.         final T nonKeplerianAz =     dCdP[5][0].multiply(aDot).
  865.                                  add(dCdP[5][1].multiply(eDot)).
  866.                                  add(dCdP[5][2].multiply(iDot)).
  867.                                  add(dCdP[5][3].multiply(paDot)).
  868.                                  add(dCdP[5][4].multiply(raanDot)).
  869.                                  add(dCdP[5][5].multiply(nonKeplerianMeanMotion));

  871.         return new FieldVector3D<>(nonKeplerianAx, nonKeplerianAy, nonKeplerianAz);

  873.     }

  875.     /** {@inheritDoc} */
  876.     protected FieldVector3D<T> initPosition() {
  877.         final FieldVector3D<T>[] axes = referenceAxes();

  879.         if (isElliptical()) {

  881.             // elliptical case

  883.             // elliptic eccentric anomaly
  884.             final T uME2             = e.negate().add(1).multiply(e.add(1));
  885.             final T s1Me2            = uME2.sqrt();
  886.             final FieldSinCos<T> scE = FastMath.sinCos(getEccentricAnomaly());
  887.             final T cosE             = scE.cos();
  888.             final T sinE             = scE.sin();

  890.             return new FieldVector3D<>(a.multiply(cosE.subtract(e)), axes[0], a.multiply(sinE).multiply(s1Me2), axes[1]);

  892.         } else {

  894.             // hyperbolic case

  896.             // compute position and velocity factors
  897.             final FieldSinCos<T> scV = FastMath.sinCos(v);
  898.             final T sinV             = scV.sin();
  899.             final T cosV             = scV.cos();
  900.             final T f                = a.multiply(e.multiply(e).negate().add(1));
  901.             final T posFactor        = f.divide(e.multiply(cosV).add(1));

  903.             return new FieldVector3D<>(posFactor.multiply(cosV), axes[0], posFactor.multiply(sinV), axes[1]);

  905.         }

  907.     }

  909.     /** {@inheritDoc} */
  910.     protected TimeStampedFieldPVCoordinates<T> initPVCoordinates() {

  912.         // position and velocity
  913.         computePVWithoutA();

  915.         // acceleration
  916.         final T r2 = partialPV.getPosition().getNormSq();
  917.         final FieldVector3D<T> keplerianAcceleration = new FieldVector3D<>(r2.multiply(FastMath.sqrt(r2)).reciprocal().multiply(getMu().negate()),
  918.                                                                            partialPV.getPosition());
  919.         final FieldVector3D<T> acceleration = hasDerivatives() ?
  920.                                               keplerianAcceleration.add(nonKeplerianAcceleration()) :
  921.                                               keplerianAcceleration;

  923.         return new TimeStampedFieldPVCoordinates<>(getDate(), partialPV.getPosition(), partialPV.getVelocity(), acceleration);

  925.     }

  927.     /** {@inheritDoc} */
  928.     public FieldKeplerianOrbit<T> shiftedBy(final double dt) {
  929.         return shiftedBy(getZero().add(dt));
  930.     }

  932.     /** {@inheritDoc} */
  933.     public FieldKeplerianOrbit<T> shiftedBy(final T dt) {

  935.         // use Keplerian-only motion
  936.         final FieldKeplerianOrbit<T> keplerianShifted = new FieldKeplerianOrbit<>(a, e, i, pa, raan,
  937.                                                                                   getKeplerianMeanMotion().multiply(dt).add(getMeanAnomaly()),
  938.                                                                                   PositionAngleType.MEAN, getFrame(), getDate().shiftedBy(dt), getMu());

  940.         if (hasDerivatives()) {

  942.             // extract non-Keplerian acceleration from first time derivatives
  943.             final FieldVector3D<T> nonKeplerianAcceleration = nonKeplerianAcceleration();

  945.             // add quadratic effect of non-Keplerian acceleration to Keplerian-only shift
  946.             keplerianShifted.computePVWithoutA();
  947.             final FieldVector3D<T> fixedP   = new FieldVector3D<>(getOne(), keplerianShifted.partialPV.getPosition(),
  948.                                                                   dt.multiply(dt).multiply(0.5), nonKeplerianAcceleration);
  949.             final T   fixedR2 = fixedP.getNormSq();
  950.             final T   fixedR  = fixedR2.sqrt();
  951.             final FieldVector3D<T> fixedV  = new FieldVector3D<>(getOne(), keplerianShifted.partialPV.getVelocity(),
  952.                                                                  dt, nonKeplerianAcceleration);
  953.             final FieldVector3D<T> fixedA  = new FieldVector3D<>(fixedR2.multiply(fixedR).reciprocal().multiply(getMu().negate()),
  954.                                                                  keplerianShifted.partialPV.getPosition(),
  955.                                                                  getOne(), nonKeplerianAcceleration);

  957.             // build a new orbit, taking non-Keplerian acceleration into account
  958.             return new FieldKeplerianOrbit<>(new TimeStampedFieldPVCoordinates<>(keplerianShifted.getDate(),
  959.                                                                                  fixedP, fixedV, fixedA),
  960.                                              keplerianShifted.getFrame(), keplerianShifted.getMu());

  962.         } else {
  963.             // Keplerian-only motion is all we can do
  964.             return keplerianShifted;
  965.         }

  967.     }

  969.     /** {@inheritDoc} */
  970.     protected T[][] computeJacobianMeanWrtCartesian() {
  971.         if (isElliptical()) {
  972.             return computeJacobianMeanWrtCartesianElliptical();
  973.         } else {
  974.             return computeJacobianMeanWrtCartesianHyperbolic();
  975.         }
  976.     }

  978.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  979.      * <p>
  980.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  981.      * respect to Cartesian coordinate j (x for j=0, y for j=1, z for j=2, xDot for j=3,
  982.      * yDot for j=4, zDot for j=5).
  983.      * </p>
  984.      * @return 6x6 Jacobian matrix
  985.      */
  986.     private T[][] computeJacobianMeanWrtCartesianElliptical() {

  988.         final T[][] jacobian = MathArrays.buildArray(getA().getField(), 6, 6);

  990.         // compute various intermediate parameters
  991.         computePVWithoutA();
  992.         final FieldVector3D<T> position = partialPV.getPosition();
  993.         final FieldVector3D<T> velocity = partialPV.getVelocity();
  994.         final FieldVector3D<T> momentum = partialPV.getMomentum();
  995.         final T v2         = velocity.getNormSq();
  996.         final T r2         = position.getNormSq();
  997.         final T r          = r2.sqrt();
  998.         final T r3         = r.multiply(r2);

  1000.         final T px         = position.getX();
  1001.         final T py         = position.getY();
  1002.         final T pz         = position.getZ();
  1003.         final T vx         = velocity.getX();
  1004.         final T vy         = velocity.getY();
  1005.         final T vz         = velocity.getZ();
  1006.         final T mx         = momentum.getX();
  1007.         final T my         = momentum.getY();
  1008.         final T mz         = momentum.getZ();

  1010.         final T mu         = getMu();
  1011.         final T sqrtMuA    = FastMath.sqrt(a.multiply(mu));
  1012.         final T sqrtAoMu   = FastMath.sqrt(a.divide(mu));
  1013.         final T a2         = a.multiply(a);
  1014.         final T twoA       = a.multiply(2);
  1015.         final T rOnA       = r.divide(a);

  1017.         final T oMe2       = e.multiply(e).negate().add(1);
  1018.         final T epsilon    = oMe2.sqrt();
  1019.         final T sqrtRec    = epsilon.reciprocal();

  1021.         final FieldSinCos<T> scI  = FastMath.sinCos(i);
  1022.         final FieldSinCos<T> scPA = FastMath.sinCos(pa);
  1023.         final T cosI       = scI.cos();
  1024.         final T sinI       = scI.sin();
  1025.         final T cosPA      = scPA.cos();
  1026.         final T sinPA      = scPA.sin();

  1028.         final T pv         = FieldVector3D.dotProduct(position, velocity);
  1029.         final T cosE       = a.subtract(r).divide(a.multiply(e));
  1030.         final T sinE       = pv.divide(e.multiply(sqrtMuA));

  1032.         // da
  1033.         final FieldVector3D<T> vectorAR = new FieldVector3D<>(a2.multiply(2).divide(r3), position);
  1034.         final FieldVector3D<T> vectorARDot = velocity.scalarMultiply(a2.multiply(mu.divide(2.).reciprocal()));
  1035.         fillHalfRow(getOne(), vectorAR,    jacobian[0], 0);
  1036.         fillHalfRow(getOne(), vectorARDot, jacobian[0], 3);

  1038.         // de
  1039.         final T factorER3 = pv.divide(twoA);
  1040.         final FieldVector3D<T> vectorER   = new FieldVector3D<>(cosE.multiply(v2).divide(r.multiply(mu)), position,
  1041.                                                                 sinE.divide(sqrtMuA), velocity,
  1042.                                                                 factorER3.negate().multiply(sinE).divide(sqrtMuA), vectorAR);
  1043.         final FieldVector3D<T> vectorERDot = new FieldVector3D<>(sinE.divide(sqrtMuA), position,
  1044.                                                                  cosE.multiply(mu.divide(2.).reciprocal()).multiply(r), velocity,
  1045.                                                                  factorER3.negate().multiply(sinE).divide(sqrtMuA), vectorARDot);
  1046.         fillHalfRow(getOne(), vectorER,    jacobian[1], 0);
  1047.         fillHalfRow(getOne(), vectorERDot, jacobian[1], 3);

  1049.         // dE / dr (Eccentric anomaly)
  1050.         final T coefE = cosE.divide(e.multiply(sqrtMuA));
  1051.         final FieldVector3D<T>  vectorEAnR =
  1052.             new FieldVector3D<>(sinE.negate().multiply(v2).divide(e.multiply(r).multiply(mu)), position, coefE, velocity,
  1053.                                 factorER3.negate().multiply(coefE), vectorAR);

  1055.         // dE / drDot
  1056.         final FieldVector3D<T>  vectorEAnRDot =
  1057.             new FieldVector3D<>(sinE.multiply(-2).multiply(r).divide(e.multiply(mu)), velocity, coefE, position,
  1058.                                 factorER3.negate().multiply(coefE), vectorARDot);

  1060.         // precomputing some more factors
  1061.         final T s1 = sinE.negate().multiply(pz).divide(r).subtract(cosE.multiply(vz).multiply(sqrtAoMu));
  1062.         final T s2 = cosE.negate().multiply(pz).divide(r3);
  1063.         final T s3 = sinE.multiply(vz).divide(sqrtMuA.multiply(-2));
  1064.         final T t1 = sqrtRec.multiply(cosE.multiply(pz).divide(r).subtract(sinE.multiply(vz).multiply(sqrtAoMu)));
  1065.         final T t2 = sqrtRec.multiply(sinE.negate().multiply(pz).divide(r3));
  1066.         final T t3 = sqrtRec.multiply(cosE.subtract(e)).multiply(vz).divide(sqrtMuA.multiply(2));
  1067.         final T t4 = sqrtRec.multiply(e.multiply(sinI).multiply(cosPA).multiply(sqrtRec).subtract(vz.multiply(sqrtAoMu)));
  1068.         final FieldVector3D<T> s = new FieldVector3D<>(cosE.divide(r), this.PLUS_K,
  1069.                                                        s1,       vectorEAnR,
  1070.                                                        s2,       position,
  1071.                                                        s3,       vectorAR);
  1072.         final FieldVector3D<T> sDot = new FieldVector3D<>(sinE.negate().multiply(sqrtAoMu), this.PLUS_K,
  1073.                                                           s1,               vectorEAnRDot,
  1074.                                                           s3,               vectorARDot);
  1075.         final FieldVector3D<T> t =
  1076.             new FieldVector3D<>(sqrtRec.multiply(sinE).divide(r), this.PLUS_K).add(new FieldVector3D<>(t1, vectorEAnR,
  1077.                                                                                                        t2, position,
  1078.                                                                                                        t3, vectorAR,
  1079.                                                                                                        t4, vectorER));
  1080.         final FieldVector3D<T> tDot = new FieldVector3D<>(sqrtRec.multiply(cosE.subtract(e)).multiply(sqrtAoMu), this.PLUS_K,
  1081.                                                           t1,                                                    vectorEAnRDot,
  1082.                                                           t3,                                                    vectorARDot,
  1083.                                                           t4,                                                    vectorERDot);

  1085.         // di
  1086.         final T factorI1 = sinI.negate().multiply(sqrtRec).divide(sqrtMuA);
  1087.         final T i1 =  factorI1;
  1088.         final T i2 =  factorI1.negate().multiply(mz).divide(twoA);
  1089.         final T i3 =  factorI1.multiply(mz).multiply(e).divide(oMe2);
  1090.         final T i4 = cosI.multiply(sinPA);
  1091.         final T i5 = cosI.multiply(cosPA);
  1092.         fillHalfRow(i1, new FieldVector3D<>(vy, vx.negate(), getZero()), i2, vectorAR, i3, vectorER, i4, s, i5, t,
  1093.                     jacobian[2], 0);
  1094.         fillHalfRow(i1, new FieldVector3D<>(py.negate(), px, getZero()), i2, vectorARDot, i3, vectorERDot, i4, sDot, i5, tDot,
  1095.                     jacobian[2], 3);

  1097.         // dpa
  1098.         fillHalfRow(cosPA.divide(sinI), s,    sinPA.negate().divide(sinI), t,    jacobian[3], 0);
  1099.         fillHalfRow(cosPA.divide(sinI), sDot, sinPA.negate().divide(sinI), tDot, jacobian[3], 3);

  1101.         // dRaan
  1102.         final T factorRaanR = (a.multiply(mu).multiply(oMe2).multiply(sinI).multiply(sinI)).reciprocal();
  1103.         fillHalfRow( factorRaanR.negate().multiply(my), new FieldVector3D<>(getZero(), vz, vy.negate()),
  1104.                      factorRaanR.multiply(mx), new FieldVector3D<>(vz.negate(), getZero(),  vx),
  1105.                      jacobian[4], 0);
  1106.         fillHalfRow(factorRaanR.negate().multiply(my), new FieldVector3D<>(getZero(), pz.negate(),  py),
  1107.                      factorRaanR.multiply(mx), new FieldVector3D<>(pz, getZero(), px.negate()),
  1108.                      jacobian[4], 3);

  1110.         // dM
  1111.         fillHalfRow(rOnA, vectorEAnR,    sinE.negate(), vectorER,    jacobian[5], 0);
  1112.         fillHalfRow(rOnA, vectorEAnRDot, sinE.negate(), vectorERDot, jacobian[5], 3);

  1114.         return jacobian;

  1116.     }

  1118.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  1119.      * <p>
  1120.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  1121.      * respect to Cartesian coordinate j (x for j=0, y for j=1, z for j=2, xDot for j=3,
  1122.      * yDot for j=4, zDot for j=5).
  1123.      * </p>
  1124.      * @return 6x6 Jacobian matrix
  1125.      */
  1126.     private T[][] computeJacobianMeanWrtCartesianHyperbolic() {

  1128.         final T[][] jacobian = MathArrays.buildArray(getA().getField(), 6, 6);

  1130.         // compute various intermediate parameters
  1131.         computePVWithoutA();
  1132.         final FieldVector3D<T> position = partialPV.getPosition();
  1133.         final FieldVector3D<T> velocity = partialPV.getVelocity();
  1134.         final FieldVector3D<T> momentum = partialPV.getMomentum();
  1135.         final T r2         = position.getNormSq();
  1136.         final T r          = r2.sqrt();
  1137.         final T r3         = r.multiply(r2);

  1139.         final T x          = position.getX();
  1140.         final T y          = position.getY();
  1141.         final T z          = position.getZ();
  1142.         final T vx         = velocity.getX();
  1143.         final T vy         = velocity.getY();
  1144.         final T vz         = velocity.getZ();
  1145.         final T mx         = momentum.getX();
  1146.         final T my         = momentum.getY();
  1147.         final T mz         = momentum.getZ();

  1149.         final T mu         = getMu();
  1150.         final T absA       = a.negate();
  1151.         final T sqrtMuA    = absA.multiply(mu).sqrt();
  1152.         final T a2         = a.multiply(a);
  1153.         final T rOa        = r.divide(absA);

  1155.         final FieldSinCos<T> scI = FastMath.sinCos(i);
  1156.         final T cosI       = scI.cos();
  1157.         final T sinI       = scI.sin();

  1159.         final T pv         = FieldVector3D.dotProduct(position, velocity);

  1161.         // da
  1162.         final FieldVector3D<T> vectorAR = new FieldVector3D<>(a2.multiply(-2).divide(r3), position);
  1163.         final FieldVector3D<T> vectorARDot = velocity.scalarMultiply(a2.multiply(-2).divide(mu));
  1164.         fillHalfRow(getOne().negate(), vectorAR,    jacobian[0], 0);
  1165.         fillHalfRow(getOne().negate(), vectorARDot, jacobian[0], 3);

  1167.         // differentials of the momentum
  1168.         final T m      = momentum.getNorm();
  1169.         final T oOm    = m.reciprocal();
  1170.         final FieldVector3D<T> dcXP = new FieldVector3D<>(getZero(), vz, vy.negate());
  1171.         final FieldVector3D<T> dcYP = new FieldVector3D<>(vz.negate(),   getZero(),  vx);
  1172.         final FieldVector3D<T> dcZP = new FieldVector3D<>( vy, vx.negate(),   getZero());
  1173.         final FieldVector3D<T> dcXV = new FieldVector3D<>(  getZero(),  z.negate(),   y);
  1174.         final FieldVector3D<T> dcYV = new FieldVector3D<>(  z,   getZero(),  x.negate());
  1175.         final FieldVector3D<T> dcZV = new FieldVector3D<>( y.negate(),   x,   getZero());
  1176.         final FieldVector3D<T> dCP  = new FieldVector3D<>(mx.multiply(oOm), dcXP, my.multiply(oOm), dcYP, mz.multiply(oOm), dcZP);
  1177.         final FieldVector3D<T> dCV  = new FieldVector3D<>(mx.multiply(oOm), dcXV, my.multiply(oOm), dcYV, mz.multiply(oOm), dcZV);

  1179.         // dp
  1180.         final T mOMu   = m.divide(mu);
  1181.         final FieldVector3D<T> dpP  = new FieldVector3D<>(mOMu.multiply(2), dCP);
  1182.         final FieldVector3D<T> dpV  = new FieldVector3D<>(mOMu.multiply(2), dCV);

  1184.         // de
  1185.         final T p      = m.multiply(mOMu);
  1186.         final T moO2ae = absA.multiply(2).multiply(e).reciprocal();
  1187.         final T m2OaMu = p.negate().divide(absA);
  1188.         fillHalfRow(moO2ae, dpP, m2OaMu.multiply(moO2ae), vectorAR,    jacobian[1], 0);
  1189.         fillHalfRow(moO2ae, dpV, m2OaMu.multiply(moO2ae), vectorARDot, jacobian[1], 3);

  1191.         // di
  1192.         final T cI1 = m.multiply(sinI).reciprocal();
  1193.         final T cI2 = cosI.multiply(cI1);
  1194.         fillHalfRow(cI2, dCP, cI1.negate(), dcZP, jacobian[2], 0);
  1195.         fillHalfRow(cI2, dCV, cI1.negate(), dcZV, jacobian[2], 3);


  1198.         // dPA
  1199.         final T cP1     =  y.multiply(oOm);
  1200.         final T cP2     =  x.negate().multiply(oOm);
  1201.         final T cP3     =  mx.multiply(cP1).add(my.multiply(cP2)).negate();
  1202.         final T cP4     =  cP3.multiply(oOm);
  1203.         final T cP5     =  r2.multiply(sinI).multiply(sinI).negate().reciprocal();
  1204.         final T cP6     = z.multiply(cP5);
  1205.         final T cP7     = cP3.multiply(cP5);
  1206.         final FieldVector3D<T> dacP  = new FieldVector3D<>(cP1, dcXP, cP2, dcYP, cP4, dCP, oOm, new FieldVector3D<>(my.negate(), mx, getZero()));
  1207.         final FieldVector3D<T> dacV  = new FieldVector3D<>(cP1, dcXV, cP2, dcYV, cP4, dCV);
  1208.         final FieldVector3D<T> dpoP  = new FieldVector3D<>(cP6, dacP, cP7, this.PLUS_K);
  1209.         final FieldVector3D<T> dpoV  = new FieldVector3D<>(cP6, dacV);

  1211.         final T re2     = r2.multiply(e).multiply(e);
  1212.         final T recOre2 = p.subtract(r).divide(re2);
  1213.         final T resOre2 = pv.multiply(mOMu).divide(re2);
  1214.         final FieldVector3D<T> dreP  = new FieldVector3D<>(mOMu, velocity, pv.divide(mu), dCP);
  1215.         final FieldVector3D<T> dreV  = new FieldVector3D<>(mOMu, position, pv.divide(mu), dCV);
  1216.         final FieldVector3D<T> davP  = new FieldVector3D<>(resOre2.negate(), dpP, recOre2, dreP, resOre2.divide(r), position);
  1217.         final FieldVector3D<T> davV  = new FieldVector3D<>(resOre2.negate(), dpV, recOre2, dreV);
  1218.         fillHalfRow(getOne(), dpoP, getOne().negate(), davP, jacobian[3], 0);
  1219.         fillHalfRow(getOne(), dpoV, getOne().negate(), davV, jacobian[3], 3);

  1221.         // dRAAN
  1222.         final T cO0 = cI1.multiply(cI1);
  1223.         final T cO1 =  mx.multiply(cO0);
  1224.         final T cO2 =  my.negate().multiply(cO0);
  1225.         fillHalfRow(cO1, dcYP, cO2, dcXP, jacobian[4], 0);
  1226.         fillHalfRow(cO1, dcYV, cO2, dcXV, jacobian[4], 3);

  1228.         // dM
  1229.         final T s2a    = pv.divide(absA.multiply(2));
  1230.         final T oObux  = m.multiply(m).add(absA.multiply(mu)).sqrt().reciprocal();
  1231.         final T scasbu = pv.multiply(oObux);
  1232.         final FieldVector3D<T> dauP = new FieldVector3D<>(sqrtMuA.reciprocal(), velocity, s2a.negate().divide(sqrtMuA), vectorAR);
  1233.         final FieldVector3D<T> dauV = new FieldVector3D<>(sqrtMuA.reciprocal(), position, s2a.negate().divide(sqrtMuA), vectorARDot);
  1234.         final FieldVector3D<T> dbuP = new FieldVector3D<>(oObux.multiply(mu.divide(2.)), vectorAR,    m.multiply(oObux), dCP);
  1235.         final FieldVector3D<T> dbuV = new FieldVector3D<>(oObux.multiply(mu.divide(2.)), vectorARDot, m.multiply(oObux), dCV);
  1236.         final FieldVector3D<T> dcuP = new FieldVector3D<>(oObux, velocity, scasbu.negate().multiply(oObux), dbuP);
  1237.         final FieldVector3D<T> dcuV = new FieldVector3D<>(oObux, position, scasbu.negate().multiply(oObux), dbuV);
  1238.         fillHalfRow(getOne(), dauP, e.negate().divide(rOa.add(1)), dcuP, jacobian[5], 0);
  1239.         fillHalfRow(getOne(), dauV, e.negate().divide(rOa.add(1)), dcuV, jacobian[5], 3);

  1241.         return jacobian;

  1243.     }

  1245.     /** {@inheritDoc} */
  1246.     protected T[][] computeJacobianEccentricWrtCartesian() {
  1247.         if (isElliptical()) {
  1248.             return computeJacobianEccentricWrtCartesianElliptical();
  1249.         } else {
  1250.             return computeJacobianEccentricWrtCartesianHyperbolic();
  1251.         }
  1252.     }

  1254.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  1255.      * <p>
  1256.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  1257.      * respect to Cartesian coordinate j (x for j=0, y for j=1, z for j=2, xDot for j=3,
  1258.      * yDot for j=4, zDot for j=5).
  1259.      * </p>
  1260.      * @return 6x6 Jacobian matrix
  1261.      */
  1262.     private T[][] computeJacobianEccentricWrtCartesianElliptical() {

  1264.         // start by computing the Jacobian with mean angle
  1265.         final T[][] jacobian = computeJacobianMeanWrtCartesianElliptical();

  1267.         // Differentiating the Kepler equation M = E - e sin E leads to:
  1268.         // dM = (1 - e cos E) dE - sin E de
  1269.         // which is inverted and rewritten as:
  1270.         // dE = a/r dM + sin E a/r de
  1271.         final FieldSinCos<T> scE = FastMath.sinCos(getEccentricAnomaly());
  1272.         final T aOr              = e.negate().multiply(scE.cos()).add(1).reciprocal();

  1274.         // update anomaly row
  1275.         final T[] eRow           = jacobian[1];
  1276.         final T[] anomalyRow     = jacobian[5];
  1277.         for (int j = 0; j < anomalyRow.length; ++j) {
  1278.             anomalyRow[j] = aOr.multiply(anomalyRow[j].add(scE.sin().multiply(eRow[j])));
  1279.         }

  1281.         return jacobian;

  1283.     }

  1285.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  1286.      * <p>
  1287.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  1288.      * respect to Cartesian coordinate j (x for j=0, y for j=1, z for j=2, xDot for j=3,
  1289.      * yDot for j=4, zDot for j=5).
  1290.      * </p>
  1291.      * @return 6x6 Jacobian matrix
  1292.      */
  1293.     private T[][] computeJacobianEccentricWrtCartesianHyperbolic() {

  1295.         // start by computing the Jacobian with mean angle
  1296.         final T[][] jacobian = computeJacobianMeanWrtCartesianHyperbolic();

  1298.         // Differentiating the Kepler equation M = e sinh H - H leads to:
  1299.         // dM = (e cosh H - 1) dH + sinh H de
  1300.         // which is inverted and rewritten as:
  1301.         // dH = 1 / (e cosh H - 1) dM - sinh H / (e cosh H - 1) de
  1302.         final T H      = getEccentricAnomaly();
  1303.         final T coshH  = H.cosh();
  1304.         final T sinhH  = H.sinh();
  1305.         final T absaOr = e.multiply(coshH).subtract(1).reciprocal();
  1306.         // update anomaly row
  1307.         final T[] eRow       = jacobian[1];
  1308.         final T[] anomalyRow = jacobian[5];

  1310.         for (int j = 0; j < anomalyRow.length; ++j) {
  1311.             anomalyRow[j] = absaOr.multiply(anomalyRow[j].subtract(sinhH.multiply(eRow[j])));

  1313.         }

  1315.         return jacobian;

  1317.     }

  1319.     /** {@inheritDoc} */
  1320.     protected T[][] computeJacobianTrueWrtCartesian() {
  1321.         if (isElliptical()) {
  1322.             return computeJacobianTrueWrtCartesianElliptical();
  1323.         } else {
  1324.             return computeJacobianTrueWrtCartesianHyperbolic();
  1325.         }
  1326.     }

  1328.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  1329.      * <p>
  1330.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  1331.      * respect to Cartesian coordinate j (x for j=0, y for j=1, z for j=2, xDot for j=3,
  1332.      * yDot for j=4, zDot for j=5).
  1333.      * </p>
  1334.      * @return 6x6 Jacobian matrix
  1335.      */
  1336.     private T[][] computeJacobianTrueWrtCartesianElliptical() {

  1338.         // start by computing the Jacobian with eccentric angle
  1339.         final T[][] jacobian = computeJacobianEccentricWrtCartesianElliptical();
  1340.         // Differentiating the eccentric anomaly equation sin E = sqrt(1-e^2) sin v / (1 + e cos v)
  1341.         // and using cos E = (e + cos v) / (1 + e cos v) to get rid of cos E leads to:
  1342.         // dE = [sqrt (1 - e^2) / (1 + e cos v)] dv - [sin E / (1 - e^2)] de
  1343.         // which is inverted and rewritten as:
  1344.         // dv = sqrt (1 - e^2) a/r dE + [sin E / sqrt (1 - e^2)] a/r de
  1345.         final T e2               = e.multiply(e);
  1346.         final T oMe2             = e2.negate().add(1);
  1347.         final T epsilon          = oMe2.sqrt();
  1348.         final FieldSinCos<T> scE = FastMath.sinCos(getEccentricAnomaly());
  1349.         final T aOr              = e.multiply(scE.cos()).negate().add(1).reciprocal();
  1350.         final T aFactor          = epsilon.multiply(aOr);
  1351.         final T eFactor          = scE.sin().multiply(aOr).divide(epsilon);

  1353.         // update anomaly row
  1354.         final T[] eRow           = jacobian[1];
  1355.         final T[] anomalyRow     = jacobian[5];
  1356.         for (int j = 0; j < anomalyRow.length; ++j) {
  1357.             anomalyRow[j] = aFactor.multiply(anomalyRow[j]).add(eFactor.multiply(eRow[j]));
  1358.         }
  1359.         return jacobian;

  1361.     }

  1363.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  1364.      * <p>
  1365.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  1366.      * respect to Cartesian coordinate j (x for j=0, y for j=1, z for j=2, xDot for j=3,
  1367.      * yDot for j=4, zDot for j=5).
  1368.      * </p>
  1369.      * @return 6x6 Jacobian matrix
  1370.      */
  1371.     private T[][] computeJacobianTrueWrtCartesianHyperbolic() {

  1373.         // start by computing the Jacobian with eccentric angle
  1374.         final T[][] jacobian = computeJacobianEccentricWrtCartesianHyperbolic();

  1376.         // Differentiating the eccentric anomaly equation sinh H = sqrt(e^2-1) sin v / (1 + e cos v)
  1377.         // and using cosh H = (e + cos v) / (1 + e cos v) to get rid of cosh H leads to:
  1378.         // dH = [sqrt (e^2 - 1) / (1 + e cos v)] dv + [sinh H / (e^2 - 1)] de
  1379.         // which is inverted and rewritten as:
  1380.         // dv = sqrt (1 - e^2) a/r dH - [sinh H / sqrt (e^2 - 1)] a/r de
  1381.         final T e2       = e.multiply(e);
  1382.         final T e2Mo     = e2.subtract(1);
  1383.         final T epsilon  = e2Mo.sqrt();
  1384.         final T H        = getEccentricAnomaly();
  1385.         final T coshH    = H.cosh();
  1386.         final T sinhH    = H.sinh();
  1387.         final T aOr      = e.multiply(coshH).subtract(1).reciprocal();
  1388.         final T aFactor  = epsilon.multiply(aOr);
  1389.         final T eFactor  = sinhH .multiply(aOr).divide(epsilon);

  1391.         // update anomaly row
  1392.         final T[] eRow           = jacobian[1];
  1393.         final T[] anomalyRow     = jacobian[5];
  1394.         for (int j = 0; j < anomalyRow.length; ++j) {
  1395.             anomalyRow[j] = aFactor.multiply(anomalyRow[j]).subtract(eFactor.multiply(eRow[j]));
  1396.         }

  1398.         return jacobian;

  1400.     }

  1402.     /** {@inheritDoc} */
  1403.     public void addKeplerContribution(final PositionAngleType type, final T gm,
  1404.                                       final T[] pDot) {
  1405.         final T oMe2;
  1406.         final T ksi;
  1407.         final T absA = a.abs();
  1408.         final T n    = absA.reciprocal().multiply(gm).sqrt().divide(absA);
  1409.         switch (type) {
  1410.             case MEAN :
  1411.                 pDot[5] = pDot[5].add(n);
  1412.                 break;
  1413.             case ECCENTRIC :
  1414.                 oMe2 = e.multiply(e).negate().add(1).abs();
  1415.                 ksi  = e.multiply(v.cos()).add(1);
  1416.                 pDot[5] = pDot[5].add( n.multiply(ksi).divide(oMe2));
  1417.                 break;
  1418.             case TRUE :
  1419.                 oMe2 = e.multiply(e).negate().add(1).abs();
  1420.                 ksi  = e.multiply(v.cos()).add(1);
  1421.                 pDot[5] = pDot[5].add(n.multiply(ksi).multiply(ksi).divide(oMe2.multiply(oMe2.sqrt())));
  1422.                 break;
  1423.             default :
  1424.                 throw new OrekitInternalError(null);
  1425.         }
  1426.     }

  1428.     /**  Returns a string representation of this Keplerian parameters object.
  1429.      * @return a string representation of this object
  1430.      */
  1431.     public String toString() {
  1432.         return new StringBuilder().append("Keplerian parameters: ").append('{').
  1433.                                   append("a: ").append(a.getReal()).
  1434.                                   append("; e: ").append(e.getReal()).
  1435.                                   append("; i: ").append(FastMath.toDegrees(i.getReal())).
  1436.                                   append("; pa: ").append(FastMath.toDegrees(pa.getReal())).
  1437.                                   append("; raan: ").append(FastMath.toDegrees(raan.getReal())).
  1438.                                   append("; v: ").append(FastMath.toDegrees(v.getReal())).
  1439.                                   append(";}").toString();
  1440.     }

  1442.     /** {@inheritDoc} */
  1443.     @Override
  1444.     public PositionAngleType getCachedPositionAngleType() {
  1445.         return PositionAngleType.TRUE;
  1446.     }

  1448.     /** {@inheritDoc} */
  1449.     @Override
  1450.     public boolean hasRates() {
  1451.         return hasDerivatives();
  1452.     }

  1454.     /** {@inheritDoc} */
  1455.     @Override
  1456.     public FieldKeplerianOrbit<T> removeRates() {
  1457.         final PositionAngleType positionAngleType = getCachedPositionAngleType();
  1458.         return new FieldKeplerianOrbit<>(getA(), getE(), getI(), getPerigeeArgument(), getRightAscensionOfAscendingNode(),
  1459.                 getAnomaly(positionAngleType), positionAngleType, getFrame(), getDate(), getMu());
  1460.     }


  1463.     /** Check if the given parameter is within an acceptable range.
  1464.      * The bounds are inclusive: an exception is raised when either of those conditions are met:
  1465.      * <ul>
  1466.      *     <li>The parameter is strictly greater than upperBound</li>
  1467.      *     <li>The parameter is strictly lower than lowerBound</li>
  1468.      * </ul>
  1469.      * <p>
  1470.      * In either of these cases, an OrekitException is raised.
  1471.      * </p>
  1472.      * @param parameterName name of the parameter
  1473.      * @param parameter value of the parameter
  1474.      * @param lowerBound lower bound of the acceptable range (inclusive)
  1475.      * @param upperBound upper bound of the acceptable range (inclusive)
  1476.      */
  1477.     private void checkParameterRangeInclusive(final String parameterName, final double parameter,
  1478.                                               final double lowerBound, final double upperBound) {
  1479.         if (parameter < lowerBound || parameter > upperBound) {
  1480.             throw new OrekitException(OrekitMessages.INVALID_PARAMETER_RANGE, parameterName,
  1481.                                       parameter, lowerBound, upperBound);
  1482.         }
  1483.     }

  1485.     /** {@inheritDoc} */
  1486.     @Override
  1487.     public KeplerianOrbit toOrbit() {
  1488.         if (hasDerivatives()) {
  1489.             return new KeplerianOrbit(a.getReal(), e.getReal(), i.getReal(),
  1490.                                       pa.getReal(), raan.getReal(), v.getReal(),
  1491.                                       aDot.getReal(), eDot.getReal(), iDot.getReal(),
  1492.                                       paDot.getReal(), raanDot.getReal(), vDot.getReal(),
  1493.                                       PositionAngleType.TRUE,
  1494.                                       getFrame(), getDate().toAbsoluteDate(), getMu().getReal());
  1495.         } else {
  1496.             return new KeplerianOrbit(a.getReal(), e.getReal(), i.getReal(),
  1497.                                       pa.getReal(), raan.getReal(), v.getReal(),
  1498.                                       PositionAngleType.TRUE,
  1499.                                       getFrame(), getDate().toAbsoluteDate(), getMu().getReal());
  1500.         }
  1501.     }


  1504. }
Created with JaCoCo 0.8.11.202310140853