CircularOrbit.java

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

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


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

  37.  * <p>
  38.  * The parameters used internally are the circular elements which can be
  39.  * related to Keplerian elements as follows:
  40.  *   <ul>
  41.  *     <li>a</li>
  42.  *     <li>e<sub>x</sub> = e cos(ω)</li>
  43.  *     <li>e<sub>y</sub> = e sin(ω)</li>
  44.  *     <li>i</li>
  45.  *     <li>Ω</li>
  46.  *     <li>α<sub>v</sub> = v + ω</li>
  47.  *   </ul>
  48.  * where Ω stands for the Right Ascension of the Ascending Node and
  49.  * α<sub>v</sub> stands for the true latitude argument
  50.  *
  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.  */

  73. public class CircularOrbit extends Orbit implements PositionAngleBased<CircularOrbit> {

  75.     /** Semi-major axis (m). */
  76.     private final double a;

  78.     /** First component of the circular eccentricity vector. */
  79.     private final double ex;

  81.     /** Second component of the circular eccentricity vector. */
  82.     private final double ey;

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

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

  90.     /** Cached latitude argument (rad). */
  91.     private final double cachedAlpha;

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

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

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

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

  105.     /** Inclination derivative (rad/s). */
  106.     private final double iDot;

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

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

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

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

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

  167.     /** Creates a new instance.
  168.      * @param a  semi-major axis (m)
  169.      * @param ex e cos(ω), first component of circular eccentricity vector
  170.      * @param ey e sin(ω), second component of circular eccentricity vector
  171.      * @param i inclination (rad)
  172.      * @param raan right ascension of ascending node (Ω, rad)
  173.      * @param alpha  an + ω, mean, eccentric or true latitude argument (rad)
  174.      * @param aDot  semi-major axis derivative (m/s)
  175.      * @param exDot d(e cos(ω))/dt, first component of circular eccentricity vector derivative
  176.      * @param eyDot d(e sin(ω))/dt, second component of circular eccentricity vector derivative
  177.      * @param iDot inclination  derivative(rad/s)
  178.      * @param raanDot right ascension of ascending node derivative (rad/s)
  179.      * @param alphaDot  d(an + ω), mean, eccentric or true latitude argument derivative (rad/s)
  180.      * @param type type of latitude argument
  181.      * @param cachedPositionAngleType type of cached 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.      * @since 12.1
  189.      */
  190.     public CircularOrbit(final double a, final double ex, final double ey,
  191.                          final double i, final double raan, final double alpha,
  192.                          final double aDot, final double exDot, final double eyDot,
  193.                          final double iDot, final double raanDot, final double alphaDot,
  194.                          final PositionAngleType type, final PositionAngleType cachedPositionAngleType,
  195.                          final Frame frame, final AbsoluteDate date, final double mu)
  196.         throws IllegalArgumentException {
  197.         super(frame, date, mu);
  198.         if (ex * ex + ey * ey >= 1.0) {
  199.             throw new OrekitIllegalArgumentException(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS,
  200.                                                      getClass().getName());
  201.         }
  202.         this.a       =  a;
  203.         this.aDot    =  aDot;
  204.         this.ex      = ex;
  205.         this.exDot   = exDot;
  206.         this.ey      = ey;
  207.         this.eyDot   = eyDot;
  208.         this.i       = i;
  209.         this.iDot    = iDot;
  210.         this.raan    = raan;
  211.         this.raanDot = raanDot;
  212.         this.cachedPositionAngleType = cachedPositionAngleType;

  214.         final UnivariateDerivative1 alphaUD = initializeCachedAlpha(alpha, alphaDot, type);
  215.         this.cachedAlpha = alphaUD.getValue();
  216.         this.cachedAlphaDot = alphaUD.getFirstDerivative();

  218.         partialPV   = null;

  220.     }

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

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

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

  318.         // compute semi-major axis
  319.         final Vector3D pvP = pvCoordinates.getPosition();
  320.         final Vector3D pvV = pvCoordinates.getVelocity();
  321.         final Vector3D pvA = pvCoordinates.getAcceleration();
  322.         final double r2 = pvP.getNormSq();
  323.         final double r  = FastMath.sqrt(r2);
  324.         final double V2 = pvV.getNormSq();
  325.         final double rV2OnMu = r * V2 / mu;
  326.         a = r / (2 - rV2OnMu);

  328.         if (!isElliptical()) {
  329.             throw new OrekitIllegalArgumentException(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS,
  330.                                                      getClass().getName());
  331.         }

  333.         // compute inclination
  334.         final Vector3D momentum = pvCoordinates.getMomentum();
  335.         i = Vector3D.angle(momentum, Vector3D.PLUS_K);

  337.         // compute right ascension of ascending node
  338.         final Vector3D node  = Vector3D.crossProduct(Vector3D.PLUS_K, momentum);
  339.         raan = FastMath.atan2(node.getY(), node.getX());

  341.         // 2D-coordinates in the canonical frame
  342.         final SinCos scRaan = FastMath.sinCos(raan);
  343.         final SinCos scI    = FastMath.sinCos(i);
  344.         final double xP     = pvP.getX();
  345.         final double yP     = pvP.getY();
  346.         final double zP     = pvP.getZ();
  347.         final double x2     = (xP * scRaan.cos() + yP * scRaan.sin()) / a;
  348.         final double y2     = ((yP * scRaan.cos() - xP * scRaan.sin()) * scI.cos() + zP * scI.sin()) / a;

  350.         // compute eccentricity vector
  351.         final double eSE    = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(mu * a);
  352.         final double eCE    = rV2OnMu - 1;
  353.         final double e2     = eCE * eCE + eSE * eSE;
  354.         final double f      = eCE - e2;
  355.         final double g      = FastMath.sqrt(1 - e2) * eSE;
  356.         final double aOnR   = a / r;
  357.         final double a2OnR2 = aOnR * aOnR;
  358.         ex = a2OnR2 * (f * x2 + g * y2);
  359.         ey = a2OnR2 * (f * y2 - g * x2);

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

  365.         partialPV   = pvCoordinates;

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

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

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

  384.             // in order to compute latitude argument derivative, we must compute
  385.             // mean latitude argument derivative including Keplerian motion and convert to true latitude argument
  386.             final double alphaMDot = getKeplerianMeanMotion() +
  387.                                      jacobian[5][3] * aX + jacobian[5][4] * aY + jacobian[5][5] * aZ;
  388.             final UnivariateDerivative1 exUD     = new UnivariateDerivative1(ex, exDot);
  389.             final UnivariateDerivative1 eyUD     = new UnivariateDerivative1(ey, eyDot);
  390.             final UnivariateDerivative1 alphaMUD = new UnivariateDerivative1(getAlphaM(), alphaMDot);
  391.             final UnivariateDerivative1 alphavUD = FieldCircularLatitudeArgumentUtility.meanToTrue(exUD, eyUD, alphaMUD);
  392.             cachedAlphaDot = alphavUD.getFirstDerivative();

  394.         } else {
  395.             // acceleration is either almost zero or NaN,
  396.             // we assume acceleration was not known
  397.             // we don't set up derivatives
  398.             aDot      = 0.;
  399.             exDot     = 0.;
  400.             eyDot     = 0.;
  401.             iDot      = 0.;
  402.             raanDot   = 0.;
  403.             cachedAlphaDot = computeKeplerianAlphaDot(cachedPositionAngleType, a, ex, ey, mu, cachedAlpha, cachedPositionAngleType);
  404.         }

  406.     }

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

  429.     /** Constructor from any kind of orbital parameters.
  430.      * @param op orbital parameters to copy
  431.      */
  432.     public CircularOrbit(final Orbit op) {

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

  436.         a    = op.getA();
  437.         i    = op.getI();
  438.         final double hx = op.getHx();
  439.         final double hy = op.getHy();
  440.         final double h2 = hx * hx + hy * hy;
  441.         final double h  = FastMath.sqrt(h2);
  442.         raan = FastMath.atan2(hy, hx);
  443.         final SinCos scRaan  = FastMath.sinCos(raan);
  444.         final double cosRaan = h == 0 ? scRaan.cos() : hx / h;
  445.         final double sinRaan = h == 0 ? scRaan.sin() : hy / h;
  446.         final double equiEx = op.getEquinoctialEx();
  447.         final double equiEy = op.getEquinoctialEy();
  448.         ex     = equiEx * cosRaan + equiEy * sinRaan;
  449.         ey     = equiEy * cosRaan - equiEx * sinRaan;
  450.         cachedPositionAngleType = PositionAngleType.TRUE;
  451.         cachedAlpha = op.getLv() - raan;

  453.         if (op.hasNonKeplerianAcceleration()) {
  454.             aDot    = op.getADot();
  455.             final double hxDot = op.getHxDot();
  456.             final double hyDot = op.getHyDot();
  457.             iDot    = 2 * (cosRaan * hxDot + sinRaan * hyDot) / (1 + h2);
  458.             raanDot = (hx * hyDot - hy * hxDot) / h2;
  459.             final double equiExDot = op.getEquinoctialExDot();
  460.             final double equiEyDot = op.getEquinoctialEyDot();
  461.             exDot   = (equiExDot + equiEy * raanDot) * cosRaan +
  462.                       (equiEyDot - equiEx * raanDot) * sinRaan;
  463.             eyDot   = (equiEyDot - equiEx * raanDot) * cosRaan -
  464.                       (equiExDot + equiEy * raanDot) * sinRaan;
  465.             cachedAlphaDot = op.getLvDot() - raanDot;
  466.         } else {
  467.             aDot      = 0.;
  468.             exDot     = 0.;
  469.             eyDot     = 0.;
  470.             iDot      = 0.;
  471.             raanDot   = 0.;
  472.             cachedAlphaDot = computeKeplerianAlphaDot(cachedPositionAngleType, a, ex, ey, getMu(), cachedAlpha, cachedPositionAngleType);
  473.         }

  475.         partialPV   = null;

  477.     }

  479.     /** {@inheritDoc} */
  480.     @Override
  481.     public boolean hasNonKeplerianAcceleration() {
  482.         return aDot != 0. || exDot != 0. || eyDot != 0. || iDot != 0. || raanDot != 0. ||
  483.                 FastMath.abs(cachedAlphaDot - computeKeplerianAlphaDot(cachedPositionAngleType, a, ex, ey, getMu(), cachedAlpha, cachedPositionAngleType)) > TOLERANCE_POSITION_ANGLE_RATE;
  484.     }

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

  492.     /** {@inheritDoc} */
  493.     @Override
  494.     public double getA() {
  495.         return a;
  496.     }

  498.     /** {@inheritDoc} */
  499.     @Override
  500.     public double getADot() {
  501.         return aDot;
  502.     }

  504.     /** {@inheritDoc} */
  505.     @Override
  506.     public double getEquinoctialEx() {
  507.         final SinCos sc = FastMath.sinCos(raan);
  508.         return ex * sc.cos() - ey * sc.sin();
  509.     }

  511.     /** {@inheritDoc} */
  512.     @Override
  513.     public double getEquinoctialExDot() {
  514.         if (!hasNonKeplerianAcceleration()) {
  515.             return 0.;
  516.         }
  517.         final SinCos sc = FastMath.sinCos(raan);
  518.         return (exDot - ey * raanDot) * sc.cos() - (eyDot + ex * raanDot) * sc.sin();
  519.     }

  521.     /** {@inheritDoc} */
  522.     @Override
  523.     public double getEquinoctialEy() {
  524.         final SinCos sc = FastMath.sinCos(raan);
  525.         return ey * sc.cos() + ex * sc.sin();
  526.     }

  528.     /** {@inheritDoc} */
  529.     @Override
  530.     public double getEquinoctialEyDot() {
  531.         if (!hasNonKeplerianAcceleration()) {
  532.             return 0.;
  533.         }
  534.         final SinCos sc = FastMath.sinCos(raan);
  535.         return (eyDot + ex * raanDot) * sc.cos() + (exDot - ey * raanDot) * sc.sin();
  536.     }

  538.     /** Get the first component of the circular eccentricity vector.
  539.      * @return ex = e cos(ω), first component of the circular eccentricity vector
  540.      */
  541.     public double getCircularEx() {
  542.         return ex;
  543.     }

  545.     /** Get the first component of the circular eccentricity vector derivative.
  546.      * @return ex = e cos(ω), first component of the circular eccentricity vector derivative
  547.      * @since 9.0
  548.      */
  549.     public double getCircularExDot() {
  550.         return exDot;
  551.     }

  553.     /** Get the second component of the circular eccentricity vector.
  554.      * @return ey = e sin(ω), second component of the circular eccentricity vector
  555.      */
  556.     public double getCircularEy() {
  557.         return ey;
  558.     }

  560.     /** Get the second component of the circular eccentricity vector derivative.
  561.      * @return ey = e sin(ω), second component of the circular eccentricity vector derivative
  562.      */
  563.     public double getCircularEyDot() {
  564.         return eyDot;
  565.     }

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

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

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

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

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

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

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

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

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

  653.             case TRUE:
  654.                 return cachedAlphaDot;

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

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

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

  677.             case ECCENTRIC:
  678.                 return cachedAlpha;

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

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

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

  705.             case ECCENTRIC:
  706.                 return cachedAlphaDot;

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

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

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

  729.             case MEAN:
  730.                 return cachedAlpha;

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

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

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

  757.             case MEAN:
  758.                 return cachedAlphaDot;

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

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

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

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

  797.     /** {@inheritDoc} */
  798.     @Override
  799.     public double getE() {
  800.         return FastMath.sqrt(ex * ex + ey * ey);
  801.     }

  803.     /** {@inheritDoc} */
  804.     @Override
  805.     public double getEDot() {
  806.         if (!hasNonKeplerianAcceleration()) {
  807.             return 0.;
  808.         }
  809.         return (ex * exDot + ey * eyDot) / getE();
  810.     }

  812.     /** {@inheritDoc} */
  813.     @Override
  814.     public double getI() {
  815.         return i;
  816.     }

  818.     /** {@inheritDoc} */
  819.     @Override
  820.     public double getIDot() {
  821.         return iDot;
  822.     }

  824.     /** Get the right ascension of the ascending node.
  825.      * @return right ascension of the ascending node (rad)
  826.      */
  827.     public double getRightAscensionOfAscendingNode() {
  828.         return raan;
  829.     }

  831.     /** Get the right ascension of the ascending node derivative.
  832.      * <p>
  833.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  834.      * </p>
  835.      * @return right ascension of the ascending node derivative (rad/s)
  836.      * @since 9.0
  837.      */
  838.     public double getRightAscensionOfAscendingNodeDot() {
  839.         return raanDot;
  840.     }

  842.     /** {@inheritDoc} */
  843.     @Override
  844.     public double getLv() {
  845.         return getAlphaV() + raan;
  846.     }

  848.     /** {@inheritDoc} */
  849.     @Override
  850.     public double getLvDot() {
  851.         return getAlphaVDot() + raanDot;
  852.     }

  854.     /** {@inheritDoc} */
  855.     @Override
  856.     public double getLE() {
  857.         return getAlphaE() + raan;
  858.     }

  860.     /** {@inheritDoc} */
  861.     @Override
  862.     public double getLEDot() {
  863.         return getAlphaEDot() + raanDot;
  864.     }

  866.     /** {@inheritDoc} */
  867.     @Override
  868.     public double getLM() {
  869.         return getAlphaM() + raan;
  870.     }

  872.     /** {@inheritDoc} */
  873.     @Override
  874.     public double getLMDot() {
  875.         return getAlphaMDot() + raanDot;
  876.     }

  878.     /** Compute position and velocity but not acceleration.
  879.      */
  880.     private void computePVWithoutA() {

  882.         if (partialPV != null) {
  883.             // already computed
  884.             return;
  885.         }

  887.         // get equinoctial parameters
  888.         final double equEx = getEquinoctialEx();
  889.         final double equEy = getEquinoctialEy();
  890.         final double hx = getHx();
  891.         final double hy = getHy();
  892.         final double lE = getLE();

  894.         // inclination-related intermediate parameters
  895.         final double hx2   = hx * hx;
  896.         final double hy2   = hy * hy;
  897.         final double factH = 1. / (1 + hx2 + hy2);

  899.         // reference axes defining the orbital plane
  900.         final double ux = (1 + hx2 - hy2) * factH;
  901.         final double uy =  2 * hx * hy * factH;
  902.         final double uz = -2 * hy * factH;

  904.         final double vx = uy;
  905.         final double vy = (1 - hx2 + hy2) * factH;
  906.         final double vz =  2 * hx * factH;

  908.         // eccentricity-related intermediate parameters
  909.         final double exey = equEx * equEy;
  910.         final double ex2  = equEx * equEx;
  911.         final double ey2  = equEy * equEy;
  912.         final double e2   = ex2 + ey2;
  913.         final double eta  = 1 + FastMath.sqrt(1 - e2);
  914.         final double beta = 1. / eta;

  916.         // eccentric latitude argument
  917.         final SinCos scLe   = FastMath.sinCos(lE);
  918.         final double cLe    = scLe.cos();
  919.         final double sLe    = scLe.sin();
  920.         final double exCeyS = equEx * cLe + equEy * sLe;

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

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

  930.         final Vector3D position =
  931.                         new Vector3D(x * ux + y * vx, x * uy + y * vy, x * uz + y * vz);
  932.         final Vector3D velocity =
  933.                         new Vector3D(xdot * ux + ydot * vx, xdot * uy + ydot * vy, xdot * uz + ydot * vz);

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

  937.     }

  939.     /** Initialize cached alpha with rate.
  940.      * @param alpha input alpha
  941.      * @param alphaDot rate of input alpha
  942.      * @param inputType position angle type passed as input
  943.      * @return alpha to cache with rate
  944.      * @since 12.1
  945.      */
  946.     private UnivariateDerivative1 initializeCachedAlpha(final double alpha, final double alphaDot,
  947.                                                         final PositionAngleType inputType) {
  948.         if (cachedPositionAngleType == inputType) {
  949.             return new UnivariateDerivative1(alpha, alphaDot);

  951.         } else {
  952.             final UnivariateDerivative1 exUD = new UnivariateDerivative1(ex, exDot);
  953.             final UnivariateDerivative1 eyUD = new UnivariateDerivative1(ey, eyDot);
  954.             final UnivariateDerivative1 alphaUD = new UnivariateDerivative1(alpha, alphaDot);

  956.             switch (cachedPositionAngleType) {

  958.                 case ECCENTRIC:
  959.                     if (inputType == PositionAngleType.MEAN) {
  960.                         return FieldCircularLatitudeArgumentUtility.meanToEccentric(exUD, eyUD, alphaUD);
  961.                     } else {
  962.                         return FieldCircularLatitudeArgumentUtility.trueToEccentric(exUD, eyUD, alphaUD);
  963.                     }

  965.                 case TRUE:
  966.                     if (inputType == PositionAngleType.MEAN) {
  967.                         return FieldCircularLatitudeArgumentUtility.meanToTrue(exUD, eyUD, alphaUD);
  968.                     } else {
  969.                         return FieldCircularLatitudeArgumentUtility.eccentricToTrue(exUD, eyUD, alphaUD);
  970.                     }

  972.                 case MEAN:
  973.                     if (inputType == PositionAngleType.TRUE) {
  974.                         return FieldCircularLatitudeArgumentUtility.trueToMean(exUD, eyUD, alphaUD);
  975.                     } else {
  976.                         return FieldCircularLatitudeArgumentUtility.eccentricToMean(exUD, eyUD, alphaUD);
  977.                     }

  979.                 default:
  980.                     throw new OrekitInternalError(null);

  982.             }

  984.         }

  986.     }

  988.     /** {@inheritDoc} */
  989.     @Override
  990.     protected Vector3D initPosition() {

  992.         // get equinoctial parameters
  993.         final double equEx = getEquinoctialEx();
  994.         final double equEy = getEquinoctialEy();
  995.         final double hx = getHx();
  996.         final double hy = getHy();
  997.         final double lE = getLE();

  999.         // inclination-related intermediate parameters
  1000.         final double hx2   = hx * hx;
  1001.         final double hy2   = hy * hy;
  1002.         final double factH = 1. / (1 + hx2 + hy2);

  1004.         // reference axes defining the orbital plane
  1005.         final double ux = (1 + hx2 - hy2) * factH;
  1006.         final double uy =  2 * hx * hy * factH;
  1007.         final double uz = -2 * hy * factH;

  1009.         final double vx = uy;
  1010.         final double vy = (1 - hx2 + hy2) * factH;
  1011.         final double vz =  2 * hx * factH;

  1013.         // eccentricity-related intermediate parameters
  1014.         final double exey = equEx * equEy;
  1015.         final double ex2  = equEx * equEx;
  1016.         final double ey2  = equEy * equEy;
  1017.         final double e2   = ex2 + ey2;
  1018.         final double eta  = 1 + FastMath.sqrt(1 - e2);
  1019.         final double beta = 1. / eta;

  1021.         // eccentric latitude argument
  1022.         final SinCos scLe   = FastMath.sinCos(lE);
  1023.         final double cLe    = scLe.cos();
  1024.         final double sLe    = scLe.sin();

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

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

  1032.     }

  1034.     /** {@inheritDoc} */
  1035.     @Override
  1036.     protected TimeStampedPVCoordinates initPVCoordinates() {

  1038.         // position and velocity
  1039.         computePVWithoutA();

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

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

  1050.     }

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

  1070.     /** {@inheritDoc} */
  1071.     @Override
  1072.     public CircularOrbit withCachedPositionAngleType(final PositionAngleType positionAngleType) {
  1073.         return new CircularOrbit(a, ex, ey, i, raan, getAlpha(positionAngleType), aDot, exDot, eyDot, iDot, raanDot,
  1074.                 getAlphaDot(positionAngleType), positionAngleType, getFrame(), getDate(), getMu());
  1075.     }

  1077.     /** {@inheritDoc} */
  1078.     @Override
  1079.     public CircularOrbit shiftedBy(final double dt) {
  1080.         return shiftedBy(new TimeOffset(dt));
  1081.     }

  1083.     /** {@inheritDoc} */
  1084.     @Override
  1085.     public CircularOrbit shiftedBy(final TimeOffset dt) {

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

  1089.         // use Keplerian-only motion
  1090.         final CircularOrbit keplerianShifted = new CircularOrbit(a, ex, ey, i, raan,
  1091.                                                                  getAlphaM() + getKeplerianMeanMotion() * dtS,
  1092.                                                                  PositionAngleType.MEAN, cachedPositionAngleType,
  1093.                                                                  getFrame(), getDate().shiftedBy(dt), getMu());

  1095.         if (dtS != 0. && hasNonKeplerianRates()) {

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

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

  1111.             // build a new orbit, taking non-Keplerian acceleration into account
  1112.             return new CircularOrbit(new TimeStampedPVCoordinates(keplerianShifted.getDate(),
  1113.                                                                   fixedP, fixedV, fixedA),
  1114.                                      keplerianShifted.getFrame(), keplerianShifted.getMu());

  1116.         } else {
  1117.             // Keplerian-only motion is all we can do
  1118.             return keplerianShifted;
  1119.         }

  1121.     }

  1123.     /** {@inheritDoc} */
  1124.     @Override
  1125.     protected double[][] computeJacobianMeanWrtCartesian() {


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

  1130.         computePVWithoutA();
  1131.         final Vector3D position = partialPV.getPosition();
  1132.         final Vector3D velocity = partialPV.getVelocity();
  1133.         final double x          = position.getX();
  1134.         final double y          = position.getY();
  1135.         final double z          = position.getZ();
  1136.         final double vx         = velocity.getX();
  1137.         final double vy         = velocity.getY();
  1138.         final double vz         = velocity.getZ();
  1139.         final double pv         = Vector3D.dotProduct(position, velocity);
  1140.         final double r2         = position.getNormSq();
  1141.         final double r          = FastMath.sqrt(r2);
  1142.         final double v2         = velocity.getNormSq();

  1144.         final double mu         = getMu();
  1145.         final double oOsqrtMuA  = 1 / FastMath.sqrt(mu * a);
  1146.         final double rOa        = r / a;
  1147.         final double aOr        = a / r;
  1148.         final double aOr2       = a / r2;
  1149.         final double a2         = a * a;

  1151.         final double ex2        = ex * ex;
  1152.         final double ey2        = ey * ey;
  1153.         final double e2         = ex2 + ey2;
  1154.         final double epsilon    = FastMath.sqrt(1 - e2);
  1155.         final double beta       = 1 / (1 + epsilon);

  1157.         final double eCosE      = 1 - rOa;
  1158.         final double eSinE      = pv * oOsqrtMuA;

  1160.         final SinCos scI    = FastMath.sinCos(i);
  1161.         final SinCos scRaan = FastMath.sinCos(raan);
  1162.         final double cosI       = scI.cos();
  1163.         final double sinI       = scI.sin();
  1164.         final double cosRaan    = scRaan.cos();
  1165.         final double sinRaan    = scRaan.sin();

  1167.         // da
  1168.         fillHalfRow(2 * aOr * aOr2, position, jacobian[0], 0);
  1169.         fillHalfRow(2 * a2 / mu, velocity, jacobian[0], 3);

  1171.         // differentials of the normalized momentum
  1172.         final Vector3D danP = new Vector3D(v2, position, -pv, velocity);
  1173.         final Vector3D danV = new Vector3D(r2, velocity, -pv, position);
  1174.         final double recip  = 1 / partialPV.getMomentum().getNorm();
  1175.         final double recip2 = recip * recip;
  1176.         final Vector3D dwXP = new Vector3D(recip, new Vector3D(  0,  vz, -vy), -recip2 * sinRaan * sinI, danP);
  1177.         final Vector3D dwYP = new Vector3D(recip, new Vector3D(-vz,   0,  vx),  recip2 * cosRaan * sinI, danP);
  1178.         final Vector3D dwZP = new Vector3D(recip, new Vector3D( vy, -vx,   0), -recip2 * cosI,           danP);
  1179.         final Vector3D dwXV = new Vector3D(recip, new Vector3D(  0,  -z,   y), -recip2 * sinRaan * sinI, danV);
  1180.         final Vector3D dwYV = new Vector3D(recip, new Vector3D(  z,   0,  -x),  recip2 * cosRaan * sinI, danV);
  1181.         final Vector3D dwZV = new Vector3D(recip, new Vector3D( -y,   x,   0), -recip2 * cosI,           danV);

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

  1187.         // dRaan
  1188.         fillHalfRow(sinRaan / sinI, dwYP, cosRaan / sinI, dwXP, jacobian[4], 0);
  1189.         fillHalfRow(sinRaan / sinI, dwYV, cosRaan / sinI, dwXV, jacobian[4], 3);

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

  1197.         // du
  1198.         final Vector3D duP = new Vector3D(cv * cosRaan / sinI, dwXP,
  1199.                                           cv * sinRaan / sinI, dwYP,
  1200.                                           1, new Vector3D(cosRaan, sinRaan, 0));
  1201.         final Vector3D duV = new Vector3D(cv * cosRaan / sinI, dwXV,
  1202.                                           cv * sinRaan / sinI, dwYV);

  1204.         // dv
  1205.         final Vector3D dvP = new Vector3D(-u * cosRaan * cosI / sinI + sinRaan * z, dwXP,
  1206.                                           -u * sinRaan * cosI / sinI - cosRaan * z, dwYP,
  1207.                                           cv, dwZP,
  1208.                                           1, new Vector3D(-sinRaan * cosI, cosRaan * cosI, sinI));
  1209.         final Vector3D dvV = new Vector3D(-u * cosRaan * cosI / sinI + sinRaan * z, dwXV,
  1210.                                           -u * sinRaan * cosI / sinI - cosRaan * z, dwYV,
  1211.                                           cv, dwZV);

  1213.         final Vector3D dc1P = new Vector3D(aOr2 * (2 * eSinE * eSinE + 1 - eCosE) / r2, position,
  1214.                                             -2 * aOr2 * eSinE * oOsqrtMuA, velocity);
  1215.         final Vector3D dc1V = new Vector3D(-2 * aOr2 * eSinE * oOsqrtMuA, position,
  1216.                                             2 / mu, velocity);
  1217.         final Vector3D dc2P = new Vector3D(aOr2 * eSinE * (eSinE * eSinE - (1 - e2)) / (r2 * epsilon), position,
  1218.                                             aOr2 * (1 - e2 - eSinE * eSinE) * oOsqrtMuA / epsilon, velocity);
  1219.         final Vector3D dc2V = new Vector3D(aOr2 * (1 - e2 - eSinE * eSinE) * oOsqrtMuA / epsilon, position,
  1220.                                             eSinE / (mu * epsilon), velocity);

  1222.         final double cof1   = aOr2 * (eCosE - e2);
  1223.         final double cof2   = aOr2 * epsilon * eSinE;
  1224.         final Vector3D dexP = new Vector3D(u, dc1P,  v, dc2P, cof1, duP,  cof2, dvP);
  1225.         final Vector3D dexV = new Vector3D(u, dc1V,  v, dc2V, cof1, duV,  cof2, dvV);
  1226.         final Vector3D deyP = new Vector3D(v, dc1P, -u, dc2P, cof1, dvP, -cof2, duP);
  1227.         final Vector3D deyV = new Vector3D(v, dc1V, -u, dc2V, cof1, dvV, -cof2, duV);
  1228.         fillHalfRow(1, dexP, jacobian[1], 0);
  1229.         fillHalfRow(1, dexV, jacobian[1], 3);
  1230.         fillHalfRow(1, deyP, jacobian[2], 0);
  1231.         fillHalfRow(1, deyV, jacobian[2], 3);

  1233.         final double cle = u / a + ex - eSinE * beta * ey;
  1234.         final double sle = v / a + ey + eSinE * beta * ex;
  1235.         final double m1  = beta * eCosE;
  1236.         final double m2  = 1 - m1 * eCosE;
  1237.         final double m3  = (u * ey - v * ex) + eSinE * beta * (u * ex + v * ey);
  1238.         final double m4  = -sle + cle * eSinE * beta;
  1239.         final double m5  = cle + sle * eSinE * beta;
  1240.         fillHalfRow((2 * m3 / r + aOr * eSinE + m1 * eSinE * (1 + m1 - (1 + aOr) * m2) / epsilon) / r2, position,
  1241.                     (m1 * m2 / epsilon - 1) * oOsqrtMuA, velocity,
  1242.                     m4, dexP, m5, deyP, -sle / a, duP, cle / a, dvP,
  1243.                     jacobian[5], 0);
  1244.         fillHalfRow((m1 * m2 / epsilon - 1) * oOsqrtMuA, position,
  1245.                     (2 * m3 + eSinE * a + m1 * eSinE * r * (eCosE * beta * 2 - aOr * m2) / epsilon) / mu, velocity,
  1246.                     m4, dexV, m5, deyV, -sle / a, duV, cle / a, dvV,
  1247.                     jacobian[5], 3);

  1249.         return jacobian;

  1251.     }

  1253.     /** {@inheritDoc} */
  1254.     @Override
  1255.     protected double[][] computeJacobianEccentricWrtCartesian() {

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

  1260.         // Differentiating the Kepler equation aM = aE - ex sin aE + ey cos aE leads to:
  1261.         // daM = (1 - ex cos aE - ey sin aE) daE - sin aE dex + cos aE dey
  1262.         // which is inverted and rewritten as:
  1263.         // daE = a/r daM + sin aE a/r dex - cos aE a/r dey
  1264.         final double alphaE = getAlphaE();
  1265.         final SinCos scAe   = FastMath.sinCos(alphaE);
  1266.         final double cosAe  = scAe.cos();
  1267.         final double sinAe  = scAe.sin();
  1268.         final double aOr    = 1 / (1 - ex * cosAe - ey * sinAe);

  1270.         // update longitude row
  1271.         final double[] rowEx = jacobian[1];
  1272.         final double[] rowEy = jacobian[2];
  1273.         final double[] rowL  = jacobian[5];
  1274.         for (int j = 0; j < 6; ++j) {
  1275.             rowL[j] = aOr * (rowL[j] + sinAe * rowEx[j] - cosAe * rowEy[j]);
  1276.         }

  1278.         return jacobian;

  1280.     }

  1282.     /** {@inheritDoc} */
  1283.     @Override
  1284.     protected double[][] computeJacobianTrueWrtCartesian() {

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

  1289.         // Differentiating the eccentric latitude equation
  1290.         // tan((aV - aE)/2) = [ex sin aE - ey cos aE] / [sqrt(1-ex^2-ey^2) + 1 - ex cos aE - ey sin aE]
  1291.         // leads to
  1292.         // cT (daV - daE) = cE daE + cX dex + cY dey
  1293.         // with
  1294.         // cT = [d^2 + (ex sin aE - ey cos aE)^2] / 2
  1295.         // d  = 1 + sqrt(1-ex^2-ey^2) - ex cos aE - ey sin aE
  1296.         // cE = (ex cos aE + ey sin aE) (sqrt(1-ex^2-ey^2) + 1) - ex^2 - ey^2
  1297.         // cX =  sin aE (sqrt(1-ex^2-ey^2) + 1) - ey + ex (ex sin aE - ey cos aE) / sqrt(1-ex^2-ey^2)
  1298.         // cY = -cos aE (sqrt(1-ex^2-ey^2) + 1) + ex + ey (ex sin aE - ey cos aE) / sqrt(1-ex^2-ey^2)
  1299.         // which can be solved to find the differential of the true latitude
  1300.         // daV = (cT + cE) / cT daE + cX / cT deX + cY / cT deX
  1301.         final double alphaE    = getAlphaE();
  1302.         final SinCos scAe      = FastMath.sinCos(alphaE);
  1303.         final double cosAe     = scAe.cos();
  1304.         final double sinAe     = scAe.sin();
  1305.         final double eSinE     = ex * sinAe - ey * cosAe;
  1306.         final double ecosE     = ex * cosAe + ey * sinAe;
  1307.         final double e2        = ex * ex + ey * ey;
  1308.         final double epsilon   = FastMath.sqrt(1 - e2);
  1309.         final double onePeps   = 1 + epsilon;
  1310.         final double d         = onePeps - ecosE;
  1311.         final double cT        = (d * d + eSinE * eSinE) / 2;
  1312.         final double cE        = ecosE * onePeps - e2;
  1313.         final double cX        = ex * eSinE / epsilon - ey + sinAe * onePeps;
  1314.         final double cY        = ey * eSinE / epsilon + ex - cosAe * onePeps;
  1315.         final double factorLe  = (cT + cE) / cT;
  1316.         final double factorEx  = cX / cT;
  1317.         final double factorEy  = cY / cT;

  1319.         // update latitude row
  1320.         final double[] rowEx = jacobian[1];
  1321.         final double[] rowEy = jacobian[2];
  1322.         final double[] rowA = jacobian[5];
  1323.         for (int j = 0; j < 6; ++j) {
  1324.             rowA[j] = factorLe * rowA[j] + factorEx * rowEx[j] + factorEy * rowEy[j];
  1325.         }

  1327.         return jacobian;

  1329.     }

  1331.     /** {@inheritDoc} */
  1332.     @Override
  1333.     public void addKeplerContribution(final PositionAngleType type, final double gm,
  1334.                                       final double[] pDot) {
  1335.         pDot[5] += computeKeplerianAlphaDot(type, a, ex, ey, gm, cachedAlpha, cachedPositionAngleType);
  1336.     }

  1338.     /**
  1339.      * Compute rate of argument of latitude.
  1340.      * @param type position angle type of rate
  1341.      * @param a semi major axis
  1342.      * @param ex ex
  1343.      * @param ey ey
  1344.      * @param mu mu
  1345.      * @param alpha argument of latitude
  1346.      * @param cachedType position angle type of passed alpha
  1347.      * @return first-order time derivative for alpha
  1348.      * @since 12.2
  1349.      */
  1350.     private static double computeKeplerianAlphaDot(final PositionAngleType type, final double a, final double ex,
  1351.                                                    final double ey, final double mu,
  1352.                                                    final double alpha, final PositionAngleType cachedType) {
  1353.         final double n  = FastMath.sqrt(mu / a) / a;
  1354.         if (type == PositionAngleType.MEAN) {
  1355.             return n;
  1356.         }
  1357.         final double ksi;
  1358.         final SinCos sc;
  1359.         if (type == PositionAngleType.ECCENTRIC) {
  1360.             sc = FastMath.sinCos(CircularLatitudeArgumentUtility.convertAlpha(cachedType, alpha, ex, ey, type));
  1361.             ksi   = 1. / (1 - ex * sc.cos() - ey * sc.sin());
  1362.             return n * ksi;
  1363.         } else { // TRUE
  1364.             sc = FastMath.sinCos(CircularLatitudeArgumentUtility.convertAlpha(cachedType, alpha, ex, ey, type));
  1365.             final double oMe2  = 1 - ex * ex - ey * ey;
  1366.             ksi   = 1 + ex * sc.cos() + ey * sc.sin();
  1367.             return n * ksi * ksi / (oMe2 * FastMath.sqrt(oMe2));
  1368.         }
  1369.     }

  1371.     /**  Returns a string representation of this Orbit object.
  1372.      * @return a string representation of this object
  1373.      */
  1374.     public String toString() {
  1375.         return new StringBuilder().append("circular parameters: ").append('{').
  1376.                                   append("a: ").append(a).
  1377.                                   append(", ex: ").append(ex).append(", ey: ").append(ey).
  1378.                                   append(", i: ").append(FastMath.toDegrees(i)).
  1379.                                   append(", raan: ").append(FastMath.toDegrees(raan)).
  1380.                                   append(", alphaV: ").append(FastMath.toDegrees(getAlphaV())).
  1381.                                   append(";}").toString();
  1382.     }

  1384.     /** {@inheritDoc} */
  1385.     @Override
  1386.     public PositionAngleType getCachedPositionAngleType() {
  1387.         return cachedPositionAngleType;
  1388.     }

  1390.     /** {@inheritDoc} */
  1391.     @Override
  1392.     public boolean hasNonKeplerianRates() {
  1393.         return hasNonKeplerianAcceleration();
  1394.     }

  1396.     /** {@inheritDoc} */
  1397.     @Override
  1398.     public CircularOrbit withKeplerianRates() {
  1399.         final PositionAngleType positionAngleType = getCachedPositionAngleType();
  1400.         return new CircularOrbit(a, ex, ey, i, raan, cachedAlpha, positionAngleType, positionAngleType,
  1401.                 getFrame(), getDate(), getMu());
  1402.     }

  1404. }
Created with JaCoCo 0.8.12.202403310830