FieldCircularOrbit.java

  1. /* Copyright 2002-2024 CS GROUP
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.orbits;

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


  35. /**
  36.  * This class handles circular orbital parameters.

  38.  * <p>
  39.  * The parameters used internally are the circular elements which can be
  40.  * related to Keplerian elements as follows:
  41.  *   <ul>
  42.  *     <li>a</li>
  43.  *     <li>e<sub>x</sub> = e cos(ω)</li>
  44.  *     <li>e<sub>y</sub> = e sin(ω)</li>
  45.  *     <li>i</li>
  46.  *     <li>Ω</li>
  47.  *     <li>α<sub>v</sub> = v + ω</li>
  48.  *   </ul>
  49.  * where Ω stands for the Right Ascension of the Ascending Node and
  50.  * α<sub>v</sub> stands for the true latitude argument
  51.  * <p>
  52.  * The conversion equations from and to Keplerian elements given above hold only
  53.  * when both sides are unambiguously defined, i.e. when orbit is neither equatorial
  54.  * nor circular. When orbit is circular (but not equatorial), the circular
  55.  * parameters are still unambiguously defined whereas some Keplerian elements
  56.  * (more precisely ω and Ω) become ambiguous. When orbit is equatorial,
  57.  * neither the Keplerian nor the circular parameters can be defined unambiguously.
  58.  * {@link EquinoctialOrbit equinoctial orbits} is the recommended way to represent
  59.  * orbits.
  60.  * </p>
  61.  * <p>
  62.  * The instance <code>CircularOrbit</code> is guaranteed to be immutable.
  63.  * </p>
  64.  * @see    Orbit
  65.  * @see    KeplerianOrbit
  66.  * @see    CartesianOrbit
  67.  * @see    EquinoctialOrbit
  68.  * @author Luc Maisonobe
  69.  * @author Fabien Maussion
  70.  * @author V&eacute;ronique Pommier-Maurussane
  71.  * @since 9.0
  72.  * @param <T> type of the field elements
  73.  */

  75. public class FieldCircularOrbit<T extends CalculusFieldElement<T>> extends FieldOrbit<T> implements PositionAngleBased {

  77.     /** Semi-major axis (m). */
  78.     private final T a;

  80.     /** First component of the circular eccentricity vector. */
  81.     private final T ex;

  83.     /** Second component of the circular eccentricity vector. */
  84.     private final T ey;

  86.     /** Inclination (rad). */
  87.     private final T i;

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

  92.     /** Cached latitude argument (rad). */
  93.     private final T cachedAlpha;

  95.     /** Type of cached position angle (latitude argument). */
  96.     private final PositionAngleType cachedPositionAngleType;

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

  101.     /** First component of the circular eccentricity vector derivative. */
  102.     private final T exDot;

  104.     /** Second component of the circular eccentricity vector derivative. */
  105.     private final T eyDot;

  107.     /** Inclination derivative (rad/s). */
  108.     private final T iDot;

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

  113.     /** True latitude argument derivative (rad/s). */
  114.     private final T cachedAlphaDot;

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

  119.     /** Creates a new instance.
  120.      * @param a  semi-major axis (m)
  121.      * @param ex e cos(ω), first component of circular eccentricity vector
  122.      * @param ey e sin(ω), second component of circular eccentricity vector
  123.      * @param i inclination (rad)
  124.      * @param raan right ascension of ascending node (Ω, rad)
  125.      * @param alpha  an + ω, mean, eccentric or true latitude argument (rad)
  126.      * @param type type of latitude argument
  127.      * @param cachedPositionAngleType type of cached latitude argument
  128.      * @param frame the frame in which are defined the parameters
  129.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  130.      * @param date date of the orbital parameters
  131.      * @param mu central attraction coefficient (m³/s²)
  132.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  133.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  134.      * @since 12.1
  135.      */
  136.     public FieldCircularOrbit(final T a, final T ex, final T ey, final T i, final T raan,
  137.                               final T alpha, final PositionAngleType type,
  138.                               final PositionAngleType cachedPositionAngleType,
  139.                               final Frame frame, final FieldAbsoluteDate<T> date, final T mu)
  140.         throws IllegalArgumentException {
  141.         this(a, ex, ey, i, raan, alpha,
  142.              null, null, null, null, null, null,
  143.              type, cachedPositionAngleType, frame, date, mu);
  144.     }

  146.     /** Creates a new instance without derivatives and with cached position angle same as value inputted.
  147.      * @param a  semi-major axis (m)
  148.      * @param ex e cos(ω), first component of circular eccentricity vector
  149.      * @param ey e sin(ω), second component of circular eccentricity vector
  150.      * @param i inclination (rad)
  151.      * @param raan right ascension of ascending node (Ω, rad)
  152.      * @param alpha  an + ω, mean, eccentric or true latitude argument (rad)
  153.      * @param type type of latitude argument
  154.      * @param frame the frame in which are defined the parameters
  155.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  156.      * @param date date of the orbital parameters
  157.      * @param mu central attraction coefficient (m³/s²)
  158.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  159.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  160.      */
  161.     public FieldCircularOrbit(final T a, final T ex, final T ey, final T i, final T raan,
  162.                               final T alpha, final PositionAngleType type,
  163.                               final Frame frame, final FieldAbsoluteDate<T> date, final T mu)
  164.             throws IllegalArgumentException {
  165.         this(a, ex, ey, i, raan, alpha, type, type, frame, date, mu);
  166.     }

  168.     /** Creates a new instance.
  169.      * @param a  semi-major axis (m)
  170.      * @param ex e cos(ω), first component of circular eccentricity vector
  171.      * @param ey e sin(ω), second component of circular eccentricity vector
  172.      * @param i inclination (rad)
  173.      * @param raan right ascension of ascending node (Ω, rad)
  174.      * @param alpha  an + ω, mean, eccentric or true latitude argument (rad)
  175.      * @param aDot  semi-major axis derivative (m/s)
  176.      * @param exDot d(e cos(ω))/dt, first component of circular eccentricity vector derivative
  177.      * @param eyDot d(e sin(ω))/dt, second component of circular eccentricity vector derivative
  178.      * @param iDot inclination  derivative(rad/s)
  179.      * @param raanDot right ascension of ascending node derivative (rad/s)
  180.      * @param alphaDot  d(an + ω), mean, eccentric or true latitude argument derivative (rad/s)
  181.      * @param type type of latitude argument
  182.      * @param frame the frame in which are defined the parameters
  183.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  184.      * @param date date of the orbital parameters
  185.      * @param mu central attraction coefficient (m³/s²)
  186.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  187.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  188.      */
  189.     public FieldCircularOrbit(final T a, final T ex, final T ey,
  190.                               final T i, final T raan, final T alpha,
  191.                               final T aDot, final T exDot, final T eyDot,
  192.                               final T iDot, final T raanDot, final T alphaDot, final PositionAngleType type,
  193.                               final Frame frame, final FieldAbsoluteDate<T> date, final T mu)
  194.             throws IllegalArgumentException {
  195.         this(a, ex, ey, i, raan, alpha, aDot, exDot, eyDot, iDot, raanDot, alphaDot, type, type, frame, date, mu);
  196.     }

  198.     /** Creates a new instance.
  199.      * @param a  semi-major axis (m)
  200.      * @param ex e cos(ω), first component of circular eccentricity vector
  201.      * @param ey e sin(ω), second component of circular eccentricity vector
  202.      * @param i inclination (rad)
  203.      * @param raan right ascension of ascending node (Ω, rad)
  204.      * @param alpha  an + ω, mean, eccentric or true latitude argument (rad)
  205.      * @param aDot  semi-major axis derivative (m/s)
  206.      * @param exDot d(e cos(ω))/dt, first component of circular eccentricity vector derivative
  207.      * @param eyDot d(e sin(ω))/dt, second component of circular eccentricity vector derivative
  208.      * @param iDot inclination  derivative(rad/s)
  209.      * @param raanDot right ascension of ascending node derivative (rad/s)
  210.      * @param alphaDot  d(an + ω), mean, eccentric or true latitude argument derivative (rad/s)
  211.      * @param type type of latitude argument
  212.      * @param cachedPositionAngleType type of cached latitude argument
  213.      * @param frame the frame in which are defined the parameters
  214.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  215.      * @param date date of the orbital parameters
  216.      * @param mu central attraction coefficient (m³/s²)
  217.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  218.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  219.      * @since 12.1
  220.      */
  221.     public FieldCircularOrbit(final T a, final T ex, final T ey,
  222.                               final T i, final T raan, final T alpha,
  223.                               final T aDot, final T exDot, final T eyDot,
  224.                               final T iDot, final T raanDot, final T alphaDot,
  225.                               final PositionAngleType type, final PositionAngleType cachedPositionAngleType,
  226.                               final Frame frame, final FieldAbsoluteDate<T> date, final T mu)
  227.         throws IllegalArgumentException {
  228.         super(frame, date, mu);
  229.         if (ex.getReal() * ex.getReal() + ey.getReal() * ey.getReal() >= 1.0) {
  230.             throw new OrekitIllegalArgumentException(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS,
  231.                                                      getClass().getName());
  232.         }

  234.         this.a       =  a;
  235.         this.aDot    =  aDot;
  236.         this.ex      = ex;
  237.         this.exDot   = exDot;
  238.         this.ey      = ey;
  239.         this.eyDot   = eyDot;
  240.         this.i       = i;
  241.         this.iDot    = iDot;
  242.         this.raan    = raan;
  243.         this.raanDot = raanDot;
  244.         this.cachedPositionAngleType = cachedPositionAngleType;

  246.         if (hasDerivatives()) {
  247.             final FieldUnivariateDerivative1<T> alphaUD = initializeCachedAlpha(alpha, alphaDot, type);
  248.             this.cachedAlpha = alphaUD.getValue();
  249.             this.cachedAlphaDot = alphaUD.getFirstDerivative();
  250.         } else {
  251.             this.cachedAlpha = initializeCachedAlpha(alpha, type);
  252.             this.cachedAlphaDot = null;
  253.         }

  255.         partialPV   = null;

  257.     }

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

  279.         // compute semi-major axis
  280.         final FieldVector3D<T> pvP = pvCoordinates.getPosition();
  281.         final FieldVector3D<T> pvV = pvCoordinates.getVelocity();
  282.         final FieldVector3D<T> pvA = pvCoordinates.getAcceleration();
  283.         final T r2 = pvP.getNormSq();
  284.         final T r  = r2.sqrt();
  285.         final T V2 = pvV.getNormSq();
  286.         final T rV2OnMu = r.multiply(V2).divide(mu);

  288.         a = r.divide(rV2OnMu.negate().add(2));

  290.         if (!isElliptical()) {
  291.             throw new OrekitIllegalArgumentException(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS,
  292.                                                      getClass().getName());
  293.         }

  295.         // compute inclination
  296.         final FieldVector3D<T> momentum = pvCoordinates.getMomentum();
  297.         final FieldVector3D<T> plusK = FieldVector3D.getPlusK(r.getField());
  298.         i = FieldVector3D.angle(momentum, plusK);

  300.         // compute right ascension of ascending node
  301.         final FieldVector3D<T> node  = FieldVector3D.crossProduct(plusK, momentum);
  302.         raan = node.getY().atan2(node.getX());

  304.         // 2D-coordinates in the canonical frame
  305.         final FieldSinCos<T> scRaan = FastMath.sinCos(raan);
  306.         final FieldSinCos<T> scI    = FastMath.sinCos(i);
  307.         final T xP      = pvP.getX();
  308.         final T yP      = pvP.getY();
  309.         final T zP      = pvP.getZ();
  310.         final T x2      = (xP.multiply(scRaan.cos()).add(yP .multiply(scRaan.sin()))).divide(a);
  311.         final T y2      = (yP.multiply(scRaan.cos()).subtract(xP.multiply(scRaan.sin()))).multiply(scI.cos()).add(zP.multiply(scI.sin())).divide(a);

  313.         // compute eccentricity vector
  314.         final T eSE    = FieldVector3D.dotProduct(pvP, pvV).divide(a.multiply(mu).sqrt());
  315.         final T eCE    = rV2OnMu.subtract(1);
  316.         final T e2     = eCE.multiply(eCE).add(eSE.multiply(eSE));
  317.         final T f      = eCE.subtract(e2);
  318.         final T g      = eSE.multiply(e2.negate().add(1).sqrt());
  319.         final T aOnR   = a.divide(r);
  320.         final T a2OnR2 = aOnR.square();
  321.         ex = a2OnR2.multiply(f.multiply(x2).add(g.multiply(y2)));
  322.         ey = a2OnR2.multiply(f.multiply(y2).subtract(g.multiply(x2)));

  324.         // compute latitude argument
  325.         final T beta = (ex.multiply(ex).add(ey.multiply(ey)).negate().add(1)).sqrt().add(1).reciprocal();
  326.         cachedAlpha = FieldCircularLatitudeArgumentUtility.eccentricToTrue(ex, ey, y2.add(ey).add(eSE.multiply(beta).multiply(ex)).atan2(x2.add(ex).subtract(eSE.multiply(beta).multiply(ey)))
  327.         );

  329.         partialPV = pvCoordinates;

  331.         if (hasNonKeplerianAcceleration(pvCoordinates, mu)) {
  332.             // we have a relevant acceleration, we can compute derivatives

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

  337.             final FieldVector3D<T> keplerianAcceleration    = new FieldVector3D<>(r.multiply(r2).reciprocal().multiply(mu.negate()), pvP);
  338.             final FieldVector3D<T> nonKeplerianAcceleration = pvA.subtract(keplerianAcceleration);
  339.             final T   aX                       = nonKeplerianAcceleration.getX();
  340.             final T   aY                       = nonKeplerianAcceleration.getY();
  341.             final T   aZ                       = nonKeplerianAcceleration.getZ();
  342.             aDot    = jacobian[0][3].multiply(aX).add(jacobian[0][4].multiply(aY)).add(jacobian[0][5].multiply(aZ));
  343.             exDot   = jacobian[1][3].multiply(aX).add(jacobian[1][4].multiply(aY)).add(jacobian[1][5].multiply(aZ));
  344.             eyDot   = jacobian[2][3].multiply(aX).add(jacobian[2][4].multiply(aY)).add(jacobian[2][5].multiply(aZ));
  345.             iDot    = jacobian[3][3].multiply(aX).add(jacobian[3][4].multiply(aY)).add(jacobian[3][5].multiply(aZ));
  346.             raanDot = jacobian[4][3].multiply(aX).add(jacobian[4][4].multiply(aY)).add(jacobian[4][5].multiply(aZ));

  348.             // in order to compute true anomaly derivative, we must compute
  349.             // mean anomaly derivative including Keplerian motion and convert to true anomaly
  350.             final T alphaMDot = getKeplerianMeanMotion().
  351.                                 add(jacobian[5][3].multiply(aX)).add(jacobian[5][4].multiply(aY)).add(jacobian[5][5].multiply(aZ));
  352.             final FieldUnivariateDerivative1<T> exUD     = new FieldUnivariateDerivative1<>(ex, exDot);
  353.             final FieldUnivariateDerivative1<T> eyUD     = new FieldUnivariateDerivative1<>(ey, eyDot);
  354.             final FieldUnivariateDerivative1<T> alphaMUD = new FieldUnivariateDerivative1<>(getAlphaM(), alphaMDot);
  355.             final FieldUnivariateDerivative1<T> alphavUD = FieldCircularLatitudeArgumentUtility.meanToTrue(exUD, eyUD, alphaMUD);
  356.             cachedAlphaDot = alphavUD.getFirstDerivative();

  358.         } else {
  359.             // acceleration is either almost zero or NaN,
  360.             // we assume acceleration was not known
  361.             // we don't set up derivatives
  362.             aDot      = null;
  363.             exDot     = null;
  364.             eyDot     = null;
  365.             iDot      = null;
  366.             raanDot   = null;
  367.             cachedAlphaDot = null;
  368.         }

  370.     }

  372.     /** Constructor from Cartesian parameters.
  373.      *
  374.      * <p> The acceleration provided in {@code FieldPVCoordinates} is accessible using
  375.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  376.      * use {@code mu} and the position to compute the acceleration, including
  377.      * {@link #shiftedBy(CalculusFieldElement)} and {@link #getPVCoordinates(FieldAbsoluteDate, Frame)}.
  378.      *
  379.      * @param PVCoordinates the {@link FieldPVCoordinates} in inertial frame
  380.      * @param frame the frame in which are defined the {@link FieldPVCoordinates}
  381.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  382.      * @param date date of the orbital parameters
  383.      * @param mu central attraction coefficient (m³/s²)
  384.      * @exception IllegalArgumentException if frame is not a {@link
  385.      * Frame#isPseudoInertial pseudo-inertial frame}
  386.      */
  387.     public FieldCircularOrbit(final FieldPVCoordinates<T> PVCoordinates, final Frame frame,
  388.                               final FieldAbsoluteDate<T> date, final T mu)
  389.         throws IllegalArgumentException {
  390.         this(new TimeStampedFieldPVCoordinates<>(date, PVCoordinates), frame, mu);
  391.     }

  393.     /** Constructor from any kind of orbital parameters.
  394.      * @param op orbital parameters to copy
  395.      */
  396.     public FieldCircularOrbit(final FieldOrbit<T> op) {
  397.         super(op.getFrame(), op.getDate(), op.getMu());
  398.         a    = op.getA();
  399.         i    = op.getI();
  400.         final T hx = op.getHx();
  401.         final T hy = op.getHy();
  402.         final T h2 = hx.square().add(hy.square());
  403.         final T h  = h2.sqrt();
  404.         raan = hy.atan2(hx);
  405.         final FieldSinCos<T> scRaan = FastMath.sinCos(raan);
  406.         final T cosRaan = h.getReal() == 0 ? scRaan.cos() : hx.divide(h);
  407.         final T sinRaan = h.getReal() == 0 ? scRaan.sin() : hy.divide(h);
  408.         final T equiEx = op.getEquinoctialEx();
  409.         final T equiEy = op.getEquinoctialEy();
  410.         ex   = equiEx.multiply(cosRaan).add(equiEy.multiply(sinRaan));
  411.         ey   = equiEy.multiply(cosRaan).subtract(equiEx.multiply(sinRaan));
  412.         cachedPositionAngleType = PositionAngleType.TRUE;
  413.         cachedAlpha = op.getLv().subtract(raan);

  415.         if (op.hasDerivatives()) {
  416.             aDot      = op.getADot();
  417.             final T      hxDot = op.getHxDot();
  418.             final T      hyDot = op.getHyDot();
  419.             iDot      = cosRaan.multiply(hxDot).add(sinRaan.multiply(hyDot)).multiply(2).divide(h2.add(1));
  420.             raanDot   = hx.multiply(hyDot).subtract(hy.multiply(hxDot)).divide(h2);
  421.             final T equiExDot = op.getEquinoctialExDot();
  422.             final T equiEyDot = op.getEquinoctialEyDot();
  423.             exDot     = equiExDot.add(equiEy.multiply(raanDot)).multiply(cosRaan).
  424.                         add(equiEyDot.subtract(equiEx.multiply(raanDot)).multiply(sinRaan));
  425.             eyDot     = equiEyDot.subtract(equiEx.multiply(raanDot)).multiply(cosRaan).
  426.                         subtract(equiExDot.add(equiEy.multiply(raanDot)).multiply(sinRaan));
  427.             cachedAlphaDot = op.getLvDot().subtract(raanDot);
  428.         } else {
  429.             aDot      = null;
  430.             exDot     = null;
  431.             eyDot     = null;
  432.             iDot      = null;
  433.             raanDot   = null;
  434.             cachedAlphaDot = null;
  435.         }

  437.         partialPV = null;

  439.     }

  441.     /** Constructor from Field and CircularOrbit.
  442.      * <p>Build a FieldCircularOrbit from non-Field CircularOrbit.</p>
  443.      * @param field CalculusField to base object on
  444.      * @param op non-field orbit with only "constant" terms
  445.      * @since 12.0
  446.      */
  447.     public FieldCircularOrbit(final Field<T> field, final CircularOrbit op) {
  448.         super(op.getFrame(), new FieldAbsoluteDate<>(field, op.getDate()), field.getZero().newInstance(op.getMu()));

  450.         a    = getZero().newInstance(op.getA());
  451.         i    = getZero().newInstance(op.getI());
  452.         raan = getZero().newInstance(op.getRightAscensionOfAscendingNode());
  453.         ex   = getZero().newInstance(op.getCircularEx());
  454.         ey   = getZero().newInstance(op.getCircularEy());
  455.         cachedPositionAngleType = op.getCachedPositionAngleType();
  456.         cachedAlpha = getZero().newInstance(op.getAlpha(cachedPositionAngleType));

  458.         if (op.hasDerivatives()) {
  459.             aDot      = getZero().newInstance(op.getADot());
  460.             iDot      = getZero().newInstance(op.getIDot());
  461.             raanDot   = getZero().newInstance(op.getRightAscensionOfAscendingNodeDot());
  462.             exDot     = getZero().newInstance(op.getCircularExDot());
  463.             eyDot     = getZero().newInstance(op.getCircularEyDot());
  464.             cachedAlphaDot = getZero().newInstance(op.getAlphaDot(cachedPositionAngleType));
  465.         } else {
  466.             aDot      = null;
  467.             exDot     = null;
  468.             eyDot     = null;
  469.             iDot      = null;
  470.             raanDot   = null;
  471.             cachedAlphaDot = null;
  472.         }

  474.         partialPV = null;

  476.     }

  478.     /** Constructor from Field and Orbit.
  479.      * <p>Build a FieldCircularOrbit from any non-Field Orbit.</p>
  480.      * @param field CalculusField to base object on
  481.      * @param op non-field orbit with only "constant" terms
  482.      * @since 12.0
  483.      */
  484.     public FieldCircularOrbit(final Field<T> field, final Orbit op) {
  485.         this(field, (CircularOrbit) OrbitType.CIRCULAR.convertType(op));
  486.     }

  488.     /** {@inheritDoc} */
  489.     @Override
  490.     public OrbitType getType() {
  491.         return OrbitType.CIRCULAR;
  492.     }

  494.     /** {@inheritDoc} */
  495.     @Override
  496.     public T getA() {
  497.         return a;
  498.     }

  500.     /** {@inheritDoc} */
  501.     @Override
  502.     public T getADot() {
  503.         return aDot;
  504.     }

  506.     /** {@inheritDoc} */
  507.     @Override
  508.     public T getEquinoctialEx() {
  509.         final FieldSinCos<T> sc = FastMath.sinCos(raan);
  510.         return ex.multiply(sc.cos()).subtract(ey.multiply(sc.sin()));
  511.     }

  513.     /** {@inheritDoc} */
  514.     @Override
  515.     public T getEquinoctialExDot() {

  517.         if (!hasDerivatives()) {
  518.             return null;
  519.         }

  521.         final FieldSinCos<T> sc = FastMath.sinCos(raan);
  522.         return exDot.subtract(ey.multiply(raanDot)).multiply(sc.cos()).
  523.                subtract(eyDot.add(ex.multiply(raanDot)).multiply(sc.sin()));

  525.     }

  527.     /** {@inheritDoc} */
  528.     @Override
  529.     public T getEquinoctialEy() {
  530.         final FieldSinCos<T> sc = FastMath.sinCos(raan);
  531.         return ey.multiply(sc.cos()).add(ex.multiply(sc.sin()));
  532.     }

  534.     /** {@inheritDoc} */
  535.     @Override
  536.     public T getEquinoctialEyDot() {

  538.         if (!hasDerivatives()) {
  539.             return null;
  540.         }

  542.         final FieldSinCos<T> sc = FastMath.sinCos(raan);
  543.         return eyDot.add(ex.multiply(raanDot)).multiply(sc.cos()).
  544.                add(exDot.subtract(ey.multiply(raanDot)).multiply(sc.sin()));

  546.     }

  548.     /** Get the first component of the circular eccentricity vector.
  549.      * @return ex = e cos(ω), first component of the circular eccentricity vector
  550.      */
  551.     public T getCircularEx() {
  552.         return ex;
  553.     }

  555.     /** Get the first component of the circular eccentricity vector derivative.
  556.      * @return d(ex)/dt = d(e cos(ω))/dt, first component of the circular eccentricity vector derivative
  557.      */
  558.     public T getCircularExDot() {
  559.         return exDot;
  560.     }

  562.     /** Get the second component of the circular eccentricity vector.
  563.      * @return ey = e sin(ω), second component of the circular eccentricity vector
  564.      */
  565.     public T getCircularEy() {
  566.         return ey;
  567.     }

  569.     /** Get the second component of the circular eccentricity vector derivative.
  570.      * @return d(ey)/dt = d(e sin(ω))/dt, second component of the circular eccentricity vector derivative
  571.      */
  572.     public T getCircularEyDot() {
  573.         return eyDot;
  574.     }

  576.     /** {@inheritDoc} */
  577.     @Override
  578.     public T getHx() {
  579.         // Check for equatorial retrograde orbit
  580.         if (FastMath.abs(i.subtract(i.getPi()).getReal()) < 1.0e-10) {
  581.             return getZero().add(Double.NaN);
  582.         }
  583.         return raan.cos().multiply(i.divide(2).tan());
  584.     }

  586.     /** {@inheritDoc} */
  587.     @Override
  588.     public T getHxDot() {

  590.         if (!hasDerivatives()) {
  591.             return null;
  592.         }

  594.         // Check for equatorial retrograde orbit
  595.         if (FastMath.abs(i.subtract(i.getPi()).getReal()) < 1.0e-10) {
  596.             return getZero().add(Double.NaN);
  597.         }

  599.         final FieldSinCos<T> sc = FastMath.sinCos(raan);
  600.         final T tan             = i.multiply(0.5).tan();
  601.         return sc.cos().multiply(0.5).multiply(tan.multiply(tan).add(1)).multiply(iDot).
  602.                subtract(sc.sin().multiply(tan).multiply(raanDot));

  604.     }

  606.     /** {@inheritDoc} */
  607.     @Override
  608.     public T getHy() {
  609.         // Check for equatorial retrograde orbit
  610.         if (FastMath.abs(i.subtract(i.getPi()).getReal()) < 1.0e-10) {
  611.             return getZero().add(Double.NaN);
  612.         }
  613.         return raan.sin().multiply(i.divide(2).tan());
  614.     }

  616.     /** {@inheritDoc} */
  617.     @Override
  618.     public T getHyDot() {

  620.         if (!hasDerivatives()) {
  621.             return null;
  622.         }

  624.         // Check for equatorial retrograde orbit
  625.         if (FastMath.abs(i.subtract(i.getPi()).getReal()) < 1.0e-10) {
  626.             return getZero().add(Double.NaN);
  627.         }

  629.         final FieldSinCos<T> sc = FastMath.sinCos(raan);
  630.         final T tan     = i.multiply(0.5).tan();
  631.         return sc.sin().multiply(0.5).multiply(tan.multiply(tan).add(1)).multiply(iDot).
  632.                add(sc.cos().multiply(tan).multiply(raanDot));

  634.     }

  636.     /** Get the true latitude argument.
  637.      * @return v + ω true latitude argument (rad)
  638.      */
  639.     public T getAlphaV() {
  640.         switch (cachedPositionAngleType) {
  641.             case TRUE:
  642.                 return cachedAlpha;

  644.             case ECCENTRIC:
  645.                 return FieldCircularLatitudeArgumentUtility.eccentricToTrue(ex, ey, cachedAlpha);

  647.             case MEAN:
  648.                 return FieldCircularLatitudeArgumentUtility.meanToTrue(ex, ey, cachedAlpha);

  650.             default:
  651.                 throw new OrekitInternalError(null);
  652.         }
  653.     }

  655.     /** Get the true latitude argument derivative.
  656.      * @return d(v + ω)/dt true latitude argument derivative (rad/s)
  657.      */
  658.     public T getAlphaVDot() {

  660.         if (!hasDerivatives()) {
  661.             return null;
  662.         }
  663.         switch (cachedPositionAngleType) {
  664.             case ECCENTRIC:
  665.                 final FieldUnivariateDerivative1<T> alphaEUD = new FieldUnivariateDerivative1<>(cachedAlpha, cachedAlphaDot);
  666.                 final FieldUnivariateDerivative1<T> exUD     = new FieldUnivariateDerivative1<>(ex,     exDot);
  667.                 final FieldUnivariateDerivative1<T> eyUD     = new FieldUnivariateDerivative1<>(ey,     eyDot);
  668.                 final FieldUnivariateDerivative1<T> alphaVUD = FieldCircularLatitudeArgumentUtility.eccentricToTrue(exUD, eyUD,
  669.                         alphaEUD);
  670.                 return alphaVUD.getFirstDerivative();

  672.             case TRUE:
  673.                 return cachedAlphaDot;

  675.             case MEAN:
  676.                 final FieldUnivariateDerivative1<T> alphaMUD = new FieldUnivariateDerivative1<>(cachedAlpha, cachedAlphaDot);
  677.                 final FieldUnivariateDerivative1<T> exUD2    = new FieldUnivariateDerivative1<>(ex,     exDot);
  678.                 final FieldUnivariateDerivative1<T> eyUD2    = new FieldUnivariateDerivative1<>(ey,     eyDot);
  679.                 final FieldUnivariateDerivative1<T> alphaVUD2 = FieldCircularLatitudeArgumentUtility.meanToTrue(exUD2,
  680.                         eyUD2, alphaMUD);
  681.                 return alphaVUD2.getFirstDerivative();

  683.             default:
  684.                 throw new OrekitInternalError(null);
  685.         }
  686.     }

  688.     /** Get the eccentric latitude argument.
  689.      * @return E + ω eccentric latitude argument (rad)
  690.      */
  691.     public T getAlphaE() {
  692.         switch (cachedPositionAngleType) {
  693.             case TRUE:
  694.                 return FieldCircularLatitudeArgumentUtility.trueToEccentric(ex, ey, cachedAlpha);

  696.             case ECCENTRIC:
  697.                 return cachedAlpha;

  699.             case MEAN:
  700.                 return FieldCircularLatitudeArgumentUtility.meanToEccentric(ex, ey, cachedAlpha);

  702.             default:
  703.                 throw new OrekitInternalError(null);
  704.         }
  705.     }

  707.     /** Get the eccentric latitude argument derivative.
  708.      * @return d(E + ω)/dt eccentric latitude argument derivative (rad/s)
  709.      */
  710.     public T getAlphaEDot() {

  712.         if (!hasDerivatives()) {
  713.             return null;
  714.         }
  715.         switch (cachedPositionAngleType) {
  716.             case TRUE:
  717.                 final FieldUnivariateDerivative1<T> alphaVUD = new FieldUnivariateDerivative1<>(cachedAlpha, cachedAlphaDot);
  718.                 final FieldUnivariateDerivative1<T> exUD = new FieldUnivariateDerivative1<>(ex, exDot);
  719.                 final FieldUnivariateDerivative1<T> eyUD = new FieldUnivariateDerivative1<>(ey, eyDot);
  720.                 final FieldUnivariateDerivative1<T> alphaEUD = FieldCircularLatitudeArgumentUtility.trueToEccentric(exUD, eyUD,
  721.                         alphaVUD);
  722.                 return alphaEUD.getFirstDerivative();

  724.             case ECCENTRIC:
  725.                 return cachedAlphaDot;

  727.             case MEAN:
  728.                 final FieldUnivariateDerivative1<T> alphaMUD = new FieldUnivariateDerivative1<>(cachedAlpha, cachedAlphaDot);
  729.                 final FieldUnivariateDerivative1<T> exUD2 = new FieldUnivariateDerivative1<>(ex, exDot);
  730.                 final FieldUnivariateDerivative1<T> eyUD2 = new FieldUnivariateDerivative1<>(ey, eyDot);
  731.                 final FieldUnivariateDerivative1<T> alphaVUD2 = FieldCircularLatitudeArgumentUtility.meanToEccentric(exUD2,
  732.                         eyUD2, alphaMUD);
  733.                 return alphaVUD2.getFirstDerivative();

  735.             default:
  736.                 throw new OrekitInternalError(null);
  737.         }

  739.     }

  741.     /** Get the mean latitude argument.
  742.      * @return M + ω mean latitude argument (rad)
  743.      */
  744.     public T getAlphaM() {
  745.         switch (cachedPositionAngleType) {
  746.             case TRUE:
  747.                 return FieldCircularLatitudeArgumentUtility.trueToMean(ex, ey, cachedAlpha);

  749.             case MEAN:
  750.                 return cachedAlpha;

  752.             case ECCENTRIC:
  753.                 return FieldCircularLatitudeArgumentUtility.eccentricToMean(ex, ey, cachedAlpha);

  755.             default:
  756.                 throw new OrekitInternalError(null);
  757.         }
  758.     }

  760.     /** Get the mean latitude argument derivative.
  761.      * @return d(M + ω)/dt mean latitude argument derivative (rad/s)
  762.      */
  763.     public T getAlphaMDot() {

  765.         if (!hasDerivatives()) {
  766.             return null;
  767.         }
  768.         switch (cachedPositionAngleType) {
  769.             case TRUE:
  770.                 final FieldUnivariateDerivative1<T> alphaVUD = new FieldUnivariateDerivative1<>(cachedAlpha, cachedAlphaDot);
  771.                 final FieldUnivariateDerivative1<T> exUD = new FieldUnivariateDerivative1<>(ex, exDot);
  772.                 final FieldUnivariateDerivative1<T> eyUD = new FieldUnivariateDerivative1<>(ey, eyDot);
  773.                 final FieldUnivariateDerivative1<T> alphaMUD = FieldCircularLatitudeArgumentUtility.trueToMean(exUD, eyUD,
  774.                         alphaVUD);
  775.                 return alphaMUD.getFirstDerivative();

  777.             case MEAN:
  778.                 return cachedAlphaDot;

  780.             case ECCENTRIC:
  781.                 final FieldUnivariateDerivative1<T> alphaEUD = new FieldUnivariateDerivative1<>(cachedAlpha, cachedAlphaDot);
  782.                 final FieldUnivariateDerivative1<T> exUD2 = new FieldUnivariateDerivative1<>(ex, exDot);
  783.                 final FieldUnivariateDerivative1<T> eyUD2 = new FieldUnivariateDerivative1<>(ey, eyDot);
  784.                 final FieldUnivariateDerivative1<T> alphaMUD2 = FieldCircularLatitudeArgumentUtility.eccentricToMean(exUD2,
  785.                         eyUD2, alphaEUD);
  786.                 return alphaMUD2.getFirstDerivative();

  788.             default:
  789.                 throw new OrekitInternalError(null);
  790.         }
  791.     }

  793.     /** Get the latitude argument.
  794.      * @param type type of the angle
  795.      * @return latitude argument (rad)
  796.      */
  797.     public T getAlpha(final PositionAngleType type) {
  798.         return (type == PositionAngleType.MEAN) ? getAlphaM() :
  799.                                               ((type == PositionAngleType.ECCENTRIC) ? getAlphaE() :
  800.                                                                                    getAlphaV());
  801.     }

  803.     /** Get the latitude argument derivative.
  804.      * @param type type of the angle
  805.      * @return latitude argument derivative (rad/s)
  806.      */
  807.     public T getAlphaDot(final PositionAngleType type) {
  808.         return (type == PositionAngleType.MEAN) ? getAlphaMDot() :
  809.                                               ((type == PositionAngleType.ECCENTRIC) ? getAlphaEDot() :
  810.                                                                                    getAlphaVDot());
  811.     }

  813.     /** Computes the true latitude argument from the eccentric latitude argument.
  814.      * @param alphaE = E + ω eccentric latitude argument (rad)
  815.      * @param ex e cos(ω), first component of circular eccentricity vector
  816.      * @param ey e sin(ω), second component of circular eccentricity vector
  817.      * @param <T> Type of the field elements
  818.      * @return the true latitude argument.
  819.      */
  820.     @Deprecated
  821.     public static <T extends CalculusFieldElement<T>> T eccentricToTrue(final T alphaE, final T ex, final T ey) {
  822.         return FieldCircularLatitudeArgumentUtility.eccentricToTrue(ex, ey, alphaE);
  823.     }

  825.     /** Computes the eccentric latitude argument from the true latitude argument.
  826.      * @param alphaV = v + ω true latitude argument (rad)
  827.      * @param ex e cos(ω), first component of circular eccentricity vector
  828.      * @param ey e sin(ω), second component of circular eccentricity vector
  829.      * @param <T> Type of the field elements
  830.      * @return the eccentric latitude argument.
  831.      */
  832.     @Deprecated
  833.     public static <T extends CalculusFieldElement<T>> T trueToEccentric(final T alphaV, final T ex, final T ey) {
  834.         return FieldCircularLatitudeArgumentUtility.trueToEccentric(ex, ey, alphaV);
  835.     }

  837.     /** Computes the eccentric latitude argument from the mean latitude argument.
  838.      * @param alphaM = M + ω  mean latitude argument (rad)
  839.      * @param ex e cos(ω), first component of circular eccentricity vector
  840.      * @param ey e sin(ω), second component of circular eccentricity vector
  841.      * @param <T> Type of the field elements
  842.      * @return the eccentric latitude argument.
  843.      */
  844.     @Deprecated
  845.     public static <T extends CalculusFieldElement<T>> T meanToEccentric(final T alphaM, final T ex, final T ey) {
  846.         return FieldCircularLatitudeArgumentUtility.meanToEccentric(ex, ey, alphaM);
  847.     }

  849.     /** Computes the mean latitude argument from the eccentric latitude argument.
  850.      * @param alphaE = E + ω  eccentric latitude argument (rad)
  851.      * @param ex e cos(ω), first component of circular eccentricity vector
  852.      * @param ey e sin(ω), second component of circular eccentricity vector
  853.      * @param <T> Type of the field elements
  854.      * @return the mean latitude argument.
  855.      */
  856.     @Deprecated
  857.     public static <T extends CalculusFieldElement<T>> T eccentricToMean(final T alphaE, final T ex, final T ey) {
  858.         return FieldCircularLatitudeArgumentUtility.eccentricToMean(ex, ey, alphaE);
  859.     }

  861.     /** {@inheritDoc} */
  862.     @Override
  863.     public T getE() {
  864.         return ex.multiply(ex).add(ey.multiply(ey)).sqrt();
  865.     }

  867.     /** {@inheritDoc} */
  868.     @Override
  869.     public T getEDot() {

  871.         if (!hasDerivatives()) {
  872.             return null;
  873.         }

  875.         return ex.multiply(exDot).add(ey.multiply(eyDot)).divide(getE());

  877.     }

  879.     /** {@inheritDoc} */
  880.     @Override
  881.     public T getI() {
  882.         return i;
  883.     }

  885.     /** {@inheritDoc} */
  886.     @Override
  887.     public T getIDot() {
  888.         return iDot;
  889.     }

  891.     /** Get the right ascension of the ascending node.
  892.      * @return right ascension of the ascending node (rad)
  893.      */
  894.     public T getRightAscensionOfAscendingNode() {
  895.         return raan;
  896.     }

  898.     /** Get the right ascension of the ascending node derivative.
  899.      * @return right ascension of the ascending node derivative (rad/s)
  900.      */
  901.     public T getRightAscensionOfAscendingNodeDot() {
  902.         return raanDot;
  903.     }

  905.     /** {@inheritDoc} */
  906.     @Override
  907.     public T getLv() {
  908.         return getAlphaV().add(raan);
  909.     }

  911.     /** {@inheritDoc} */
  912.     @Override
  913.     public T getLvDot() {
  914.         return hasDerivatives() ? getAlphaVDot().add(raanDot) : null;
  915.     }

  917.     /** {@inheritDoc} */
  918.     @Override
  919.     public T getLE() {
  920.         return getAlphaE().add(raan);
  921.     }

  923.     /** {@inheritDoc} */
  924.     @Override
  925.     public T getLEDot() {
  926.         return hasDerivatives() ? getAlphaEDot().add(raanDot) : null;
  927.     }

  929.     /** {@inheritDoc} */
  930.     @Override
  931.     public T getLM() {
  932.         return getAlphaM().add(raan);
  933.     }

  935.     /** {@inheritDoc} */
  936.     @Override
  937.     public T getLMDot() {
  938.         return hasDerivatives() ? getAlphaMDot().add(raanDot) : null;
  939.     }

  941.     /** {@inheritDoc} */
  942.     @Override
  943.     public boolean hasDerivatives() {
  944.         return aDot != null;
  945.     }

  947.     /** Compute position and velocity but not acceleration.
  948.      */
  949.     private void computePVWithoutA() {

  951.         if (partialPV != null) {
  952.             // already computed
  953.             return;
  954.         }

  956.         // get equinoctial parameters
  957.         final T equEx = getEquinoctialEx();
  958.         final T equEy = getEquinoctialEy();
  959.         final T hx = getHx();
  960.         final T hy = getHy();
  961.         final T lE = getLE();
  962.         // inclination-related intermediate parameters
  963.         final T hx2   = hx.multiply(hx);
  964.         final T hy2   = hy.multiply(hy);
  965.         final T factH = (hx2.add(1).add(hy2)).reciprocal();

  967.         // reference axes defining the orbital plane
  968.         final T ux = (hx2.add(1).subtract(hy2)).multiply(factH);
  969.         final T uy =  hx.multiply(2).multiply(hy).multiply(factH);
  970.         final T uz = hy.multiply(-2).multiply(factH);

  972.         final T vx = uy;
  973.         final T vy = (hy2.subtract(hx2).add(1)).multiply(factH);
  974.         final T vz =  hx.multiply(factH).multiply(2);

  976.         // eccentricity-related intermediate parameters
  977.         final T exey = equEx.multiply(equEy);
  978.         final T ex2  = equEx.square();
  979.         final T ey2  = equEy.square();
  980.         final T e2   = ex2.add(ey2);
  981.         final T eta  = e2.negate().add(1).sqrt().add(1);
  982.         final T beta = eta.reciprocal();

  984.         // eccentric latitude argument
  985.         final FieldSinCos<T> scLe = FastMath.sinCos(lE);
  986.         final T cLe    = scLe.cos();
  987.         final T sLe    = scLe.sin();
  988.         final T exCeyS = equEx.multiply(cLe).add(equEy.multiply(sLe));
  989.         // coordinates of position and velocity in the orbital plane
  990.         final T x      = a.multiply(beta.negate().multiply(ey2).add(1).multiply(cLe).add(beta.multiply(exey).multiply(sLe)).subtract(equEx));
  991.         final T y      = a.multiply(beta.negate().multiply(ex2).add(1).multiply(sLe).add(beta.multiply(exey).multiply(cLe)).subtract(equEy));

  993.         final T factor = getOne().add(getMu()).divide(a).sqrt().divide(exCeyS.negate().add(1));
  994.         final T xdot   = factor.multiply( beta.multiply(equEy).multiply(exCeyS).subtract(sLe ));
  995.         final T ydot   = factor.multiply( cLe.subtract(beta.multiply(equEx).multiply(exCeyS)));

  997.         final FieldVector3D<T> position = new FieldVector3D<>(x.multiply(ux).add(y.multiply(vx)),
  998.                                                               x.multiply(uy).add(y.multiply(vy)),
  999.                                                               x.multiply(uz).add(y.multiply(vz)));
  1000.         final FieldVector3D<T> velocity = new FieldVector3D<>(xdot.multiply(ux).add(ydot.multiply(vx)),
  1001.                                                               xdot.multiply(uy).add(ydot.multiply(vy)),
  1002.                                                               xdot.multiply(uz).add(ydot.multiply(vz)));

  1004.         partialPV = new FieldPVCoordinates<>(position, velocity);

  1006.     }


  1009.     /** Initialize cached alpha with rate.
  1010.      * @param alpha input alpha
  1011.      * @param alphaDot rate of input alpha
  1012.      * @param inputType position angle type passed as input
  1013.      * @return alpha to cache with rate
  1014.      * @since 12.1
  1015.      */
  1016.     private FieldUnivariateDerivative1<T> initializeCachedAlpha(final T alpha, final T alphaDot,
  1017.                                                                 final PositionAngleType inputType) {
  1018.         if (cachedPositionAngleType == inputType) {
  1019.             return new FieldUnivariateDerivative1<>(alpha, alphaDot);

  1021.         } else {
  1022.             final FieldUnivariateDerivative1<T> exUD = new FieldUnivariateDerivative1<>(ex, exDot);
  1023.             final FieldUnivariateDerivative1<T> eyUD = new FieldUnivariateDerivative1<>(ey, eyDot);
  1024.             final FieldUnivariateDerivative1<T> alphaUD = new FieldUnivariateDerivative1<>(alpha, alphaDot);

  1026.             switch (cachedPositionAngleType) {

  1028.                 case ECCENTRIC:
  1029.                     if (inputType == PositionAngleType.MEAN) {
  1030.                         return FieldCircularLatitudeArgumentUtility.meanToEccentric(exUD, eyUD, alphaUD);
  1031.                     } else {
  1032.                         return FieldCircularLatitudeArgumentUtility.trueToEccentric(exUD, eyUD, alphaUD);
  1033.                     }

  1035.                 case TRUE:
  1036.                     if (inputType == PositionAngleType.MEAN) {
  1037.                         return FieldCircularLatitudeArgumentUtility.meanToTrue(exUD, eyUD, alphaUD);
  1038.                     } else {
  1039.                         return FieldCircularLatitudeArgumentUtility.eccentricToTrue(exUD, eyUD, alphaUD);
  1040.                     }

  1042.                 case MEAN:
  1043.                     if (inputType == PositionAngleType.TRUE) {
  1044.                         return FieldCircularLatitudeArgumentUtility.trueToMean(exUD, eyUD, alphaUD);
  1045.                     } else {
  1046.                         return FieldCircularLatitudeArgumentUtility.eccentricToMean(exUD, eyUD, alphaUD);
  1047.                     }

  1049.                 default:
  1050.                     throw new OrekitInternalError(null);

  1052.             }

  1054.         }

  1056.     }

  1058.     /** Initialize cached alpha.
  1059.      * @param alpha input alpha
  1060.      * @param positionAngleType position angle type passed as input
  1061.      * @return alpha to cache
  1062.      * @since 12.1
  1063.      */
  1064.     private T initializeCachedAlpha(final T alpha, final PositionAngleType positionAngleType) {
  1065.         return FieldCircularLatitudeArgumentUtility.convertAlpha(positionAngleType, alpha, ex, ey, cachedPositionAngleType);
  1066.     }

  1068.     /** Compute non-Keplerian part of the acceleration from first time derivatives.
  1069.      * <p>
  1070.      * This method should be called only when {@link #hasDerivatives()} returns true.
  1071.      * </p>
  1072.      * @return non-Keplerian part of the acceleration
  1073.      */
  1074.     private FieldVector3D<T> nonKeplerianAcceleration() {

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

  1079.         final T nonKeplerianMeanMotion = getAlphaMDot().subtract(getKeplerianMeanMotion());
  1080.         final T nonKeplerianAx =     dCdP[3][0].multiply(aDot).
  1081.                                  add(dCdP[3][1].multiply(exDot)).
  1082.                                  add(dCdP[3][2].multiply(eyDot)).
  1083.                                  add(dCdP[3][3].multiply(iDot)).
  1084.                                  add(dCdP[3][4].multiply(raanDot)).
  1085.                                  add(dCdP[3][5].multiply(nonKeplerianMeanMotion));
  1086.         final T nonKeplerianAy =     dCdP[4][0].multiply(aDot).
  1087.                                  add(dCdP[4][1].multiply(exDot)).
  1088.                                  add(dCdP[4][2].multiply(eyDot)).
  1089.                                  add(dCdP[4][3].multiply(iDot)).
  1090.                                  add(dCdP[4][4].multiply(raanDot)).
  1091.                                  add(dCdP[4][5].multiply(nonKeplerianMeanMotion));
  1092.         final T nonKeplerianAz =     dCdP[5][0].multiply(aDot).
  1093.                                  add(dCdP[5][1].multiply(exDot)).
  1094.                                  add(dCdP[5][2].multiply(eyDot)).
  1095.                                  add(dCdP[5][3].multiply(iDot)).
  1096.                                  add(dCdP[5][4].multiply(raanDot)).
  1097.                                  add(dCdP[5][5].multiply(nonKeplerianMeanMotion));

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

  1101.     }

  1103.     /** {@inheritDoc} */
  1104.     @Override
  1105.     protected FieldVector3D<T> initPosition() {
  1106.         // get equinoctial parameters
  1107.         final T equEx = getEquinoctialEx();
  1108.         final T equEy = getEquinoctialEy();
  1109.         final T hx = getHx();
  1110.         final T hy = getHy();
  1111.         final T lE = getLE();
  1112.         // inclination-related intermediate parameters
  1113.         final T hx2   = hx.multiply(hx);
  1114.         final T hy2   = hy.multiply(hy);
  1115.         final T factH = (hx2.add(1).add(hy2)).reciprocal();

  1117.         // reference axes defining the orbital plane
  1118.         final T ux = (hx2.add(1).subtract(hy2)).multiply(factH);
  1119.         final T uy =  hx.multiply(2).multiply(hy).multiply(factH);
  1120.         final T uz = hy.multiply(-2).multiply(factH);

  1122.         final T vx = uy;
  1123.         final T vy = (hy2.subtract(hx2).add(1)).multiply(factH);
  1124.         final T vz =  hx.multiply(factH).multiply(2);

  1126.         // eccentricity-related intermediate parameters
  1127.         final T exey = equEx.multiply(equEy);
  1128.         final T ex2  = equEx.square();
  1129.         final T ey2  = equEy.square();
  1130.         final T e2   = ex2.add(ey2);
  1131.         final T eta  = e2.negate().add(1).sqrt().add(1);
  1132.         final T beta = eta.reciprocal();

  1134.         // eccentric latitude argument
  1135.         final FieldSinCos<T> scLe = FastMath.sinCos(lE);
  1136.         final T cLe    = scLe.cos();
  1137.         final T sLe    = scLe.sin();

  1139.         // coordinates of position and velocity in the orbital plane
  1140.         final T x      = a.multiply(beta.negate().multiply(ey2).add(1).multiply(cLe).add(beta.multiply(exey).multiply(sLe)).subtract(equEx));
  1141.         final T y      = a.multiply(beta.negate().multiply(ex2).add(1).multiply(sLe).add(beta.multiply(exey).multiply(cLe)).subtract(equEy));

  1143.         return new FieldVector3D<>(x.multiply(ux).add(y.multiply(vx)),
  1144.                                    x.multiply(uy).add(y.multiply(vy)),
  1145.                                    x.multiply(uz).add(y.multiply(vz)));

  1147.     }

  1149.     /** {@inheritDoc} */
  1150.     @Override
  1151.     protected TimeStampedFieldPVCoordinates<T> initPVCoordinates() {

  1153.         // position and velocity
  1154.         computePVWithoutA();

  1156.         // acceleration
  1157.         final T r2 = partialPV.getPosition().getNormSq();
  1158.         final FieldVector3D<T> keplerianAcceleration = new FieldVector3D<>(r2.multiply(r2.sqrt()).reciprocal().multiply(getMu().negate()),
  1159.                                                                            partialPV.getPosition());
  1160.         final FieldVector3D<T> acceleration = hasDerivatives() ?
  1161.                                               keplerianAcceleration.add(nonKeplerianAcceleration()) :
  1162.                                               keplerianAcceleration;

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

  1166.     }

  1168.     /** {@inheritDoc} */
  1169.     @Override
  1170.     public FieldCircularOrbit<T> shiftedBy(final double dt) {
  1171.         return shiftedBy(getZero().newInstance(dt));
  1172.     }

  1174.     /** {@inheritDoc} */
  1175.     @Override
  1176.     public FieldCircularOrbit<T> shiftedBy(final T dt) {

  1178.         // use Keplerian-only motion
  1179.         final FieldCircularOrbit<T> keplerianShifted = new FieldCircularOrbit<>(a, ex, ey, i, raan,
  1180.                                                                                 getAlphaM().add(getKeplerianMeanMotion().multiply(dt)),
  1181.                                                                                 PositionAngleType.MEAN, cachedPositionAngleType, getFrame(),
  1182.                                                                                 getDate().shiftedBy(dt), getMu());

  1184.         if (hasDerivatives()) {

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

  1189.             // add quadratic effect of non-Keplerian acceleration to Keplerian-only shift
  1190.             keplerianShifted.computePVWithoutA();
  1191.             final FieldVector3D<T> fixedP   = new FieldVector3D<>(getOne(), keplerianShifted.partialPV.getPosition(),
  1192.                                                                   dt.square().multiply(0.5), nonKeplerianAcceleration);
  1193.             final T   fixedR2 = fixedP.getNormSq();
  1194.             final T   fixedR  = fixedR2.sqrt();
  1195.             final FieldVector3D<T> fixedV  = new FieldVector3D<>(getOne(), keplerianShifted.partialPV.getVelocity(),
  1196.                                                                  dt, nonKeplerianAcceleration);
  1197.             final FieldVector3D<T> fixedA  = new FieldVector3D<>(fixedR2.multiply(fixedR).reciprocal().multiply(getMu().negate()),
  1198.                                                                  keplerianShifted.partialPV.getPosition(),
  1199.                                                                  getOne(), nonKeplerianAcceleration);

  1201.             // build a new orbit, taking non-Keplerian acceleration into account
  1202.             return new FieldCircularOrbit<>(new TimeStampedFieldPVCoordinates<>(keplerianShifted.getDate(),
  1203.                                                                                 fixedP, fixedV, fixedA),
  1204.                                             keplerianShifted.getFrame(), keplerianShifted.getMu());

  1206.         } else {
  1207.             // Keplerian-only motion is all we can do
  1208.             return keplerianShifted;
  1209.         }

  1211.     }

  1213.     /** {@inheritDoc} */
  1214.     @Override
  1215.     protected T[][] computeJacobianMeanWrtCartesian() {

  1217.         final T[][] jacobian = MathArrays.buildArray(getOne().getField(), 6, 6);

  1219.         // compute various intermediate parameters
  1220.         computePVWithoutA();
  1221.         final FieldVector3D<T> position = partialPV.getPosition();
  1222.         final FieldVector3D<T> velocity = partialPV.getVelocity();

  1224.         final T x          = position.getX();
  1225.         final T y          = position.getY();
  1226.         final T z          = position.getZ();
  1227.         final T vx         = velocity.getX();
  1228.         final T vy         = velocity.getY();
  1229.         final T vz         = velocity.getZ();
  1230.         final T pv         = FieldVector3D.dotProduct(position, velocity);
  1231.         final T r2         = position.getNormSq();
  1232.         final T r          = r2.sqrt();
  1233.         final T v2         = velocity.getNormSq();

  1235.         final T mu         = getMu();
  1236.         final T oOsqrtMuA  = getOne().divide(a.multiply(mu).sqrt());
  1237.         final T rOa        = r.divide(a);
  1238.         final T aOr        = a.divide(r);
  1239.         final T aOr2       = a.divide(r2);
  1240.         final T a2         = a.square();

  1242.         final T ex2        = ex.square();
  1243.         final T ey2        = ey.square();
  1244.         final T e2         = ex2.add(ey2);
  1245.         final T epsilon    = e2.negate().add(1.0).sqrt();
  1246.         final T beta       = epsilon.add(1).reciprocal();

  1248.         final T eCosE      = rOa.negate().add(1);
  1249.         final T eSinE      = pv.multiply(oOsqrtMuA);

  1251.         final FieldSinCos<T> scI    = FastMath.sinCos(i);
  1252.         final FieldSinCos<T> scRaan = FastMath.sinCos(raan);
  1253.         final T cosI       = scI.cos();
  1254.         final T sinI       = scI.sin();
  1255.         final T cosRaan    = scRaan.cos();
  1256.         final T sinRaan    = scRaan.sin();

  1258.         // da
  1259.         fillHalfRow(aOr.multiply(2.0).multiply(aOr2), position, jacobian[0], 0);
  1260.         fillHalfRow(a2.multiply(mu.divide(2.).reciprocal()), velocity, jacobian[0], 3);

  1262.         // differentials of the normalized momentum
  1263.         final FieldVector3D<T> danP = new FieldVector3D<>(v2, position, pv.negate(), velocity);
  1264.         final FieldVector3D<T> danV = new FieldVector3D<>(r2, velocity, pv.negate(), position);
  1265.         final T recip   = partialPV.getMomentum().getNorm().reciprocal();
  1266.         final T recip2  = recip.multiply(recip);
  1267.         final T recip2N = recip2.negate();
  1268.         final FieldVector3D<T> dwXP = new FieldVector3D<>(recip,
  1269.                                                           new FieldVector3D<>(getZero(), vz, vy.negate()),
  1270.                                                           recip2N.multiply(sinRaan).multiply(sinI),
  1271.                                                           danP);
  1272.         final FieldVector3D<T> dwYP = new FieldVector3D<>(recip,
  1273.                                                           new FieldVector3D<>(vz.negate(), getZero(), vx),
  1274.                                                           recip2.multiply(cosRaan).multiply(sinI),
  1275.                                                           danP);
  1276.         final FieldVector3D<T> dwZP = new FieldVector3D<>(recip,
  1277.                                                           new FieldVector3D<>(vy, vx.negate(), getZero()),
  1278.                                                           recip2N.multiply(cosI),
  1279.                                                           danP);
  1280.         final FieldVector3D<T> dwXV = new FieldVector3D<>(recip,
  1281.                                                           new FieldVector3D<>(getZero(), z.negate(), y),
  1282.                                                           recip2N.multiply(sinRaan).multiply(sinI),
  1283.                                                           danV);
  1284.         final FieldVector3D<T> dwYV = new FieldVector3D<>(recip,
  1285.                                                           new FieldVector3D<>(z, getZero(), x.negate()),
  1286.                                                           recip2.multiply(cosRaan).multiply(sinI),
  1287.                                                           danV);
  1288.         final FieldVector3D<T> dwZV = new FieldVector3D<>(recip,
  1289.                                                           new FieldVector3D<>(y.negate(), x, getZero()),
  1290.                                                           recip2N.multiply(cosI),
  1291.                                                           danV);

  1293.         // di
  1294.         fillHalfRow(sinRaan.multiply(cosI), dwXP, cosRaan.negate().multiply(cosI), dwYP, sinI.negate(), dwZP, jacobian[3], 0);
  1295.         fillHalfRow(sinRaan.multiply(cosI), dwXV, cosRaan.negate().multiply(cosI), dwYV, sinI.negate(), dwZV, jacobian[3], 3);

  1297.         // dRaan
  1298.         fillHalfRow(sinRaan.divide(sinI), dwYP, cosRaan.divide(sinI), dwXP, jacobian[4], 0);
  1299.         fillHalfRow(sinRaan.divide(sinI), dwYV, cosRaan.divide(sinI), dwXV, jacobian[4], 3);

  1301.         // orbital frame: (p, q, w) p along ascending node, w along momentum
  1302.         // the coordinates of the spacecraft in this frame are: (u, v, 0)
  1303.         final T u     =  x.multiply(cosRaan).add(y.multiply(sinRaan));
  1304.         final T cv    =  x.negate().multiply(sinRaan).add(y.multiply(cosRaan));
  1305.         final T v     = cv.multiply(cosI).add(z.multiply(sinI));

  1307.         // du
  1308.         final FieldVector3D<T> duP = new FieldVector3D<>(cv.multiply(cosRaan).divide(sinI), dwXP,
  1309.                                                          cv.multiply(sinRaan).divide(sinI), dwYP,
  1310.                                                          getOne(), new FieldVector3D<>(cosRaan, sinRaan, getZero()));
  1311.         final FieldVector3D<T> duV = new FieldVector3D<>(cv.multiply(cosRaan).divide(sinI), dwXV,
  1312.                                                          cv.multiply(sinRaan).divide(sinI), dwYV);

  1314.         // dv
  1315.         final FieldVector3D<T> dvP = new FieldVector3D<>(u.negate().multiply(cosRaan).multiply(cosI).divide(sinI).add(sinRaan.multiply(z)), dwXP,
  1316.                                                          u.negate().multiply(sinRaan).multiply(cosI).divide(sinI).subtract(cosRaan.multiply(z)), dwYP,
  1317.                                                          cv, dwZP,
  1318.                                                          getOne(), new FieldVector3D<>(sinRaan.negate().multiply(cosI), cosRaan.multiply(cosI), sinI));
  1319.         final FieldVector3D<T> dvV = new FieldVector3D<>(u.negate().multiply(cosRaan).multiply(cosI).divide(sinI).add(sinRaan.multiply(z)), dwXV,
  1320.                                                          u.negate().multiply(sinRaan).multiply(cosI).divide(sinI).subtract(cosRaan.multiply(z)), dwYV,
  1321.                                                          cv, dwZV);

  1323.         final FieldVector3D<T> dc1P = new FieldVector3D<>(aOr2.multiply(eSinE.multiply(eSinE).multiply(2).add(1).subtract(eCosE)).divide(r2), position,
  1324.                                                           aOr2.multiply(-2).multiply(eSinE).multiply(oOsqrtMuA), velocity);
  1325.         final FieldVector3D<T> dc1V = new FieldVector3D<>(aOr2.multiply(-2).multiply(eSinE).multiply(oOsqrtMuA), position,
  1326.                                                           getZero().newInstance(2).divide(mu), velocity);
  1327.         final FieldVector3D<T> dc2P = new FieldVector3D<>(aOr2.multiply(eSinE).multiply(eSinE.multiply(eSinE).subtract(e2.negate().add(1))).divide(r2.multiply(epsilon)), position,
  1328.                                                           aOr2.multiply(e2.negate().add(1).subtract(eSinE.multiply(eSinE))).multiply(oOsqrtMuA).divide(epsilon), velocity);
  1329.         final FieldVector3D<T> dc2V = new FieldVector3D<>(aOr2.multiply(e2.negate().add(1).subtract(eSinE.multiply(eSinE))).multiply(oOsqrtMuA).divide(epsilon), position,
  1330.                                                           eSinE.divide(epsilon.multiply(mu)), velocity);

  1332.         final T cof1   = aOr2.multiply(eCosE.subtract(e2));
  1333.         final T cof2   = aOr2.multiply(epsilon).multiply(eSinE);
  1334.         final FieldVector3D<T> dexP = new FieldVector3D<>(u, dc1P,  v, dc2P, cof1, duP,  cof2, dvP);
  1335.         final FieldVector3D<T> dexV = new FieldVector3D<>(u, dc1V,  v, dc2V, cof1, duV,  cof2, dvV);
  1336.         final FieldVector3D<T> deyP = new FieldVector3D<>(v, dc1P, u.negate(), dc2P, cof1, dvP, cof2.negate(), duP);
  1337.         final FieldVector3D<T> deyV = new FieldVector3D<>(v, dc1V, u.negate(), dc2V, cof1, dvV, cof2.negate(), duV);
  1338.         fillHalfRow(getOne(), dexP, jacobian[1], 0);
  1339.         fillHalfRow(getOne(), dexV, jacobian[1], 3);
  1340.         fillHalfRow(getOne(), deyP, jacobian[2], 0);
  1341.         fillHalfRow(getOne(), deyV, jacobian[2], 3);

  1343.         final T cle = u.divide(a).add(ex).subtract(eSinE.multiply(beta).multiply(ey));
  1344.         final T sle = v.divide(a).add(ey).add(eSinE.multiply(beta).multiply(ex));
  1345.         final T m1  = beta.multiply(eCosE);
  1346.         final T m2  = m1.multiply(eCosE).negate().add(1);
  1347.         final T m3  = (u.multiply(ey).subtract(v.multiply(ex))).add(eSinE.multiply(beta).multiply(u.multiply(ex).add(v.multiply(ey))));
  1348.         final T m4  = sle.negate().add(cle.multiply(eSinE).multiply(beta));
  1349.         final T m5  = cle.add(sle.multiply(eSinE).multiply(beta));
  1350.         final T kk = m3.multiply(2).divide(r).add(aOr.multiply(eSinE)).add(m1.multiply(eSinE).multiply(m1.add(1).subtract(aOr.add(1).multiply(m2))).divide(epsilon)).divide(r2);
  1351.         final T jj = (m1.multiply(m2).divide(epsilon).subtract(1)).multiply(oOsqrtMuA);
  1352.         fillHalfRow(kk, position,
  1353.                     jj, velocity,
  1354.                     m4, dexP,
  1355.                     m5, deyP,
  1356.                     sle.negate().divide(a), duP,
  1357.                     cle.divide(a), dvP,
  1358.                     jacobian[5], 0);
  1359.         final T ll = (m1.multiply(m2).divide(epsilon ).subtract(1)).multiply(oOsqrtMuA);
  1360.         final T mm = m3.multiply(2).add(eSinE.multiply(a)).add(m1.multiply(eSinE).multiply(r).multiply(eCosE.multiply(beta).multiply(2).subtract(aOr.multiply(m2))).divide(epsilon)).divide(mu);

  1362.         fillHalfRow(ll, position,
  1363.                     mm, velocity,
  1364.                     m4, dexV,
  1365.                     m5, deyV,
  1366.                     sle.negate().divide(a), duV,
  1367.                     cle.divide(a), dvV,
  1368.                     jacobian[5], 3);
  1369.         return jacobian;

  1371.     }

  1373.     /** {@inheritDoc} */
  1374.     @Override
  1375.     protected T[][] computeJacobianEccentricWrtCartesian() {

  1377.         // start by computing the Jacobian with mean angle
  1378.         final T[][] jacobian = computeJacobianMeanWrtCartesian();

  1380.         // Differentiating the Kepler equation aM = aE - ex sin aE + ey cos aE leads to:
  1381.         // daM = (1 - ex cos aE - ey sin aE) dE - sin aE dex + cos aE dey
  1382.         // which is inverted and rewritten as:
  1383.         // daE = a/r daM + sin aE a/r dex - cos aE a/r dey
  1384.         final T alphaE            = getAlphaE();
  1385.         final FieldSinCos<T> scAe = FastMath.sinCos(alphaE);
  1386.         final T cosAe             = scAe.cos();
  1387.         final T sinAe             = scAe.sin();
  1388.         final T aOr               = getOne().divide(getOne().subtract(ex.multiply(cosAe)).subtract(ey.multiply(sinAe)));

  1390.         // update longitude row
  1391.         final T[] rowEx = jacobian[1];
  1392.         final T[] rowEy = jacobian[2];
  1393.         final T[] rowL  = jacobian[5];
  1394.         for (int j = 0; j < 6; ++j) {
  1395.          // rowL[j] = aOr * (      rowL[j] +   sinAe *        rowEx[j]   -         cosAe *        rowEy[j]);
  1396.             rowL[j] = aOr.multiply(rowL[j].add(sinAe.multiply(rowEx[j])).subtract( cosAe.multiply(rowEy[j])));
  1397.         }
  1398.         jacobian[5] = rowL;
  1399.         return jacobian;

  1401.     }
  1402.     /** {@inheritDoc} */
  1403.     @Override
  1404.     protected T[][] computeJacobianTrueWrtCartesian() {

  1406.         // start by computing the Jacobian with eccentric angle
  1407.         final T[][] jacobian = computeJacobianEccentricWrtCartesian();
  1408.         // Differentiating the eccentric latitude equation
  1409.         // tan((aV - aE)/2) = [ex sin aE - ey cos aE] / [sqrt(1-ex^2-ey^2) + 1 - ex cos aE - ey sin aE]
  1410.         // leads to
  1411.         // cT (daV - daE) = cE daE + cX dex + cY dey
  1412.         // with
  1413.         // cT = [d^2 + (ex sin aE - ey cos aE)^2] / 2
  1414.         // d  = 1 + sqrt(1-ex^2-ey^2) - ex cos aE - ey sin aE
  1415.         // cE = (ex cos aE + ey sin aE) (sqrt(1-ex^2-ey^2) + 1) - ex^2 - ey^2
  1416.         // cX =  sin aE (sqrt(1-ex^2-ey^2) + 1) - ey + ex (ex sin aE - ey cos aE) / sqrt(1-ex^2-ey^2)
  1417.         // cY = -cos aE (sqrt(1-ex^2-ey^2) + 1) + ex + ey (ex sin aE - ey cos aE) / sqrt(1-ex^2-ey^2)
  1418.         // which can be solved to find the differential of the true latitude
  1419.         // daV = (cT + cE) / cT daE + cX / cT deX + cY / cT deX
  1420.         final T alphaE            = getAlphaE();
  1421.         final FieldSinCos<T> scAe = FastMath.sinCos(alphaE);
  1422.         final T cosAe             = scAe.cos();
  1423.         final T sinAe             = scAe.sin();
  1424.         final T eSinE             = ex.multiply(sinAe).subtract(ey.multiply(cosAe));
  1425.         final T ecosE             = ex.multiply(cosAe).add(ey.multiply(sinAe));
  1426.         final T e2                = ex.multiply(ex).add(ey.multiply(ey));
  1427.         final T epsilon           = (getOne().subtract(e2)).sqrt();
  1428.         final T onePeps           = getOne().add(epsilon);
  1429.         final T d                 = onePeps.subtract(ecosE);
  1430.         final T cT                = (d.multiply(d).add(eSinE.multiply(eSinE))).divide(2);
  1431.         final T cE                = ecosE.multiply(onePeps).subtract(e2);
  1432.         final T cX                = ex.multiply(eSinE).divide(epsilon).subtract(ey).add(sinAe.multiply(onePeps));
  1433.         final T cY                = ey.multiply(eSinE).divide(epsilon).add(ex).subtract(cosAe.multiply(onePeps));
  1434.         final T factorLe          = (cT.add(cE)).divide(cT);
  1435.         final T factorEx          = cX.divide(cT);
  1436.         final T factorEy          = cY.divide(cT);

  1438.         // update latitude row
  1439.         final T[] rowEx = jacobian[1];
  1440.         final T[] rowEy = jacobian[2];
  1441.         final T[] rowA = jacobian[5];
  1442.         for (int j = 0; j < 6; ++j) {
  1443.             rowA[j] = factorLe.multiply(rowA[j]).add(factorEx.multiply(rowEx[j])).add(factorEy.multiply(rowEy[j]));
  1444.         }
  1445.         return jacobian;

  1447.     }

  1449.     /** {@inheritDoc} */
  1450.     @Override
  1451.     public void addKeplerContribution(final PositionAngleType type, final T gm, final T[] pDot) {
  1452.         pDot[5] = pDot[5].add(computeKeplerianAlphaDot(type, a, ex, ey, gm, cachedAlpha, cachedPositionAngleType));
  1453.     }

  1455.     /**
  1456.      * Compute rate of argument of latitude.
  1457.      * @param type position angle type of rate
  1458.      * @param a semi major axis
  1459.      * @param ex ex
  1460.      * @param ey ey
  1461.      * @param mu mu
  1462.      * @param alpha argument of latitude
  1463.      * @param cachedType position angle type of passed alpha
  1464.      * @param <T> field type
  1465.      * @return first-order time derivative for alpha
  1466.      * @since 12.2
  1467.      */
  1468.     private static <T extends CalculusFieldElement<T>> T computeKeplerianAlphaDot(final PositionAngleType type, final T a,
  1469.                                                                                   final T ex, final T ey, final T mu,
  1470.                                                                                   final T alpha, final PositionAngleType cachedType) {
  1471.         final T n = a.reciprocal().multiply(mu).sqrt().divide(a);
  1472.         if (type == PositionAngleType.MEAN) {
  1473.             return n;
  1474.         }
  1475.         final FieldSinCos<T> sc;
  1476.         final T ksi;
  1477.         if (type == PositionAngleType.ECCENTRIC) {
  1478.             sc = FastMath.sinCos(FieldCircularLatitudeArgumentUtility.convertAlpha(cachedType, alpha, ex, ey, type));
  1479.             ksi  = ((ex.multiply(sc.cos())).add(ey.multiply(sc.sin()))).negate().add(1).reciprocal();
  1480.             return n.multiply(ksi);
  1481.         } else {  // TRUE
  1482.             sc = FastMath.sinCos(FieldCircularLatitudeArgumentUtility.convertAlpha(cachedType, alpha, ex, ey, type));
  1483.             final T one = n.getField().getOne();
  1484.             final T oMe2  = one.subtract(ex.square()).subtract(ey.square());
  1485.             ksi   = one.add(ex.multiply(sc.cos())).add(ey.multiply(sc.sin()));
  1486.             return n.multiply(ksi.square()).divide(oMe2.multiply(oMe2.sqrt()));
  1487.         }
  1488.     }

  1490.     /**  Returns a string representation of this Orbit object.
  1491.      * @return a string representation of this object
  1492.      */
  1493.     public String toString() {
  1494.         return new StringBuilder().append("circular parameters: ").append('{').
  1495.                                   append("a: ").append(a.getReal()).
  1496.                                   append(", ex: ").append(ex.getReal()).append(", ey: ").append(ey.getReal()).
  1497.                                   append(", i: ").append(FastMath.toDegrees(i.getReal())).
  1498.                                   append(", raan: ").append(FastMath.toDegrees(raan.getReal())).
  1499.                                   append(", alphaV: ").append(FastMath.toDegrees(getAlphaV().getReal())).
  1500.                                   append(";}").toString();
  1501.     }

  1503.     /** {@inheritDoc} */
  1504.     @Override
  1505.     public PositionAngleType getCachedPositionAngleType() {
  1506.         return cachedPositionAngleType;
  1507.     }

  1509.     /** {@inheritDoc} */
  1510.     @Override
  1511.     public boolean hasRates() {
  1512.         return hasDerivatives();
  1513.     }

  1515.     /** {@inheritDoc} */
  1516.     @Override
  1517.     public FieldCircularOrbit<T> removeRates() {
  1518.         return new FieldCircularOrbit<>(getA(), getCircularEx(), getCircularEy(),
  1519.                 getI(), getRightAscensionOfAscendingNode(), cachedAlpha,
  1520.                 cachedPositionAngleType, getFrame(), getDate(), getMu());
  1521.     }

  1523.     /** {@inheritDoc} */
  1524.     @Override
  1525.     public CircularOrbit toOrbit() {
  1526.         final double cachedPositionAngle = cachedAlpha.getReal();
  1527.         if (hasDerivatives()) {
  1528.             return new CircularOrbit(a.getReal(), ex.getReal(), ey.getReal(),
  1529.                                      i.getReal(), raan.getReal(), cachedPositionAngle,
  1530.                                      aDot.getReal(), exDot.getReal(), eyDot.getReal(),
  1531.                                      iDot.getReal(), raanDot.getReal(), cachedAlphaDot.getReal(),
  1532.                                      cachedPositionAngleType, getFrame(),
  1533.                                      getDate().toAbsoluteDate(), getMu().getReal());
  1534.         } else {
  1535.             return new CircularOrbit(a.getReal(), ex.getReal(), ey.getReal(),
  1536.                                      i.getReal(), raan.getReal(), cachedPositionAngle,
  1537.                                      cachedPositionAngleType, getFrame(),
  1538.                                      getDate().toAbsoluteDate(), getMu().getReal());
  1539.         }
  1540.     }


  1543. }
Created with JaCoCo 0.8.12.202403310830