FieldCircularOrbit.java

  1. /* Copyright 2002-2018 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;

  19. import java.util.HashMap;
  20. import java.util.List;
  21. import java.util.Map;
  22. import java.util.stream.Collectors;
  23. import java.util.stream.Stream;

  25. import org.hipparchus.Field;
  26. import org.hipparchus.RealFieldElement;
  27. import org.hipparchus.analysis.differentiation.FDSFactory;
  28. import org.hipparchus.analysis.differentiation.FieldDerivativeStructure;
  29. import org.hipparchus.analysis.interpolation.FieldHermiteInterpolator;
  30. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  31. import org.hipparchus.util.FastMath;
  32. import org.hipparchus.util.MathArrays;
  33. import org.orekit.errors.OrekitIllegalArgumentException;
  34. import org.orekit.errors.OrekitInternalError;
  35. import org.orekit.errors.OrekitMessages;
  36. import org.orekit.frames.Frame;
  37. import org.orekit.time.FieldAbsoluteDate;
  38. import org.orekit.utils.FieldPVCoordinates;
  39. import org.orekit.utils.TimeStampedFieldPVCoordinates;


  42. /**
  43.  * This class handles circular orbital parameters.

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

  82. public  class FieldCircularOrbit<T extends RealFieldElement<T>>
  83.     extends FieldOrbit<T> {

  85.     /** Factory for first time derivatives. */
  86.     private static final Map<Field<? extends RealFieldElement<?>>, FDSFactory<? extends RealFieldElement<?>>> FACTORIES =
  87.                     new HashMap<>();

  89.     /** Semi-major axis (m). */
  90.     private final T a;

  92.     /** First component of the circular eccentricity vector. */
  93.     private final T ex;

  95.     /** Second component of the circular eccentricity vector. */
  96.     private final T ey;

  98.     /** Inclination (rad). */
  99.     private final T i;

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

  104.     /** True latitude argument (rad). */
  105.     private final T alphaV;

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

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

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

  116.     /** Inclination derivative (rad/s). */
  117.     private final T iDot;

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

  122.     /** True latitude argument derivative (rad/s). */
  123.     private final T alphaVDot;

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

  128.     /** one. */
  129.     private final T one;

  131.     /** zero. */
  132.     private final T zero;

  134.     /** Creates a new instance.
  135.      * @param a  semi-major axis (m)
  136.      * @param ex e cos(ω), first component of circular eccentricity vector
  137.      * @param ey e sin(ω), second component of circular eccentricity vector
  138.      * @param i inclination (rad)
  139.      * @param raan right ascension of ascending node (Ω, rad)
  140.      * @param alpha  an + ω, mean, eccentric or true latitude argument (rad)
  141.      * @param type type of latitude argument
  142.      * @param frame the frame in which are defined the parameters
  143.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  144.      * @param date date of the orbital parameters
  145.      * @param mu central attraction coefficient (m³/s²)
  146.      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
  147.      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
  148.      */
  149.     public FieldCircularOrbit(final T a, final T ex, final T ey,
  150.                               final T i, final T raan,
  151.                               final T alpha, final PositionAngle type,
  152.                               final Frame frame, final FieldAbsoluteDate<T> date, final double mu)
  153.         throws IllegalArgumentException {
  154.         this(a, ex, ey, i, raan, alpha,
  155.              null, null, null, null, null, null,
  156.              type, frame, date, mu);
  157.     }

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

  193.         if (!FACTORIES.containsKey(a.getField())) {
  194.             FACTORIES.put(a.getField(), new FDSFactory<>(a.getField(), 1, 1));
  195.         }

  197.         this.a       =  a;
  198.         this.aDot    =  aDot;
  199.         this.ex      = ex;
  200.         this.exDot   = exDot;
  201.         this.ey      = ey;
  202.         this.eyDot   = eyDot;
  203.         this.i       = i;
  204.         this.iDot    = iDot;
  205.         this.raan    = raan;
  206.         this.raanDot = raanDot;

  208.         one = a.getField().getOne();
  209.         zero = a.getField().getZero();

  211.         if (hasDerivatives()) {
  212.             @SuppressWarnings("unchecked")
  213.             final FDSFactory<T> factory = (FDSFactory<T>) FACTORIES.get(a.getField());
  214.             final FieldDerivativeStructure<T> exDS    = factory.build(ex,    exDot);
  215.             final FieldDerivativeStructure<T> eyDS    = factory.build(ey,    eyDot);
  216.             final FieldDerivativeStructure<T> alphaDS = factory.build(alpha, alphaDot);
  217.             final FieldDerivativeStructure<T> alphavDS;
  218.             switch (type) {
  219.                 case MEAN :
  220.                     alphavDS = eccentricToTrue(meanToEccentric(alphaDS, exDS, eyDS), exDS, eyDS);
  221.                     break;
  222.                 case ECCENTRIC :
  223.                     alphavDS = eccentricToTrue(alphaDS, exDS, eyDS);
  224.                     break;
  225.                 case TRUE :
  226.                     alphavDS = alphaDS;
  227.                     break;
  228.                 default :
  229.                     throw new OrekitInternalError(null);
  230.             }
  231.             this.alphaV    = alphavDS.getValue();
  232.             this.alphaVDot = alphavDS.getPartialDerivative(1);
  233.         } else {
  234.             switch (type) {
  235.                 case MEAN :
  236.                     this.alphaV = eccentricToTrue(meanToEccentric(alpha, ex, ey), ex, ey);
  237.                     break;
  238.                 case ECCENTRIC :
  239.                     this.alphaV = eccentricToTrue(alpha, ex, ey);
  240.                     break;
  241.                 case TRUE :
  242.                     this.alphaV = alpha;
  243.                     break;
  244.                 default :
  245.                     throw new OrekitInternalError(null);
  246.             }
  247.             this.alphaVDot = null;
  248.         }

  250.         this.partialPV = null;

  252.     }

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

  273.         // compute semi-major axis
  274.         final FieldVector3D<T> pvP = pvCoordinates.getPosition();
  275.         final FieldVector3D<T> pvV = pvCoordinates.getVelocity();
  276.         final FieldVector3D<T> pvA = pvCoordinates.getAcceleration();
  277.         final T r2 = pvP.getNormSq();
  278.         final T r  = r2.sqrt();
  279.         final T V2 = pvV.getNormSq();
  280.         final T rV2OnMu = r.multiply(V2).divide(mu);

  282.         zero = r.getField().getZero();
  283.         one = r.getField().getOne();

  285.         if (rV2OnMu.getReal() > 2) {
  286.             throw new OrekitIllegalArgumentException(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS,
  287.                                                      getClass().getName());
  288.         }

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

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

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

  301.         // 2D-coordinates in the canonical frame
  302.         final T cosRaan = raan.cos();
  303.         final T sinRaan = raan.sin();
  304.         final T cosI    = i.cos();
  305.         final T sinI    = i.sin();
  306.         final T xP      = pvP.getX();
  307.         final T yP      = pvP.getY();
  308.         final T zP      = pvP.getZ();
  309.         final T x2      = (xP.multiply(cosRaan).add(yP .multiply(sinRaan))).divide(a);
  310.         final T y2      = (yP.multiply(cosRaan).subtract(xP.multiply(sinRaan))).multiply(cosI).add(zP.multiply(sinI)).divide(a);

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

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

  328.         partialPV = pvCoordinates;

  330.         if (!FACTORIES.containsKey(a.getField())) {
  331.             FACTORIES.put(a.getField(), new FDSFactory<>(a.getField(), 1, 1));
  332.         }

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

  337.             final T[][] jacobian = MathArrays.buildArray(a.getField(), 6, 6);
  338.             getJacobianWrtCartesian(PositionAngle.MEAN, jacobian);

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

  351.             // in order to compute true anomaly derivative, we must compute
  352.             // mean anomaly derivative including Keplerian motion and convert to true anomaly
  353.             final T alphaMDot = getKeplerianMeanMotion().
  354.                                 add(jacobian[5][3].multiply(aX)).add(jacobian[5][4].multiply(aY)).add(jacobian[5][5].multiply(aZ));
  355.             @SuppressWarnings("unchecked")
  356.             final FDSFactory<T> factory = (FDSFactory<T>) FACTORIES.get(a.getField());
  357.             final FieldDerivativeStructure<T> exDS     = factory.build(ex, exDot);
  358.             final FieldDerivativeStructure<T> eyDS     = factory.build(ey, eyDot);
  359.             final FieldDerivativeStructure<T> alphaMDS = factory.build(getAlphaM(), alphaMDot);
  360.             final FieldDerivativeStructure<T> alphavDS = eccentricToTrue(meanToEccentric(alphaMDS, exDS, eyDS), exDS, eyDS);
  361.             alphaVDot = alphavDS.getPartialDerivative(1);

  363.         } else {
  364.             // acceleration is either almost zero or NaN,
  365.             // we assume acceleration was not known
  366.             // we don't set up derivatives
  367.             aDot      = null;
  368.             exDot     = null;
  369.             eyDot     = null;
  370.             iDot      = null;
  371.             raanDot   = null;
  372.             alphaVDot = null;
  373.         }

  375.     }

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

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

  418.         if (!FACTORIES.containsKey(a.getField())) {
  419.             FACTORIES.put(a.getField(), new FDSFactory<>(a.getField(), 1, 1));
  420.         }

  422.         if (op.hasDerivatives()) {
  423.             aDot      = op.getADot();
  424.             final T      hxDot = op.getHxDot();
  425.             final T      hyDot = op.getHyDot();
  426.             iDot      = cosRaan.multiply(hxDot).add(sinRaan.multiply(hyDot)).multiply(2).divide(h2.add(1));
  427.             raanDot   = hx.multiply(hyDot).subtract(hy.multiply(hxDot)).divide(h2);
  428.             final T equiExDot = op.getEquinoctialExDot();
  429.             final T equiEyDot = op.getEquinoctialEyDot();
  430.             exDot     = equiExDot.add(equiEy.multiply(raanDot)).multiply(cosRaan).
  431.                         add(equiEyDot.subtract(equiEx.multiply(raanDot)).multiply(sinRaan));
  432.             eyDot     = equiEyDot.subtract(equiEx.multiply(raanDot)).multiply(cosRaan).
  433.                         subtract(equiExDot.add(equiEy.multiply(raanDot)).multiply(sinRaan));
  434.             alphaVDot = op.getLvDot().subtract(raanDot);
  435.         } else {
  436.             aDot      = null;
  437.             exDot     = null;
  438.             eyDot     = null;
  439.             iDot      = null;
  440.             raanDot   = null;
  441.             alphaVDot = null;
  442.         }

  444.         partialPV = null;

  446.         one = a.getField().getOne();
  447.         zero = a.getField().getZero();

  449.     }

  451.     /** {@inheritDoc} */
  452.     public OrbitType getType() {
  453.         return OrbitType.CIRCULAR;
  454.     }

  456.     /** {@inheritDoc} */
  457.     public T getA() {
  458.         return a;
  459.     }

  461.     /** {@inheritDoc} */
  462.     public T getADot() {
  463.         return aDot;
  464.     }

  466.     /** {@inheritDoc} */
  467.     public T getEquinoctialEx() {
  468.         return ex.multiply(raan.cos()).subtract(ey.multiply(raan.sin()));
  469.     }

  471.     /** {@inheritDoc} */
  472.     public T getEquinoctialExDot() {

  474.         if (!hasDerivatives()) {
  475.             return null;
  476.         }

  478.         final T cosRaan = raan.cos();
  479.         final T sinRaan = raan.sin();
  480.         return exDot.subtract(ey.multiply(raanDot)).multiply(cosRaan).
  481.                subtract(eyDot.add(ex.multiply(raanDot)).multiply(sinRaan));

  483.     }

  485.     /** {@inheritDoc} */
  486.     public T getEquinoctialEy() {
  487.         return ey.multiply(raan.cos()).add(ex.multiply(raan.sin()));
  488.     }

  490.     /** {@inheritDoc} */
  491.     public T getEquinoctialEyDot() {

  493.         if (!hasDerivatives()) {
  494.             return null;
  495.         }

  497.         final T cosRaan = raan.cos();
  498.         final T sinRaan = raan.sin();
  499.         return eyDot.add(ex.multiply(raanDot)).multiply(cosRaan).
  500.                add(exDot.subtract(ey.multiply(raanDot)).multiply(sinRaan));

  502.     }

  504.     /** Get the first component of the circular eccentricity vector.
  505.      * @return ex = e cos(ω), first component of the circular eccentricity vector
  506.      */
  507.     public T getCircularEx() {
  508.         return ex;
  509.     }

  511.     /** Get the first component of the circular eccentricity vector derivative.
  512.      * @return d(ex)/dt = d(e cos(ω))/dt, first component of the circular eccentricity vector derivative
  513.      */
  514.     public T getCircularExDot() {
  515.         return exDot;
  516.     }

  518.     /** Get the second component of the circular eccentricity vector.
  519.      * @return ey = e sin(ω), second component of the circular eccentricity vector
  520.      */
  521.     public T getCircularEy() {
  522.         return ey;
  523.     }

  525.     /** Get the second component of the circular eccentricity vector derivative.
  526.      * @return d(ey)/dt = d(e sin(ω))/dt, second component of the circular eccentricity vector derivative
  527.      */
  528.     public T getCircularEyDot() {
  529.         return eyDot;
  530.     }

  532.     /** {@inheritDoc} */
  533.     public T getHx() {
  534.         // Check for equatorial retrograde orbit
  535.         if (FastMath.abs(i.getReal() - FastMath.PI) < 1.0e-10) {
  536.             return zero.add(Double.NaN);
  537.         }
  538.         return  raan.cos().multiply(i.divide(2).tan());
  539.     }

  541.     /** {@inheritDoc} */
  542.     public T getHxDot() {

  544.         if (!hasDerivatives()) {
  545.             return null;
  546.         }

  548.         // Check for equatorial retrograde orbit
  549.         if (FastMath.abs(i.getReal() - FastMath.PI) < 1.0e-10) {
  550.             return zero.add(Double.NaN);
  551.         }

  553.         final T cosRaan = raan.cos();
  554.         final T sinRaan = raan.sin();
  555.         final T tan     = i.multiply(0.5).tan();
  556.         return cosRaan.multiply(0.5).multiply(tan.multiply(tan).add(1)).multiply(iDot).
  557.                subtract(sinRaan.multiply(tan).multiply(raanDot));

  559.     }

  561.     /** {@inheritDoc} */
  562.     public T getHy() {
  563.         // Check for equatorial retrograde orbit
  564.         if (FastMath.abs(i.getReal() - FastMath.PI) < 1.0e-10) {
  565.             return zero.add(Double.NaN);
  566.         }
  567.         return raan.sin().multiply(i.divide(2).tan());
  568.     }

  570.     /** {@inheritDoc} */
  571.     public T getHyDot() {

  573.         if (!hasDerivatives()) {
  574.             return null;
  575.         }

  577.         // Check for equatorial retrograde orbit
  578.         if (FastMath.abs(i.getReal() - FastMath.PI) < 1.0e-10) {
  579.             return zero.add(Double.NaN);
  580.         }

  582.         final T cosRaan = raan.cos();
  583.         final T sinRaan = raan.sin();
  584.         final T tan     = i.multiply(0.5).tan();
  585.         return sinRaan.multiply(0.5).multiply(tan.multiply(tan).add(1)).multiply(iDot).
  586.                add(cosRaan.multiply(tan).multiply(raanDot));

  588.     }

  590.     /** Get the true latitude argument.
  591.      * @return v + ω true latitude argument (rad)
  592.      */
  593.     public T getAlphaV() {
  594.         return alphaV;
  595.     }

  597.     /** Get the true latitude argument derivative.
  598.      * @return d(v + ω)/dt true latitude argument derivative (rad/s)
  599.      */
  600.     public T getAlphaVDot() {
  601.         return alphaVDot;
  602.     }

  604.     /** Get the eccentric latitude argument.
  605.      * @return E + ω eccentric latitude argument (rad)
  606.      */
  607.     public T getAlphaE() {
  608.         return trueToEccentric(alphaV, ex, ey);
  609.     }

  611.     /** Get the eccentric latitude argument derivative.
  612.      * @return d(E + ω)/dt eccentric latitude argument derivative (rad/s)
  613.      */
  614.     public T getAlphaEDot() {

  616.         if (!hasDerivatives()) {
  617.             return null;
  618.         }

  620.         @SuppressWarnings("unchecked")
  621.         final FDSFactory<T> factory = (FDSFactory<T>) FACTORIES.get(a.getField());
  622.         final FieldDerivativeStructure<T> alphaVDS = factory.build(alphaV, alphaVDot);
  623.         final FieldDerivativeStructure<T> exDS     = factory.build(ex, exDot);
  624.         final FieldDerivativeStructure<T> eyDS     = factory.build(ey, eyDot);
  625.         final FieldDerivativeStructure<T> alphaEDS = trueToEccentric(alphaVDS, exDS, eyDS);
  626.         return alphaEDS.getPartialDerivative(1);

  628.     }

  630.     /** Get the mean latitude argument.
  631.      * @return M + ω mean latitude argument (rad)
  632.      */
  633.     public T getAlphaM() {
  634.         return eccentricToMean(trueToEccentric(alphaV, ex, ey), ex, ey);
  635.     }

  637.     /** Get the mean latitude argument derivative.
  638.      * @return d(M + ω)/dt mean latitude argument derivative (rad/s)
  639.      */
  640.     public T getAlphaMDot() {

  642.         if (!hasDerivatives()) {
  643.             return null;
  644.         }

  646.         @SuppressWarnings("unchecked")
  647.         final FDSFactory<T> factory = (FDSFactory<T>) FACTORIES.get(a.getField());
  648.         final FieldDerivativeStructure<T> alphaVDS = factory.build(alphaV, alphaVDot);
  649.         final FieldDerivativeStructure<T> exDS     = factory.build(ex, exDot);
  650.         final FieldDerivativeStructure<T> eyDS     = factory.build(ey, eyDot);
  651.         final FieldDerivativeStructure<T> alphaMDS = eccentricToMean(trueToEccentric(alphaVDS, exDS, eyDS), exDS, eyDS);
  652.         return alphaMDS.getPartialDerivative(1);

  654.     }

  656.     /** Get the latitude argument.
  657.      * @param type type of the angle
  658.      * @return latitude argument (rad)
  659.      */
  660.     public T getAlpha(final PositionAngle type) {
  661.         return (type == PositionAngle.MEAN) ? getAlphaM() :
  662.                                               ((type == PositionAngle.ECCENTRIC) ? getAlphaE() :
  663.                                                                                    getAlphaV());
  664.     }

  666.     /** Get the latitude argument derivative.
  667.      * @param type type of the angle
  668.      * @return latitude argument derivative (rad/s)
  669.      */
  670.     public T getAlphaDot(final PositionAngle type) {
  671.         return (type == PositionAngle.MEAN) ? getAlphaMDot() :
  672.                                               ((type == PositionAngle.ECCENTRIC) ? getAlphaEDot() :
  673.                                                                                    getAlphaVDot());
  674.     }

  676.     /** Computes the true latitude argument from the eccentric latitude argument.
  677.      * @param alphaE = E + ω eccentric latitude argument (rad)
  678.      * @param ex e cos(ω), first component of circular eccentricity vector
  679.      * @param ey e sin(ω), second component of circular eccentricity vector
  680.      * @param <T> Type of the field elements
  681.      * @return the true latitude argument.
  682.      */
  683.     public static <T extends RealFieldElement<T>> T eccentricToTrue(final T alphaE, final T ex, final T ey) {
  684.         final T epsilon   = ex.multiply(ex).add(ey.multiply(ey)).negate().add(1).sqrt();
  685.         final T cosAlphaE = alphaE.cos();
  686.         final T sinAlphaE = alphaE.sin();
  687.         return alphaE.add(ex.multiply(sinAlphaE).subtract(ey.multiply(cosAlphaE)).divide(
  688.                                       epsilon.add(1).subtract(ex.multiply(cosAlphaE)).subtract(
  689.                                       ey.multiply(sinAlphaE))).atan().multiply(2));
  690.     }

  692.     /** Computes the eccentric latitude argument from the true latitude argument.
  693.      * @param alphaV = v + ω true latitude argument (rad)
  694.      * @param ex e cos(ω), first component of circular eccentricity vector
  695.      * @param ey e sin(ω), second component of circular eccentricity vector
  696.      * @param <T> Type of the field elements
  697.      * @return the eccentric latitude argument.
  698.      */
  699.     public static <T extends RealFieldElement<T>> T trueToEccentric(final T alphaV, final T ex, final T ey) {
  700.         final T epsilon   = ex.multiply(ex).add(ey.multiply(ey)).negate().add(1).sqrt();
  701.         final T cosAlphaV = alphaV.cos();
  702.         final T sinAlphaV = alphaV.sin();
  703.         return alphaV.add(ey.multiply(cosAlphaV).subtract(ex.multiply(sinAlphaV)).divide
  704.                                       (epsilon.add(1).add(ex.multiply(cosAlphaV).add(ey.multiply(sinAlphaV)))).atan().multiply(2));
  705.     }

  707.     /** Computes the eccentric latitude argument from the mean latitude argument.
  708.      * @param alphaM = M + ω  mean latitude argument (rad)
  709.      * @param ex e cos(ω), first component of circular eccentricity vector
  710.      * @param ey e sin(ω), second component of circular eccentricity vector
  711.      * @param <T> Type of the field elements
  712.      * @return the eccentric latitude argument.
  713.      */
  714.     public static <T extends RealFieldElement<T>> T meanToEccentric(final T alphaM, final T ex, final T ey) {
  715.         // Generalization of Kepler equation to circular parameters
  716.         // with alphaE = PA + E and
  717.         //      alphaM = PA + M = alphaE - ex.sin(alphaE) + ey.cos(alphaE)

  719.         T alphaE        = alphaM;
  720.         T shift         = alphaM.getField().getZero();
  721.         T alphaEMalphaM = alphaM.getField().getZero();
  722.         T cosAlphaE     = alphaE.cos();
  723.         T sinAlphaE     = alphaE.sin();
  724.         int    iter     = 0;
  725.         do {
  726.             final T f2 = ex.multiply(sinAlphaE).subtract(ey.multiply(cosAlphaE));
  727.             final T f1 = ex.negate().multiply(cosAlphaE).subtract(ey.multiply(sinAlphaE)).add(1);
  728.             final T f0 = alphaEMalphaM.subtract(f2);

  730.             final T f12 = f1.multiply(2);
  731.             shift = f0.multiply(f12).divide(f1.multiply(f12).subtract(f0.multiply(f2)));

  733.             alphaEMalphaM  = alphaEMalphaM.subtract(shift);
  734.             alphaE         = alphaM.add(alphaEMalphaM);
  735.             cosAlphaE      = alphaE.cos();
  736.             sinAlphaE      = alphaE.sin();
  737.         } while ((++iter < 50) && (FastMath.abs(shift.getReal()) > 1.0e-12));
  738.         return alphaE;

  740.     }

  742.     /** Computes the mean latitude argument from the eccentric latitude argument.
  743.      * @param alphaE = E + ω  eccentric latitude argument (rad)
  744.      * @param ex e cos(ω), first component of circular eccentricity vector
  745.      * @param ey e sin(ω), second component of circular eccentricity vector
  746.      * @param <T> Type of the field elements
  747.      * @return the mean latitude argument.
  748.      */
  749.     public static <T extends RealFieldElement<T>> T eccentricToMean(final T alphaE, final T ex, final T ey) {
  750.         return alphaE.subtract(ex.multiply(alphaE.sin()).subtract(ey.multiply(alphaE.cos())));
  751.     }

  753.     /** Compute position from circular parameters.
  754.      * @param a  semi-major axis (m)
  755.      * @param ex e cos(ω), first component of circular eccentricity vector
  756.      * @param ey e sin(ω), second component of circular eccentricity vector
  757.      * @param i inclination (rad)
  758.      * @param raan right ascension of ascending node (Ω, rad)
  759.      * @param alphaV  v + ω true latitude argument (rad)
  760.      * @param mu central attraction coefficient (m³/s²)
  761.      * @param <T> type of the fiels elements
  762.      * @return position vector
  763.      */
  764.     public static <T extends RealFieldElement<T>> FieldVector3D<T> circularToPosition(final T a, final T ex, final T ey,
  765.                                                                                       final T i, final T raan, final T alphaV,
  766.                                                                                       final double mu) {

  768.         final T zero = a.getField().getZero();

  770.         // get equinoctial parameters
  771.         final T equEx = ex.multiply(raan.cos()).subtract(ey.multiply(raan.sin()));
  772.         final T equEy = ey.multiply(raan.cos()).add(ex.multiply(raan.sin()));
  773.         final T hx;
  774.         final T hy;
  775.         if (FastMath.abs(i.getReal() - FastMath.PI) < 1.0e-10) {
  776.             hx = zero.add(Double.NaN);
  777.             hy = zero.add(Double.NaN);
  778.         } else {
  779.             final T tan = i.divide(2).tan();
  780.             hx = raan.cos().multiply(tan);
  781.             hy = raan.sin().multiply(tan);
  782.         }
  783.         final T lE = trueToEccentric(alphaV, ex, ey).add(raan);

  785.         // inclination-related intermediate parameters
  786.         final T hx2   = hx.multiply(hx);
  787.         final T hy2   = hy.multiply(hy);
  788.         final T factH = (hx2.add(1).add(hy2)).reciprocal();

  790.         // reference axes defining the orbital plane
  791.         final T ux = (hx2.add(1).subtract(hy2)).multiply(factH);
  792.         final T uy =  hx.multiply(2).multiply(hy).multiply(factH);
  793.         final T uz = hy.multiply(-2).multiply(factH);

  795.         final T vx = uy;
  796.         final T vy = (hy2.subtract(hx2).add(1)).multiply(factH);
  797.         final T vz =  hx.multiply(factH).multiply(2);

  799.         // eccentricity-related intermediate parameters
  800.         final T exey = equEx.multiply(equEy);
  801.         final T ex2  = equEx.multiply(equEx);
  802.         final T ey2  = equEy.multiply(equEy);
  803.         final T e2   = ex2.add(ey2);
  804.         final T eta  = e2.negate().add(1).sqrt().add(1);
  805.         final T beta = eta.reciprocal();

  807.         // eccentric latitude argument
  808.         final T cLe    = lE.cos();
  809.         final T sLe    = lE.sin();

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

  815.         return new FieldVector3D<>(x.multiply(ux).add(y.multiply(vx)),
  816.                                    x.multiply(uy).add(y.multiply(vy)),
  817.                                    x.multiply(uz).add(y.multiply(vz)));

  819.     }

  821.     /** {@inheritDoc} */
  822.     public T getE() {
  823.         return ex.multiply(ex).add(ey.multiply(ey)).sqrt();
  824.     }

  826.     /** {@inheritDoc} */
  827.     public T getEDot() {

  829.         if (!hasDerivatives()) {
  830.             return null;
  831.         }

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

  835.     }

  837.     /** {@inheritDoc} */
  838.     public T getI() {
  839.         return i;
  840.     }

  842.     /** {@inheritDoc} */
  843.     public T getIDot() {
  844.         return iDot;
  845.     }

  847.     /** Get the right ascension of the ascending node.
  848.      * @return right ascension of the ascending node (rad)
  849.      */
  850.     public T getRightAscensionOfAscendingNode() {
  851.         return raan;
  852.     }

  854.     /** Get the right ascension of the ascending node derivative.
  855.      * @return right ascension of the ascending node derivative (rad/s)
  856.      */
  857.     public T getRightAscensionOfAscendingNodeDot() {
  858.         return raanDot;
  859.     }

  861.     /** {@inheritDoc} */
  862.     public T getLv() {
  863.         return alphaV.add(raan);
  864.     }

  866.     /** {@inheritDoc} */
  867.     public T getLvDot() {
  868.         return hasDerivatives() ? alphaVDot.add(raanDot) : null;
  869.     }

  871.     /** {@inheritDoc} */
  872.     public T getLE() {
  873.         return getAlphaE().add(raan);
  874.     }

  876.     /** {@inheritDoc} */
  877.     public T getLEDot() {
  878.         return hasDerivatives() ? getAlphaEDot().add(raanDot) : null;
  879.     }

  881.     /** {@inheritDoc} */
  882.     public T getLM() {
  883.         return getAlphaM().add(raan);
  884.     }

  886.     /** {@inheritDoc} */
  887.     public T getLMDot() {
  888.         return hasDerivatives() ? getAlphaMDot().add(raanDot) : null;
  889.     }

  891.     /** {@inheritDoc} */
  892.     @Override
  893.     public boolean hasDerivatives() {
  894.         return aDot != null;
  895.     }

  897.     /** Compute position and velocity but not acceleration.
  898.      */
  899.     private void computePVWithoutA() {

  901.         if (partialPV != null) {
  902.             // already computed
  903.             return;
  904.         }

  906.         // get equinoctial parameters
  907.         final T equEx = getEquinoctialEx();
  908.         final T equEy = getEquinoctialEy();
  909.         final T hx = getHx();
  910.         final T hy = getHy();
  911.         final T lE = getLE();
  912.         // inclination-related intermediate parameters
  913.         final T hx2   = hx.multiply(hx);
  914.         final T hy2   = hy.multiply(hy);
  915.         final T factH = (hx2.add(1).add(hy2)).reciprocal();

  917.         // reference axes defining the orbital plane
  918.         final T ux = (hx2.add(1).subtract(hy2)).multiply(factH);
  919.         final T uy =  hx.multiply(2).multiply(hy).multiply(factH);
  920.         final T uz = hy.multiply(-2).multiply(factH);

  922.         final T vx = uy;
  923.         final T vy = (hy2.subtract(hx2).add(1)).multiply(factH);
  924.         final T vz =  hx.multiply(factH).multiply(2);

  926.         // eccentricity-related intermediate parameters
  927.         final T exey = equEx.multiply(equEy);
  928.         final T ex2  = equEx.multiply(equEx);
  929.         final T ey2  = equEy.multiply(equEy);
  930.         final T e2   = ex2.add(ey2);
  931.         final T eta  = e2.negate().add(1).sqrt().add(1);
  932.         final T beta = eta.reciprocal();

  934.         // eccentric latitude argument
  935.         final T cLe    = lE.cos();
  936.         final T sLe    = lE.sin();
  937.         final T exCeyS = equEx.multiply(cLe).add(equEy.multiply(sLe));
  938.         // coordinates of position and velocity in the orbital plane
  939.         final T x      = a.multiply(beta.negate().multiply(ey2).add(1).multiply(cLe).add(beta.multiply(exey).multiply(sLe)).subtract(equEx));
  940.         final T y      = a.multiply(beta.negate().multiply(ex2).add(1).multiply(sLe).add(beta.multiply(exey).multiply(cLe)).subtract(equEy));

  942.         final T factor = one.add(getMu()).divide(a).sqrt().divide(exCeyS.negate().add(1));
  943.         final T xdot   = factor.multiply( beta.multiply(equEy).multiply(exCeyS).subtract(sLe ));
  944.         final T ydot   = factor.multiply( cLe.subtract(beta.multiply(equEx).multiply(exCeyS)));

  946.         final FieldVector3D<T> position = new FieldVector3D<>(x.multiply(ux).add(y.multiply(vx)),
  947.                                                               x.multiply(uy).add(y.multiply(vy)),
  948.                                                               x.multiply(uz).add(y.multiply(vz)));
  949.         final FieldVector3D<T> velocity = new FieldVector3D<>(xdot.multiply(ux).add(ydot.multiply(vx)),
  950.                                                               xdot.multiply(uy).add(ydot.multiply(vy)),
  951.                                                               xdot.multiply(uz).add(ydot.multiply(vz)));

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

  955.     }

  957.     /** Compute non-Keplerian part of the acceleration from first time derivatives.
  958.      * <p>
  959.      * This method should be called only when {@link #hasDerivatives()} returns true.
  960.      * </p>
  961.      * @return non-Keplerian part of the acceleration
  962.      */
  963.     private FieldVector3D<T> nonKeplerianAcceleration() {

  965.         final T[][] dCdP = MathArrays.buildArray(a.getField(), 6, 6);
  966.         getJacobianWrtParameters(PositionAngle.MEAN, dCdP);

  968.         final T nonKeplerianMeanMotion = getAlphaMDot().subtract(getKeplerianMeanMotion());
  969.         final T nonKeplerianAx =     dCdP[3][0].multiply(aDot).
  970.                                  add(dCdP[3][1].multiply(exDot)).
  971.                                  add(dCdP[3][2].multiply(eyDot)).
  972.                                  add(dCdP[3][3].multiply(iDot)).
  973.                                  add(dCdP[3][4].multiply(raanDot)).
  974.                                  add(dCdP[3][5].multiply(nonKeplerianMeanMotion));
  975.         final T nonKeplerianAy =     dCdP[4][0].multiply(aDot).
  976.                                  add(dCdP[4][1].multiply(exDot)).
  977.                                  add(dCdP[4][2].multiply(eyDot)).
  978.                                  add(dCdP[4][3].multiply(iDot)).
  979.                                  add(dCdP[4][4].multiply(raanDot)).
  980.                                  add(dCdP[4][5].multiply(nonKeplerianMeanMotion));
  981.         final T nonKeplerianAz =     dCdP[5][0].multiply(aDot).
  982.                                  add(dCdP[5][1].multiply(exDot)).
  983.                                  add(dCdP[5][2].multiply(eyDot)).
  984.                                  add(dCdP[5][3].multiply(iDot)).
  985.                                  add(dCdP[5][4].multiply(raanDot)).
  986.                                  add(dCdP[5][5].multiply(nonKeplerianMeanMotion));

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

  990.     }

  992.     /** {@inheritDoc} */
  993.     protected TimeStampedFieldPVCoordinates<T> initPVCoordinates() {

  995.         // position and velocity
  996.         computePVWithoutA();

  998.         // acceleration
  999.         final T r2 = partialPV.getPosition().getNormSq();
  1000.         final FieldVector3D<T> keplerianAcceleration = new FieldVector3D<>(r2.multiply(r2.sqrt()).reciprocal().multiply(-getMu()),
  1001.                                                                            partialPV.getPosition());
  1002.         final FieldVector3D<T> acceleration = hasDerivatives() ?
  1003.                                               keplerianAcceleration.add(nonKeplerianAcceleration()) :
  1004.                                               keplerianAcceleration;

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

  1008.     }

  1010.     /** {@inheritDoc} */
  1011.     public FieldCircularOrbit<T> shiftedBy(final double dt) {
  1012.         return shiftedBy(getDate().getField().getZero().add(dt));
  1013.     }

  1015.     /** {@inheritDoc} */
  1016.     public FieldCircularOrbit<T> shiftedBy(final T dt) {

  1018.         // use Keplerian-only motion
  1019.         final FieldCircularOrbit<T> keplerianShifted = new FieldCircularOrbit<>(a, ex, ey, i, raan,
  1020.                                                                                 getAlphaM().add(getKeplerianMeanMotion().multiply(dt)),
  1021.                                                                                 PositionAngle.MEAN, getFrame(),
  1022.                                                                                 getDate().shiftedBy(dt), getMu());

  1024.         if (hasDerivatives()) {

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

  1029.             // add quadratic effect of non-Keplerian acceleration to Keplerian-only shift
  1030.             keplerianShifted.computePVWithoutA();
  1031.             final FieldVector3D<T> fixedP   = new FieldVector3D<>(one, keplerianShifted.partialPV.getPosition(),
  1032.                                                                   dt.multiply(dt).multiply(0.5), nonKeplerianAcceleration);
  1033.             final T   fixedR2 = fixedP.getNormSq();
  1034.             final T   fixedR  = fixedR2.sqrt();
  1035.             final FieldVector3D<T> fixedV  = new FieldVector3D<>(one, keplerianShifted.partialPV.getVelocity(),
  1036.                                                                  dt, nonKeplerianAcceleration);
  1037.             final FieldVector3D<T> fixedA  = new FieldVector3D<>(fixedR2.multiply(fixedR).reciprocal().multiply(-getMu()),
  1038.                                                                  keplerianShifted.partialPV.getPosition(),
  1039.                                                                  one, nonKeplerianAcceleration);

  1041.             // build a new orbit, taking non-Keplerian acceleration into account
  1042.             return new FieldCircularOrbit<>(new TimeStampedFieldPVCoordinates<>(keplerianShifted.getDate(),
  1043.                                                                                 fixedP, fixedV, fixedA),
  1044.                                             keplerianShifted.getFrame(), keplerianShifted.getMu());

  1046.         } else {
  1047.             // Keplerian-only motion is all we can do
  1048.             return keplerianShifted;
  1049.         }

  1051.     }

  1053.     /** {@inheritDoc}
  1054.      * <p>
  1055.      * The interpolated instance is created by polynomial Hermite interpolation
  1056.      * on circular elements, without derivatives (which means the interpolation
  1057.      * falls back to Lagrange interpolation only).
  1058.      * </p>
  1059.      * <p>
  1060.      * As this implementation of interpolation is polynomial, it should be used only
  1061.      * with small samples (about 10-20 points) in order to avoid <a
  1062.      * href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
  1063.      * and numerical problems (including NaN appearing).
  1064.      * </p>
  1065.      * <p>
  1066.      * If orbit interpolation on large samples is needed, using the {@link
  1067.      * org.orekit.propagation.analytical.Ephemeris} class is a better way than using this
  1068.      * low-level interpolation. The Ephemeris class automatically handles selection of
  1069.      * a neighboring sub-sample with a predefined number of point from a large global sample
  1070.      * in a thread-safe way.
  1071.      * </p>
  1072.      */
  1073.     public FieldCircularOrbit<T> interpolate(final FieldAbsoluteDate<T> date, final Stream<FieldOrbit<T>> sample) {

  1075.         // first pass to check if derivatives are available throughout the sample
  1076.         final List<FieldOrbit<T>> list = sample.collect(Collectors.toList());
  1077.         boolean useDerivatives = true;
  1078.         for (final FieldOrbit<T> orbit : list) {
  1079.             useDerivatives = useDerivatives && orbit.hasDerivatives();
  1080.         }

  1082.         // set up an interpolator
  1083.         final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();

  1085.         // second pass to feed interpolator
  1086.         FieldAbsoluteDate<T> previousDate   = null;
  1087.         T                    previousRAAN   = zero.add(Double.NaN);
  1088.         T                    previousAlphaM = zero.add(Double.NaN);
  1089.         for (final FieldOrbit<T> orbit : list) {
  1090.             final FieldCircularOrbit<T> circ = (FieldCircularOrbit<T>) OrbitType.CIRCULAR.convertType(orbit);
  1091.             final T continuousRAAN;
  1092.             final T continuousAlphaM;
  1093.             if (previousDate == null) {
  1094.                 continuousRAAN   = circ.getRightAscensionOfAscendingNode();
  1095.                 continuousAlphaM = circ.getAlphaM();
  1096.             } else {
  1097.                 final T dt       = circ.getDate().durationFrom(previousDate);
  1098.                 final T keplerAM = previousAlphaM .add(circ.getKeplerianMeanMotion().multiply(dt));
  1099.                 continuousRAAN   = normalizeAngle(circ.getRightAscensionOfAscendingNode(), previousRAAN);
  1100.                 continuousAlphaM = normalizeAngle(circ.getAlphaM(), keplerAM);
  1101.             }
  1102.             previousDate   = circ.getDate();
  1103.             previousRAAN   = continuousRAAN;
  1104.             previousAlphaM = continuousAlphaM;
  1105.             final T[] toAdd = MathArrays.buildArray(one.getField(), 6);
  1106.             toAdd[0] = circ.getA();
  1107.             toAdd[1] = circ.getCircularEx();
  1108.             toAdd[2] = circ.getCircularEy();
  1109.             toAdd[3] = circ.getI();
  1110.             toAdd[4] = continuousRAAN;
  1111.             toAdd[5] = continuousAlphaM;
  1112.             if (useDerivatives) {
  1113.                 final T[] toAddDot = MathArrays.buildArray(one.getField(), 6);
  1114.                 toAddDot[0] = circ.getADot();
  1115.                 toAddDot[1] = circ.getCircularExDot();
  1116.                 toAddDot[2] = circ.getCircularEyDot();
  1117.                 toAddDot[3] = circ.getIDot();
  1118.                 toAddDot[4] = circ.getRightAscensionOfAscendingNodeDot();
  1119.                 toAddDot[5] = circ.getAlphaMDot();
  1120.                 interpolator.addSamplePoint(circ.getDate().durationFrom(date),
  1121.                                             toAdd, toAddDot);
  1122.             } else {
  1123.                 interpolator.addSamplePoint(circ.getDate().durationFrom(date),
  1124.                                             toAdd);
  1125.             }
  1126.         }

  1128.         // interpolate
  1129.         final T[][] interpolated = interpolator.derivatives(zero, 1);

  1131.         // build a new interpolated instance
  1132.         return new FieldCircularOrbit<>(interpolated[0][0], interpolated[0][1], interpolated[0][2],
  1133.                                         interpolated[0][3], interpolated[0][4], interpolated[0][5],
  1134.                                         interpolated[1][0], interpolated[1][1], interpolated[1][2],
  1135.                                         interpolated[1][3], interpolated[1][4], interpolated[1][5],
  1136.                                         PositionAngle.MEAN, getFrame(), date, getMu());

  1138.     }

  1140.     /** {@inheritDoc} */
  1141.     protected T[][] computeJacobianMeanWrtCartesian() {

  1143.         final T[][] jacobian = MathArrays.buildArray(one.getField(), 6, 6);

  1145.         // compute various intermediate parameters
  1146.         computePVWithoutA();
  1147.         final FieldVector3D<T> position = partialPV.getPosition();
  1148.         final FieldVector3D<T> velocity = partialPV.getVelocity();

  1150.         final T x          = position.getX();
  1151.         final T y          = position.getY();
  1152.         final T z          = position.getZ();
  1153.         final T vx         = velocity.getX();
  1154.         final T vy         = velocity.getY();
  1155.         final T vz         = velocity.getZ();
  1156.         final T pv         = FieldVector3D.dotProduct(position, velocity);
  1157.         final T r2         = position.getNormSq();
  1158.         final T r          = r2.sqrt();
  1159.         final T v2         = velocity.getNormSq();

  1161.         final double mu         = getMu();
  1162.         final T oOsqrtMuA  = one.divide(a.multiply(mu).sqrt());
  1163.         final T rOa        = r.divide(a);
  1164.         final T aOr        = a.divide(r);
  1165.         final T aOr2       = a.divide(r2);
  1166.         final T a2         = a.multiply(a);

  1168.         final T ex2        = ex.multiply(ex);
  1169.         final T ey2        = ey.multiply(ey);
  1170.         final T e2         = ex2.add(ey2);
  1171.         final T epsilon    = e2.negate().add(1.0).sqrt();
  1172.         final T beta       = epsilon.add(1).reciprocal();

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

  1177.         final T cosI       = i.cos();
  1178.         final T sinI       = i.sin();
  1179.         final T cosRaan    = raan.cos();
  1180.         final T sinRaan    = raan.sin();

  1182.         // da
  1183.         fillHalfRow(aOr.multiply(2.0).multiply(aOr2), position, jacobian[0], 0);
  1184.         fillHalfRow(a2.multiply(2.0 / mu), velocity, jacobian[0], 3);

  1186.         // differentials of the normalized momentum
  1187.         final FieldVector3D<T> danP = new FieldVector3D<>(v2, position, pv.negate(), velocity);
  1188.         final FieldVector3D<T> danV = new FieldVector3D<>(r2, velocity, pv.negate(), position);
  1189.         final T recip   = partialPV.getMomentum().getNorm().reciprocal();
  1190.         final T recip2  = recip.multiply(recip);
  1191.         final T recip2N = recip2.negate();
  1192.         final FieldVector3D<T> dwXP = new FieldVector3D<>(recip,
  1193.                                                           new FieldVector3D<>(zero, vz, vy.negate()),
  1194.                                                           recip2N.multiply(sinRaan).multiply(sinI),
  1195.                                                           danP);
  1196.         final FieldVector3D<T> dwYP = new FieldVector3D<>(recip,
  1197.                                                           new FieldVector3D<>(vz.negate(), zero, vx),
  1198.                                                           recip2.multiply(cosRaan).multiply(sinI),
  1199.                                                           danP);
  1200.         final FieldVector3D<T> dwZP = new FieldVector3D<>(recip,
  1201.                                                           new FieldVector3D<>(vy, vx.negate(), zero),
  1202.                                                           recip2N.multiply(cosI),
  1203.                                                           danP);
  1204.         final FieldVector3D<T> dwXV = new FieldVector3D<>(recip,
  1205.                                                           new FieldVector3D<>(zero, z.negate(), y),
  1206.                                                           recip2N.multiply(sinRaan).multiply(sinI),
  1207.                                                           danV);
  1208.         final FieldVector3D<T> dwYV = new FieldVector3D<>(recip,
  1209.                                                           new FieldVector3D<>(z, zero, x.negate()),
  1210.                                                           recip2.multiply(cosRaan).multiply(sinI),
  1211.                                                           danV);
  1212.         final FieldVector3D<T> dwZV = new FieldVector3D<>(recip,
  1213.                                                           new FieldVector3D<>(y.negate(), x, zero),
  1214.                                                           recip2N.multiply(cosI),
  1215.                                                           danV);

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

  1221.         // dRaan
  1222.         fillHalfRow(sinRaan.divide(sinI), dwYP, cosRaan.divide(sinI), dwXP, jacobian[4], 0);
  1223.         fillHalfRow(sinRaan.divide(sinI), dwYV, cosRaan.divide(sinI), dwXV, jacobian[4], 3);

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

  1231.         // du
  1232.         final FieldVector3D<T> duP = new FieldVector3D<>(cv.multiply(cosRaan).divide(sinI), dwXP,
  1233.                                                          cv.multiply(sinRaan).divide(sinI), dwYP,
  1234.                                                          one, new FieldVector3D<>(cosRaan, sinRaan, zero));
  1235.         final FieldVector3D<T> duV = new FieldVector3D<>(cv.multiply(cosRaan).divide(sinI), dwXV,
  1236.                                                          cv.multiply(sinRaan).divide(sinI), dwYV);

  1238.         // dv
  1239.         final FieldVector3D<T> dvP = new FieldVector3D<>(u.negate().multiply(cosRaan).multiply(cosI).divide(sinI).add(sinRaan.multiply(z)), dwXP,
  1240.                                                          u.negate().multiply(sinRaan).multiply(cosI).divide(sinI).subtract(cosRaan.multiply(z)), dwYP,
  1241.                                                          cv, dwZP,
  1242.                                                          one, new FieldVector3D<>(sinRaan.negate().multiply(cosI), cosRaan.multiply(cosI), sinI));
  1243.         final FieldVector3D<T> dvV = new FieldVector3D<>(u.negate().multiply(cosRaan).multiply(cosI).divide(sinI).add(sinRaan.multiply(z)), dwXV,
  1244.                                                          u.negate().multiply(sinRaan).multiply(cosI).divide(sinI).subtract(cosRaan.multiply(z)), dwYV,
  1245.                                                          cv, dwZV);

  1247.         final FieldVector3D<T> dc1P = new FieldVector3D<>(aOr2.multiply(eSinE.multiply(eSinE).multiply(2).add(1).subtract(eCosE)).divide(r2), position,
  1248.                                                           aOr2.multiply(-2).multiply(eSinE).multiply(oOsqrtMuA), velocity);
  1249.         final FieldVector3D<T> dc1V = new FieldVector3D<>(aOr2.multiply(-2).multiply(eSinE).multiply(oOsqrtMuA), position,
  1250.                                                           zero.add(2).divide(mu), velocity);
  1251.         final FieldVector3D<T> dc2P = new FieldVector3D<>(aOr2.multiply(eSinE).multiply(eSinE.multiply(eSinE).subtract(e2.negate().add(1))).divide(r2.multiply(epsilon)), position,
  1252.                                                           aOr2.multiply(e2.negate().add(1).subtract(eSinE.multiply(eSinE))).multiply(oOsqrtMuA).divide(epsilon), velocity);
  1253.         final FieldVector3D<T> dc2V = new FieldVector3D<>(aOr2.multiply(e2.negate().add(1).subtract(eSinE.multiply(eSinE))).multiply(oOsqrtMuA).divide(epsilon), position,
  1254.                                                           eSinE.divide(epsilon.multiply(mu)), velocity);

  1256.         final T cof1   = aOr2.multiply(eCosE.subtract(e2));
  1257.         final T cof2   = aOr2.multiply(epsilon).multiply(eSinE);
  1258.         final FieldVector3D<T> dexP = new FieldVector3D<>(u, dc1P,  v, dc2P, cof1, duP,  cof2, dvP);
  1259.         final FieldVector3D<T> dexV = new FieldVector3D<>(u, dc1V,  v, dc2V, cof1, duV,  cof2, dvV);
  1260.         final FieldVector3D<T> deyP = new FieldVector3D<>(v, dc1P, u.negate(), dc2P, cof1, dvP, cof2.negate(), duP);
  1261.         final FieldVector3D<T> deyV = new FieldVector3D<>(v, dc1V, u.negate(), dc2V, cof1, dvV, cof2.negate(), duV);
  1262.         fillHalfRow(one, dexP, jacobian[1], 0);
  1263.         fillHalfRow(one, dexV, jacobian[1], 3);
  1264.         fillHalfRow(one, deyP, jacobian[2], 0);
  1265.         fillHalfRow(one, deyV, jacobian[2], 3);

  1267.         final T cle = u.divide(a).add(ex).subtract(eSinE.multiply(beta).multiply(ey));
  1268.         final T sle = v.divide(a).add(ey).add(eSinE.multiply(beta).multiply(ex));
  1269.         final T m1  = beta.multiply(eCosE);
  1270.         final T m2  = m1.multiply(eCosE).negate().add(1);
  1271.         final T m3  = (u.multiply(ey).subtract(v.multiply(ex))).add(eSinE.multiply(beta).multiply(u.multiply(ex).add(v.multiply(ey))));
  1272.         final T m4  = sle.negate().add(cle.multiply(eSinE).multiply(beta));
  1273.         final T m5  = cle.add(sle.multiply(eSinE).multiply(beta));
  1274.         final T kk = m3.multiply(2).divide(r).add(aOr.multiply(eSinE)).add(m1.multiply(eSinE).multiply(m1.add(1).subtract(aOr.add(1).multiply(m2))).divide(epsilon)).divide(r2);
  1275.         final T jj = (m1.multiply(m2).divide(epsilon).subtract(1)).multiply(oOsqrtMuA);
  1276.         fillHalfRow(kk, position,
  1277.                     jj, velocity,
  1278.                     m4, dexP,
  1279.                     m5, deyP,
  1280.                     sle.negate().divide(a), duP,
  1281.                     cle.divide(a), dvP,
  1282.                     jacobian[5], 0);
  1283.         final T ll = (m1.multiply(m2).divide(epsilon ).subtract(1)).multiply(oOsqrtMuA);
  1284.         final T mm = m3.multiply(2).add(eSinE.multiply(a)).add(m1.multiply(eSinE).multiply(r).multiply(eCosE.multiply(beta).multiply(2).subtract(aOr.multiply(m2))).divide(epsilon)).divide(mu);

  1286.         fillHalfRow(ll, position,
  1287.                     mm, velocity,
  1288.                     m4, dexV,
  1289.                     m5, deyV,
  1290.                     sle.negate().divide(a), duV,
  1291.                     cle.divide(a), dvV,
  1292.                     jacobian[5], 3);
  1293.         return jacobian;

  1295.     }

  1297.     /** {@inheritDoc} */
  1298.     protected T[][] computeJacobianEccentricWrtCartesian() {

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

  1303.         // Differentiating the Kepler equation aM = aE - ex sin aE + ey cos aE leads to:
  1304.         // daM = (1 - ex cos aE - ey sin aE) dE - sin aE dex + cos aE dey
  1305.         // which is inverted and rewritten as:
  1306.         // daE = a/r daM + sin aE a/r dex - cos aE a/r dey
  1307.         final T alphaE = getAlphaE();
  1308.         final T cosAe  = alphaE.cos();
  1309.         final T sinAe  = alphaE.sin();
  1310.         final T aOr    = one.divide(one.subtract(ex.multiply(cosAe)).subtract(ey.multiply(sinAe)));

  1312.         // update longitude row
  1313.         final T[] rowEx = jacobian[1];
  1314.         final T[] rowEy = jacobian[2];
  1315.         final T[] rowL  = jacobian[5];
  1316.         for (int j = 0; j < 6; ++j) {
  1317.          // rowL[j] = aOr * (      rowL[j] +   sinAe *        rowEx[j]   -         cosAe *        rowEy[j]);
  1318.             rowL[j] = aOr.multiply(rowL[j].add(sinAe.multiply(rowEx[j])).subtract( cosAe.multiply(rowEy[j])));
  1319.         }
  1320.         jacobian[5] = rowL;
  1321.         return jacobian;

  1323.     }
  1324.     /** {@inheritDoc} */
  1325.     protected T[][] computeJacobianTrueWrtCartesian() {

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

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

  1367.     }

  1369.     /** {@inheritDoc} */
  1370.     public void addKeplerContribution(final PositionAngle type, final double gm,
  1371.                                       final T[] pDot) {
  1372.         final T oMe2;
  1373.         final T ksi;
  1374.         final T n = a.reciprocal().multiply(gm).sqrt().divide(a);
  1375.         switch (type) {
  1376.             case MEAN :
  1377.                 pDot[5] = pDot[5].add(n);
  1378.                 break;
  1379.             case ECCENTRIC :
  1380.                 oMe2  = one.subtract(ex.multiply(ex)).subtract(ey.multiply(ey));
  1381.                 ksi   = one.add(ex.multiply(alphaV.cos())).add(ey.multiply(alphaV.sin()));
  1382.                 pDot[5] = pDot[5].add(n.multiply(ksi).divide(oMe2));
  1383.                 break;
  1384.             case TRUE :
  1385.                 oMe2  = one.subtract(ex.multiply(ex)).subtract(ey.multiply(ey));
  1386.                 ksi   = one.add(ex.multiply(alphaV.cos())).add(ey.multiply(alphaV.sin()));
  1387.                 pDot[5] = pDot[5].add(n.multiply(ksi).multiply(ksi).divide(oMe2.multiply(oMe2.sqrt())));
  1388.                 break;
  1389.             default :
  1390.                 throw new OrekitInternalError(null);
  1391.         }
  1392.     }

  1394.     /**  Returns a string representation of this Orbit object.
  1395.      * @return a string representation of this object
  1396.      */
  1397.     public String toString() {
  1398.         return new StringBuffer().append("circular parameters: ").append('{').
  1399.                                   append("a: ").append(a.getReal()).
  1400.                                   append(", ex: ").append(ex.getReal()).append(", ey: ").append(ey.getReal()).
  1401.                                   append(", i: ").append(FastMath.toDegrees(i.getReal())).
  1402.                                   append(", raan: ").append(FastMath.toDegrees(raan.getReal())).
  1403.                                   append(", alphaV: ").append(FastMath.toDegrees(alphaV.getReal())).
  1404.                                   append(";}").toString();
  1405.     }



  1409.     /**
  1410.      * Normalize an angle in a 2&pi; wide interval around a center value.
  1411.      * <p>This method has three main uses:</p>
  1412.      * <ul>
  1413.      *   <li>normalize an angle between 0 and 2&pi;:<br/>
  1414.      *       {@code a = MathUtils.normalizeAngle(a, FastMath.PI);}</li>
  1415.      *   <li>normalize an angle between -&pi; and +&pi;<br/>
  1416.      *       {@code a = MathUtils.normalizeAngle(a, 0.0);}</li>
  1417.      *   <li>compute the angle between two defining angular positions:<br>
  1418.      *       {@code angle = MathUtils.normalizeAngle(end, start) - start;}</li>
  1419.      * </ul>
  1420.      * <p>Note that due to numerical accuracy and since &pi; cannot be represented
  1421.      * exactly, the result interval is <em>closed</em>, it cannot be half-closed
  1422.      * as would be more satisfactory in a purely mathematical view.</p>
  1423.      * @param a angle to normalize
  1424.      * @param center center of the desired 2&pi; interval for the result
  1425.      * @param <T> the type of the field elements
  1426.      * @return a-2k&pi; with integer k and center-&pi; &lt;= a-2k&pi; &lt;= center+&pi;
  1427.      */
  1428.     public static <T extends RealFieldElement<T>> T normalizeAngle(final T a, final T center) {
  1429.         return a.subtract(2 * FastMath.PI * FastMath.floor((a.getReal() + FastMath.PI - center.getReal()) / (2 * FastMath.PI)));
  1430.     }

  1432.     @Override
  1433.     public CircularOrbit toOrbit() {
  1434.         if (hasDerivatives()) {
  1435.             return new CircularOrbit(a.getReal(), ex.getReal(), ey.getReal(),
  1436.                                      i.getReal(), raan.getReal(), alphaV.getReal(),
  1437.                                      aDot.getReal(), exDot.getReal(), eyDot.getReal(),
  1438.                                      iDot.getReal(), raanDot.getReal(), alphaVDot.getReal(),
  1439.                                      PositionAngle.TRUE, getFrame(),
  1440.                                      getDate().toAbsoluteDate(), getMu());
  1441.         } else {
  1442.             return new CircularOrbit(a.getReal(), ex.getReal(), ey.getReal(),
  1443.                                      i.getReal(), raan.getReal(), alphaV.getReal(),
  1444.                                      PositionAngle.TRUE, getFrame(),
  1445.                                      getDate().toAbsoluteDate(), getMu());
  1446.         }
  1447.     }


  1450. }
Created with JaCoCo 0.7.9.201702052155