CircularOrbit.java

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

  19. import java.io.Serializable;

  21. import org.hipparchus.analysis.differentiation.UnivariateDerivative1;
  22. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  23. import org.hipparchus.util.FastMath;
  24. import org.hipparchus.util.SinCos;
  25. import org.orekit.annotation.DefaultDataContext;
  26. import org.orekit.data.DataContext;
  27. import org.orekit.errors.OrekitIllegalArgumentException;
  28. import org.orekit.errors.OrekitInternalError;
  29. import org.orekit.errors.OrekitMessages;
  30. import org.orekit.frames.Frame;
  31. import org.orekit.time.AbsoluteDate;
  32. import org.orekit.utils.PVCoordinates;
  33. import org.orekit.utils.TimeStampedPVCoordinates;


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

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

  75. public class CircularOrbit extends Orbit implements PositionAngleBased {

  77.     /** Serializable UID. */
  78.     private static final long serialVersionUID = 20231217L;

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

  83.     /** First component of the circular eccentricity vector. */
  84.     private final double ex;

  86.     /** Second component of the circular eccentricity vector. */
  87.     private final double ey;

  89.     /** Inclination (rad). */
  90.     private final double i;

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

  95.     /** Cached latitude argument (rad). */
  96.     private final double cachedAlpha;

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

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

  104.     /** First component of the circular eccentricity vector derivative. */
  105.     private final double exDot;

  107.     /** Second component of the circular eccentricity vector derivative. */
  108.     private final double eyDot;

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

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

  116.     /** True latitude argument derivative (rad/s). */
  117.     private final double cachedAlphaDot;

  119.     /** Indicator for {@link PVCoordinates} serialization. */
  120.     private final boolean serializePV;

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

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

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

  175.     /** Creates a new instance.
  176.      * @param a  semi-major axis (m)
  177.      * @param ex e cos(ω), first component of circular eccentricity vector
  178.      * @param ey e sin(ω), second component of circular eccentricity vector
  179.      * @param i inclination (rad)
  180.      * @param raan right ascension of ascending node (Ω, rad)
  181.      * @param alpha  an + ω, mean, eccentric or true latitude argument (rad)
  182.      * @param aDot  semi-major axis derivative (m/s)
  183.      * @param exDot d(e cos(ω))/dt, first component of circular eccentricity vector derivative
  184.      * @param eyDot d(e sin(ω))/dt, second component of circular eccentricity vector derivative
  185.      * @param iDot inclination  derivative(rad/s)
  186.      * @param raanDot right ascension of ascending node derivative (rad/s)
  187.      * @param alphaDot  d(an + ω), mean, eccentric or true latitude argument derivative (rad/s)
  188.      * @param type type of latitude argument
  189.      * @param cachedPositionAngleType type of cached latitude argument
  190.      * @param frame the frame in which are defined the parameters
  191.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  192.      * @param date date of the orbital parameters
  193.      * @param mu central attraction coefficient (m³/s²)
  194.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  195.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  196.      * @since 12.1
  197.      */
  198.     public CircularOrbit(final double a, final double ex, final double ey,
  199.                          final double i, final double raan, final double alpha,
  200.                          final double aDot, final double exDot, final double eyDot,
  201.                          final double iDot, final double raanDot, final double alphaDot,
  202.                          final PositionAngleType type, final PositionAngleType cachedPositionAngleType,
  203.                          final Frame frame, final AbsoluteDate date, final double mu)
  204.         throws IllegalArgumentException {
  205.         super(frame, date, mu);
  206.         if (ex * ex + ey * ey >= 1.0) {
  207.             throw new OrekitIllegalArgumentException(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS,
  208.                                                      getClass().getName());
  209.         }
  210.         this.a       =  a;
  211.         this.aDot    =  aDot;
  212.         this.ex      = ex;
  213.         this.exDot   = exDot;
  214.         this.ey      = ey;
  215.         this.eyDot   = eyDot;
  216.         this.i       = i;
  217.         this.iDot    = iDot;
  218.         this.raan    = raan;
  219.         this.raanDot = raanDot;
  220.         this.cachedPositionAngleType = cachedPositionAngleType;

  222.         if (hasDerivatives()) {
  223.             final UnivariateDerivative1 alphaUD = initializeCachedAlpha(alpha, alphaDot, type);
  224.             this.cachedAlpha = alphaUD.getValue();
  225.             this.cachedAlphaDot = alphaUD.getFirstDerivative();
  226.         } else {
  227.             this.cachedAlpha = initializeCachedAlpha(alpha, type);
  228.             this.cachedAlphaDot = Double.NaN;
  229.         }

  231.         serializePV = false;
  232.         partialPV   = null;

  234.     }

  236.     /** Creates a new instance with derivatives and with cached position angle same as value inputted.
  237.      * @param a  semi-major axis (m)
  238.      * @param ex e cos(ω), first component of circular eccentricity vector
  239.      * @param ey e sin(ω), second component of circular eccentricity vector
  240.      * @param i inclination (rad)
  241.      * @param raan right ascension of ascending node (Ω, rad)
  242.      * @param alpha  an + ω, mean, eccentric or true latitude argument (rad)
  243.      * @param aDot  semi-major axis derivative (m/s)
  244.      * @param exDot d(e cos(ω))/dt, first component of circular eccentricity vector derivative
  245.      * @param eyDot d(e sin(ω))/dt, second component of circular eccentricity vector derivative
  246.      * @param iDot inclination  derivative(rad/s)
  247.      * @param raanDot right ascension of ascending node derivative (rad/s)
  248.      * @param alphaDot  d(an + ω), mean, eccentric or true latitude argument derivative (rad/s)
  249.      * @param type type of latitude argument
  250.      * @param frame the frame in which are defined the parameters
  251.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  252.      * @param date date of the orbital parameters
  253.      * @param mu central attraction coefficient (m³/s²)
  254.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  255.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  256.      */
  257.     public CircularOrbit(final double a, final double ex, final double ey,
  258.                          final double i, final double raan, final double alpha,
  259.                          final double aDot, final double exDot, final double eyDot,
  260.                          final double iDot, final double raanDot, final double alphaDot,
  261.                          final PositionAngleType type,
  262.                          final Frame frame, final AbsoluteDate date, final double mu)
  263.             throws IllegalArgumentException {
  264.         this(a, ex, ey, i, raan, alpha, aDot, exDot, eyDot, iDot, raanDot, alphaDot, type, type,
  265.                 frame, date, mu);
  266.     }

  268.     /** Creates a new instance.
  269.      * @param a  semi-major axis (m)
  270.      * @param ex e cos(ω), first component of circular eccentricity vector
  271.      * @param ey e sin(ω), second component of circular eccentricity vector
  272.      * @param i inclination (rad)
  273.      * @param raan right ascension of ascending node (Ω, rad)
  274.      * @param alpha  input latitude argument (rad)
  275.      * @param aDot  semi-major axis derivative (m/s)
  276.      * @param exDot d(e cos(ω))/dt, first component of circular eccentricity vector derivative
  277.      * @param eyDot d(e sin(ω))/dt, second component of circular eccentricity vector derivative
  278.      * @param iDot inclination  derivative(rad/s)
  279.      * @param raanDot right ascension of ascending node derivative (rad/s)
  280.      * @param alphaDot  input latitude argument derivative (rad/s)
  281.      * @param pvCoordinates the {@link PVCoordinates} in inertial frame
  282.      * @param positionAngleType type of position angle
  283.      * @param frame the frame in which are defined the parameters
  284.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  285.      * @param mu central attraction coefficient (m³/s²)
  286.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  287.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  288.      */
  289.     private CircularOrbit(final double a, final double ex, final double ey,
  290.                           final double i, final double raan, final double alpha,
  291.                           final double aDot, final double exDot, final double eyDot,
  292.                           final double iDot, final double raanDot, final double alphaDot,
  293.                           final TimeStampedPVCoordinates pvCoordinates,
  294.                           final PositionAngleType positionAngleType, final Frame frame, final double mu)
  295.         throws IllegalArgumentException {
  296.         super(pvCoordinates, frame, mu);
  297.         this.a           =  a;
  298.         this.aDot        =  aDot;
  299.         this.ex          = ex;
  300.         this.exDot       = exDot;
  301.         this.ey          = ey;
  302.         this.eyDot       = eyDot;
  303.         this.i           = i;
  304.         this.iDot        = iDot;
  305.         this.raan        = raan;
  306.         this.raanDot     = raanDot;
  307.         this.cachedAlpha = alpha;
  308.         this.cachedAlphaDot = alphaDot;
  309.         this.cachedPositionAngleType = positionAngleType;
  310.         this.serializePV = true;
  311.         this.partialPV   = null;
  312.     }

  314.     /** Constructor from Cartesian parameters.
  315.      *
  316.      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
  317.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  318.      * use {@code mu} and the position to compute the acceleration, including
  319.      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
  320.      *
  321.      * @param pvCoordinates the {@link PVCoordinates} in inertial frame
  322.      * @param frame the frame in which are defined the {@link PVCoordinates}
  323.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  324.      * @param mu central attraction coefficient (m³/s²)
  325.      * @exception IllegalArgumentException if frame is not a {@link
  326.      * Frame#isPseudoInertial pseudo-inertial frame}
  327.      */
  328.     public CircularOrbit(final TimeStampedPVCoordinates pvCoordinates, final Frame frame, final double mu)
  329.         throws IllegalArgumentException {
  330.         super(pvCoordinates, frame, mu);
  331.         this.cachedPositionAngleType = PositionAngleType.TRUE;

  333.         // compute semi-major axis
  334.         final Vector3D pvP = pvCoordinates.getPosition();
  335.         final Vector3D pvV = pvCoordinates.getVelocity();
  336.         final Vector3D pvA = pvCoordinates.getAcceleration();
  337.         final double r2 = pvP.getNormSq();
  338.         final double r  = FastMath.sqrt(r2);
  339.         final double V2 = pvV.getNormSq();
  340.         final double rV2OnMu = r * V2 / mu;
  341.         a = r / (2 - rV2OnMu);

  343.         if (!isElliptical()) {
  344.             throw new OrekitIllegalArgumentException(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS,
  345.                                                      getClass().getName());
  346.         }

  348.         // compute inclination
  349.         final Vector3D momentum = pvCoordinates.getMomentum();
  350.         i = Vector3D.angle(momentum, Vector3D.PLUS_K);

  352.         // compute right ascension of ascending node
  353.         final Vector3D node  = Vector3D.crossProduct(Vector3D.PLUS_K, momentum);
  354.         raan = FastMath.atan2(node.getY(), node.getX());

  356.         // 2D-coordinates in the canonical frame
  357.         final SinCos scRaan = FastMath.sinCos(raan);
  358.         final SinCos scI    = FastMath.sinCos(i);
  359.         final double xP     = pvP.getX();
  360.         final double yP     = pvP.getY();
  361.         final double zP     = pvP.getZ();
  362.         final double x2     = (xP * scRaan.cos() + yP * scRaan.sin()) / a;
  363.         final double y2     = ((yP * scRaan.cos() - xP * scRaan.sin()) * scI.cos() + zP * scI.sin()) / a;

  365.         // compute eccentricity vector
  366.         final double eSE    = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(mu * a);
  367.         final double eCE    = rV2OnMu - 1;
  368.         final double e2     = eCE * eCE + eSE * eSE;
  369.         final double f      = eCE - e2;
  370.         final double g      = FastMath.sqrt(1 - e2) * eSE;
  371.         final double aOnR   = a / r;
  372.         final double a2OnR2 = aOnR * aOnR;
  373.         ex = a2OnR2 * (f * x2 + g * y2);
  374.         ey = a2OnR2 * (f * y2 - g * x2);

  376.         // compute latitude argument
  377.         final double beta = 1 / (1 + FastMath.sqrt(1 - ex * ex - ey * ey));
  378.         cachedAlpha = CircularLatitudeArgumentUtility.eccentricToTrue(ex, ey, FastMath.atan2(y2 + ey + eSE * beta * ex, x2 + ex - eSE * beta * ey));

  380.         partialPV   = pvCoordinates;

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

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

  388.             final Vector3D keplerianAcceleration    = new Vector3D(-mu / (r * r2), pvP);
  389.             final Vector3D nonKeplerianAcceleration = pvA.subtract(keplerianAcceleration);
  390.             final double   aX                       = nonKeplerianAcceleration.getX();
  391.             final double   aY                       = nonKeplerianAcceleration.getY();
  392.             final double   aZ                       = nonKeplerianAcceleration.getZ();
  393.             aDot    = jacobian[0][3] * aX + jacobian[0][4] * aY + jacobian[0][5] * aZ;
  394.             exDot   = jacobian[1][3] * aX + jacobian[1][4] * aY + jacobian[1][5] * aZ;
  395.             eyDot   = jacobian[2][3] * aX + jacobian[2][4] * aY + jacobian[2][5] * aZ;
  396.             iDot    = jacobian[3][3] * aX + jacobian[3][4] * aY + jacobian[3][5] * aZ;
  397.             raanDot = jacobian[4][3] * aX + jacobian[4][4] * aY + jacobian[4][5] * aZ;

  399.             // in order to compute latitude argument derivative, we must compute
  400.             // mean latitude argument derivative including Keplerian motion and convert to true latitude argument
  401.             final double alphaMDot = getKeplerianMeanMotion() +
  402.                                      jacobian[5][3] * aX + jacobian[5][4] * aY + jacobian[5][5] * aZ;
  403.             final UnivariateDerivative1 exUD     = new UnivariateDerivative1(ex, exDot);
  404.             final UnivariateDerivative1 eyUD     = new UnivariateDerivative1(ey, eyDot);
  405.             final UnivariateDerivative1 alphaMUD = new UnivariateDerivative1(getAlphaM(), alphaMDot);
  406.             final UnivariateDerivative1 alphavUD = FieldCircularLatitudeArgumentUtility.meanToTrue(exUD, eyUD, alphaMUD);
  407.             cachedAlphaDot = alphavUD.getFirstDerivative();

  409.         } else {
  410.             // acceleration is either almost zero or NaN,
  411.             // we assume acceleration was not known
  412.             // we don't set up derivatives
  413.             aDot      = Double.NaN;
  414.             exDot     = Double.NaN;
  415.             eyDot     = Double.NaN;
  416.             iDot      = Double.NaN;
  417.             raanDot   = Double.NaN;
  418.             cachedAlphaDot = Double.NaN;
  419.         }

  421.         serializePV = true;

  423.     }

  425.     /** Constructor from Cartesian parameters.
  426.      *
  427.      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
  428.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  429.      * use {@code mu} and the position to compute the acceleration, including
  430.      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
  431.      *
  432.      * @param pvCoordinates the {@link PVCoordinates} in inertial frame
  433.      * @param frame the frame in which are defined the {@link PVCoordinates}
  434.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  435.      * @param date date of the orbital parameters
  436.      * @param mu central attraction coefficient (m³/s²)
  437.      * @exception IllegalArgumentException if frame is not a {@link
  438.      * Frame#isPseudoInertial pseudo-inertial frame}
  439.      */
  440.     public CircularOrbit(final PVCoordinates pvCoordinates, final Frame frame,
  441.                          final AbsoluteDate date, final double mu)
  442.         throws IllegalArgumentException {
  443.         this(new TimeStampedPVCoordinates(date, pvCoordinates), frame, mu);
  444.     }

  446.     /** Constructor from any kind of orbital parameters.
  447.      * @param op orbital parameters to copy
  448.      */
  449.     public CircularOrbit(final Orbit op) {

  451.         super(op.getFrame(), op.getDate(), op.getMu());

  453.         a    = op.getA();
  454.         i    = op.getI();
  455.         final double hx = op.getHx();
  456.         final double hy = op.getHy();
  457.         final double h2 = hx * hx + hy * hy;
  458.         final double h  = FastMath.sqrt(h2);
  459.         raan = FastMath.atan2(hy, hx);
  460.         final SinCos scRaan  = FastMath.sinCos(raan);
  461.         final double cosRaan = h == 0 ? scRaan.cos() : hx / h;
  462.         final double sinRaan = h == 0 ? scRaan.sin() : hy / h;
  463.         final double equiEx = op.getEquinoctialEx();
  464.         final double equiEy = op.getEquinoctialEy();
  465.         ex     = equiEx * cosRaan + equiEy * sinRaan;
  466.         ey     = equiEy * cosRaan - equiEx * sinRaan;
  467.         cachedPositionAngleType = PositionAngleType.TRUE;
  468.         cachedAlpha = op.getLv() - raan;

  470.         if (op.hasDerivatives()) {
  471.             aDot    = op.getADot();
  472.             final double hxDot = op.getHxDot();
  473.             final double hyDot = op.getHyDot();
  474.             iDot    = 2 * (cosRaan * hxDot + sinRaan * hyDot) / (1 + h2);
  475.             raanDot = (hx * hyDot - hy * hxDot) / h2;
  476.             final double equiExDot = op.getEquinoctialExDot();
  477.             final double equiEyDot = op.getEquinoctialEyDot();
  478.             exDot   = (equiExDot + equiEy * raanDot) * cosRaan +
  479.                       (equiEyDot - equiEx * raanDot) * sinRaan;
  480.             eyDot   = (equiEyDot - equiEx * raanDot) * cosRaan -
  481.                       (equiExDot + equiEy * raanDot) * sinRaan;
  482.             cachedAlphaDot = op.getLvDot() - raanDot;
  483.         } else {
  484.             aDot      = Double.NaN;
  485.             exDot     = Double.NaN;
  486.             eyDot     = Double.NaN;
  487.             iDot      = Double.NaN;
  488.             raanDot   = Double.NaN;
  489.             cachedAlphaDot = Double.NaN;
  490.         }

  492.         serializePV = false;
  493.         partialPV   = null;

  495.     }

  497.     /** {@inheritDoc} */
  498.     @Override
  499.     public OrbitType getType() {
  500.         return OrbitType.CIRCULAR;
  501.     }

  503.     /** {@inheritDoc} */
  504.     @Override
  505.     public double getA() {
  506.         return a;
  507.     }

  509.     /** {@inheritDoc} */
  510.     @Override
  511.     public double getADot() {
  512.         return aDot;
  513.     }

  515.     /** {@inheritDoc} */
  516.     @Override
  517.     public double getEquinoctialEx() {
  518.         final SinCos sc = FastMath.sinCos(raan);
  519.         return ex * sc.cos() - ey * sc.sin();
  520.     }

  522.     /** {@inheritDoc} */
  523.     @Override
  524.     public double getEquinoctialExDot() {
  525.         final SinCos sc = FastMath.sinCos(raan);
  526.         return (exDot - ey * raanDot) * sc.cos() - (eyDot + ex * raanDot) * sc.sin();
  527.     }

  529.     /** {@inheritDoc} */
  530.     @Override
  531.     public double getEquinoctialEy() {
  532.         final SinCos sc = FastMath.sinCos(raan);
  533.         return ey * sc.cos() + ex * sc.sin();
  534.     }

  536.     /** {@inheritDoc} */
  537.     @Override
  538.     public double getEquinoctialEyDot() {
  539.         final SinCos sc = FastMath.sinCos(raan);
  540.         return (eyDot + ex * raanDot) * sc.cos() + (exDot - ey * raanDot) * sc.sin();
  541.     }

  543.     /** Get the first component of the circular eccentricity vector.
  544.      * @return ex = e cos(ω), first component of the circular eccentricity vector
  545.      */
  546.     public double getCircularEx() {
  547.         return ex;
  548.     }

  550.     /** Get the first component of the circular eccentricity vector derivative.
  551.      * @return ex = e cos(ω), first component of the circular eccentricity vector derivative
  552.      * @since 9.0
  553.      */
  554.     public double getCircularExDot() {
  555.         return exDot;
  556.     }

  558.     /** Get the second component of the circular eccentricity vector.
  559.      * @return ey = e sin(ω), second component of the circular eccentricity vector
  560.      */
  561.     public double getCircularEy() {
  562.         return ey;
  563.     }

  565.     /** Get the second component of the circular eccentricity vector derivative.
  566.      * @return ey = e sin(ω), second component of the circular eccentricity vector derivative
  567.      */
  568.     public double getCircularEyDot() {
  569.         return eyDot;
  570.     }

  572.     /** {@inheritDoc} */
  573.     @Override
  574.     public double getHx() {
  575.         // Check for equatorial retrograde orbit
  576.         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
  577.             return Double.NaN;
  578.         }
  579.         return  FastMath.cos(raan) * FastMath.tan(i / 2);
  580.     }

  582.     /** {@inheritDoc} */
  583.     @Override
  584.     public double getHxDot() {
  585.         // Check for equatorial retrograde orbit
  586.         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
  587.             return Double.NaN;
  588.         }
  589.         final SinCos sc  = FastMath.sinCos(raan);
  590.         final double tan = FastMath.tan(0.5 * i);
  591.         return 0.5 * sc.cos() * (1 + tan * tan) * iDot - sc.sin() * tan * raanDot;
  592.     }

  594.     /** {@inheritDoc} */
  595.     @Override
  596.     public double getHy() {
  597.         // Check for equatorial retrograde orbit
  598.         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
  599.             return Double.NaN;
  600.         }
  601.         return  FastMath.sin(raan) * FastMath.tan(i / 2);
  602.     }

  604.     /** {@inheritDoc} */
  605.     @Override
  606.     public double getHyDot() {
  607.         // Check for equatorial retrograde orbit
  608.         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
  609.             return Double.NaN;
  610.         }
  611.         final SinCos sc  = FastMath.sinCos(raan);
  612.         final double tan = FastMath.tan(0.5 * i);
  613.         return 0.5 * sc.sin() * (1 + tan * tan) * iDot + sc.cos() * tan * raanDot;
  614.     }

  616.     /** Get the true latitude argument.
  617.      * @return v + ω true latitude argument (rad)
  618.      */
  619.     public double getAlphaV() {
  620.         switch (cachedPositionAngleType) {
  621.             case TRUE:
  622.                 return cachedAlpha;

  624.             case ECCENTRIC:
  625.                 return CircularLatitudeArgumentUtility.eccentricToTrue(ex, ey, cachedAlpha);

  627.             case MEAN:
  628.                 return CircularLatitudeArgumentUtility.meanToTrue(ex, ey, cachedAlpha);

  630.             default:
  631.                 throw new OrekitInternalError(null);
  632.         }
  633.     }

  635.     /** Get the true latitude argument derivative.
  636.      * <p>
  637.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  638.      * </p>
  639.      * @return v + ω true latitude argument derivative (rad/s)
  640.      * @since 9.0
  641.      */
  642.     public double getAlphaVDot() {
  643.         switch (cachedPositionAngleType) {
  644.             case ECCENTRIC:
  645.                 final UnivariateDerivative1 alphaEUD = new UnivariateDerivative1(cachedAlpha, cachedAlphaDot);
  646.                 final UnivariateDerivative1 exUD     = new UnivariateDerivative1(ex,     exDot);
  647.                 final UnivariateDerivative1 eyUD     = new UnivariateDerivative1(ey,     eyDot);
  648.                 final UnivariateDerivative1 alphaVUD = FieldCircularLatitudeArgumentUtility.eccentricToTrue(exUD, eyUD,
  649.                         alphaEUD);
  650.                 return alphaVUD.getFirstDerivative();

  652.             case TRUE:
  653.                 return cachedAlphaDot;

  655.             case MEAN:
  656.                 final UnivariateDerivative1 alphaMUD = new UnivariateDerivative1(cachedAlpha, cachedAlphaDot);
  657.                 final UnivariateDerivative1 exUD2    = new UnivariateDerivative1(ex,     exDot);
  658.                 final UnivariateDerivative1 eyUD2    = new UnivariateDerivative1(ey,     eyDot);
  659.                 final UnivariateDerivative1 alphaVUD2 = FieldCircularLatitudeArgumentUtility.meanToTrue(exUD2,
  660.                         eyUD2, alphaMUD);
  661.                 return alphaVUD2.getFirstDerivative();

  663.             default:
  664.                 throw new OrekitInternalError(null);
  665.         }
  666.     }

  668.     /** Get the eccentric latitude argument.
  669.      * @return E + ω eccentric latitude argument (rad)
  670.      */
  671.     public double getAlphaE() {
  672.         switch (cachedPositionAngleType) {
  673.             case TRUE:
  674.                 return CircularLatitudeArgumentUtility.trueToEccentric(ex, ey, cachedAlpha);

  676.             case ECCENTRIC:
  677.                 return cachedAlpha;

  679.             case MEAN:
  680.                 return CircularLatitudeArgumentUtility.meanToEccentric(ex, ey, cachedAlpha);

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

  687.     /** Get the eccentric latitude argument derivative.
  688.      * <p>
  689.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  690.      * </p>
  691.      * @return d(E + ω)/dt eccentric latitude argument derivative (rad/s)
  692.      * @since 9.0
  693.      */
  694.     public double getAlphaEDot() {
  695.         switch (cachedPositionAngleType) {
  696.             case TRUE:
  697.                 final UnivariateDerivative1 alphaVUD = new UnivariateDerivative1(cachedAlpha, cachedAlphaDot);
  698.                 final UnivariateDerivative1 exUD     = new UnivariateDerivative1(ex,     exDot);
  699.                 final UnivariateDerivative1 eyUD     = new UnivariateDerivative1(ey,     eyDot);
  700.                 final UnivariateDerivative1 alphaEUD = FieldCircularLatitudeArgumentUtility.trueToEccentric(exUD, eyUD,
  701.                         alphaVUD);
  702.                 return alphaEUD.getFirstDerivative();

  704.             case ECCENTRIC:
  705.                 return cachedAlphaDot;

  707.             case MEAN:
  708.                 final UnivariateDerivative1 alphaMUD = new UnivariateDerivative1(cachedAlpha, cachedAlphaDot);
  709.                 final UnivariateDerivative1 exUD2    = new UnivariateDerivative1(ex,     exDot);
  710.                 final UnivariateDerivative1 eyUD2    = new UnivariateDerivative1(ey,     eyDot);
  711.                 final UnivariateDerivative1 alphaVUD2 = FieldCircularLatitudeArgumentUtility.meanToEccentric(exUD2,
  712.                         eyUD2, alphaMUD);
  713.                 return alphaVUD2.getFirstDerivative();

  715.             default:
  716.                 throw new OrekitInternalError(null);
  717.         }
  718.     }

  720.     /** Get the mean latitude argument.
  721.      * @return M + ω mean latitude argument (rad)
  722.      */
  723.     public double getAlphaM() {
  724.         switch (cachedPositionAngleType) {
  725.             case TRUE:
  726.                 return CircularLatitudeArgumentUtility.trueToMean(ex, ey, cachedAlpha);

  728.             case MEAN:
  729.                 return cachedAlpha;

  731.             case ECCENTRIC:
  732.                 return CircularLatitudeArgumentUtility.eccentricToMean(ex, ey, cachedAlpha);

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

  739.     /** Get the mean latitude argument derivative.
  740.      * <p>
  741.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  742.      * </p>
  743.      * @return d(M + ω)/dt mean latitude argument derivative (rad/s)
  744.      * @since 9.0
  745.      */
  746.     public double getAlphaMDot() {
  747.         switch (cachedPositionAngleType) {
  748.             case TRUE:
  749.                 final UnivariateDerivative1 alphaVUD = new UnivariateDerivative1(cachedAlpha, cachedAlphaDot);
  750.                 final UnivariateDerivative1 exUD     = new UnivariateDerivative1(ex,     exDot);
  751.                 final UnivariateDerivative1 eyUD     = new UnivariateDerivative1(ey,     eyDot);
  752.                 final UnivariateDerivative1 alphaMUD = FieldCircularLatitudeArgumentUtility.trueToMean(exUD, eyUD,
  753.                         alphaVUD);
  754.                 return alphaMUD.getFirstDerivative();

  756.             case MEAN:
  757.                 return cachedAlphaDot;

  759.             case ECCENTRIC:
  760.                 final UnivariateDerivative1 alphaEUD = new UnivariateDerivative1(cachedAlpha, cachedAlphaDot);
  761.                 final UnivariateDerivative1 exUD2    = new UnivariateDerivative1(ex,     exDot);
  762.                 final UnivariateDerivative1 eyUD2    = new UnivariateDerivative1(ey,     eyDot);
  763.                 final UnivariateDerivative1 alphaMUD2 = FieldCircularLatitudeArgumentUtility.eccentricToMean(exUD2,
  764.                         eyUD2, alphaEUD);
  765.                 return alphaMUD2.getFirstDerivative();

  767.             default:
  768.                 throw new OrekitInternalError(null);
  769.         }
  770.     }

  772.     /** Get the latitude argument.
  773.      * @param type type of the angle
  774.      * @return latitude argument (rad)
  775.      */
  776.     public double getAlpha(final PositionAngleType type) {
  777.         return (type == PositionAngleType.MEAN) ? getAlphaM() :
  778.                                               ((type == PositionAngleType.ECCENTRIC) ? getAlphaE() :
  779.                                                                                    getAlphaV());
  780.     }

  782.     /** Get the latitude argument derivative.
  783.      * <p>
  784.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  785.      * </p>
  786.      * @param type type of the angle
  787.      * @return latitude argument derivative (rad/s)
  788.      * @since 9.0
  789.      */
  790.     public double getAlphaDot(final PositionAngleType type) {
  791.         return (type == PositionAngleType.MEAN) ? getAlphaMDot() :
  792.                                               ((type == PositionAngleType.ECCENTRIC) ? getAlphaEDot() :
  793.                                                                                    getAlphaVDot());
  794.     }

  796.     /** Computes the true latitude argument from the eccentric latitude argument.
  797.      * @param alphaE = E + ω eccentric latitude argument (rad)
  798.      * @param ex e cos(ω), first component of circular eccentricity vector
  799.      * @param ey e sin(ω), second component of circular eccentricity vector
  800.      * @return the true latitude argument.
  801.      */
  802.     @Deprecated
  803.     public static double eccentricToTrue(final double alphaE, final double ex, final double ey) {
  804.         return CircularLatitudeArgumentUtility.eccentricToTrue(ex, ey, alphaE);
  805.     }

  807.     /** Computes the eccentric latitude argument from the true latitude argument.
  808.      * @param alphaV = V + ω true latitude argument (rad)
  809.      * @param ex e cos(ω), first component of circular eccentricity vector
  810.      * @param ey e sin(ω), second component of circular eccentricity vector
  811.      * @return the eccentric latitude argument.
  812.      */
  813.     @Deprecated
  814.     public static double trueToEccentric(final double alphaV, final double ex, final double ey) {
  815.         return CircularLatitudeArgumentUtility.trueToEccentric(ex, ey, alphaV);
  816.     }

  818.     /** Computes the eccentric latitude argument from the mean latitude argument.
  819.      * @param alphaM = M + ω  mean latitude argument (rad)
  820.      * @param ex e cos(ω), first component of circular eccentricity vector
  821.      * @param ey e sin(ω), second component of circular eccentricity vector
  822.      * @return the eccentric latitude argument.
  823.      */
  824.     @Deprecated
  825.     public static double meanToEccentric(final double alphaM, final double ex, final double ey) {
  826.         return CircularLatitudeArgumentUtility.meanToEccentric(ex, ey, alphaM);
  827.     }

  829.     /** Computes the mean latitude argument from the eccentric latitude argument.
  830.      * @param alphaE = E + ω  mean latitude argument (rad)
  831.      * @param ex e cos(ω), first component of circular eccentricity vector
  832.      * @param ey e sin(ω), second component of circular eccentricity vector
  833.      * @return the mean latitude argument.
  834.      */
  835.     @Deprecated
  836.     public static double eccentricToMean(final double alphaE, final double ex, final double ey) {
  837.         return CircularLatitudeArgumentUtility.eccentricToMean(ex, ey, alphaE);
  838.     }

  840.     /** {@inheritDoc} */
  841.     @Override
  842.     public double getE() {
  843.         return FastMath.sqrt(ex * ex + ey * ey);
  844.     }

  846.     /** {@inheritDoc} */
  847.     @Override
  848.     public double getEDot() {
  849.         return (ex * exDot + ey * eyDot) / getE();
  850.     }

  852.     /** {@inheritDoc} */
  853.     @Override
  854.     public double getI() {
  855.         return i;
  856.     }

  858.     /** {@inheritDoc} */
  859.     @Override
  860.     public double getIDot() {
  861.         return iDot;
  862.     }

  864.     /** Get the right ascension of the ascending node.
  865.      * @return right ascension of the ascending node (rad)
  866.      */
  867.     public double getRightAscensionOfAscendingNode() {
  868.         return raan;
  869.     }

  871.     /** Get the right ascension of the ascending node derivative.
  872.      * <p>
  873.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  874.      * </p>
  875.      * @return right ascension of the ascending node derivative (rad/s)
  876.      * @since 9.0
  877.      */
  878.     public double getRightAscensionOfAscendingNodeDot() {
  879.         return raanDot;
  880.     }

  882.     /** {@inheritDoc} */
  883.     @Override
  884.     public double getLv() {
  885.         return getAlphaV() + raan;
  886.     }

  888.     /** {@inheritDoc} */
  889.     @Override
  890.     public double getLvDot() {
  891.         return getAlphaVDot() + raanDot;
  892.     }

  894.     /** {@inheritDoc} */
  895.     @Override
  896.     public double getLE() {
  897.         return getAlphaE() + raan;
  898.     }

  900.     /** {@inheritDoc} */
  901.     @Override
  902.     public double getLEDot() {
  903.         return getAlphaEDot() + raanDot;
  904.     }

  906.     /** {@inheritDoc} */
  907.     @Override
  908.     public double getLM() {
  909.         return getAlphaM() + raan;
  910.     }

  912.     /** {@inheritDoc} */
  913.     @Override
  914.     public double getLMDot() {
  915.         return getAlphaMDot() + raanDot;
  916.     }

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

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

  927.         // get equinoctial parameters
  928.         final double equEx = getEquinoctialEx();
  929.         final double equEy = getEquinoctialEy();
  930.         final double hx = getHx();
  931.         final double hy = getHy();
  932.         final double lE = getLE();

  934.         // inclination-related intermediate parameters
  935.         final double hx2   = hx * hx;
  936.         final double hy2   = hy * hy;
  937.         final double factH = 1. / (1 + hx2 + hy2);

  939.         // reference axes defining the orbital plane
  940.         final double ux = (1 + hx2 - hy2) * factH;
  941.         final double uy =  2 * hx * hy * factH;
  942.         final double uz = -2 * hy * factH;

  944.         final double vx = uy;
  945.         final double vy = (1 - hx2 + hy2) * factH;
  946.         final double vz =  2 * hx * factH;

  948.         // eccentricity-related intermediate parameters
  949.         final double exey = equEx * equEy;
  950.         final double ex2  = equEx * equEx;
  951.         final double ey2  = equEy * equEy;
  952.         final double e2   = ex2 + ey2;
  953.         final double eta  = 1 + FastMath.sqrt(1 - e2);
  954.         final double beta = 1. / eta;

  956.         // eccentric latitude argument
  957.         final SinCos scLe   = FastMath.sinCos(lE);
  958.         final double cLe    = scLe.cos();
  959.         final double sLe    = scLe.sin();
  960.         final double exCeyS = equEx * cLe + equEy * sLe;

  962.         // coordinates of position and velocity in the orbital plane
  963.         final double x      = a * ((1 - beta * ey2) * cLe + beta * exey * sLe - equEx);
  964.         final double y      = a * ((1 - beta * ex2) * sLe + beta * exey * cLe - equEy);

  966.         final double factor = FastMath.sqrt(getMu() / a) / (1 - exCeyS);
  967.         final double xdot   = factor * (-sLe + beta * equEy * exCeyS);
  968.         final double ydot   = factor * ( cLe - beta * equEx * exCeyS);

  970.         final Vector3D position =
  971.                         new Vector3D(x * ux + y * vx, x * uy + y * vy, x * uz + y * vz);
  972.         final Vector3D velocity =
  973.                         new Vector3D(xdot * ux + ydot * vx, xdot * uy + ydot * vy, xdot * uz + ydot * vz);

  975.         partialPV = new PVCoordinates(position, velocity);

  977.     }

  979.     /** Initialize cached alpha with rate.
  980.      * @param alpha input alpha
  981.      * @param alphaDot rate of input alpha
  982.      * @param inputType position angle type passed as input
  983.      * @return alpha to cache with rate
  984.      * @since 12.1
  985.      */
  986.     private UnivariateDerivative1 initializeCachedAlpha(final double alpha, final double alphaDot,
  987.                                                         final PositionAngleType inputType) {
  988.         if (cachedPositionAngleType == inputType) {
  989.             return new UnivariateDerivative1(alpha, alphaDot);

  991.         } else {
  992.             final UnivariateDerivative1 exUD = new UnivariateDerivative1(ex, exDot);
  993.             final UnivariateDerivative1 eyUD = new UnivariateDerivative1(ey, eyDot);
  994.             final UnivariateDerivative1 alphaUD = new UnivariateDerivative1(alpha, alphaDot);

  996.             switch (cachedPositionAngleType) {

  998.                 case ECCENTRIC:
  999.                     if (inputType == PositionAngleType.MEAN) {
  1000.                         return FieldCircularLatitudeArgumentUtility.meanToEccentric(exUD, eyUD, alphaUD);
  1001.                     } else {
  1002.                         return FieldCircularLatitudeArgumentUtility.trueToEccentric(exUD, eyUD, alphaUD);
  1003.                     }

  1005.                 case TRUE:
  1006.                     if (inputType == PositionAngleType.MEAN) {
  1007.                         return FieldCircularLatitudeArgumentUtility.meanToTrue(exUD, eyUD, alphaUD);
  1008.                     } else {
  1009.                         return FieldCircularLatitudeArgumentUtility.eccentricToTrue(exUD, eyUD, alphaUD);
  1010.                     }

  1012.                 case MEAN:
  1013.                     if (inputType == PositionAngleType.TRUE) {
  1014.                         return FieldCircularLatitudeArgumentUtility.trueToMean(exUD, eyUD, alphaUD);
  1015.                     } else {
  1016.                         return FieldCircularLatitudeArgumentUtility.eccentricToMean(exUD, eyUD, alphaUD);
  1017.                     }

  1019.                 default:
  1020.                     throw new OrekitInternalError(null);

  1022.             }

  1024.         }

  1026.     }

  1028.     /** Initialize cached alpha.
  1029.      * @param alpha input alpha
  1030.      * @param positionAngleType position angle type passed as input
  1031.      * @return alpha to cache
  1032.      * @since 12.1
  1033.      */
  1034.     private double initializeCachedAlpha(final double alpha, final PositionAngleType positionAngleType) {
  1035.         return CircularLatitudeArgumentUtility.convertAlpha(positionAngleType, alpha, ex, ey, cachedPositionAngleType);
  1036.     }

  1038.     /** Compute non-Keplerian part of the acceleration from first time derivatives.
  1039.      * <p>
  1040.      * This method should be called only when {@link #hasDerivatives()} returns true.
  1041.      * </p>
  1042.      * @return non-Keplerian part of the acceleration
  1043.      */
  1044.     private Vector3D nonKeplerianAcceleration() {

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

  1049.         final double nonKeplerianMeanMotion = getAlphaMDot() - getKeplerianMeanMotion();
  1050.         final double nonKeplerianAx = dCdP[3][0] * aDot    + dCdP[3][1] * exDot   + dCdP[3][2] * eyDot   +
  1051.                                       dCdP[3][3] * iDot    + dCdP[3][4] * raanDot + dCdP[3][5] * nonKeplerianMeanMotion;
  1052.         final double nonKeplerianAy = dCdP[4][0] * aDot    + dCdP[4][1] * exDot   + dCdP[4][2] * eyDot   +
  1053.                                       dCdP[4][3] * iDot    + dCdP[4][4] * raanDot + dCdP[4][5] * nonKeplerianMeanMotion;
  1054.         final double nonKeplerianAz = dCdP[5][0] * aDot    + dCdP[5][1] * exDot   + dCdP[5][2] * eyDot   +
  1055.                                       dCdP[5][3] * iDot    + dCdP[5][4] * raanDot + dCdP[5][5] * nonKeplerianMeanMotion;

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

  1059.     }

  1061.     /** {@inheritDoc} */
  1062.     @Override
  1063.     protected Vector3D initPosition() {

  1065.         // get equinoctial parameters
  1066.         final double equEx = getEquinoctialEx();
  1067.         final double equEy = getEquinoctialEy();
  1068.         final double hx = getHx();
  1069.         final double hy = getHy();
  1070.         final double lE = getLE();

  1072.         // inclination-related intermediate parameters
  1073.         final double hx2   = hx * hx;
  1074.         final double hy2   = hy * hy;
  1075.         final double factH = 1. / (1 + hx2 + hy2);

  1077.         // reference axes defining the orbital plane
  1078.         final double ux = (1 + hx2 - hy2) * factH;
  1079.         final double uy =  2 * hx * hy * factH;
  1080.         final double uz = -2 * hy * factH;

  1082.         final double vx = uy;
  1083.         final double vy = (1 - hx2 + hy2) * factH;
  1084.         final double vz =  2 * hx * factH;

  1086.         // eccentricity-related intermediate parameters
  1087.         final double exey = equEx * equEy;
  1088.         final double ex2  = equEx * equEx;
  1089.         final double ey2  = equEy * equEy;
  1090.         final double e2   = ex2 + ey2;
  1091.         final double eta  = 1 + FastMath.sqrt(1 - e2);
  1092.         final double beta = 1. / eta;

  1094.         // eccentric latitude argument
  1095.         final SinCos scLe   = FastMath.sinCos(lE);
  1096.         final double cLe    = scLe.cos();
  1097.         final double sLe    = scLe.sin();

  1099.         // coordinates of position and velocity in the orbital plane
  1100.         final double x      = a * ((1 - beta * ey2) * cLe + beta * exey * sLe - equEx);
  1101.         final double y      = a * ((1 - beta * ex2) * sLe + beta * exey * cLe - equEy);

  1103.         return new Vector3D(x * ux + y * vx, x * uy + y * vy, x * uz + y * vz);

  1105.     }

  1107.     /** {@inheritDoc} */
  1108.     @Override
  1109.     protected TimeStampedPVCoordinates initPVCoordinates() {

  1111.         // position and velocity
  1112.         computePVWithoutA();

  1114.         // acceleration
  1115.         final double r2 = partialPV.getPosition().getNormSq();
  1116.         final Vector3D keplerianAcceleration = new Vector3D(-getMu() / (r2 * FastMath.sqrt(r2)), partialPV.getPosition());
  1117.         final Vector3D acceleration = hasDerivatives() ?
  1118.                                       keplerianAcceleration.add(nonKeplerianAcceleration()) :
  1119.                                       keplerianAcceleration;

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

  1123.     }

  1125.     /** {@inheritDoc} */
  1126.     @Override
  1127.     public CircularOrbit shiftedBy(final double dt) {

  1129.         // use Keplerian-only motion
  1130.         final CircularOrbit keplerianShifted = new CircularOrbit(a, ex, ey, i, raan,
  1131.                                                                  getAlphaM() + getKeplerianMeanMotion() * dt,
  1132.                                                                  PositionAngleType.MEAN, cachedPositionAngleType,
  1133.                                                                  getFrame(), getDate().shiftedBy(dt), getMu());

  1135.         if (hasDerivatives()) {

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

  1140.             // add quadratic effect of non-Keplerian acceleration to Keplerian-only shift
  1141.             keplerianShifted.computePVWithoutA();
  1142.             final Vector3D fixedP   = new Vector3D(1, keplerianShifted.partialPV.getPosition(),
  1143.                                                    0.5 * dt * dt, nonKeplerianAcceleration);
  1144.             final double   fixedR2 = fixedP.getNormSq();
  1145.             final double   fixedR  = FastMath.sqrt(fixedR2);
  1146.             final Vector3D fixedV  = new Vector3D(1, keplerianShifted.partialPV.getVelocity(),
  1147.                                                   dt, nonKeplerianAcceleration);
  1148.             final Vector3D fixedA  = new Vector3D(-getMu() / (fixedR2 * fixedR), keplerianShifted.partialPV.getPosition(),
  1149.                                                   1, nonKeplerianAcceleration);

  1151.             // build a new orbit, taking non-Keplerian acceleration into account
  1152.             return new CircularOrbit(new TimeStampedPVCoordinates(keplerianShifted.getDate(),
  1153.                                                                   fixedP, fixedV, fixedA),
  1154.                                      keplerianShifted.getFrame(), keplerianShifted.getMu());

  1156.         } else {
  1157.             // Keplerian-only motion is all we can do
  1158.             return keplerianShifted;
  1159.         }

  1161.     }

  1163.     /** {@inheritDoc} */
  1164.     @Override
  1165.     protected double[][] computeJacobianMeanWrtCartesian() {


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

  1170.         computePVWithoutA();
  1171.         final Vector3D position = partialPV.getPosition();
  1172.         final Vector3D velocity = partialPV.getVelocity();
  1173.         final double x          = position.getX();
  1174.         final double y          = position.getY();
  1175.         final double z          = position.getZ();
  1176.         final double vx         = velocity.getX();
  1177.         final double vy         = velocity.getY();
  1178.         final double vz         = velocity.getZ();
  1179.         final double pv         = Vector3D.dotProduct(position, velocity);
  1180.         final double r2         = position.getNormSq();
  1181.         final double r          = FastMath.sqrt(r2);
  1182.         final double v2         = velocity.getNormSq();

  1184.         final double mu         = getMu();
  1185.         final double oOsqrtMuA  = 1 / FastMath.sqrt(mu * a);
  1186.         final double rOa        = r / a;
  1187.         final double aOr        = a / r;
  1188.         final double aOr2       = a / r2;
  1189.         final double a2         = a * a;

  1191.         final double ex2        = ex * ex;
  1192.         final double ey2        = ey * ey;
  1193.         final double e2         = ex2 + ey2;
  1194.         final double epsilon    = FastMath.sqrt(1 - e2);
  1195.         final double beta       = 1 / (1 + epsilon);

  1197.         final double eCosE      = 1 - rOa;
  1198.         final double eSinE      = pv * oOsqrtMuA;

  1200.         final SinCos scI    = FastMath.sinCos(i);
  1201.         final SinCos scRaan = FastMath.sinCos(raan);
  1202.         final double cosI       = scI.cos();
  1203.         final double sinI       = scI.sin();
  1204.         final double cosRaan    = scRaan.cos();
  1205.         final double sinRaan    = scRaan.sin();

  1207.         // da
  1208.         fillHalfRow(2 * aOr * aOr2, position, jacobian[0], 0);
  1209.         fillHalfRow(2 * a2 / mu, velocity, jacobian[0], 3);

  1211.         // differentials of the normalized momentum
  1212.         final Vector3D danP = new Vector3D(v2, position, -pv, velocity);
  1213.         final Vector3D danV = new Vector3D(r2, velocity, -pv, position);
  1214.         final double recip  = 1 / partialPV.getMomentum().getNorm();
  1215.         final double recip2 = recip * recip;
  1216.         final Vector3D dwXP = new Vector3D(recip, new Vector3D(  0,  vz, -vy), -recip2 * sinRaan * sinI, danP);
  1217.         final Vector3D dwYP = new Vector3D(recip, new Vector3D(-vz,   0,  vx),  recip2 * cosRaan * sinI, danP);
  1218.         final Vector3D dwZP = new Vector3D(recip, new Vector3D( vy, -vx,   0), -recip2 * cosI,           danP);
  1219.         final Vector3D dwXV = new Vector3D(recip, new Vector3D(  0,  -z,   y), -recip2 * sinRaan * sinI, danV);
  1220.         final Vector3D dwYV = new Vector3D(recip, new Vector3D(  z,   0,  -x),  recip2 * cosRaan * sinI, danV);
  1221.         final Vector3D dwZV = new Vector3D(recip, new Vector3D( -y,   x,   0), -recip2 * cosI,           danV);

  1223.         // di
  1224.         fillHalfRow(sinRaan * cosI, dwXP, -cosRaan * cosI, dwYP, -sinI, dwZP, jacobian[3], 0);
  1225.         fillHalfRow(sinRaan * cosI, dwXV, -cosRaan * cosI, dwYV, -sinI, dwZV, jacobian[3], 3);

  1227.         // dRaan
  1228.         fillHalfRow(sinRaan / sinI, dwYP, cosRaan / sinI, dwXP, jacobian[4], 0);
  1229.         fillHalfRow(sinRaan / sinI, dwYV, cosRaan / sinI, dwXV, jacobian[4], 3);

  1231.         // orbital frame: (p, q, w) p along ascending node, w along momentum
  1232.         // the coordinates of the spacecraft in this frame are: (u, v, 0)
  1233.         final double u     =  x * cosRaan + y * sinRaan;
  1234.         final double cv    = -x * sinRaan + y * cosRaan;
  1235.         final double v     = cv * cosI + z * sinI;

  1237.         // du
  1238.         final Vector3D duP = new Vector3D(cv * cosRaan / sinI, dwXP,
  1239.                                           cv * sinRaan / sinI, dwYP,
  1240.                                           1, new Vector3D(cosRaan, sinRaan, 0));
  1241.         final Vector3D duV = new Vector3D(cv * cosRaan / sinI, dwXV,
  1242.                                           cv * sinRaan / sinI, dwYV);

  1244.         // dv
  1245.         final Vector3D dvP = new Vector3D(-u * cosRaan * cosI / sinI + sinRaan * z, dwXP,
  1246.                                           -u * sinRaan * cosI / sinI - cosRaan * z, dwYP,
  1247.                                           cv, dwZP,
  1248.                                           1, new Vector3D(-sinRaan * cosI, cosRaan * cosI, sinI));
  1249.         final Vector3D dvV = new Vector3D(-u * cosRaan * cosI / sinI + sinRaan * z, dwXV,
  1250.                                           -u * sinRaan * cosI / sinI - cosRaan * z, dwYV,
  1251.                                           cv, dwZV);

  1253.         final Vector3D dc1P = new Vector3D(aOr2 * (2 * eSinE * eSinE + 1 - eCosE) / r2, position,
  1254.                                             -2 * aOr2 * eSinE * oOsqrtMuA, velocity);
  1255.         final Vector3D dc1V = new Vector3D(-2 * aOr2 * eSinE * oOsqrtMuA, position,
  1256.                                             2 / mu, velocity);
  1257.         final Vector3D dc2P = new Vector3D(aOr2 * eSinE * (eSinE * eSinE - (1 - e2)) / (r2 * epsilon), position,
  1258.                                             aOr2 * (1 - e2 - eSinE * eSinE) * oOsqrtMuA / epsilon, velocity);
  1259.         final Vector3D dc2V = new Vector3D(aOr2 * (1 - e2 - eSinE * eSinE) * oOsqrtMuA / epsilon, position,
  1260.                                             eSinE / (mu * epsilon), velocity);

  1262.         final double cof1   = aOr2 * (eCosE - e2);
  1263.         final double cof2   = aOr2 * epsilon * eSinE;
  1264.         final Vector3D dexP = new Vector3D(u, dc1P,  v, dc2P, cof1, duP,  cof2, dvP);
  1265.         final Vector3D dexV = new Vector3D(u, dc1V,  v, dc2V, cof1, duV,  cof2, dvV);
  1266.         final Vector3D deyP = new Vector3D(v, dc1P, -u, dc2P, cof1, dvP, -cof2, duP);
  1267.         final Vector3D deyV = new Vector3D(v, dc1V, -u, dc2V, cof1, dvV, -cof2, duV);
  1268.         fillHalfRow(1, dexP, jacobian[1], 0);
  1269.         fillHalfRow(1, dexV, jacobian[1], 3);
  1270.         fillHalfRow(1, deyP, jacobian[2], 0);
  1271.         fillHalfRow(1, deyV, jacobian[2], 3);

  1273.         final double cle = u / a + ex - eSinE * beta * ey;
  1274.         final double sle = v / a + ey + eSinE * beta * ex;
  1275.         final double m1  = beta * eCosE;
  1276.         final double m2  = 1 - m1 * eCosE;
  1277.         final double m3  = (u * ey - v * ex) + eSinE * beta * (u * ex + v * ey);
  1278.         final double m4  = -sle + cle * eSinE * beta;
  1279.         final double m5  = cle + sle * eSinE * beta;
  1280.         fillHalfRow((2 * m3 / r + aOr * eSinE + m1 * eSinE * (1 + m1 - (1 + aOr) * m2) / epsilon) / r2, position,
  1281.                     (m1 * m2 / epsilon - 1) * oOsqrtMuA, velocity,
  1282.                     m4, dexP, m5, deyP, -sle / a, duP, cle / a, dvP,
  1283.                     jacobian[5], 0);
  1284.         fillHalfRow((m1 * m2 / epsilon - 1) * oOsqrtMuA, position,
  1285.                     (2 * m3 + eSinE * a + m1 * eSinE * r * (eCosE * beta * 2 - aOr * m2) / epsilon) / mu, velocity,
  1286.                     m4, dexV, m5, deyV, -sle / a, duV, cle / a, dvV,
  1287.                     jacobian[5], 3);

  1289.         return jacobian;

  1291.     }

  1293.     /** {@inheritDoc} */
  1294.     @Override
  1295.     protected double[][] computeJacobianEccentricWrtCartesian() {

  1297.         // start by computing the Jacobian with mean angle
  1298.         final double[][] jacobian = computeJacobianMeanWrtCartesian();

  1300.         // Differentiating the Kepler equation aM = aE - ex sin aE + ey cos aE leads to:
  1301.         // daM = (1 - ex cos aE - ey sin aE) daE - sin aE dex + cos aE dey
  1302.         // which is inverted and rewritten as:
  1303.         // daE = a/r daM + sin aE a/r dex - cos aE a/r dey
  1304.         final double alphaE = getAlphaE();
  1305.         final SinCos scAe   = FastMath.sinCos(alphaE);
  1306.         final double cosAe  = scAe.cos();
  1307.         final double sinAe  = scAe.sin();
  1308.         final double aOr    = 1 / (1 - ex * cosAe - ey * sinAe);

  1310.         // update longitude row
  1311.         final double[] rowEx = jacobian[1];
  1312.         final double[] rowEy = jacobian[2];
  1313.         final double[] rowL  = jacobian[5];
  1314.         for (int j = 0; j < 6; ++j) {
  1315.             rowL[j] = aOr * (rowL[j] + sinAe * rowEx[j] - cosAe * rowEy[j]);
  1316.         }

  1318.         return jacobian;

  1320.     }

  1322.     /** {@inheritDoc} */
  1323.     @Override
  1324.     protected double[][] computeJacobianTrueWrtCartesian() {

  1326.         // start by computing the Jacobian with eccentric angle
  1327.         final double[][] jacobian = computeJacobianEccentricWrtCartesian();

  1329.         // Differentiating the eccentric latitude equation
  1330.         // tan((aV - aE)/2) = [ex sin aE - ey cos aE] / [sqrt(1-ex^2-ey^2) + 1 - ex cos aE - ey sin aE]
  1331.         // leads to
  1332.         // cT (daV - daE) = cE daE + cX dex + cY dey
  1333.         // with
  1334.         // cT = [d^2 + (ex sin aE - ey cos aE)^2] / 2
  1335.         // d  = 1 + sqrt(1-ex^2-ey^2) - ex cos aE - ey sin aE
  1336.         // cE = (ex cos aE + ey sin aE) (sqrt(1-ex^2-ey^2) + 1) - ex^2 - ey^2
  1337.         // cX =  sin aE (sqrt(1-ex^2-ey^2) + 1) - ey + ex (ex sin aE - ey cos aE) / sqrt(1-ex^2-ey^2)
  1338.         // cY = -cos aE (sqrt(1-ex^2-ey^2) + 1) + ex + ey (ex sin aE - ey cos aE) / sqrt(1-ex^2-ey^2)
  1339.         // which can be solved to find the differential of the true latitude
  1340.         // daV = (cT + cE) / cT daE + cX / cT deX + cY / cT deX
  1341.         final double alphaE    = getAlphaE();
  1342.         final SinCos scAe      = FastMath.sinCos(alphaE);
  1343.         final double cosAe     = scAe.cos();
  1344.         final double sinAe     = scAe.sin();
  1345.         final double eSinE     = ex * sinAe - ey * cosAe;
  1346.         final double ecosE     = ex * cosAe + ey * sinAe;
  1347.         final double e2        = ex * ex + ey * ey;
  1348.         final double epsilon   = FastMath.sqrt(1 - e2);
  1349.         final double onePeps   = 1 + epsilon;
  1350.         final double d         = onePeps - ecosE;
  1351.         final double cT        = (d * d + eSinE * eSinE) / 2;
  1352.         final double cE        = ecosE * onePeps - e2;
  1353.         final double cX        = ex * eSinE / epsilon - ey + sinAe * onePeps;
  1354.         final double cY        = ey * eSinE / epsilon + ex - cosAe * onePeps;
  1355.         final double factorLe  = (cT + cE) / cT;
  1356.         final double factorEx  = cX / cT;
  1357.         final double factorEy  = cY / cT;

  1359.         // update latitude row
  1360.         final double[] rowEx = jacobian[1];
  1361.         final double[] rowEy = jacobian[2];
  1362.         final double[] rowA = jacobian[5];
  1363.         for (int j = 0; j < 6; ++j) {
  1364.             rowA[j] = factorLe * rowA[j] + factorEx * rowEx[j] + factorEy * rowEy[j];
  1365.         }

  1367.         return jacobian;

  1369.     }

  1371.     /** {@inheritDoc} */
  1372.     @Override
  1373.     public void addKeplerContribution(final PositionAngleType type, final double gm,
  1374.                                       final double[] pDot) {
  1375.         pDot[5] += computeKeplerianAlphaDot(type, a, ex, ey, gm, cachedAlpha, cachedPositionAngleType);
  1376.     }

  1378.     /**
  1379.      * Compute rate of argument of latitude.
  1380.      * @param type position angle type of rate
  1381.      * @param a semi major axis
  1382.      * @param ex ex
  1383.      * @param ey ey
  1384.      * @param mu mu
  1385.      * @param alpha argument of latitude
  1386.      * @param cachedType position angle type of passed alpha
  1387.      * @return first-order time derivative for alpha
  1388.      * @since 12.2
  1389.      */
  1390.     private static double computeKeplerianAlphaDot(final PositionAngleType type, final double a, final double ex,
  1391.                                                    final double ey, final double mu,
  1392.                                                    final double alpha, final PositionAngleType cachedType) {
  1393.         final double n  = FastMath.sqrt(mu / a) / a;
  1394.         if (type == PositionAngleType.MEAN) {
  1395.             return n;
  1396.         }
  1397.         final double ksi;
  1398.         final SinCos sc;
  1399.         if (type == PositionAngleType.ECCENTRIC) {
  1400.             sc = FastMath.sinCos(CircularLatitudeArgumentUtility.convertAlpha(cachedType, alpha, ex, ey, type));
  1401.             ksi   = 1. / (1 - ex * sc.cos() - ey * sc.sin());
  1402.             return n * ksi;
  1403.         } else { // TRUE
  1404.             sc = FastMath.sinCos(CircularLatitudeArgumentUtility.convertAlpha(cachedType, alpha, ex, ey, type));
  1405.             final double oMe2  = 1 - ex * ex - ey * ey;
  1406.             ksi   = 1 + ex * sc.cos() + ey * sc.sin();
  1407.             return n * ksi * ksi / (oMe2 * FastMath.sqrt(oMe2));
  1408.         }
  1409.     }

  1411.     /**  Returns a string representation of this Orbit object.
  1412.      * @return a string representation of this object
  1413.      */
  1414.     public String toString() {
  1415.         return new StringBuilder().append("circular parameters: ").append('{').
  1416.                                   append("a: ").append(a).
  1417.                                   append(", ex: ").append(ex).append(", ey: ").append(ey).
  1418.                                   append(", i: ").append(FastMath.toDegrees(i)).
  1419.                                   append(", raan: ").append(FastMath.toDegrees(raan)).
  1420.                                   append(", alphaV: ").append(FastMath.toDegrees(getAlphaV())).
  1421.                                   append(";}").toString();
  1422.     }

  1424.     /** {@inheritDoc} */
  1425.     @Override
  1426.     public PositionAngleType getCachedPositionAngleType() {
  1427.         return cachedPositionAngleType;
  1428.     }

  1430.     /** {@inheritDoc} */
  1431.     @Override
  1432.     public boolean hasRates() {
  1433.         return hasDerivatives();
  1434.     }

  1436.     /** {@inheritDoc} */
  1437.     @Override
  1438.     public CircularOrbit removeRates() {
  1439.         final PositionAngleType positionAngleType = getCachedPositionAngleType();
  1440.         return new CircularOrbit(a, ex, ey, i, raan, cachedAlpha, positionAngleType, positionAngleType,
  1441.                 getFrame(), getDate(), getMu());
  1442.     }

  1444.     /** Replace the instance with a data transfer object for serialization.
  1445.      * @return data transfer object that will be serialized
  1446.      */
  1447.     @DefaultDataContext
  1448.     private Object writeReplace() {
  1449.         return new DTO(this);
  1450.     }

  1452.     /** Internal class used only for serialization. */
  1453.     @DefaultDataContext
  1454.     private static class DTO implements Serializable {

  1456.         /** Serializable UID. */
  1457.         private static final long serialVersionUID = 20231217L;

  1459.         /** Double values. */
  1460.         private final double[] d;

  1462.         /** Frame in which are defined the orbital parameters. */
  1463.         private final Frame frame;

  1465.         /** Type of cached position angle. */
  1466.         private final PositionAngleType positionAngleType;

  1468.         /** Simple constructor.
  1469.          * @param orbit instance to serialize
  1470.          */
  1471.         private DTO(final CircularOrbit orbit) {

  1473.             final AbsoluteDate date = orbit.getDate();
  1474.             positionAngleType = orbit.getCachedPositionAngleType();

  1476.             // decompose date
  1477.             final AbsoluteDate j2000Epoch =
  1478.                     DataContext.getDefault().getTimeScales().getJ2000Epoch();
  1479.             final double epoch  = FastMath.floor(date.durationFrom(j2000Epoch));
  1480.             final double offset = date.durationFrom(j2000Epoch.shiftedBy(epoch));

  1482.             if (orbit.serializePV) {
  1483.                 final TimeStampedPVCoordinates pv = orbit.getPVCoordinates();
  1484.                 if (orbit.hasDerivatives()) {
  1485.                     this.d = new double[] {
  1486.                         // date + mu + orbit + derivatives + Cartesian : 24 parameters
  1487.                         epoch, offset, orbit.getMu(),
  1488.                         orbit.a, orbit.ex, orbit.ey,
  1489.                         orbit.i, orbit.raan, orbit.cachedAlpha,
  1490.                         orbit.aDot, orbit.exDot, orbit.eyDot,
  1491.                         orbit.iDot, orbit.raanDot, orbit.cachedAlphaDot,
  1492.                         pv.getPosition().getX(),     pv.getPosition().getY(),     pv.getPosition().getZ(),
  1493.                         pv.getVelocity().getX(),     pv.getVelocity().getY(),     pv.getVelocity().getZ(),
  1494.                         pv.getAcceleration().getX(), pv.getAcceleration().getY(), pv.getAcceleration().getZ(),
  1495.                     };
  1496.                 } else {
  1497.                     this.d = new double[] {
  1498.                         // date + mu + orbit + Cartesian : 18 parameters
  1499.                         epoch, offset, orbit.getMu(),
  1500.                         orbit.a, orbit.ex, orbit.ey,
  1501.                         orbit.i, orbit.raan, orbit.cachedAlpha,
  1502.                         pv.getPosition().getX(),     pv.getPosition().getY(),     pv.getPosition().getZ(),
  1503.                         pv.getVelocity().getX(),     pv.getVelocity().getY(),     pv.getVelocity().getZ(),
  1504.                         pv.getAcceleration().getX(), pv.getAcceleration().getY(), pv.getAcceleration().getZ(),
  1505.                     };
  1506.                 }
  1507.             } else {
  1508.                 if (orbit.hasDerivatives()) {
  1509.                     // date + mu + orbit + derivatives: 15 parameters
  1510.                     this.d = new double[] {
  1511.                         epoch, offset, orbit.getMu(),
  1512.                         orbit.a, orbit.ex, orbit.ey,
  1513.                         orbit.i, orbit.raan, orbit.cachedAlpha,
  1514.                         orbit.aDot, orbit.exDot, orbit.eyDot,
  1515.                         orbit.iDot, orbit.raanDot, orbit.cachedAlphaDot
  1516.                     };
  1517.                 } else {
  1518.                     // date + mu + orbit: 9 parameters
  1519.                     this.d = new double[] {
  1520.                         epoch, offset, orbit.getMu(),
  1521.                         orbit.a, orbit.ex, orbit.ey,
  1522.                         orbit.i, orbit.raan, orbit.cachedAlpha
  1523.                     };
  1524.                 }
  1525.             }

  1527.             this.frame = orbit.getFrame();

  1529.         }

  1531.         /** Replace the deserialized data transfer object with a {@link CircularOrbit}.
  1532.          * @return replacement {@link CircularOrbit}
  1533.          */
  1534.         private Object readResolve() {
  1535.             final AbsoluteDate j2000Epoch =
  1536.                     DataContext.getDefault().getTimeScales().getJ2000Epoch();
  1537.             switch (d.length) {
  1538.                 case 24 : // date + mu + orbit + derivatives + Cartesian
  1539.                     return new CircularOrbit(d[ 3], d[ 4], d[ 5], d[ 6], d[ 7], d[ 8],
  1540.                                              d[ 9], d[10], d[11], d[12], d[13], d[14],
  1541.                                              new TimeStampedPVCoordinates(j2000Epoch.shiftedBy(d[0]).shiftedBy(d[1]),
  1542.                                                                           new Vector3D(d[15], d[16], d[17]),
  1543.                                                                           new Vector3D(d[18], d[19], d[20]),
  1544.                                                                           new Vector3D(d[21], d[22], d[23])),
  1545.                                              positionAngleType, frame,
  1546.                                              d[2]);
  1547.                 case 18 : // date + mu + orbit + Cartesian
  1548.                     return new CircularOrbit(d[3], d[4], d[5], d[6], d[7], d[8],
  1549.                                              Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN,
  1550.                                              new TimeStampedPVCoordinates(j2000Epoch.shiftedBy(d[0]).shiftedBy(d[1]),
  1551.                                                                           new Vector3D(d[ 9], d[10], d[11]),
  1552.                                                                           new Vector3D(d[12], d[13], d[14]),
  1553.                                                                           new Vector3D(d[15], d[16], d[17])),
  1554.                                              positionAngleType, frame,
  1555.                                              d[2]);
  1556.                 case 15 : // date + mu + orbit + derivatives
  1557.                     return new CircularOrbit(d[ 3], d[ 4], d[ 5], d[ 6], d[ 7], d[ 8],
  1558.                                              d[ 9], d[10], d[11], d[12], d[13], d[14],
  1559.                                              positionAngleType, positionAngleType,
  1560.                                              frame, j2000Epoch.shiftedBy(d[0]).shiftedBy(d[1]),
  1561.                                              d[2]);
  1562.                 default : // date + mu + orbit
  1563.                     return new CircularOrbit(d[3], d[4], d[5], d[6], d[7], d[8], positionAngleType, positionAngleType,
  1564.                                              frame, j2000Epoch.shiftedBy(d[0]).shiftedBy(d[1]),
  1565.                                              d[2]);

  1567.             }
  1568.         }

  1570.     }

  1572. }
Created with JaCoCo 0.8.12.202403310830