EquinoctialOrbit.java

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

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

  24. import org.hipparchus.analysis.differentiation.UnivariateDerivative1;
  25. import org.hipparchus.analysis.interpolation.HermiteInterpolator;
  26. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  27. import org.hipparchus.util.FastMath;
  28. import org.hipparchus.util.MathUtils;
  29. import org.hipparchus.util.SinCos;
  30. import org.orekit.annotation.DefaultDataContext;
  31. import org.orekit.data.DataContext;
  32. import org.orekit.errors.OrekitIllegalArgumentException;
  33. import org.orekit.errors.OrekitInternalError;
  34. import org.orekit.errors.OrekitMessages;
  35. import org.orekit.frames.Frame;
  36. import org.orekit.time.AbsoluteDate;
  37. import org.orekit.utils.PVCoordinates;
  38. import org.orekit.utils.TimeStampedPVCoordinates;


  41. /**
  42.  * This class handles equinoctial orbital parameters, which can support both
  43.  * circular and equatorial orbits.
  44.  * <p>
  45.  * The parameters used internally are the equinoctial elements which can be
  46.  * related to Keplerian elements as follows:
  47.  *   <pre>
  48.  *     a
  49.  *     ex = e cos(ω + Ω)
  50.  *     ey = e sin(ω + Ω)
  51.  *     hx = tan(i/2) cos(Ω)
  52.  *     hy = tan(i/2) sin(Ω)
  53.  *     lv = v + ω + Ω
  54.  *   </pre>
  55.  * where ω stands for the Perigee Argument and Ω stands for the
  56.  * Right Ascension of the Ascending Node.
  57.  *
  58.  * <p>
  59.  * The conversion equations from and to Keplerian elements given above hold only
  60.  * when both sides are unambiguously defined, i.e. when orbit is neither equatorial
  61.  * nor circular. When orbit is either equatorial or circular, the equinoctial
  62.  * parameters are still unambiguously defined whereas some Keplerian elements
  63.  * (more precisely ω and Ω) become ambiguous. For this reason, equinoctial
  64.  * parameters are the recommended way to represent orbits.
  65.  * </p>
  66.  * <p>
  67.  * The instance <code>EquinoctialOrbit</code> is guaranteed to be immutable.
  68.  * </p>
  69.  * @see    Orbit
  70.  * @see    KeplerianOrbit
  71.  * @see    CircularOrbit
  72.  * @see    CartesianOrbit
  73.  * @author Mathieu Rom&eacute;ro
  74.  * @author Luc Maisonobe
  75.  * @author Guylaine Prat
  76.  * @author Fabien Maussion
  77.  * @author V&eacute;ronique Pommier-Maurussane
  78.  */
  79. public class EquinoctialOrbit extends Orbit {

  81.     /** Serializable UID. */
  82.     private static final long serialVersionUID = 20170414L;

  84.     /** Semi-major axis (m). */
  85.     private final double a;

  87.     /** First component of the eccentricity vector. */
  88.     private final double ex;

  90.     /** Second component of the eccentricity vector. */
  91.     private final double ey;

  93.     /** First component of the inclination vector. */
  94.     private final double hx;

  96.     /** Second component of the inclination vector. */
  97.     private final double hy;

  99.     /** True longitude argument (rad). */
  100.     private final double lv;

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

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

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

  111.     /** First component of the inclination vector derivative. */
  112.     private final double hxDot;

  114.     /** Second component of the inclination vector derivative. */
  115.     private final double hyDot;

  117.     /** True longitude argument derivative (rad/s). */
  118.     private final double lvDot;

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

  123.     /** Creates a new instance.
  124.      * @param a  semi-major axis (m)
  125.      * @param ex e cos(ω + Ω), first component of eccentricity vector
  126.      * @param ey e sin(ω + Ω), second component of eccentricity vector
  127.      * @param hx tan(i/2) cos(Ω), first component of inclination vector
  128.      * @param hy tan(i/2) sin(Ω), second component of inclination vector
  129.      * @param l  (M or E or v) + ω + Ω, mean, eccentric or true longitude argument (rad)
  130.      * @param type type of longitude argument
  131.      * @param frame the frame in which the parameters are defined
  132.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  133.      * @param date date of the orbital parameters
  134.      * @param mu central attraction coefficient (m³/s²)
  135.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  136.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  137.      */
  138.     public EquinoctialOrbit(final double a, final double ex, final double ey,
  139.                             final double hx, final double hy, final double l,
  140.                             final PositionAngle type,
  141.                             final Frame frame, final AbsoluteDate date, final double mu)
  142.         throws IllegalArgumentException {
  143.         this(a, ex, ey, hx, hy, l,
  144.              Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN,
  145.              type, frame, date, mu);
  146.     }

  148.     /** Creates a new instance.
  149.      * @param a  semi-major axis (m)
  150.      * @param ex e cos(ω + Ω), first component of eccentricity vector
  151.      * @param ey e sin(ω + Ω), second component of eccentricity vector
  152.      * @param hx tan(i/2) cos(Ω), first component of inclination vector
  153.      * @param hy tan(i/2) sin(Ω), second component of inclination vector
  154.      * @param l  (M or E or v) + ω + Ω, mean, eccentric or true longitude argument (rad)
  155.      * @param aDot  semi-major axis derivative (m/s)
  156.      * @param exDot d(e cos(ω + Ω))/dt, first component of eccentricity vector derivative
  157.      * @param eyDot d(e sin(ω + Ω))/dt, second component of eccentricity vector derivative
  158.      * @param hxDot d(tan(i/2) cos(Ω))/dt, first component of inclination vector derivative
  159.      * @param hyDot d(tan(i/2) sin(Ω))/dt, second component of inclination vector derivative
  160.      * @param lDot  d(M or E or v) + ω + Ω)/dr, mean, eccentric or true longitude argument  derivative (rad/s)
  161.      * @param type type of longitude argument
  162.      * @param frame the frame in which the parameters are defined
  163.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  164.      * @param date date of the orbital parameters
  165.      * @param mu central attraction coefficient (m³/s²)
  166.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  167.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  168.      */
  169.     public EquinoctialOrbit(final double a, final double ex, final double ey,
  170.                             final double hx, final double hy, final double l,
  171.                             final double aDot, final double exDot, final double eyDot,
  172.                             final double hxDot, final double hyDot, final double lDot,
  173.                             final PositionAngle type,
  174.                             final Frame frame, final AbsoluteDate date, final double mu)
  175.         throws IllegalArgumentException {
  176.         super(frame, date, mu);
  177.         if (ex * ex + ey * ey >= 1.0) {
  178.             throw new OrekitIllegalArgumentException(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS,
  179.                                                      getClass().getName());
  180.         }
  181.         this.a     = a;
  182.         this.aDot  = aDot;
  183.         this.ex    = ex;
  184.         this.exDot = exDot;
  185.         this.ey    = ey;
  186.         this.eyDot = eyDot;
  187.         this.hx    = hx;
  188.         this.hxDot = hxDot;
  189.         this.hy    = hy;
  190.         this.hyDot = hyDot;

  192.         if (hasDerivatives()) {
  193.             final UnivariateDerivative1 exUD = new UnivariateDerivative1(ex, exDot);
  194.             final UnivariateDerivative1 eyUD = new UnivariateDerivative1(ey, eyDot);
  195.             final UnivariateDerivative1 lUD  = new UnivariateDerivative1(l,  lDot);
  196.             final UnivariateDerivative1 lvUD;
  197.             switch (type) {
  198.                 case MEAN :
  199.                     lvUD = FieldEquinoctialOrbit.eccentricToTrue(FieldEquinoctialOrbit.meanToEccentric(lUD, exUD, eyUD), exUD, eyUD);
  200.                     break;
  201.                 case ECCENTRIC :
  202.                     lvUD = FieldEquinoctialOrbit.eccentricToTrue(lUD, exUD, eyUD);
  203.                     break;
  204.                 case TRUE :
  205.                     lvUD = lUD;
  206.                     break;
  207.                 default : // this should never happen
  208.                     throw new OrekitInternalError(null);
  209.             }
  210.             this.lv    = lvUD.getValue();
  211.             this.lvDot = lvUD.getDerivative(1);
  212.         } else {
  213.             switch (type) {
  214.                 case MEAN :
  215.                     this.lv = eccentricToTrue(meanToEccentric(l, ex, ey), ex, ey);
  216.                     break;
  217.                 case ECCENTRIC :
  218.                     this.lv = eccentricToTrue(l, ex, ey);
  219.                     break;
  220.                 case TRUE :
  221.                     this.lv = l;
  222.                     break;
  223.                 default : // this should never happen
  224.                     throw new OrekitInternalError(null);
  225.             }
  226.             this.lvDot = Double.NaN;
  227.         }

  229.         this.partialPV = null;

  231.     }

  233.     /** Constructor from Cartesian parameters.
  234.      *
  235.      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
  236.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  237.      * use {@code mu} and the position to compute the acceleration, including
  238.      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
  239.      *
  240.      * @param pvCoordinates the position, velocity and acceleration
  241.      * @param frame the frame in which are defined the {@link PVCoordinates}
  242.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  243.      * @param mu central attraction coefficient (m³/s²)
  244.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  245.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  246.      */
  247.     public EquinoctialOrbit(final TimeStampedPVCoordinates pvCoordinates,
  248.                             final Frame frame, final double mu)
  249.         throws IllegalArgumentException {
  250.         super(pvCoordinates, frame, mu);

  252.         //  compute semi-major axis
  253.         final Vector3D pvP   = pvCoordinates.getPosition();
  254.         final Vector3D pvV   = pvCoordinates.getVelocity();
  255.         final Vector3D pvA   = pvCoordinates.getAcceleration();
  256.         final double r2      = pvP.getNormSq();
  257.         final double r       = FastMath.sqrt(r2);
  258.         final double V2      = pvV.getNormSq();
  259.         final double rV2OnMu = r * V2 / mu;

  261.         if (rV2OnMu > 2) {
  262.             throw new OrekitIllegalArgumentException(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS,
  263.                                                      getClass().getName());
  264.         }

  266.         // compute inclination vector
  267.         final Vector3D w = pvCoordinates.getMomentum().normalize();
  268.         final double d = 1.0 / (1 + w.getZ());
  269.         hx = -d * w.getY();
  270.         hy =  d * w.getX();

  272.         // compute true longitude argument
  273.         final double cLv = (pvP.getX() - d * pvP.getZ() * w.getX()) / r;
  274.         final double sLv = (pvP.getY() - d * pvP.getZ() * w.getY()) / r;
  275.         lv = FastMath.atan2(sLv, cLv);


  278.         // compute semi-major axis
  279.         a = r / (2 - rV2OnMu);

  281.         // compute eccentricity vector
  282.         final double eSE = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(mu * a);
  283.         final double eCE = rV2OnMu - 1;
  284.         final double e2  = eCE * eCE + eSE * eSE;
  285.         final double f   = eCE - e2;
  286.         final double g   = FastMath.sqrt(1 - e2) * eSE;
  287.         ex = a * (f * cLv + g * sLv) / r;
  288.         ey = a * (f * sLv - g * cLv) / r;

  290.         partialPV = pvCoordinates;

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

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

  298.             final Vector3D keplerianAcceleration    = new Vector3D(-mu / (r * r2), pvP);
  299.             final Vector3D nonKeplerianAcceleration = pvA.subtract(keplerianAcceleration);
  300.             final double   aX                       = nonKeplerianAcceleration.getX();
  301.             final double   aY                       = nonKeplerianAcceleration.getY();
  302.             final double   aZ                       = nonKeplerianAcceleration.getZ();
  303.             aDot  = jacobian[0][3] * aX + jacobian[0][4] * aY + jacobian[0][5] * aZ;
  304.             exDot = jacobian[1][3] * aX + jacobian[1][4] * aY + jacobian[1][5] * aZ;
  305.             eyDot = jacobian[2][3] * aX + jacobian[2][4] * aY + jacobian[2][5] * aZ;
  306.             hxDot = jacobian[3][3] * aX + jacobian[3][4] * aY + jacobian[3][5] * aZ;
  307.             hyDot = jacobian[4][3] * aX + jacobian[4][4] * aY + jacobian[4][5] * aZ;

  309.             // in order to compute true anomaly derivative, we must compute
  310.             // mean anomaly derivative including Keplerian motion and convert to true anomaly
  311.             final double lMDot = getKeplerianMeanMotion() +
  312.                                  jacobian[5][3] * aX + jacobian[5][4] * aY + jacobian[5][5] * aZ;
  313.             final UnivariateDerivative1 exUD = new UnivariateDerivative1(ex, exDot);
  314.             final UnivariateDerivative1 eyUD = new UnivariateDerivative1(ey, eyDot);
  315.             final UnivariateDerivative1 lMUD = new UnivariateDerivative1(getLM(), lMDot);
  316.             final UnivariateDerivative1 lvUD = FieldEquinoctialOrbit.eccentricToTrue(FieldEquinoctialOrbit.meanToEccentric(lMUD, exUD, eyUD), exUD, eyUD);
  317.             lvDot = lvUD.getDerivative(1);

  319.         } else {
  320.             // acceleration is either almost zero or NaN,
  321.             // we assume acceleration was not known
  322.             // we don't set up derivatives
  323.             aDot  = Double.NaN;
  324.             exDot = Double.NaN;
  325.             eyDot = Double.NaN;
  326.             hxDot = Double.NaN;
  327.             hyDot = Double.NaN;
  328.             lvDot = Double.NaN;
  329.         }

  331.     }

  333.     /** Constructor from Cartesian parameters.
  334.      *
  335.      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
  336.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  337.      * use {@code mu} and the position to compute the acceleration, including
  338.      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
  339.      *
  340.      * @param pvCoordinates the position end velocity
  341.      * @param frame the frame in which are defined the {@link PVCoordinates}
  342.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  343.      * @param date date of the orbital parameters
  344.      * @param mu central attraction coefficient (m³/s²)
  345.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  346.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  347.      */
  348.     public EquinoctialOrbit(final PVCoordinates pvCoordinates, final Frame frame,
  349.                             final AbsoluteDate date, final double mu)
  350.         throws IllegalArgumentException {
  351.         this(new TimeStampedPVCoordinates(date, pvCoordinates), frame, mu);
  352.     }

  354.     /** Constructor from any kind of orbital parameters.
  355.      * @param op orbital parameters to copy
  356.      */
  357.     public EquinoctialOrbit(final Orbit op) {
  358.         super(op.getFrame(), op.getDate(), op.getMu());
  359.         a         = op.getA();
  360.         aDot      = op.getADot();
  361.         ex        = op.getEquinoctialEx();
  362.         exDot     = op.getEquinoctialExDot();
  363.         ey        = op.getEquinoctialEy();
  364.         eyDot     = op.getEquinoctialEyDot();
  365.         hx        = op.getHx();
  366.         hxDot     = op.getHxDot();
  367.         hy        = op.getHy();
  368.         hyDot     = op.getHyDot();
  369.         lv        = op.getLv();
  370.         lvDot     = op.getLvDot();
  371.         partialPV = null;
  372.     }

  374.     /** {@inheritDoc} */
  375.     public OrbitType getType() {
  376.         return OrbitType.EQUINOCTIAL;
  377.     }

  379.     /** {@inheritDoc} */
  380.     public double getA() {
  381.         return a;
  382.     }

  384.     /** {@inheritDoc} */
  385.     public double getADot() {
  386.         return aDot;
  387.     }

  389.     /** {@inheritDoc} */
  390.     public double getEquinoctialEx() {
  391.         return ex;
  392.     }

  394.     /** {@inheritDoc} */
  395.     public double getEquinoctialExDot() {
  396.         return exDot;
  397.     }

  399.     /** {@inheritDoc} */
  400.     public double getEquinoctialEy() {
  401.         return ey;
  402.     }

  404.     /** {@inheritDoc} */
  405.     public double getEquinoctialEyDot() {
  406.         return eyDot;
  407.     }

  409.     /** {@inheritDoc} */
  410.     public double getHx() {
  411.         return hx;
  412.     }

  414.     /** {@inheritDoc} */
  415.     public double getHxDot() {
  416.         return hxDot;
  417.     }

  419.     /** {@inheritDoc} */
  420.     public double getHy() {
  421.         return hy;
  422.     }

  424.     /** {@inheritDoc} */
  425.     public double getHyDot() {
  426.         return hyDot;
  427.     }

  429.     /** {@inheritDoc} */
  430.     public double getLv() {
  431.         return lv;
  432.     }

  434.     /** {@inheritDoc} */
  435.     public double getLvDot() {
  436.         return lvDot;
  437.     }

  439.     /** {@inheritDoc} */
  440.     public double getLE() {
  441.         return trueToEccentric(lv, ex, ey);
  442.     }

  444.     /** {@inheritDoc} */
  445.     public double getLEDot() {
  446.         final UnivariateDerivative1 lVUD = new UnivariateDerivative1(lv, lvDot);
  447.         final UnivariateDerivative1 exUD = new UnivariateDerivative1(ex, exDot);
  448.         final UnivariateDerivative1 eyUD = new UnivariateDerivative1(ey, eyDot);
  449.         final UnivariateDerivative1 lEUD = FieldEquinoctialOrbit.trueToEccentric(lVUD, exUD, eyUD);
  450.         return lEUD.getDerivative(1);
  451.     }

  453.     /** {@inheritDoc} */
  454.     public double getLM() {
  455.         return eccentricToMean(trueToEccentric(lv, ex, ey), ex, ey);
  456.     }

  458.     /** {@inheritDoc} */
  459.     public double getLMDot() {
  460.         final UnivariateDerivative1 lVUD = new UnivariateDerivative1(lv, lvDot);
  461.         final UnivariateDerivative1 exUD = new UnivariateDerivative1(ex, exDot);
  462.         final UnivariateDerivative1 eyUD = new UnivariateDerivative1(ey, eyDot);
  463.         final UnivariateDerivative1 lMUD = FieldEquinoctialOrbit.eccentricToMean(FieldEquinoctialOrbit.trueToEccentric(lVUD, exUD, eyUD), exUD, eyUD);
  464.         return lMUD.getDerivative(1);
  465.     }

  467.     /** Get the longitude argument.
  468.      * @param type type of the angle
  469.      * @return longitude argument (rad)
  470.      */
  471.     public double getL(final PositionAngle type) {
  472.         return (type == PositionAngle.MEAN) ? getLM() :
  473.                                               ((type == PositionAngle.ECCENTRIC) ? getLE() :
  474.                                                                                    getLv());
  475.     }

  477.     /** Get the longitude argument derivative.
  478.      * @param type type of the angle
  479.      * @return longitude argument derivative (rad/s)
  480.      */
  481.     public double getLDot(final PositionAngle type) {
  482.         return (type == PositionAngle.MEAN) ? getLMDot() :
  483.                                               ((type == PositionAngle.ECCENTRIC) ? getLEDot() :
  484.                                                                                    getLvDot());
  485.     }

  487.     /** Computes the true longitude argument from the eccentric longitude argument.
  488.      * @param lE = E + ω + Ω eccentric longitude argument (rad)
  489.      * @param ex first component of the eccentricity vector
  490.      * @param ey second component of the eccentricity vector
  491.      * @return the true longitude argument
  492.      */
  493.     public static double eccentricToTrue(final double lE, final double ex, final double ey) {
  494.         final double epsilon = FastMath.sqrt(1 - ex * ex - ey * ey);
  495.         final SinCos scLE    = FastMath.sinCos(lE);
  496.         final double num     = ex * scLE.sin() - ey * scLE.cos();
  497.         final double den     = epsilon + 1 - ex * scLE.cos() - ey * scLE.sin();
  498.         return lE + 2 * FastMath.atan(num / den);
  499.     }

  501.     /** Computes the eccentric longitude argument from the true longitude argument.
  502.      * @param lv = v + ω + Ω true longitude argument (rad)
  503.      * @param ex first component of the eccentricity vector
  504.      * @param ey second component of the eccentricity vector
  505.      * @return the eccentric longitude argument
  506.      */
  507.     public static double trueToEccentric(final double lv, final double ex, final double ey) {
  508.         final double epsilon = FastMath.sqrt(1 - ex * ex - ey * ey);
  509.         final SinCos scLv    = FastMath.sinCos(lv);
  510.         final double num     = ey * scLv.cos() - ex * scLv.sin();
  511.         final double den     = epsilon + 1 + ex * scLv.cos() + ey * scLv.sin();
  512.         return lv + 2 * FastMath.atan(num / den);
  513.     }

  515.     /** Computes the eccentric longitude argument from the mean longitude argument.
  516.      * @param lM = M + ω + Ω mean longitude argument (rad)
  517.      * @param ex first component of the eccentricity vector
  518.      * @param ey second component of the eccentricity vector
  519.      * @return the eccentric longitude argument
  520.      */
  521.     public static double meanToEccentric(final double lM, final double ex, final double ey) {
  522.         // Generalization of Kepler equation to equinoctial parameters
  523.         // with lE = PA + RAAN + E and
  524.         //      lM = PA + RAAN + M = lE - ex.sin(lE) + ey.cos(lE)
  525.         double lE = lM;
  526.         double shift = 0.0;
  527.         double lEmlM = 0.0;
  528.         SinCos scLE  = FastMath.sinCos(lE);
  529.         int iter = 0;
  530.         do {
  531.             final double f2 = ex * scLE.sin() - ey * scLE.cos();
  532.             final double f1 = 1.0 - ex * scLE.cos() - ey * scLE.sin();
  533.             final double f0 = lEmlM - f2;

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

  538.             lEmlM -= shift;
  539.             lE     = lM + lEmlM;
  540.             scLE   = FastMath.sinCos(lE);

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

  544.         return lE;

  546.     }

  548.     /** Computes the mean longitude argument from the eccentric longitude argument.
  549.      * @param lE = E + ω + Ω mean longitude argument (rad)
  550.      * @param ex first component of the eccentricity vector
  551.      * @param ey second component of the eccentricity vector
  552.      * @return the mean longitude argument
  553.      */
  554.     public static double eccentricToMean(final double lE, final double ex, final double ey) {
  555.         final SinCos scLE = FastMath.sinCos(lE);
  556.         return lE - ex * scLE.sin() + ey * scLE.cos();
  557.     }

  559.     /** {@inheritDoc} */
  560.     public double getE() {
  561.         return FastMath.sqrt(ex * ex + ey * ey);
  562.     }

  564.     /** {@inheritDoc} */
  565.     public double getEDot() {
  566.         return (ex * exDot + ey * eyDot) / FastMath.sqrt(ex * ex + ey * ey);
  567.     }

  569.     /** {@inheritDoc} */
  570.     public double getI() {
  571.         return 2 * FastMath.atan(FastMath.sqrt(hx * hx + hy * hy));
  572.     }

  574.     /** {@inheritDoc} */
  575.     public double getIDot() {
  576.         final double h2 = hx * hx + hy * hy;
  577.         final double h  = FastMath.sqrt(h2);
  578.         return 2 * (hx * hxDot + hy * hyDot) / (h * (1 + h2));
  579.     }

  581.     /** Compute position and velocity but not acceleration.
  582.      */
  583.     private void computePVWithoutA() {

  585.         if (partialPV != null) {
  586.             // already computed
  587.             return;
  588.         }

  590.         // get equinoctial parameters
  591.         final double lE = getLE();

  593.         // inclination-related intermediate parameters
  594.         final double hx2   = hx * hx;
  595.         final double hy2   = hy * hy;
  596.         final double factH = 1. / (1 + hx2 + hy2);

  598.         // reference axes defining the orbital plane
  599.         final double ux = (1 + hx2 - hy2) * factH;
  600.         final double uy =  2 * hx * hy * factH;
  601.         final double uz = -2 * hy * factH;

  603.         final double vx = uy;
  604.         final double vy = (1 - hx2 + hy2) * factH;
  605.         final double vz =  2 * hx * factH;

  607.         // eccentricity-related intermediate parameters
  608.         final double exey = ex * ey;
  609.         final double ex2  = ex * ex;
  610.         final double ey2  = ey * ey;
  611.         final double e2   = ex2 + ey2;
  612.         final double eta  = 1 + FastMath.sqrt(1 - e2);
  613.         final double beta = 1. / eta;

  615.         // eccentric longitude argument
  616.         final SinCos scLe   = FastMath.sinCos(lE);
  617.         final double cLe    = scLe.cos();
  618.         final double sLe    = scLe.sin();
  619.         final double exCeyS = ex * cLe + ey * sLe;

  621.         // coordinates of position and velocity in the orbital plane
  622.         final double x      = a * ((1 - beta * ey2) * cLe + beta * exey * sLe - ex);
  623.         final double y      = a * ((1 - beta * ex2) * sLe + beta * exey * cLe - ey);

  625.         final double factor = FastMath.sqrt(getMu() / a) / (1 - exCeyS);
  626.         final double xdot   = factor * (-sLe + beta * ey * exCeyS);
  627.         final double ydot   = factor * ( cLe - beta * ex * exCeyS);

  629.         final Vector3D position =
  630.                         new Vector3D(x * ux + y * vx, x * uy + y * vy, x * uz + y * vz);
  631.         final Vector3D velocity =
  632.                         new Vector3D(xdot * ux + ydot * vx, xdot * uy + ydot * vy, xdot * uz + ydot * vz);
  633.         partialPV = new PVCoordinates(position, velocity);

  635.     }

  637.     /** Compute non-Keplerian part of the acceleration from first time derivatives.
  638.      * <p>
  639.      * This method should be called only when {@link #hasDerivatives()} returns true.
  640.      * </p>
  641.      * @return non-Keplerian part of the acceleration
  642.      */
  643.     private Vector3D nonKeplerianAcceleration() {

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

  648.         final double nonKeplerianMeanMotion = getLMDot() - getKeplerianMeanMotion();
  649.         final double nonKeplerianAx = dCdP[3][0] * aDot  + dCdP[3][1] * exDot + dCdP[3][2] * eyDot +
  650.                                       dCdP[3][3] * hxDot + dCdP[3][4] * hyDot + dCdP[3][5] * nonKeplerianMeanMotion;
  651.         final double nonKeplerianAy = dCdP[4][0] * aDot  + dCdP[4][1] * exDot + dCdP[4][2] * eyDot +
  652.                                       dCdP[4][3] * hxDot + dCdP[4][4] * hyDot + dCdP[4][5] * nonKeplerianMeanMotion;
  653.         final double nonKeplerianAz = dCdP[5][0] * aDot  + dCdP[5][1] * exDot + dCdP[5][2] * eyDot +
  654.                                       dCdP[5][3] * hxDot + dCdP[5][4] * hyDot + dCdP[5][5] * nonKeplerianMeanMotion;

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

  658.     }

  660.     /** {@inheritDoc} */
  661.     protected TimeStampedPVCoordinates initPVCoordinates() {

  663.         // position and velocity
  664.         computePVWithoutA();

  666.         // acceleration
  667.         final double r2 = partialPV.getPosition().getNormSq();
  668.         final Vector3D keplerianAcceleration = new Vector3D(-getMu() / (r2 * FastMath.sqrt(r2)), partialPV.getPosition());
  669.         final Vector3D acceleration = hasDerivatives() ?
  670.                                       keplerianAcceleration.add(nonKeplerianAcceleration()) :
  671.                                       keplerianAcceleration;

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

  675.     }

  677.     /** {@inheritDoc} */
  678.     public EquinoctialOrbit shiftedBy(final double dt) {

  680.         // use Keplerian-only motion
  681.         final EquinoctialOrbit keplerianShifted = new EquinoctialOrbit(a, ex, ey, hx, hy,
  682.                                                                        getLM() + getKeplerianMeanMotion() * dt,
  683.                                                                        PositionAngle.MEAN, getFrame(),
  684.                                                                        getDate().shiftedBy(dt), getMu());

  686.         if (hasDerivatives()) {

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

  691.             // add quadratic effect of non-Keplerian acceleration to Keplerian-only shift
  692.             keplerianShifted.computePVWithoutA();
  693.             final Vector3D fixedP   = new Vector3D(1, keplerianShifted.partialPV.getPosition(),
  694.                                                    0.5 * dt * dt, nonKeplerianAcceleration);
  695.             final double   fixedR2 = fixedP.getNormSq();
  696.             final double   fixedR  = FastMath.sqrt(fixedR2);
  697.             final Vector3D fixedV  = new Vector3D(1, keplerianShifted.partialPV.getVelocity(),
  698.                                                   dt, nonKeplerianAcceleration);
  699.             final Vector3D fixedA  = new Vector3D(-getMu() / (fixedR2 * fixedR), keplerianShifted.partialPV.getPosition(),
  700.                                                   1, nonKeplerianAcceleration);

  702.             // build a new orbit, taking non-Keplerian acceleration into account
  703.             return new EquinoctialOrbit(new TimeStampedPVCoordinates(keplerianShifted.getDate(),
  704.                                                                      fixedP, fixedV, fixedA),
  705.                                         keplerianShifted.getFrame(), keplerianShifted.getMu());

  707.         } else {
  708.             // Keplerian-only motion is all we can do
  709.             return keplerianShifted;
  710.         }

  712.     }

  714.     /** {@inheritDoc}
  715.      * <p>
  716.      * The interpolated instance is created by polynomial Hermite interpolation
  717.      * on equinoctial elements, without derivatives (which means the interpolation
  718.      * falls back to Lagrange interpolation only).
  719.      * </p>
  720.      * <p>
  721.      * As this implementation of interpolation is polynomial, it should be used only
  722.      * with small samples (about 10-20 points) in order to avoid <a
  723.      * href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
  724.      * and numerical problems (including NaN appearing).
  725.      * </p>
  726.      * <p>
  727.      * If orbit interpolation on large samples is needed, using the {@link
  728.      * org.orekit.propagation.analytical.Ephemeris} class is a better way than using this
  729.      * low-level interpolation. The Ephemeris class automatically handles selection of
  730.      * a neighboring sub-sample with a predefined number of point from a large global sample
  731.      * in a thread-safe way.
  732.      * </p>
  733.      */
  734.     public EquinoctialOrbit interpolate(final AbsoluteDate date, final Stream<Orbit> sample) {

  736.         // first pass to check if derivatives are available throughout the sample
  737.         final List<Orbit> list = sample.collect(Collectors.toList());
  738.         boolean useDerivatives = true;
  739.         for (final Orbit orbit : list) {
  740.             useDerivatives = useDerivatives && orbit.hasDerivatives();
  741.         }

  743.         // set up an interpolator
  744.         final HermiteInterpolator interpolator = new HermiteInterpolator();

  746.         // second pass to feed interpolator
  747.         AbsoluteDate previousDate = null;
  748.         double       previousLm   = Double.NaN;
  749.         for (final Orbit orbit : list) {
  750.             final EquinoctialOrbit equi = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(orbit);
  751.             final double continuousLm;
  752.             if (previousDate == null) {
  753.                 continuousLm = equi.getLM();
  754.             } else {
  755.                 final double dt       = equi.getDate().durationFrom(previousDate);
  756.                 final double keplerLm = previousLm + equi.getKeplerianMeanMotion() * dt;
  757.                 continuousLm = MathUtils.normalizeAngle(equi.getLM(), keplerLm);
  758.             }
  759.             previousDate = equi.getDate();
  760.             previousLm   = continuousLm;
  761.             if (useDerivatives) {
  762.                 interpolator.addSamplePoint(equi.getDate().durationFrom(date),
  763.                                             new double[] {
  764.                                                 equi.getA(),
  765.                                                 equi.getEquinoctialEx(),
  766.                                                 equi.getEquinoctialEy(),
  767.                                                 equi.getHx(),
  768.                                                 equi.getHy(),
  769.                                                 continuousLm
  770.                                             },
  771.                                             new double[] {
  772.                                                 equi.getADot(),
  773.                                                 equi.getEquinoctialExDot(),
  774.                                                 equi.getEquinoctialEyDot(),
  775.                                                 equi.getHxDot(),
  776.                                                 equi.getHyDot(),
  777.                                                 equi.getLMDot()
  778.                                             });
  779.             } else {
  780.                 interpolator.addSamplePoint(equi.getDate().durationFrom(date),
  781.                                             new double[] {
  782.                                                 equi.getA(),
  783.                                                 equi.getEquinoctialEx(),
  784.                                                 equi.getEquinoctialEy(),
  785.                                                 equi.getHx(),
  786.                                                 equi.getHy(),
  787.                                                 continuousLm
  788.                                             });
  789.             }
  790.         }

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

  795.         // build a new interpolated instance
  796.         return new EquinoctialOrbit(interpolated[0][0], interpolated[0][1], interpolated[0][2],
  797.                                     interpolated[0][3], interpolated[0][4], interpolated[0][5],
  798.                                     interpolated[1][0], interpolated[1][1], interpolated[1][2],
  799.                                     interpolated[1][3], interpolated[1][4], interpolated[1][5],
  800.                                     PositionAngle.MEAN, getFrame(), date, getMu());

  802.     }

  804.     /** {@inheritDoc} */
  805.     protected double[][] computeJacobianMeanWrtCartesian() {

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

  809.         // compute various intermediate parameters
  810.         computePVWithoutA();
  811.         final Vector3D position = partialPV.getPosition();
  812.         final Vector3D velocity = partialPV.getVelocity();
  813.         final double r2         = position.getNormSq();
  814.         final double r          = FastMath.sqrt(r2);
  815.         final double r3         = r * r2;

  817.         final double mu         = getMu();
  818.         final double sqrtMuA    = FastMath.sqrt(a * mu);
  819.         final double a2         = a * a;

  821.         final double e2         = ex * ex + ey * ey;
  822.         final double oMe2       = 1 - e2;
  823.         final double epsilon    = FastMath.sqrt(oMe2);
  824.         final double beta       = 1 / (1 + epsilon);
  825.         final double ratio      = epsilon * beta;

  827.         final double hx2        = hx * hx;
  828.         final double hy2        = hy * hy;
  829.         final double hxhy       = hx * hy;

  831.         // precomputing equinoctial frame unit vectors (f, g, w)
  832.         final Vector3D f  = new Vector3D(1 - hy2 + hx2, 2 * hxhy, -2 * hy).normalize();
  833.         final Vector3D g  = new Vector3D(2 * hxhy, 1 + hy2 - hx2, 2 * hx).normalize();
  834.         final Vector3D w  = Vector3D.crossProduct(position, velocity).normalize();

  836.         // coordinates of the spacecraft in the equinoctial frame
  837.         final double x    = Vector3D.dotProduct(position, f);
  838.         final double y    = Vector3D.dotProduct(position, g);
  839.         final double xDot = Vector3D.dotProduct(velocity, f);
  840.         final double yDot = Vector3D.dotProduct(velocity, g);

  842.         // drDot / dEx = dXDot / dEx * f + dYDot / dEx * g
  843.         final double c1 = a / (sqrtMuA * epsilon);
  844.         final double c2 = a * sqrtMuA * beta / r3;
  845.         final double c3 = sqrtMuA / (r3 * epsilon);
  846.         final Vector3D drDotSdEx = new Vector3D( c1 * xDot * yDot - c2 * ey * x - c3 * x * y, f,
  847.                                                 -c1 * xDot * xDot - c2 * ey * y + c3 * x * x, g);

  849.         // drDot / dEy = dXDot / dEy * f + dYDot / dEy * g
  850.         final Vector3D drDotSdEy = new Vector3D( c1 * yDot * yDot + c2 * ex * x - c3 * y * y, f,
  851.                                                 -c1 * xDot * yDot + c2 * ex * y + c3 * x * y, g);

  853.         // da
  854.         final Vector3D vectorAR = new Vector3D(2 * a2 / r3, position);
  855.         final Vector3D vectorARDot = new Vector3D(2 * a2 / mu, velocity);
  856.         fillHalfRow(1, vectorAR,    jacobian[0], 0);
  857.         fillHalfRow(1, vectorARDot, jacobian[0], 3);

  859.         // dEx
  860.         final double d1 = -a * ratio / r3;
  861.         final double d2 = (hy * xDot - hx * yDot) / (sqrtMuA * epsilon);
  862.         final double d3 = (hx * y - hy * x) / sqrtMuA;
  863.         final Vector3D vectorExRDot =
  864.             new Vector3D((2 * x * yDot - xDot * y) / mu, g, -y * yDot / mu, f, -ey * d3 / epsilon, w);
  865.         fillHalfRow(ex * d1, position, -ey * d2, w, epsilon / sqrtMuA, drDotSdEy, jacobian[1], 0);
  866.         fillHalfRow(1, vectorExRDot, jacobian[1], 3);

  868.         // dEy
  869.         final Vector3D vectorEyRDot =
  870.             new Vector3D((2 * xDot * y - x * yDot) / mu, f, -x * xDot / mu, g, ex * d3 / epsilon, w);
  871.         fillHalfRow(ey * d1, position, ex * d2, w, -epsilon / sqrtMuA, drDotSdEx, jacobian[2], 0);
  872.         fillHalfRow(1, vectorEyRDot, jacobian[2], 3);

  874.         // dHx
  875.         final double h = (1 + hx2 + hy2) / (2 * sqrtMuA * epsilon);
  876.         fillHalfRow(-h * xDot, w, jacobian[3], 0);
  877.         fillHalfRow( h * x,    w, jacobian[3], 3);

  879.         // dHy
  880.         fillHalfRow(-h * yDot, w, jacobian[4], 0);
  881.         fillHalfRow( h * y,    w, jacobian[4], 3);

  883.         // dLambdaM
  884.         final double l = -ratio / sqrtMuA;
  885.         fillHalfRow(-1 / sqrtMuA, velocity, d2, w, l * ex, drDotSdEx, l * ey, drDotSdEy, jacobian[5], 0);
  886.         fillHalfRow(-2 / sqrtMuA, position, ex * beta, vectorEyRDot, -ey * beta, vectorExRDot, d3, w, jacobian[5], 3);

  888.         return jacobian;

  890.     }

  892.     /** {@inheritDoc} */
  893.     protected double[][] computeJacobianEccentricWrtCartesian() {

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

  898.         // Differentiating the Kepler equation lM = lE - ex sin lE + ey cos lE leads to:
  899.         // dlM = (1 - ex cos lE - ey sin lE) dE - sin lE dex + cos lE dey
  900.         // which is inverted and rewritten as:
  901.         // dlE = a/r dlM + sin lE a/r dex - cos lE a/r dey
  902.         final SinCos scLe  = FastMath.sinCos(getLE());
  903.         final double cosLe = scLe.cos();
  904.         final double sinLe = scLe.sin();
  905.         final double aOr   = 1 / (1 - ex * cosLe - ey * sinLe);

  907.         // update longitude row
  908.         final double[] rowEx = jacobian[1];
  909.         final double[] rowEy = jacobian[2];
  910.         final double[] rowL  = jacobian[5];
  911.         for (int j = 0; j < 6; ++j) {
  912.             rowL[j] = aOr * (rowL[j] + sinLe * rowEx[j] - cosLe * rowEy[j]);
  913.         }

  915.         return jacobian;

  917.     }

  919.     /** {@inheritDoc} */
  920.     protected double[][] computeJacobianTrueWrtCartesian() {

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

  925.         // Differentiating the eccentric longitude equation
  926.         // tan((lV - lE)/2) = [ex sin lE - ey cos lE] / [sqrt(1-ex^2-ey^2) + 1 - ex cos lE - ey sin lE]
  927.         // leads to
  928.         // cT (dlV - dlE) = cE dlE + cX dex + cY dey
  929.         // with
  930.         // cT = [d^2 + (ex sin lE - ey cos lE)^2] / 2
  931.         // d  = 1 + sqrt(1-ex^2-ey^2) - ex cos lE - ey sin lE
  932.         // cE = (ex cos lE + ey sin lE) (sqrt(1-ex^2-ey^2) + 1) - ex^2 - ey^2
  933.         // cX =  sin lE (sqrt(1-ex^2-ey^2) + 1) - ey + ex (ex sin lE - ey cos lE) / sqrt(1-ex^2-ey^2)
  934.         // cY = -cos lE (sqrt(1-ex^2-ey^2) + 1) + ex + ey (ex sin lE - ey cos lE) / sqrt(1-ex^2-ey^2)
  935.         // which can be solved to find the differential of the true longitude
  936.         // dlV = (cT + cE) / cT dlE + cX / cT deX + cY / cT deX
  937.         final SinCos scLe      = FastMath.sinCos(getLE());
  938.         final double cosLe     = scLe.cos();
  939.         final double sinLe     = scLe.sin();
  940.         final double eSinE     = ex * sinLe - ey * cosLe;
  941.         final double ecosE     = ex * cosLe + ey * sinLe;
  942.         final double e2        = ex * ex + ey * ey;
  943.         final double epsilon   = FastMath.sqrt(1 - e2);
  944.         final double onePeps   = 1 + epsilon;
  945.         final double d         = onePeps - ecosE;
  946.         final double cT        = (d * d + eSinE * eSinE) / 2;
  947.         final double cE        = ecosE * onePeps - e2;
  948.         final double cX        = ex * eSinE / epsilon - ey + sinLe * onePeps;
  949.         final double cY        = ey * eSinE / epsilon + ex - cosLe * onePeps;
  950.         final double factorLe  = (cT + cE) / cT;
  951.         final double factorEx  = cX / cT;
  952.         final double factorEy  = cY / cT;

  954.         // update longitude row
  955.         final double[] rowEx = jacobian[1];
  956.         final double[] rowEy = jacobian[2];
  957.         final double[] rowL = jacobian[5];
  958.         for (int j = 0; j < 6; ++j) {
  959.             rowL[j] = factorLe * rowL[j] + factorEx * rowEx[j] + factorEy * rowEy[j];
  960.         }

  962.         return jacobian;

  964.     }

  966.     /** {@inheritDoc} */
  967.     public void addKeplerContribution(final PositionAngle type, final double gm,
  968.                                       final double[] pDot) {
  969.         final double oMe2;
  970.         final double ksi;
  971.         final double n  = FastMath.sqrt(gm / a) / a;
  972.         final SinCos sc = FastMath.sinCos(lv);
  973.         switch (type) {
  974.             case MEAN :
  975.                 pDot[5] += n;
  976.                 break;
  977.             case ECCENTRIC :
  978.                 oMe2 = 1 - ex * ex - ey * ey;
  979.                 ksi  = 1 + ex * sc.cos() + ey * sc.sin();
  980.                 pDot[5] += n * ksi / oMe2;
  981.                 break;
  982.             case TRUE :
  983.                 oMe2 = 1 - ex * ex - ey * ey;
  984.                 ksi  = 1 + ex * sc.cos() + ey * sc.sin();
  985.                 pDot[5] += n * ksi * ksi / (oMe2 * FastMath.sqrt(oMe2));
  986.                 break;
  987.             default :
  988.                 throw new OrekitInternalError(null);
  989.         }
  990.     }

  992.     /**  Returns a string representation of this equinoctial parameters object.
  993.      * @return a string representation of this object
  994.      */
  995.     public String toString() {
  996.         return new StringBuilder().append("equinoctial parameters: ").append('{').
  997.                                   append("a: ").append(a).
  998.                                   append("; ex: ").append(ex).append("; ey: ").append(ey).
  999.                                   append("; hx: ").append(hx).append("; hy: ").append(hy).
  1000.                                   append("; lv: ").append(FastMath.toDegrees(lv)).
  1001.                                   append(";}").toString();
  1002.     }

  1004.     /** Replace the instance with a data transfer object for serialization.
  1005.      * @return data transfer object that will be serialized
  1006.      */
  1007.     @DefaultDataContext
  1008.     private Object writeReplace() {
  1009.         return new DTO(this);
  1010.     }

  1012.     /** Internal class used only for serialization. */
  1013.     @DefaultDataContext
  1014.     private static class DTO implements Serializable {

  1016.         /** Serializable UID. */
  1017.         private static final long serialVersionUID = 20170414L;

  1019.         /** Double values. */
  1020.         private double[] d;

  1022.         /** Frame in which are defined the orbital parameters. */
  1023.         private final Frame frame;

  1025.         /** Simple constructor.
  1026.          * @param orbit instance to serialize
  1027.          */
  1028.         private DTO(final EquinoctialOrbit orbit) {

  1030.             final TimeStampedPVCoordinates pv = orbit.getPVCoordinates();

  1032.             // decompose date
  1033.             final AbsoluteDate j2000Epoch =
  1034.                     DataContext.getDefault().getTimeScales().getJ2000Epoch();
  1035.             final double epoch  = FastMath.floor(pv.getDate().durationFrom(j2000Epoch));
  1036.             final double offset = pv.getDate().durationFrom(j2000Epoch.shiftedBy(epoch));

  1038.             if (orbit.hasDerivatives()) {
  1039.                 // we have derivatives
  1040.                 this.d = new double[] {
  1041.                     epoch, offset, orbit.getMu(),
  1042.                     orbit.a, orbit.ex, orbit.ey,
  1043.                     orbit.hx, orbit.hy, orbit.lv,
  1044.                     orbit.aDot, orbit.exDot, orbit.eyDot,
  1045.                     orbit.hxDot, orbit.hyDot, orbit.lvDot
  1046.                 };
  1047.             } else {
  1048.                 // we don't have derivatives
  1049.                 this.d = new double[] {
  1050.                     epoch, offset, orbit.getMu(),
  1051.                     orbit.a, orbit.ex, orbit.ey,
  1052.                     orbit.hx, orbit.hy, orbit.lv
  1053.                 };
  1054.             }

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

  1058.         }

  1060.         /** Replace the deserialized data transfer object with a {@link EquinoctialOrbit}.
  1061.          * @return replacement {@link EquinoctialOrbit}
  1062.          */
  1063.         private Object readResolve() {
  1064.             final AbsoluteDate j2000Epoch =
  1065.                     DataContext.getDefault().getTimeScales().getJ2000Epoch();
  1066.             if (d.length >= 15) {
  1067.                 // we have derivatives
  1068.                 return new EquinoctialOrbit(d[ 3], d[ 4], d[ 5], d[ 6], d[ 7], d[ 8],
  1069.                                             d[ 9], d[10], d[11], d[12], d[13], d[14],
  1070.                                             PositionAngle.TRUE,
  1071.                                             frame, j2000Epoch.shiftedBy(d[0]).shiftedBy(d[1]),
  1072.                                             d[2]);
  1073.             } else {
  1074.                 // we don't have derivatives
  1075.                 return new EquinoctialOrbit(d[ 3], d[ 4], d[ 5], d[ 6], d[ 7], d[ 8], PositionAngle.TRUE,
  1076.                                             frame, j2000Epoch.shiftedBy(d[0]).shiftedBy(d[1]),
  1077.                                             d[2]);
  1078.             }
  1079.         }

  1081.     }

  1083. }
Created with JaCoCo 0.8.7.202105040129