CircularOrbit.java

  1. /* Copyright 2002-2017 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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;

  18. import java.io.Serializable;
  19. import java.util.List;
  20. import java.util.stream.Collectors;
  21. import java.util.stream.Stream;

  22. import org.hipparchus.analysis.differentiation.DSFactory;
  23. import org.hipparchus.analysis.differentiation.DerivativeStructure;
  24. import org.hipparchus.analysis.interpolation.HermiteInterpolator;
  25. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  26. import org.hipparchus.util.FastMath;
  27. import org.hipparchus.util.MathUtils;
  28. import org.orekit.errors.OrekitIllegalArgumentException;
  29. import org.orekit.errors.OrekitInternalError;
  30. import org.orekit.errors.OrekitMessages;
  31. import org.orekit.frames.Frame;
  32. import org.orekit.time.AbsoluteDate;
  33. import org.orekit.utils.PVCoordinates;
  34. import org.orekit.utils.TimeStampedPVCoordinates;


  35. /**
  36.  * 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.  */

  72. public class CircularOrbit
  73.     extends Orbit {

  74.     /** Serializable UID. */
  75.     private static final long serialVersionUID = 20170414L;

  76.     /** Factory for first time derivatives. */
  77.     private static final DSFactory FACTORY = new DSFactory(1, 1);

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

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

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

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

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

  88.     /** True latitude argument (rad). */
  89.     private final double alphaV;

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

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

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

  96.     /** Inclination derivative (rad/s). */
  97.     private final double iDot;

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

  100.     /** True latitude argument derivative (rad/s). */
  101.     private final double alphaVDot;

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

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

  106.     /** Creates a new instance.
  107.      * @param a  semi-major axis (m)
  108.      * @param ex e cos(ω), first component of circular eccentricity vector
  109.      * @param ey e sin(ω), second component of circular eccentricity vector
  110.      * @param i inclination (rad)
  111.      * @param raan right ascension of ascending node (Ω, rad)
  112.      * @param alpha  an + ω, mean, eccentric or true latitude argument (rad)
  113.      * @param type type of latitude argument
  114.      * @param frame the frame in which are defined the parameters
  115.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  116.      * @param date date of the orbital parameters
  117.      * @param mu central attraction coefficient (m³/s²)
  118.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  119.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  120.      */
  121.     public CircularOrbit(final double a, final double ex, final double ey,
  122.                          final double i, final double raan, final double alpha,
  123.                          final PositionAngle type,
  124.                          final Frame frame, final AbsoluteDate date, final double mu)
  125.         throws IllegalArgumentException {
  126.         this(a, ex, ey, i, raan, alpha,
  127.              Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN,
  128.              type, frame, date, mu);
  129.     }

  130.     /** Creates a new instance.
  131.      * @param a  semi-major axis (m)
  132.      * @param ex e cos(ω), first component of circular eccentricity vector
  133.      * @param ey e sin(ω), second component of circular eccentricity vector
  134.      * @param i inclination (rad)
  135.      * @param raan right ascension of ascending node (Ω, rad)
  136.      * @param alpha  an + ω, mean, eccentric or true latitude argument (rad)
  137.      * @param aDot  semi-major axis derivative (m/s)
  138.      * @param exDot d(e cos(ω))/dt, first component of circular eccentricity vector derivative
  139.      * @param eyDot d(e sin(ω))/dt, second component of circular eccentricity vector derivative
  140.      * @param iDot inclination  derivative(rad/s)
  141.      * @param raanDot right ascension of ascending node derivative (rad/s)
  142.      * @param alphaDot  d(an + ω), mean, eccentric or true latitude argument derivative (rad/s)
  143.      * @param type type of latitude argument
  144.      * @param frame the frame in which are defined the parameters
  145.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  146.      * @param date date of the orbital parameters
  147.      * @param mu central attraction coefficient (m³/s²)
  148.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  149.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  150.      */
  151.     public CircularOrbit(final double a, final double ex, final double ey,
  152.                          final double i, final double raan, final double alpha,
  153.                          final double aDot, final double exDot, final double eyDot,
  154.                          final double iDot, final double raanDot, final double alphaDot,
  155.                          final PositionAngle type,
  156.                          final Frame frame, final AbsoluteDate date, final double mu)
  157.         throws IllegalArgumentException {
  158.         super(frame, date, mu);
  159.         if (ex * ex + ey * ey >= 1.0) {
  160.             throw new OrekitIllegalArgumentException(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS,
  161.                                                      getClass().getName());
  162.         }
  163.         this.a       =  a;
  164.         this.aDot    =  aDot;
  165.         this.ex      = ex;
  166.         this.exDot   = exDot;
  167.         this.ey      = ey;
  168.         this.eyDot   = eyDot;
  169.         this.i       = i;
  170.         this.iDot    = iDot;
  171.         this.raan    = raan;
  172.         this.raanDot = raanDot;

  173.         if (hasDerivatives()) {
  174.             final DerivativeStructure exDS    = FACTORY.build(ex,    exDot);
  175.             final DerivativeStructure eyDS    = FACTORY.build(ey,    eyDot);
  176.             final DerivativeStructure alphaDS = FACTORY.build(alpha, alphaDot);
  177.             final DerivativeStructure alphavDS;
  178.             switch (type) {
  179.                 case MEAN :
  180.                     alphavDS = FieldCircularOrbit.eccentricToTrue(FieldCircularOrbit.meanToEccentric(alphaDS, exDS, eyDS), exDS, eyDS);
  181.                     break;
  182.                 case ECCENTRIC :
  183.                     alphavDS = FieldCircularOrbit.eccentricToTrue(alphaDS, exDS, eyDS);
  184.                     break;
  185.                 case TRUE :
  186.                     alphavDS = alphaDS;
  187.                     break;
  188.                 default :
  189.                     throw new OrekitInternalError(null);
  190.             }
  191.             this.alphaV    = alphavDS.getValue();
  192.             this.alphaVDot = alphavDS.getPartialDerivative(1);
  193.         } else {
  194.             switch (type) {
  195.                 case MEAN :
  196.                     this.alphaV = eccentricToTrue(meanToEccentric(alpha, ex, ey), ex, ey);
  197.                     break;
  198.                 case ECCENTRIC :
  199.                     this.alphaV = eccentricToTrue(alpha, ex, ey);
  200.                     break;
  201.                 case TRUE :
  202.                     this.alphaV = alpha;
  203.                     break;
  204.                 default :
  205.                     throw new OrekitInternalError(null);
  206.             }
  207.             this.alphaVDot = Double.NaN;
  208.         }

  209.         serializePV = false;
  210.         partialPV   = null;

  211.     }

  212.     /** Creates a new instance.
  213.      * @param a  semi-major axis (m)
  214.      * @param ex e cos(ω), first component of circular eccentricity vector
  215.      * @param ey e sin(ω), second component of circular eccentricity vector
  216.      * @param i inclination (rad)
  217.      * @param raan right ascension of ascending node (Ω, rad)
  218.      * @param alphaV  v + ω, true latitude argument (rad)
  219.      * @param aDot  semi-major axis derivative (m/s)
  220.      * @param exDot d(e cos(ω))/dt, first component of circular eccentricity vector derivative
  221.      * @param eyDot d(e sin(ω))/dt, second component of circular eccentricity vector derivative
  222.      * @param iDot inclination  derivative(rad/s)
  223.      * @param raanDot right ascension of ascending node derivative (rad/s)
  224.      * @param alphaVDot  d(v + ω), true latitude argument derivative (rad/s)
  225.      * @param pvCoordinates the {@link PVCoordinates} in inertial frame
  226.      * @param frame the frame in which are defined the parameters
  227.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  228.      * @param mu central attraction coefficient (m³/s²)
  229.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  230.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  231.      */
  232.     private CircularOrbit(final double a, final double ex, final double ey,
  233.                           final double i, final double raan, final double alphaV,
  234.                           final double aDot, final double exDot, final double eyDot,
  235.                           final double iDot, final double raanDot, final double alphaVDot,
  236.                           final TimeStampedPVCoordinates pvCoordinates, final Frame frame,
  237.                           final double mu)
  238.         throws IllegalArgumentException {
  239.         super(pvCoordinates, frame, mu);
  240.         this.a           =  a;
  241.         this.aDot        =  aDot;
  242.         this.ex          = ex;
  243.         this.exDot       = exDot;
  244.         this.ey          = ey;
  245.         this.eyDot       = eyDot;
  246.         this.i           = i;
  247.         this.iDot        = iDot;
  248.         this.raan        = raan;
  249.         this.raanDot     = raanDot;
  250.         this.alphaV      = alphaV;
  251.         this.alphaVDot   = alphaVDot;
  252.         this.serializePV = true;
  253.         this.partialPV   = null;
  254.     }

  255.     /** Constructor from Cartesian parameters.
  256.      *
  257.      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
  258.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  259.      * use {@code mu} and the position to compute the acceleration, including
  260.      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
  261.      *
  262.      * @param pvCoordinates the {@link PVCoordinates} in inertial frame
  263.      * @param frame the frame in which are defined the {@link PVCoordinates}
  264.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  265.      * @param mu central attraction coefficient (m³/s²)
  266.      * @exception IllegalArgumentException if frame is not a {@link
  267.      * Frame#isPseudoInertial pseudo-inertial frame}
  268.      */
  269.     public CircularOrbit(final TimeStampedPVCoordinates pvCoordinates, final Frame frame, final double mu)
  270.         throws IllegalArgumentException {
  271.         super(pvCoordinates, frame, mu);

  272.         // compute semi-major axis
  273.         final Vector3D pvP = pvCoordinates.getPosition();
  274.         final Vector3D pvV = pvCoordinates.getVelocity();
  275.         final Vector3D pvA = pvCoordinates.getAcceleration();
  276.         final double r2 = pvP.getNormSq();
  277.         final double r  = FastMath.sqrt(r2);
  278.         final double V2 = pvV.getNormSq();
  279.         final double rV2OnMu = r * V2 / mu;

  280.         if (rV2OnMu > 2) {
  281.             throw new OrekitIllegalArgumentException(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS,
  282.                                                      getClass().getName());
  283.         }

  284.         a = r / (2 - rV2OnMu);

  285.         // compute inclination
  286.         final Vector3D momentum = pvCoordinates.getMomentum();
  287.         i = Vector3D.angle(momentum, Vector3D.PLUS_K);

  288.         // compute right ascension of ascending node
  289.         final Vector3D node  = Vector3D.crossProduct(Vector3D.PLUS_K, momentum);
  290.         raan = FastMath.atan2(node.getY(), node.getX());

  291.         // 2D-coordinates in the canonical frame
  292.         final double cosRaan = FastMath.cos(raan);
  293.         final double sinRaan = FastMath.sin(raan);
  294.         final double cosI    = FastMath.cos(i);
  295.         final double sinI    = FastMath.sin(i);
  296.         final double xP      = pvP.getX();
  297.         final double yP      = pvP.getY();
  298.         final double zP      = pvP.getZ();
  299.         final double x2      = (xP * cosRaan + yP * sinRaan) / a;
  300.         final double y2      = ((yP * cosRaan - xP * sinRaan) * cosI + zP * sinI) / a;

  301.         // compute eccentricity vector
  302.         final double eSE    = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(mu * a);
  303.         final double eCE    = rV2OnMu - 1;
  304.         final double e2     = eCE * eCE + eSE * eSE;
  305.         final double f      = eCE - e2;
  306.         final double g      = FastMath.sqrt(1 - e2) * eSE;
  307.         final double aOnR   = a / r;
  308.         final double a2OnR2 = aOnR * aOnR;
  309.         ex = a2OnR2 * (f * x2 + g * y2);
  310.         ey = a2OnR2 * (f * y2 - g * x2);

  311.         // compute latitude argument
  312.         final double beta = 1 / (1 + FastMath.sqrt(1 - ex * ex - ey * ey));
  313.         alphaV = eccentricToTrue(FastMath.atan2(y2 + ey + eSE * beta * ex, x2 + ex - eSE * beta * ey), ex, ey);

  314.         partialPV   = pvCoordinates;

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

  317.             final double[][] jacobian = new double[6][6];
  318.             getJacobianWrtCartesian(PositionAngle.MEAN, jacobian);

  319.             final Vector3D keplerianAcceleration    = new Vector3D(-mu / (r * r2), pvP);
  320.             final Vector3D nonKeplerianAcceleration = pvA.subtract(keplerianAcceleration);
  321.             final double   aX                       = nonKeplerianAcceleration.getX();
  322.             final double   aY                       = nonKeplerianAcceleration.getY();
  323.             final double   aZ                       = nonKeplerianAcceleration.getZ();
  324.             aDot    = jacobian[0][3] * aX + jacobian[0][4] * aY + jacobian[0][5] * aZ;
  325.             exDot   = jacobian[1][3] * aX + jacobian[1][4] * aY + jacobian[1][5] * aZ;
  326.             eyDot   = jacobian[2][3] * aX + jacobian[2][4] * aY + jacobian[2][5] * aZ;
  327.             iDot    = jacobian[3][3] * aX + jacobian[3][4] * aY + jacobian[3][5] * aZ;
  328.             raanDot = jacobian[4][3] * aX + jacobian[4][4] * aY + jacobian[4][5] * aZ;

  329.             // in order to compute true anomaly derivative, we must compute
  330.             // mean anomaly derivative including Keplerian motion and convert to true anomaly
  331.             final double alphaMDot = getKeplerianMeanMotion() +
  332.                                      jacobian[5][3] * aX + jacobian[5][4] * aY + jacobian[5][5] * aZ;
  333.             final DerivativeStructure exDS     = FACTORY.build(ex, exDot);
  334.             final DerivativeStructure eyDS     = FACTORY.build(ey, eyDot);
  335.             final DerivativeStructure alphaMDS = FACTORY.build(getAlphaM(), alphaMDot);
  336.             final DerivativeStructure alphavDS = FieldCircularOrbit.eccentricToTrue(FieldCircularOrbit.meanToEccentric(alphaMDS, exDS, eyDS), exDS, eyDS);
  337.             alphaVDot = alphavDS.getPartialDerivative(1);

  338.         } else {
  339.             // acceleration is either almost zero or NaN,
  340.             // we assume acceleration was not known
  341.             // we don't set up derivatives
  342.             aDot      = Double.NaN;
  343.             exDot     = Double.NaN;
  344.             eyDot     = Double.NaN;
  345.             iDot      = Double.NaN;
  346.             raanDot   = Double.NaN;
  347.             alphaVDot = Double.NaN;
  348.         }

  349.         serializePV = true;

  350.     }

  351.     /** Constructor from Cartesian parameters.
  352.      *
  353.      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
  354.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  355.      * use {@code mu} and the position to compute the acceleration, including
  356.      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
  357.      *
  358.      * @param pvCoordinates the {@link PVCoordinates} in inertial frame
  359.      * @param frame the frame in which are defined the {@link PVCoordinates}
  360.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  361.      * @param date date of the orbital parameters
  362.      * @param mu central attraction coefficient (m³/s²)
  363.      * @exception IllegalArgumentException if frame is not a {@link
  364.      * Frame#isPseudoInertial pseudo-inertial frame}
  365.      */
  366.     public CircularOrbit(final PVCoordinates pvCoordinates, final Frame frame,
  367.                          final AbsoluteDate date, final double mu)
  368.         throws IllegalArgumentException {
  369.         this(new TimeStampedPVCoordinates(date, pvCoordinates), frame, mu);
  370.     }

  371.     /** Constructor from any kind of orbital parameters.
  372.      * @param op orbital parameters to copy
  373.      */
  374.     public CircularOrbit(final Orbit op) {

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

  376.         a    = op.getA();
  377.         i    = op.getI();
  378.         final double hx = op.getHx();
  379.         final double hy = op.getHy();
  380.         final double h2 = hx * hx + hy * hy;
  381.         final double h  = FastMath.sqrt(h2);
  382.         raan = FastMath.atan2(hy, hx);
  383.         final double cosRaan = h == 0 ? FastMath.cos(raan) : hx / h;
  384.         final double sinRaan = h == 0 ? FastMath.sin(raan) : hy / h;
  385.         final double equiEx = op.getEquinoctialEx();
  386.         final double equiEy = op.getEquinoctialEy();
  387.         ex     = equiEx * cosRaan + equiEy * sinRaan;
  388.         ey     = equiEy * cosRaan - equiEx * sinRaan;
  389.         alphaV = op.getLv() - raan;

  390.         if (op.hasDerivatives()) {
  391.             aDot    = op.getADot();
  392.             final double hxDot = op.getHxDot();
  393.             final double hyDot = op.getHyDot();
  394.             iDot    = 2 * (cosRaan * hxDot + sinRaan * hyDot) / (1 + h2);
  395.             raanDot = (hx * hyDot - hy * hxDot) / h2;
  396.             final double equiExDot = op.getEquinoctialExDot();
  397.             final double equiEyDot = op.getEquinoctialEyDot();
  398.             exDot   = (equiExDot + equiEy * raanDot) * cosRaan +
  399.                       (equiEyDot - equiEx * raanDot) * sinRaan;
  400.             eyDot   = (equiEyDot - equiEx * raanDot) * cosRaan -
  401.                       (equiExDot + equiEy * raanDot) * sinRaan;
  402.             alphaVDot = op.getLvDot() - raanDot;
  403.         } else {
  404.             aDot      = Double.NaN;
  405.             exDot     = Double.NaN;
  406.             eyDot     = Double.NaN;
  407.             iDot      = Double.NaN;
  408.             raanDot   = Double.NaN;
  409.             alphaVDot = Double.NaN;
  410.         }

  411.         serializePV = false;
  412.         partialPV   = null;

  413.     }

  414.     /** {@inheritDoc} */
  415.     public OrbitType getType() {
  416.         return OrbitType.CIRCULAR;
  417.     }

  418.     /** {@inheritDoc} */
  419.     public double getA() {
  420.         return a;
  421.     }

  422.     /** {@inheritDoc} */
  423.     public double getADot() {
  424.         return aDot;
  425.     }

  426.     /** {@inheritDoc} */
  427.     public double getEquinoctialEx() {
  428.         final double cosRaan = FastMath.cos(raan);
  429.         final double sinRaan = FastMath.sin(raan);
  430.         return ex * cosRaan - ey * sinRaan;
  431.     }

  432.     /** {@inheritDoc} */
  433.     public double getEquinoctialExDot() {
  434.         final double cosRaan = FastMath.cos(raan);
  435.         final double sinRaan = FastMath.sin(raan);
  436.         return (exDot - ey * raanDot) * cosRaan - (eyDot + ex * raanDot) * sinRaan;
  437.     }

  438.     /** {@inheritDoc} */
  439.     public double getEquinoctialEy() {
  440.         final double cosRaan = FastMath.cos(raan);
  441.         final double sinRaan = FastMath.sin(raan);
  442.         return ey * cosRaan + ex * sinRaan;
  443.     }

  444.     /** {@inheritDoc} */
  445.     public double getEquinoctialEyDot() {
  446.         final double cosRaan = FastMath.cos(raan);
  447.         final double sinRaan = FastMath.sin(raan);
  448.         return (eyDot + ex * raanDot) * cosRaan + (exDot - ey * raanDot) * sinRaan;
  449.     }

  450.     /** Get the first component of the circular eccentricity vector.
  451.      * @return ex = e cos(ω), first component of the circular eccentricity vector
  452.      */
  453.     public double getCircularEx() {
  454.         return ex;
  455.     }

  456.     /** Get the first component of the circular eccentricity vector derivative.
  457.      * @return ex = e cos(ω), first component of the circular eccentricity vector derivative
  458.      * @since 9.0
  459.      */
  460.     public double getCircularExDot() {
  461.         return exDot;
  462.     }

  463.     /** Get the second component of the circular eccentricity vector.
  464.      * @return ey = e sin(ω), second component of the circular eccentricity vector
  465.      */
  466.     public double getCircularEy() {
  467.         return ey;
  468.     }

  469.     /** Get the second component of the circular eccentricity vector derivative.
  470.      * @return ey = e sin(ω), second component of the circular eccentricity vector derivative
  471.      */
  472.     public double getCircularEyDot() {
  473.         return eyDot;
  474.     }

  475.     /** {@inheritDoc} */
  476.     public double getHx() {
  477.         // Check for equatorial retrograde orbit
  478.         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
  479.             return Double.NaN;
  480.         }
  481.         return  FastMath.cos(raan) * FastMath.tan(i / 2);
  482.     }

  483.     /** {@inheritDoc} */
  484.     public double getHxDot() {
  485.         // Check for equatorial retrograde orbit
  486.         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
  487.             return Double.NaN;
  488.         }
  489.         final double cosRaan = FastMath.cos(raan);
  490.         final double sinRaan = FastMath.sin(raan);
  491.         final double tan     = FastMath.tan(0.5 * i);
  492.         return 0.5 * cosRaan * (1 + tan * tan) * iDot - sinRaan * tan * raanDot;
  493.     }

  494.     /** {@inheritDoc} */
  495.     public double getHy() {
  496.         // Check for equatorial retrograde orbit
  497.         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
  498.             return Double.NaN;
  499.         }
  500.         return  FastMath.sin(raan) * FastMath.tan(i / 2);
  501.     }

  502.     /** {@inheritDoc} */
  503.     public double getHyDot() {
  504.         // Check for equatorial retrograde orbit
  505.         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
  506.             return Double.NaN;
  507.         }
  508.         final double cosRaan = FastMath.cos(raan);
  509.         final double sinRaan = FastMath.sin(raan);
  510.         final double tan     = FastMath.tan(0.5 * i);
  511.         return 0.5 * sinRaan * (1 + tan * tan) * iDot + cosRaan * tan * raanDot;
  512.     }

  513.     /** Get the true latitude argument.
  514.      * @return v + ω true latitude argument (rad)
  515.      */
  516.     public double getAlphaV() {
  517.         return alphaV;
  518.     }

  519.     /** Get the true latitude argument derivative.
  520.      * <p>
  521.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  522.      * </p>
  523.      * @return v + ω true latitude argument derivative (rad/s)
  524.      * @since 9.0
  525.      */
  526.     public double getAlphaVDot() {
  527.         return alphaVDot;
  528.     }

  529.     /** Get the eccentric latitude argument.
  530.      * @return E + ω eccentric latitude argument (rad)
  531.      */
  532.     public double getAlphaE() {
  533.         return trueToEccentric(alphaV, ex, ey);
  534.     }

  535.     /** Get the eccentric latitude argument derivative.
  536.      * <p>
  537.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  538.      * </p>
  539.      * @return d(E + ω)/dt eccentric latitude argument derivative (rad/s)
  540.      * @since 9.0
  541.      */
  542.     public double getAlphaEDot() {
  543.         final DerivativeStructure alphaVDS = FACTORY.build(alphaV, alphaVDot);
  544.         final DerivativeStructure exDS     = FACTORY.build(ex,     exDot);
  545.         final DerivativeStructure eyDS     = FACTORY.build(ey,     eyDot);
  546.         final DerivativeStructure alphaEDS = FieldCircularOrbit.trueToEccentric(alphaVDS, exDS, eyDS);
  547.         return alphaEDS.getPartialDerivative(1);
  548.     }

  549.     /** Get the mean latitude argument.
  550.      * @return M + ω mean latitude argument (rad)
  551.      */
  552.     public double getAlphaM() {
  553.         return eccentricToMean(trueToEccentric(alphaV, ex, ey), ex, ey);
  554.     }

  555.     /** Get the mean latitude argument derivative.
  556.      * <p>
  557.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  558.      * </p>
  559.      * @return d(M + ω)/dt mean latitude argument derivative (rad/s)
  560.      * @since 9.0
  561.      */
  562.     public double getAlphaMDot() {
  563.         final DerivativeStructure alphaVDS = FACTORY.build(alphaV, alphaVDot);
  564.         final DerivativeStructure exDS     = FACTORY.build(ex,     exDot);
  565.         final DerivativeStructure eyDS     = FACTORY.build(ey,     eyDot);
  566.         final DerivativeStructure alphaMDS = FieldCircularOrbit.eccentricToMean(FieldCircularOrbit.trueToEccentric(alphaVDS, exDS, eyDS), exDS, eyDS);
  567.         return alphaMDS.getPartialDerivative(1);
  568.     }

  569.     /** Get the latitude argument.
  570.      * @param type type of the angle
  571.      * @return latitude argument (rad)
  572.      */
  573.     public double getAlpha(final PositionAngle type) {
  574.         return (type == PositionAngle.MEAN) ? getAlphaM() :
  575.                                               ((type == PositionAngle.ECCENTRIC) ? getAlphaE() :
  576.                                                                                    getAlphaV());
  577.     }

  578.     /** Get the latitude argument derivative.
  579.      * <p>
  580.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  581.      * </p>
  582.      * @param type type of the angle
  583.      * @return latitude argument derivative (rad/s)
  584.      * @since 9.0
  585.      */
  586.     public double getAlphaDot(final PositionAngle type) {
  587.         return (type == PositionAngle.MEAN) ? getAlphaMDot() :
  588.                                               ((type == PositionAngle.ECCENTRIC) ? getAlphaEDot() :
  589.                                                                                    getAlphaVDot());
  590.     }

  591.     /** Computes the true latitude argument from the eccentric latitude argument.
  592.      * @param alphaE = E + ω eccentric latitude argument (rad)
  593.      * @param ex e cos(ω), first component of circular eccentricity vector
  594.      * @param ey e sin(ω), second component of circular eccentricity vector
  595.      * @return the true latitude argument.
  596.      */
  597.     public static double eccentricToTrue(final double alphaE, final double ex, final double ey) {
  598.         final double epsilon   = FastMath.sqrt(1 - ex * ex - ey * ey);
  599.         final double cosAlphaE = FastMath.cos(alphaE);
  600.         final double sinAlphaE = FastMath.sin(alphaE);
  601.         return alphaE + 2 * FastMath.atan((ex * sinAlphaE - ey * cosAlphaE) /
  602.                                       (epsilon + 1 - ex * cosAlphaE - ey * sinAlphaE));
  603.     }

  604.     /** Computes the eccentric latitude argument from the true latitude argument.
  605.      * @param alphaV = V + ω true latitude argument (rad)
  606.      * @param ex e cos(ω), first component of circular eccentricity vector
  607.      * @param ey e sin(ω), second component of circular eccentricity vector
  608.      * @return the eccentric latitude argument.
  609.      */
  610.     public static double trueToEccentric(final double alphaV, final double ex, final double ey) {
  611.         final double epsilon   = FastMath.sqrt(1 - ex * ex - ey * ey);
  612.         final double cosAlphaV = FastMath.cos(alphaV);
  613.         final double sinAlphaV = FastMath.sin(alphaV);
  614.         return alphaV + 2 * FastMath.atan((ey * cosAlphaV - ex * sinAlphaV) /
  615.                                       (epsilon + 1 + ex * cosAlphaV + ey * sinAlphaV));
  616.     }

  617.     /** Computes the eccentric latitude argument from the mean latitude argument.
  618.      * @param alphaM = M + ω  mean latitude argument (rad)
  619.      * @param ex e cos(ω), first component of circular eccentricity vector
  620.      * @param ey e sin(ω), second component of circular eccentricity vector
  621.      * @return the eccentric latitude argument.
  622.      */
  623.     public static double meanToEccentric(final double alphaM, final double ex, final double ey) {
  624.         // Generalization of Kepler equation to circular parameters
  625.         // with alphaE = PA + E and
  626.         //      alphaM = PA + M = alphaE - ex.sin(alphaE) + ey.cos(alphaE)
  627.         double alphaE        = alphaM;
  628.         double shift         = 0.0;
  629.         double alphaEMalphaM = 0.0;
  630.         double cosAlphaE     = FastMath.cos(alphaE);
  631.         double sinAlphaE     = FastMath.sin(alphaE);
  632.         int    iter          = 0;
  633.         do {
  634.             final double f2 = ex * sinAlphaE - ey * cosAlphaE;
  635.             final double f1 = 1.0 - ex * cosAlphaE - ey * sinAlphaE;
  636.             final double f0 = alphaEMalphaM - f2;

  637.             final double f12 = 2.0 * f1;
  638.             shift = f0 * f12 / (f1 * f12 - f0 * f2);

  639.             alphaEMalphaM -= shift;
  640.             alphaE         = alphaM + alphaEMalphaM;
  641.             cosAlphaE      = FastMath.cos(alphaE);
  642.             sinAlphaE      = FastMath.sin(alphaE);

  643.         } while ((++iter < 50) && (FastMath.abs(shift) > 1.0e-12));

  644.         return alphaE;

  645.     }

  646.     /** Computes the mean latitude argument from the eccentric latitude argument.
  647.      * @param alphaE = E + ω  mean latitude argument (rad)
  648.      * @param ex e cos(ω), first component of circular eccentricity vector
  649.      * @param ey e sin(ω), second component of circular eccentricity vector
  650.      * @return the mean latitude argument.
  651.      */
  652.     public static double eccentricToMean(final double alphaE, final double ex, final double ey) {
  653.         return alphaE + (ey * FastMath.cos(alphaE) - ex * FastMath.sin(alphaE));
  654.     }

  655.     /** {@inheritDoc} */
  656.     public double getE() {
  657.         return FastMath.sqrt(ex * ex + ey * ey);
  658.     }

  659.     /** {@inheritDoc} */
  660.     public double getEDot() {
  661.         return (ex * exDot + ey * eyDot) / getE();
  662.     }

  663.     /** {@inheritDoc} */
  664.     public double getI() {
  665.         return i;
  666.     }

  667.     /** {@inheritDoc} */
  668.     public double getIDot() {
  669.         return iDot;
  670.     }

  671.     /** Get the right ascension of the ascending node.
  672.      * @return right ascension of the ascending node (rad)
  673.      */
  674.     public double getRightAscensionOfAscendingNode() {
  675.         return raan;
  676.     }

  677.     /** Get the right ascension of the ascending node derivative.
  678.      * <p>
  679.      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
  680.      * </p>
  681.      * @return right ascension of the ascending node derivative (rad/s)
  682.      * @since 9.0
  683.      */
  684.     public double getRightAscensionOfAscendingNodeDot() {
  685.         return raanDot;
  686.     }

  687.     /** {@inheritDoc} */
  688.     public double getLv() {
  689.         return alphaV + raan;
  690.     }

  691.     /** {@inheritDoc} */
  692.     public double getLvDot() {
  693.         return alphaVDot + raanDot;
  694.     }

  695.     /** {@inheritDoc} */
  696.     public double getLE() {
  697.         return getAlphaE() + raan;
  698.     }

  699.     /** {@inheritDoc} */
  700.     public double getLEDot() {
  701.         return getAlphaEDot() + raanDot;
  702.     }

  703.     /** {@inheritDoc} */
  704.     public double getLM() {
  705.         return getAlphaM() + raan;
  706.     }

  707.     /** {@inheritDoc} */
  708.     public double getLMDot() {
  709.         return getAlphaMDot() + raanDot;
  710.     }

  711.     /** Compute position and velocity but not acceleration.
  712.      */
  713.     private void computePVWithoutA() {

  714.         if (partialPV != null) {
  715.             // already computed
  716.             return;
  717.         }

  718.         // get equinoctial parameters
  719.         final double equEx = getEquinoctialEx();
  720.         final double equEy = getEquinoctialEy();
  721.         final double hx = getHx();
  722.         final double hy = getHy();
  723.         final double lE = getLE();

  724.         // inclination-related intermediate parameters
  725.         final double hx2   = hx * hx;
  726.         final double hy2   = hy * hy;
  727.         final double factH = 1. / (1 + hx2 + hy2);

  728.         // reference axes defining the orbital plane
  729.         final double ux = (1 + hx2 - hy2) * factH;
  730.         final double uy =  2 * hx * hy * factH;
  731.         final double uz = -2 * hy * factH;

  732.         final double vx = uy;
  733.         final double vy = (1 - hx2 + hy2) * factH;
  734.         final double vz =  2 * hx * factH;

  735.         // eccentricity-related intermediate parameters
  736.         final double exey = equEx * equEy;
  737.         final double ex2  = equEx * equEx;
  738.         final double ey2  = equEy * equEy;
  739.         final double e2   = ex2 + ey2;
  740.         final double eta  = 1 + FastMath.sqrt(1 - e2);
  741.         final double beta = 1. / eta;

  742.         // eccentric latitude argument
  743.         final double cLe    = FastMath.cos(lE);
  744.         final double sLe    = FastMath.sin(lE);
  745.         final double exCeyS = equEx * cLe + equEy * sLe;

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

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

  752.         final Vector3D position =
  753.                         new Vector3D(x * ux + y * vx, x * uy + y * vy, x * uz + y * vz);
  754.         final Vector3D velocity =
  755.                         new Vector3D(xdot * ux + ydot * vx, xdot * uy + ydot * vy, xdot * uz + ydot * vz);

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

  757.     }

  758.     /** Compute non-Keplerian part of the acceleration from first time derivatives.
  759.      * <p>
  760.      * This method should be called only when {@link #hasDerivatives()} returns true.
  761.      * </p>
  762.      * @return non-Keplerian part of the acceleration
  763.      */
  764.     private Vector3D nonKeplerianAcceleration() {

  765.         final double[][] dCdP = new double[6][6];
  766.         getJacobianWrtParameters(PositionAngle.MEAN, dCdP);

  767.         final double nonKeplerianMeanMotion = getAlphaMDot() - getKeplerianMeanMotion();
  768.         final double nonKeplerianAx = dCdP[3][0] * aDot    + dCdP[3][1] * exDot   + dCdP[3][2] * eyDot   +
  769.                                       dCdP[3][3] * iDot    + dCdP[3][4] * raanDot + dCdP[3][5] * nonKeplerianMeanMotion;
  770.         final double nonKeplerianAy = dCdP[4][0] * aDot    + dCdP[4][1] * exDot   + dCdP[4][2] * eyDot   +
  771.                                       dCdP[4][3] * iDot    + dCdP[4][4] * raanDot + dCdP[4][5] * nonKeplerianMeanMotion;
  772.         final double nonKeplerianAz = dCdP[5][0] * aDot    + dCdP[5][1] * exDot   + dCdP[5][2] * eyDot   +
  773.                                       dCdP[5][3] * iDot    + dCdP[5][4] * raanDot + dCdP[5][5] * nonKeplerianMeanMotion;

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

  775.     }

  776.     /** {@inheritDoc} */
  777.     protected TimeStampedPVCoordinates initPVCoordinates() {

  778.         // position and velocity
  779.         computePVWithoutA();

  780.         // acceleration
  781.         final double r2 = partialPV.getPosition().getNormSq();
  782.         final Vector3D keplerianAcceleration = new Vector3D(-getMu() / (r2 * FastMath.sqrt(r2)), partialPV.getPosition());
  783.         final Vector3D acceleration = hasDerivatives() ?
  784.                                       keplerianAcceleration.add(nonKeplerianAcceleration()) :
  785.                                       keplerianAcceleration;

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

  787.     }

  788.     /** {@inheritDoc} */
  789.     public CircularOrbit shiftedBy(final double dt) {

  790.         // use Keplerian-only motion
  791.         final CircularOrbit keplerianShifted = new CircularOrbit(a, ex, ey, i, raan,
  792.                                                                  getAlphaM() + getKeplerianMeanMotion() * dt,
  793.                                                                  PositionAngle.MEAN, getFrame(),
  794.                                                                  getDate().shiftedBy(dt), getMu());

  795.         if (hasDerivatives()) {

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

  798.             // add quadratic effect of non-Keplerian acceleration to Keplerian-only shift
  799.             keplerianShifted.computePVWithoutA();
  800.             final Vector3D fixedP   = new Vector3D(1, keplerianShifted.partialPV.getPosition(),
  801.                                                    0.5 * dt * dt, nonKeplerianAcceleration);
  802.             final double   fixedR2 = fixedP.getNormSq();
  803.             final double   fixedR  = FastMath.sqrt(fixedR2);
  804.             final Vector3D fixedV  = new Vector3D(1, keplerianShifted.partialPV.getVelocity(),
  805.                                                   dt, nonKeplerianAcceleration);
  806.             final Vector3D fixedA  = new Vector3D(-getMu() / (fixedR2 * fixedR), keplerianShifted.partialPV.getPosition(),
  807.                                                   1, nonKeplerianAcceleration);

  808.             // build a new orbit, taking non-Keplerian acceleration into account
  809.             return new CircularOrbit(new TimeStampedPVCoordinates(keplerianShifted.getDate(),
  810.                                                                   fixedP, fixedV, fixedA),
  811.                                      keplerianShifted.getFrame(), keplerianShifted.getMu());

  812.         } else {
  813.             // Keplerian-only motion is all we can do
  814.             return keplerianShifted;
  815.         }

  816.     }

  817.     /** {@inheritDoc}
  818.      * <p>
  819.      * The interpolated instance is created by polynomial Hermite interpolation
  820.      * on circular elements, without derivatives (which means the interpolation
  821.      * falls back to Lagrange interpolation only).
  822.      * </p>
  823.      * <p>
  824.      * As this implementation of interpolation is polynomial, it should be used only
  825.      * with small samples (about 10-20 points) in order to avoid <a
  826.      * href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
  827.      * and numerical problems (including NaN appearing).
  828.      * </p>
  829.      * <p>
  830.      * If orbit interpolation on large samples is needed, using the {@link
  831.      * org.orekit.propagation.analytical.Ephemeris} class is a better way than using this
  832.      * low-level interpolation. The Ephemeris class automatically handles selection of
  833.      * a neighboring sub-sample with a predefined number of point from a large global sample
  834.      * in a thread-safe way.
  835.      * </p>
  836.      */
  837.     public CircularOrbit interpolate(final AbsoluteDate date, final Stream<Orbit> sample) {

  838.         // first pass to check if derivatives are available throughout the sample
  839.         final List<Orbit> list = sample.collect(Collectors.toList());
  840.         boolean useDerivatives = true;
  841.         for (final Orbit orbit : list) {
  842.             useDerivatives = useDerivatives && orbit.hasDerivatives();
  843.         }

  844.         // set up an interpolator
  845.         final HermiteInterpolator interpolator = new HermiteInterpolator();

  846.         // second pass to feed interpolator
  847.         AbsoluteDate previousDate   = null;
  848.         double       previousRAAN   = Double.NaN;
  849.         double       previousAlphaM = Double.NaN;
  850.         for (final Orbit orbit : list) {
  851.             final CircularOrbit circ = (CircularOrbit) OrbitType.CIRCULAR.convertType(orbit);
  852.             final double continuousRAAN;
  853.             final double continuousAlphaM;
  854.             if (previousDate == null) {
  855.                 continuousRAAN   = circ.getRightAscensionOfAscendingNode();
  856.                 continuousAlphaM = circ.getAlphaM();
  857.             } else {
  858.                 final double dt       = circ.getDate().durationFrom(previousDate);
  859.                 final double keplerAM = previousAlphaM + circ.getKeplerianMeanMotion() * dt;
  860.                 continuousRAAN   = MathUtils.normalizeAngle(circ.getRightAscensionOfAscendingNode(), previousRAAN);
  861.                 continuousAlphaM = MathUtils.normalizeAngle(circ.getAlphaM(), keplerAM);
  862.             }
  863.             previousDate   = circ.getDate();
  864.             previousRAAN   = continuousRAAN;
  865.             previousAlphaM = continuousAlphaM;
  866.             if (useDerivatives) {
  867.                 interpolator.addSamplePoint(circ.getDate().durationFrom(date),
  868.                                             new double[] {
  869.                                                 circ.getA(),
  870.                                                 circ.getCircularEx(),
  871.                                                 circ.getCircularEy(),
  872.                                                 circ.getI(),
  873.                                                 continuousRAAN,
  874.                                                 continuousAlphaM
  875.                                             }, new double[] {
  876.                                                 circ.getADot(),
  877.                                                 circ.getCircularExDot(),
  878.                                                 circ.getCircularEyDot(),
  879.                                                 circ.getIDot(),
  880.                                                 circ.getRightAscensionOfAscendingNodeDot(),
  881.                                                 circ.getAlphaMDot()
  882.                                             });
  883.             } else {
  884.                 interpolator.addSamplePoint(circ.getDate().durationFrom(date),
  885.                                             new double[] {
  886.                                                 circ.getA(),
  887.                                                 circ.getCircularEx(),
  888.                                                 circ.getCircularEy(),
  889.                                                 circ.getI(),
  890.                                                 continuousRAAN,
  891.                                                 continuousAlphaM
  892.                                             });
  893.             }
  894.         }

  895.         // interpolate
  896.         final double[][] interpolated = interpolator.derivatives(0.0, 1);

  897.         // build a new interpolated instance
  898.         return new CircularOrbit(interpolated[0][0], interpolated[0][1], interpolated[0][2],
  899.                                  interpolated[0][3], interpolated[0][4], interpolated[0][5],
  900.                                  interpolated[1][0], interpolated[1][1], interpolated[1][2],
  901.                                  interpolated[1][3], interpolated[1][4], interpolated[1][5],
  902.                                  PositionAngle.MEAN, getFrame(), date, getMu());

  903.     }

  904.     /** {@inheritDoc} */
  905.     protected double[][] computeJacobianMeanWrtCartesian() {


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

  907.         computePVWithoutA();
  908.         final Vector3D position = partialPV.getPosition();
  909.         final Vector3D velocity = partialPV.getVelocity();
  910.         final double x          = position.getX();
  911.         final double y          = position.getY();
  912.         final double z          = position.getZ();
  913.         final double vx         = velocity.getX();
  914.         final double vy         = velocity.getY();
  915.         final double vz         = velocity.getZ();
  916.         final double pv         = Vector3D.dotProduct(position, velocity);
  917.         final double r2         = position.getNormSq();
  918.         final double r          = FastMath.sqrt(r2);
  919.         final double v2         = velocity.getNormSq();

  920.         final double mu         = getMu();
  921.         final double oOsqrtMuA  = 1 / FastMath.sqrt(mu * a);
  922.         final double rOa        = r / a;
  923.         final double aOr        = a / r;
  924.         final double aOr2       = a / r2;
  925.         final double a2         = a * a;

  926.         final double ex2        = ex * ex;
  927.         final double ey2        = ey * ey;
  928.         final double e2         = ex2 + ey2;
  929.         final double epsilon    = FastMath.sqrt(1 - e2);
  930.         final double beta       = 1 / (1 + epsilon);

  931.         final double eCosE      = 1 - rOa;
  932.         final double eSinE      = pv * oOsqrtMuA;

  933.         final double cosI       = FastMath.cos(i);
  934.         final double sinI       = FastMath.sin(i);
  935.         final double cosRaan    = FastMath.cos(raan);
  936.         final double sinRaan    = FastMath.sin(raan);

  937.         // da
  938.         fillHalfRow(2 * aOr * aOr2, position, jacobian[0], 0);
  939.         fillHalfRow(2 * a2 / mu, velocity, jacobian[0], 3);

  940.         // differentials of the normalized momentum
  941.         final Vector3D danP = new Vector3D(v2, position, -pv, velocity);
  942.         final Vector3D danV = new Vector3D(r2, velocity, -pv, position);
  943.         final double recip  = 1 / partialPV.getMomentum().getNorm();
  944.         final double recip2 = recip * recip;
  945.         final Vector3D dwXP = new Vector3D(recip, new Vector3D(  0,  vz, -vy), -recip2 * sinRaan * sinI, danP);
  946.         final Vector3D dwYP = new Vector3D(recip, new Vector3D(-vz,   0,  vx),  recip2 * cosRaan * sinI, danP);
  947.         final Vector3D dwZP = new Vector3D(recip, new Vector3D( vy, -vx,   0), -recip2 * cosI,           danP);
  948.         final Vector3D dwXV = new Vector3D(recip, new Vector3D(  0,  -z,   y), -recip2 * sinRaan * sinI, danV);
  949.         final Vector3D dwYV = new Vector3D(recip, new Vector3D(  z,   0,  -x),  recip2 * cosRaan * sinI, danV);
  950.         final Vector3D dwZV = new Vector3D(recip, new Vector3D( -y,   x,   0), -recip2 * cosI,           danV);

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

  954.         // dRaan
  955.         fillHalfRow(sinRaan / sinI, dwYP, cosRaan / sinI, dwXP, jacobian[4], 0);
  956.         fillHalfRow(sinRaan / sinI, dwYV, cosRaan / sinI, dwXV, jacobian[4], 3);

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

  962.         // du
  963.         final Vector3D duP = new Vector3D(cv * cosRaan / sinI, dwXP,
  964.                                           cv * sinRaan / sinI, dwYP,
  965.                                           1, new Vector3D(cosRaan, sinRaan, 0));
  966.         final Vector3D duV = new Vector3D(cv * cosRaan / sinI, dwXV,
  967.                                           cv * sinRaan / sinI, dwYV);

  968.         // dv
  969.         final Vector3D dvP = new Vector3D(-u * cosRaan * cosI / sinI + sinRaan * z, dwXP,
  970.                                           -u * sinRaan * cosI / sinI - cosRaan * z, dwYP,
  971.                                           cv, dwZP,
  972.                                           1, new Vector3D(-sinRaan * cosI, cosRaan * cosI, sinI));
  973.         final Vector3D dvV = new Vector3D(-u * cosRaan * cosI / sinI + sinRaan * z, dwXV,
  974.                                           -u * sinRaan * cosI / sinI - cosRaan * z, dwYV,
  975.                                           cv, dwZV);

  976.         final Vector3D dc1P = new Vector3D(aOr2 * (2 * eSinE * eSinE + 1 - eCosE) / r2, position,
  977.                                             -2 * aOr2 * eSinE * oOsqrtMuA, velocity);
  978.         final Vector3D dc1V = new Vector3D(-2 * aOr2 * eSinE * oOsqrtMuA, position,
  979.                                             2 / mu, velocity);
  980.         final Vector3D dc2P = new Vector3D(aOr2 * eSinE * (eSinE * eSinE - (1 - e2)) / (r2 * epsilon), position,
  981.                                             aOr2 * (1 - e2 - eSinE * eSinE) * oOsqrtMuA / epsilon, velocity);
  982.         final Vector3D dc2V = new Vector3D(aOr2 * (1 - e2 - eSinE * eSinE) * oOsqrtMuA / epsilon, position,
  983.                                             eSinE / (mu * epsilon), velocity);

  984.         final double cof1   = aOr2 * (eCosE - e2);
  985.         final double cof2   = aOr2 * epsilon * eSinE;
  986.         final Vector3D dexP = new Vector3D(u, dc1P,  v, dc2P, cof1, duP,  cof2, dvP);
  987.         final Vector3D dexV = new Vector3D(u, dc1V,  v, dc2V, cof1, duV,  cof2, dvV);
  988.         final Vector3D deyP = new Vector3D(v, dc1P, -u, dc2P, cof1, dvP, -cof2, duP);
  989.         final Vector3D deyV = new Vector3D(v, dc1V, -u, dc2V, cof1, dvV, -cof2, duV);
  990.         fillHalfRow(1, dexP, jacobian[1], 0);
  991.         fillHalfRow(1, dexV, jacobian[1], 3);
  992.         fillHalfRow(1, deyP, jacobian[2], 0);
  993.         fillHalfRow(1, deyV, jacobian[2], 3);

  994.         final double cle = u / a + ex - eSinE * beta * ey;
  995.         final double sle = v / a + ey + eSinE * beta * ex;
  996.         final double m1  = beta * eCosE;
  997.         final double m2  = 1 - m1 * eCosE;
  998.         final double m3  = (u * ey - v * ex) + eSinE * beta * (u * ex + v * ey);
  999.         final double m4  = -sle + cle * eSinE * beta;
  1000.         final double m5  = cle + sle * eSinE * beta;
  1001.         fillHalfRow((2 * m3 / r + aOr * eSinE + m1 * eSinE * (1 + m1 - (1 + aOr) * m2) / epsilon) / r2, position,
  1002.                     (m1 * m2 / epsilon - 1) * oOsqrtMuA, velocity,
  1003.                     m4, dexP, m5, deyP, -sle / a, duP, cle / a, dvP,
  1004.                     jacobian[5], 0);
  1005.         fillHalfRow((m1 * m2 / epsilon - 1) * oOsqrtMuA, position,
  1006.                     (2 * m3 + eSinE * a + m1 * eSinE * r * (eCosE * beta * 2 - aOr * m2) / epsilon) / mu, velocity,
  1007.                     m4, dexV, m5, deyV, -sle / a, duV, cle / a, dvV,
  1008.                     jacobian[5], 3);

  1009.         return jacobian;

  1010.     }

  1011.     /** {@inheritDoc} */
  1012.     protected double[][] computeJacobianEccentricWrtCartesian() {

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

  1015.         // Differentiating the Kepler equation aM = aE - ex sin aE + ey cos aE leads to:
  1016.         // daM = (1 - ex cos aE - ey sin aE) daE - sin aE dex + cos aE dey
  1017.         // which is inverted and rewritten as:
  1018.         // daE = a/r daM + sin aE a/r dex - cos aE a/r dey
  1019.         final double alphaE = getAlphaE();
  1020.         final double cosAe  = FastMath.cos(alphaE);
  1021.         final double sinAe  = FastMath.sin(alphaE);
  1022.         final double aOr    = 1 / (1 - ex * cosAe - ey * sinAe);

  1023.         // update longitude row
  1024.         final double[] rowEx = jacobian[1];
  1025.         final double[] rowEy = jacobian[2];
  1026.         final double[] rowL  = jacobian[5];
  1027.         for (int j = 0; j < 6; ++j) {
  1028.             rowL[j] = aOr * (rowL[j] + sinAe * rowEx[j] - cosAe * rowEy[j]);
  1029.         }

  1030.         return jacobian;

  1031.     }

  1032.     /** {@inheritDoc} */
  1033.     protected double[][] computeJacobianTrueWrtCartesian() {

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

  1036.         // Differentiating the eccentric latitude equation
  1037.         // tan((aV - aE)/2) = [ex sin aE - ey cos aE] / [sqrt(1-ex^2-ey^2) + 1 - ex cos aE - ey sin aE]
  1038.         // leads to
  1039.         // cT (daV - daE) = cE daE + cX dex + cY dey
  1040.         // with
  1041.         // cT = [d^2 + (ex sin aE - ey cos aE)^2] / 2
  1042.         // d  = 1 + sqrt(1-ex^2-ey^2) - ex cos aE - ey sin aE
  1043.         // cE = (ex cos aE + ey sin aE) (sqrt(1-ex^2-ey^2) + 1) - ex^2 - ey^2
  1044.         // cX =  sin aE (sqrt(1-ex^2-ey^2) + 1) - ey + ex (ex sin aE - ey cos aE) / sqrt(1-ex^2-ey^2)
  1045.         // cY = -cos aE (sqrt(1-ex^2-ey^2) + 1) + ex + ey (ex sin aE - ey cos aE) / sqrt(1-ex^2-ey^2)
  1046.         // which can be solved to find the differential of the true latitude
  1047.         // daV = (cT + cE) / cT daE + cX / cT deX + cY / cT deX
  1048.         final double alphaE    = getAlphaE();
  1049.         final double cosAe     = FastMath.cos(alphaE);
  1050.         final double sinAe     = FastMath.sin(alphaE);
  1051.         final double eSinE     = ex * sinAe - ey * cosAe;
  1052.         final double ecosE     = ex * cosAe + ey * sinAe;
  1053.         final double e2        = ex * ex + ey * ey;
  1054.         final double epsilon   = FastMath.sqrt(1 - e2);
  1055.         final double onePeps   = 1 + epsilon;
  1056.         final double d         = onePeps - ecosE;
  1057.         final double cT        = (d * d + eSinE * eSinE) / 2;
  1058.         final double cE        = ecosE * onePeps - e2;
  1059.         final double cX        = ex * eSinE / epsilon - ey + sinAe * onePeps;
  1060.         final double cY        = ey * eSinE / epsilon + ex - cosAe * onePeps;
  1061.         final double factorLe  = (cT + cE) / cT;
  1062.         final double factorEx  = cX / cT;
  1063.         final double factorEy  = cY / cT;

  1064.         // update latitude row
  1065.         final double[] rowEx = jacobian[1];
  1066.         final double[] rowEy = jacobian[2];
  1067.         final double[] rowA = jacobian[5];
  1068.         for (int j = 0; j < 6; ++j) {
  1069.             rowA[j] = factorLe * rowA[j] + factorEx * rowEx[j] + factorEy * rowEy[j];
  1070.         }

  1071.         return jacobian;

  1072.     }

  1073.     /** {@inheritDoc} */
  1074.     public void addKeplerContribution(final PositionAngle type, final double gm,
  1075.                                       final double[] pDot) {
  1076.         final double oMe2;
  1077.         final double ksi;
  1078.         final double n = FastMath.sqrt(gm / a) / a;
  1079.         switch (type) {
  1080.             case MEAN :
  1081.                 pDot[5] += n;
  1082.                 break;
  1083.             case ECCENTRIC :
  1084.                 oMe2  = 1 - ex * ex - ey * ey;
  1085.                 ksi   = 1 + ex * FastMath.cos(alphaV) + ey * FastMath.sin(alphaV);
  1086.                 pDot[5] += n * ksi / oMe2;
  1087.                 break;
  1088.             case TRUE :
  1089.                 oMe2  = 1 - ex * ex - ey * ey;
  1090.                 ksi   = 1 + ex * FastMath.cos(alphaV) + ey * FastMath.sin(alphaV);
  1091.                 pDot[5] += n * ksi * ksi / (oMe2 * FastMath.sqrt(oMe2));
  1092.                 break;
  1093.             default :
  1094.                 throw new OrekitInternalError(null);
  1095.         }
  1096.     }

  1097.     /**  Returns a string representation of this Orbit object.
  1098.      * @return a string representation of this object
  1099.      */
  1100.     public String toString() {
  1101.         return new StringBuffer().append("circular parameters: ").append('{').
  1102.                                   append("a: ").append(a).
  1103.                                   append(", ex: ").append(ex).append(", ey: ").append(ey).
  1104.                                   append(", i: ").append(FastMath.toDegrees(i)).
  1105.                                   append(", raan: ").append(FastMath.toDegrees(raan)).
  1106.                                   append(", alphaV: ").append(FastMath.toDegrees(alphaV)).
  1107.                                   append(";}").toString();
  1108.     }

  1109.     /** Replace the instance with a data transfer object for serialization.
  1110.      * @return data transfer object that will be serialized
  1111.      */
  1112.     private Object writeReplace() {
  1113.         return new DTO(this);
  1114.     }

  1115.     /** Internal class used only for serialization. */
  1116.     private static class DTO implements Serializable {

  1117.         /** Serializable UID. */
  1118.         private static final long serialVersionUID = 20170414L;

  1119.         /** Double values. */
  1120.         private double[] d;

  1121.         /** Frame in which are defined the orbital parameters. */
  1122.         private final Frame frame;

  1123.         /** Simple constructor.
  1124.          * @param orbit instance to serialize
  1125.          */
  1126.         private DTO(final CircularOrbit orbit) {

  1127.             final AbsoluteDate date = orbit.getDate();

  1128.             // decompose date
  1129.             final double epoch  = FastMath.floor(date.durationFrom(AbsoluteDate.J2000_EPOCH));
  1130.             final double offset = date.durationFrom(AbsoluteDate.J2000_EPOCH.shiftedBy(epoch));

  1131.             if (orbit.serializePV) {
  1132.                 final TimeStampedPVCoordinates pv = orbit.getPVCoordinates();
  1133.                 if (orbit.hasDerivatives()) {
  1134.                     this.d = new double[] {
  1135.                         // date + mu + orbit + derivatives + Cartesian : 24 parameters
  1136.                         epoch, offset, orbit.getMu(),
  1137.                         orbit.a, orbit.ex, orbit.ey,
  1138.                         orbit.i, orbit.raan, orbit.alphaV,
  1139.                         orbit.aDot, orbit.exDot, orbit.eyDot,
  1140.                         orbit.iDot, orbit.raanDot, orbit.alphaVDot,
  1141.                         pv.getPosition().getX(),     pv.getPosition().getY(),     pv.getPosition().getZ(),
  1142.                         pv.getVelocity().getX(),     pv.getVelocity().getY(),     pv.getVelocity().getZ(),
  1143.                         pv.getAcceleration().getX(), pv.getAcceleration().getY(), pv.getAcceleration().getZ(),
  1144.                     };
  1145.                 } else {
  1146.                     this.d = new double[] {
  1147.                         // date + mu + orbit + Cartesian : 18 parameters
  1148.                         epoch, offset, orbit.getMu(),
  1149.                         orbit.a, orbit.ex, orbit.ey,
  1150.                         orbit.i, orbit.raan, orbit.alphaV,
  1151.                         pv.getPosition().getX(),     pv.getPosition().getY(),     pv.getPosition().getZ(),
  1152.                         pv.getVelocity().getX(),     pv.getVelocity().getY(),     pv.getVelocity().getZ(),
  1153.                         pv.getAcceleration().getX(), pv.getAcceleration().getY(), pv.getAcceleration().getZ(),
  1154.                     };
  1155.                 }
  1156.             } else {
  1157.                 if (orbit.hasDerivatives()) {
  1158.                     // date + mu + orbit + derivatives: 15 parameters
  1159.                     this.d = new double[] {
  1160.                         epoch, offset, orbit.getMu(),
  1161.                         orbit.a, orbit.ex, orbit.ey,
  1162.                         orbit.i, orbit.raan, orbit.alphaV,
  1163.                         orbit.aDot, orbit.exDot, orbit.eyDot,
  1164.                         orbit.iDot, orbit.raanDot, orbit.alphaVDot
  1165.                     };
  1166.                 } else {
  1167.                     // date + mu + orbit: 9 parameters
  1168.                     this.d = new double[] {
  1169.                         epoch, offset, orbit.getMu(),
  1170.                         orbit.a, orbit.ex, orbit.ey,
  1171.                         orbit.i, orbit.raan, orbit.alphaV
  1172.                     };
  1173.                 }
  1174.             }

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

  1176.         }

  1177.         /** Replace the deserialized data transfer object with a {@link CircularOrbit}.
  1178.          * @return replacement {@link CircularOrbit}
  1179.          */
  1180.         private Object readResolve() {
  1181.             switch (d.length) {
  1182.                 case 24 : // date + mu + orbit + derivatives + Cartesian
  1183.                     return new CircularOrbit(d[ 3], d[ 4], d[ 5], d[ 6], d[ 7], d[ 8],
  1184.                                              d[ 9], d[10], d[11], d[12], d[13], d[14],
  1185.                                              new TimeStampedPVCoordinates(AbsoluteDate.J2000_EPOCH.shiftedBy(d[0]).shiftedBy(d[1]),
  1186.                                                                           new Vector3D(d[15], d[16], d[17]),
  1187.                                                                           new Vector3D(d[18], d[19], d[20]),
  1188.                                                                           new Vector3D(d[21], d[22], d[23])),
  1189.                                              frame,
  1190.                                              d[2]);
  1191.                 case 18 : // date + mu + orbit + Cartesian
  1192.                     return new CircularOrbit(d[3], d[4], d[5], d[6], d[7], d[8],
  1193.                                              Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN,
  1194.                                              new TimeStampedPVCoordinates(AbsoluteDate.J2000_EPOCH.shiftedBy(d[0]).shiftedBy(d[1]),
  1195.                                                                           new Vector3D(d[ 9], d[10], d[11]),
  1196.                                                                           new Vector3D(d[12], d[13], d[14]),
  1197.                                                                           new Vector3D(d[15], d[16], d[17])),
  1198.                                              frame,
  1199.                                              d[2]);
  1200.                 case 15 : // date + mu + orbit + derivatives
  1201.                     return new CircularOrbit(d[ 3], d[ 4], d[ 5], d[ 6], d[ 7], d[ 8],
  1202.                                              d[ 9], d[10], d[11], d[12], d[13], d[14],
  1203.                                              PositionAngle.TRUE,
  1204.                                              frame, AbsoluteDate.J2000_EPOCH.shiftedBy(d[0]).shiftedBy(d[1]),
  1205.                                              d[2]);
  1206.                 default : // date + mu + orbit
  1207.                     return new CircularOrbit(d[3], d[4], d[5], d[6], d[7], d[8], PositionAngle.TRUE,
  1208.                                              frame, AbsoluteDate.J2000_EPOCH.shiftedBy(d[0]).shiftedBy(d[1]),
  1209.                                              d[2]);

  1210.             }
  1211.         }

  1212.     }

  1213. }