CartesianOrbit.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;

  18. import java.io.Serializable;
  19. import java.util.stream.Stream;

  20. import org.hipparchus.analysis.differentiation.UnivariateDerivative2;
  21. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  22. import org.hipparchus.geometry.euclidean.threed.Rotation;
  23. import org.hipparchus.geometry.euclidean.threed.RotationConvention;
  24. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  25. import org.hipparchus.util.FastMath;
  26. import org.hipparchus.util.SinCos;
  27. import org.orekit.annotation.DefaultDataContext;
  28. import org.orekit.data.DataContext;
  29. import org.orekit.frames.Frame;
  30. import org.orekit.time.AbsoluteDate;
  31. import org.orekit.utils.CartesianDerivativesFilter;
  32. import org.orekit.utils.FieldPVCoordinates;
  33. import org.orekit.utils.PVCoordinates;
  34. import org.orekit.utils.TimeStampedPVCoordinates;


  35. /** This class holds Cartesian orbital parameters.

  36.  * <p>
  37.  * The parameters used internally are the Cartesian coordinates:
  38.  *   <ul>
  39.  *     <li>x</li>
  40.  *     <li>y</li>
  41.  *     <li>z</li>
  42.  *     <li>xDot</li>
  43.  *     <li>yDot</li>
  44.  *     <li>zDot</li>
  45.  *   </ul>
  46.  * contained in {@link PVCoordinates}.
  47.  *

  48.  * <p>
  49.  * Note that the implementation of this class delegates all non-Cartesian related
  50.  * computations ({@link #getA()}, {@link #getEquinoctialEx()}, ...) to an underlying
  51.  * instance of the {@link EquinoctialOrbit} class. This implies that using this class
  52.  * only for analytical computations which are always based on non-Cartesian
  53.  * parameters is perfectly possible but somewhat sub-optimal.
  54.  * </p>
  55.  * <p>
  56.  * The instance <code>CartesianOrbit</code> is guaranteed to be immutable.
  57.  * </p>
  58.  * @see    Orbit
  59.  * @see    KeplerianOrbit
  60.  * @see    CircularOrbit
  61.  * @see    EquinoctialOrbit
  62.  * @author Luc Maisonobe
  63.  * @author Guylaine Prat
  64.  * @author Fabien Maussion
  65.  * @author V&eacute;ronique Pommier-Maurussane
  66.  * @author Andrew Goetz
  67.  */
  68. public class CartesianOrbit extends Orbit {

  69.     /** Serializable UID. */
  70.     private static final long serialVersionUID = 20170414L;

  71.     /** Indicator for non-Keplerian derivatives. */
  72.     private final transient boolean hasNonKeplerianAcceleration;

  73.     /** Underlying equinoctial orbit to which high-level methods are delegated. */
  74.     private transient EquinoctialOrbit equinoctial;

  75.     /** Constructor from Cartesian parameters.
  76.      *
  77.      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
  78.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  79.      * use {@code mu} and the position to compute the acceleration, including
  80.      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
  81.      *
  82.      * @param pvaCoordinates the position, velocity and acceleration of the satellite.
  83.      * @param frame the frame in which the {@link PVCoordinates} are defined
  84.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  85.      * @param mu central attraction coefficient (m³/s²)
  86.      * @exception IllegalArgumentException if frame is not a {@link
  87.      * Frame#isPseudoInertial pseudo-inertial frame}
  88.      */
  89.     public CartesianOrbit(final TimeStampedPVCoordinates pvaCoordinates,
  90.                           final Frame frame, final double mu)
  91.         throws IllegalArgumentException {
  92.         super(pvaCoordinates, frame, mu);
  93.         hasNonKeplerianAcceleration = hasNonKeplerianAcceleration(pvaCoordinates, mu);
  94.         equinoctial = null;
  95.     }

  96.     /** Constructor from Cartesian parameters.
  97.      *
  98.      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
  99.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  100.      * use {@code mu} and the position to compute the acceleration, including
  101.      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
  102.      *
  103.      * @param pvaCoordinates the position and velocity of the satellite.
  104.      * @param frame the frame in which the {@link PVCoordinates} are defined
  105.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  106.      * @param date date of the orbital parameters
  107.      * @param mu central attraction coefficient (m³/s²)
  108.      * @exception IllegalArgumentException if frame is not a {@link
  109.      * Frame#isPseudoInertial pseudo-inertial frame}
  110.      */
  111.     public CartesianOrbit(final PVCoordinates pvaCoordinates, final Frame frame,
  112.                           final AbsoluteDate date, final double mu)
  113.         throws IllegalArgumentException {
  114.         this(new TimeStampedPVCoordinates(date, pvaCoordinates), frame, mu);
  115.     }

  116.     /** Constructor from any kind of orbital parameters.
  117.      * @param op orbital parameters to copy
  118.      */
  119.     public CartesianOrbit(final Orbit op) {
  120.         super(op.getPVCoordinates(), op.getFrame(), op.getMu());
  121.         hasNonKeplerianAcceleration = op.hasDerivatives();
  122.         if (op instanceof EquinoctialOrbit) {
  123.             equinoctial = (EquinoctialOrbit) op;
  124.         } else if (op instanceof CartesianOrbit) {
  125.             equinoctial = ((CartesianOrbit) op).equinoctial;
  126.         } else {
  127.             equinoctial = null;
  128.         }
  129.     }

  130.     /** {@inheritDoc} */
  131.     public OrbitType getType() {
  132.         return OrbitType.CARTESIAN;
  133.     }

  134.     /** Lazy evaluation of equinoctial parameters. */
  135.     private void initEquinoctial() {
  136.         if (equinoctial == null) {
  137.             if (hasDerivatives()) {
  138.                 // getPVCoordinates includes accelerations that will be interpreted as derivatives
  139.                 equinoctial = new EquinoctialOrbit(getPVCoordinates(), getFrame(), getDate(), getMu());
  140.             } else {
  141.                 // get rid of Keplerian acceleration so we don't assume
  142.                 // we have derivatives when in fact we don't have them
  143.                 equinoctial = new EquinoctialOrbit(new PVCoordinates(getPVCoordinates().getPosition(),
  144.                                                                      getPVCoordinates().getVelocity()),
  145.                                                    getFrame(), getDate(), getMu());
  146.             }
  147.         }
  148.     }

  149.     /** Get the position/velocity with derivatives.
  150.      * @return position/velocity with derivatives
  151.      * @since 10.2
  152.      */
  153.     private FieldPVCoordinates<UnivariateDerivative2> getPVDerivatives() {
  154.         // PVA coordinates
  155.         final PVCoordinates pva = getPVCoordinates();
  156.         final Vector3D      p   = pva.getPosition();
  157.         final Vector3D      v   = pva.getVelocity();
  158.         final Vector3D      a   = pva.getAcceleration();
  159.         // Field coordinates
  160.         final FieldVector3D<UnivariateDerivative2> pG = new FieldVector3D<>(new UnivariateDerivative2(p.getX(), v.getX(), a.getX()),
  161.                                                                new UnivariateDerivative2(p.getY(), v.getY(), a.getY()),
  162.                                                                new UnivariateDerivative2(p.getZ(), v.getZ(), a.getZ()));
  163.         final FieldVector3D<UnivariateDerivative2> vG = new FieldVector3D<>(new UnivariateDerivative2(v.getX(), a.getX(), 0.0),
  164.                                                                new UnivariateDerivative2(v.getY(), a.getY(), 0.0),
  165.                                                                new UnivariateDerivative2(v.getZ(), a.getZ(), 0.0));
  166.         return new FieldPVCoordinates<>(pG, vG);
  167.     }

  168.     /** {@inheritDoc} */
  169.     public double getA() {
  170.         final double r  = getPVCoordinates().getPosition().getNorm();
  171.         final double V2 = getPVCoordinates().getVelocity().getNormSq();
  172.         return r / (2 - r * V2 / getMu());
  173.     }

  174.     /** {@inheritDoc} */
  175.     public double getADot() {
  176.         if (hasDerivatives()) {
  177.             final FieldPVCoordinates<UnivariateDerivative2> pv = getPVDerivatives();
  178.             final UnivariateDerivative2 r  = pv.getPosition().getNorm();
  179.             final UnivariateDerivative2 V2 = pv.getVelocity().getNormSq();
  180.             final UnivariateDerivative2 a  = r.divide(r.multiply(V2).divide(getMu()).subtract(2).negate());
  181.             return a.getDerivative(1);
  182.         } else {
  183.             return Double.NaN;
  184.         }
  185.     }

  186.     /** {@inheritDoc} */
  187.     public double getE() {
  188.         final double muA = getMu() * getA();
  189.         if (muA > 0) {
  190.             // elliptic or circular orbit
  191.             final Vector3D pvP   = getPVCoordinates().getPosition();
  192.             final Vector3D pvV   = getPVCoordinates().getVelocity();
  193.             final double rV2OnMu = pvP.getNorm() * pvV.getNormSq() / getMu();
  194.             final double eSE     = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(muA);
  195.             final double eCE     = rV2OnMu - 1;
  196.             return FastMath.sqrt(eCE * eCE + eSE * eSE);
  197.         } else {
  198.             // hyperbolic orbit
  199.             final Vector3D pvM = getPVCoordinates().getMomentum();
  200.             return FastMath.sqrt(1 - pvM.getNormSq() / muA);
  201.         }
  202.     }

  203.     /** {@inheritDoc} */
  204.     public double getEDot() {
  205.         if (hasDerivatives()) {
  206.             final FieldPVCoordinates<UnivariateDerivative2> pv = getPVDerivatives();
  207.             final FieldVector3D<UnivariateDerivative2> pvP   = pv.getPosition();
  208.             final FieldVector3D<UnivariateDerivative2> pvV   = pv.getVelocity();
  209.             final UnivariateDerivative2 r       = pvP.getNorm();
  210.             final UnivariateDerivative2 V2      = pvV.getNormSq();
  211.             final UnivariateDerivative2 rV2OnMu = r.multiply(V2).divide(getMu());
  212.             final UnivariateDerivative2 a       = r.divide(rV2OnMu.negate().add(2));
  213.             final UnivariateDerivative2 eSE     = FieldVector3D.dotProduct(pvP, pvV).divide(a.multiply(getMu()).sqrt());
  214.             final UnivariateDerivative2 eCE     = rV2OnMu.subtract(1);
  215.             final UnivariateDerivative2 e       = eCE.multiply(eCE).add(eSE.multiply(eSE)).sqrt();
  216.             return e.getDerivative(1);
  217.         } else {
  218.             return Double.NaN;
  219.         }
  220.     }

  221.     /** {@inheritDoc} */
  222.     public double getI() {
  223.         return Vector3D.angle(Vector3D.PLUS_K, getPVCoordinates().getMomentum());
  224.     }

  225.     /** {@inheritDoc} */
  226.     public double getIDot() {
  227.         if (hasDerivatives()) {
  228.             final FieldPVCoordinates<UnivariateDerivative2> pv = getPVDerivatives();
  229.             final FieldVector3D<UnivariateDerivative2> momentum =
  230.                             FieldVector3D.crossProduct(pv.getPosition(), pv.getVelocity());
  231.             final UnivariateDerivative2 i = FieldVector3D.angle(Vector3D.PLUS_K, momentum);
  232.             return i.getDerivative(1);
  233.         } else {
  234.             return Double.NaN;
  235.         }
  236.     }

  237.     /** {@inheritDoc} */
  238.     public double getEquinoctialEx() {
  239.         initEquinoctial();
  240.         return equinoctial.getEquinoctialEx();
  241.     }

  242.     /** {@inheritDoc} */
  243.     public double getEquinoctialExDot() {
  244.         initEquinoctial();
  245.         return equinoctial.getEquinoctialExDot();
  246.     }

  247.     /** {@inheritDoc} */
  248.     public double getEquinoctialEy() {
  249.         initEquinoctial();
  250.         return equinoctial.getEquinoctialEy();
  251.     }

  252.     /** {@inheritDoc} */
  253.     public double getEquinoctialEyDot() {
  254.         initEquinoctial();
  255.         return equinoctial.getEquinoctialEyDot();
  256.     }

  257.     /** {@inheritDoc} */
  258.     public double getHx() {
  259.         final Vector3D w = getPVCoordinates().getMomentum().normalize();
  260.         // Check for equatorial retrograde orbit
  261.         if ((w.getX() * w.getX() + w.getY() * w.getY()) == 0 && w.getZ() < 0) {
  262.             return Double.NaN;
  263.         }
  264.         return -w.getY() / (1 + w.getZ());
  265.     }

  266.     /** {@inheritDoc} */
  267.     public double getHxDot() {
  268.         if (hasDerivatives()) {
  269.             final FieldPVCoordinates<UnivariateDerivative2> pv = getPVDerivatives();
  270.             final FieldVector3D<UnivariateDerivative2> w =
  271.                             FieldVector3D.crossProduct(pv.getPosition(), pv.getVelocity()).normalize();
  272.             // Check for equatorial retrograde orbit
  273.             final double x = w.getX().getValue();
  274.             final double y = w.getY().getValue();
  275.             final double z = w.getZ().getValue();
  276.             if ((x * x + y * y) == 0 && z < 0) {
  277.                 return Double.NaN;
  278.             }
  279.             final UnivariateDerivative2 hx = w.getY().negate().divide(w.getZ().add(1));
  280.             return hx.getDerivative(1);
  281.         } else {
  282.             return Double.NaN;
  283.         }
  284.     }

  285.     /** {@inheritDoc} */
  286.     public double getHy() {
  287.         final Vector3D w = getPVCoordinates().getMomentum().normalize();
  288.         // Check for equatorial retrograde orbit
  289.         if ((w.getX() * w.getX() + w.getY() * w.getY()) == 0 && w.getZ() < 0) {
  290.             return Double.NaN;
  291.         }
  292.         return  w.getX() / (1 + w.getZ());
  293.     }

  294.     /** {@inheritDoc} */
  295.     public double getHyDot() {
  296.         if (hasDerivatives()) {
  297.             final FieldPVCoordinates<UnivariateDerivative2> pv = getPVDerivatives();
  298.             final FieldVector3D<UnivariateDerivative2> w =
  299.                             FieldVector3D.crossProduct(pv.getPosition(), pv.getVelocity()).normalize();
  300.             // Check for equatorial retrograde orbit
  301.             final double x = w.getX().getValue();
  302.             final double y = w.getY().getValue();
  303.             final double z = w.getZ().getValue();
  304.             if ((x * x + y * y) == 0 && z < 0) {
  305.                 return Double.NaN;
  306.             }
  307.             final UnivariateDerivative2 hy = w.getX().divide(w.getZ().add(1));
  308.             return hy.getDerivative(1);
  309.         } else {
  310.             return Double.NaN;
  311.         }
  312.     }

  313.     /** {@inheritDoc} */
  314.     public double getLv() {
  315.         initEquinoctial();
  316.         return equinoctial.getLv();
  317.     }

  318.     /** {@inheritDoc} */
  319.     public double getLvDot() {
  320.         initEquinoctial();
  321.         return equinoctial.getLvDot();
  322.     }

  323.     /** {@inheritDoc} */
  324.     public double getLE() {
  325.         initEquinoctial();
  326.         return equinoctial.getLE();
  327.     }

  328.     /** {@inheritDoc} */
  329.     public double getLEDot() {
  330.         initEquinoctial();
  331.         return equinoctial.getLEDot();
  332.     }

  333.     /** {@inheritDoc} */
  334.     public double getLM() {
  335.         initEquinoctial();
  336.         return equinoctial.getLM();
  337.     }

  338.     /** {@inheritDoc} */
  339.     public double getLMDot() {
  340.         initEquinoctial();
  341.         return equinoctial.getLMDot();
  342.     }

  343.     /** {@inheritDoc} */
  344.     public boolean hasDerivatives() {
  345.         return hasNonKeplerianAcceleration;
  346.     }

  347.     /** {@inheritDoc} */
  348.     protected TimeStampedPVCoordinates initPVCoordinates() {
  349.         // nothing to do here, as the canonical elements are already the Cartesian ones
  350.         return getPVCoordinates();
  351.     }

  352.     /** {@inheritDoc} */
  353.     public CartesianOrbit shiftedBy(final double dt) {
  354.         final PVCoordinates shiftedPV = (getA() < 0) ? shiftPVHyperbolic(dt) : shiftPVElliptic(dt);
  355.         return new CartesianOrbit(shiftedPV, getFrame(), getDate().shiftedBy(dt), getMu());
  356.     }

  357.     /** {@inheritDoc}
  358.      * <p>
  359.      * The interpolated instance is created by polynomial Hermite interpolation
  360.      * ensuring velocity remains the exact derivative of position.
  361.      * </p>
  362.      * <p>
  363.      * As this implementation of interpolation is polynomial, it should be used only
  364.      * with small samples (about 10-20 points) in order to avoid <a
  365.      * href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
  366.      * and numerical problems (including NaN appearing).
  367.      * </p>
  368.      * <p>
  369.      * If orbit interpolation on large samples is needed, using the {@link
  370.      * org.orekit.propagation.analytical.Ephemeris} class is a better way than using this
  371.      * low-level interpolation. The Ephemeris class automatically handles selection of
  372.      * a neighboring sub-sample with a predefined number of point from a large global sample
  373.      * in a thread-safe way.
  374.      * </p>
  375.      */
  376.     public CartesianOrbit interpolate(final AbsoluteDate date, final Stream<Orbit> sample) {
  377.         final TimeStampedPVCoordinates interpolated =
  378.                 TimeStampedPVCoordinates.interpolate(date, CartesianDerivativesFilter.USE_PVA,
  379.                                                      sample.map(orbit -> orbit.getPVCoordinates()));
  380.         return new CartesianOrbit(interpolated, getFrame(), date, getMu());
  381.     }

  382.     /** Compute shifted position and velocity in elliptic case.
  383.      * @param dt time shift
  384.      * @return shifted position and velocity
  385.      */
  386.     private PVCoordinates shiftPVElliptic(final double dt) {

  387.         final PVCoordinates pv = getPVCoordinates();
  388.         final Vector3D pvP     = pv.getPosition();
  389.         final Vector3D pvV     = pv.getVelocity();
  390.         final Vector3D pvM     = pv.getMomentum();
  391.         final double r2        = pvP.getNormSq();
  392.         final double r         = FastMath.sqrt(r2);
  393.         final double rV2OnMu   = r * pvV.getNormSq() / getMu();
  394.         final double a         = r / (2 - rV2OnMu);
  395.         final double muA       = getMu() * a;

  396.         // compute mean anomaly
  397.         final double eSE    = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(muA);
  398.         final double eCE    = rV2OnMu - 1;
  399.         final double E0     = FastMath.atan2(eSE, eCE);
  400.         final double M0     = E0 - eSE;

  401.         final double e         = FastMath.sqrt(eCE * eCE + eSE * eSE);
  402.         final double sqrt      = FastMath.sqrt((1 + e) / (1 - e));

  403.         // find canonical 2D frame with p pointing to perigee
  404.         final double v0     = 2 * FastMath.atan(sqrt * FastMath.tan(E0 / 2));
  405.         final Vector3D p    = new Rotation(pvM, v0, RotationConvention.FRAME_TRANSFORM).applyTo(pvP).normalize();
  406.         final Vector3D q    = Vector3D.crossProduct(pvM, p).normalize();

  407.         // compute shifted eccentric anomaly
  408.         final double M1     = M0 + getKeplerianMeanMotion() * dt;
  409.         final double E1     = KeplerianAnomalyUtility.ellipticMeanToEccentric(e, M1);

  410.         // compute shifted in-plane Cartesian coordinates
  411.         final SinCos scE    = FastMath.sinCos(E1);
  412.         final double cE     = scE.cos();
  413.         final double sE     = scE.sin();
  414.         final double sE2m1  = FastMath.sqrt((1 - e) * (1 + e));

  415.         // coordinates of position and velocity in the orbital plane
  416.         final double x      = a * (cE - e);
  417.         final double y      = a * sE2m1 * sE;
  418.         final double factor = FastMath.sqrt(getMu() / a) / (1 - e * cE);
  419.         final double xDot   = -factor * sE;
  420.         final double yDot   =  factor * sE2m1 * cE;

  421.         final Vector3D shiftedP = new Vector3D(x, p, y, q);
  422.         final Vector3D shiftedV = new Vector3D(xDot, p, yDot, q);
  423.         if (hasNonKeplerianAcceleration) {

  424.             // extract non-Keplerian part of the initial acceleration
  425.             final Vector3D nonKeplerianAcceleration = new Vector3D(1, getPVCoordinates().getAcceleration(),
  426.                                                                    getMu() / (r2 * r), pvP);

  427.             // add the quadratic motion due to the non-Keplerian acceleration to the Keplerian motion
  428.             final Vector3D fixedP   = new Vector3D(1, shiftedP,
  429.                                                    0.5 * dt * dt, nonKeplerianAcceleration);
  430.             final double   fixedR2 = fixedP.getNormSq();
  431.             final double   fixedR  = FastMath.sqrt(fixedR2);
  432.             final Vector3D fixedV  = new Vector3D(1, shiftedV,
  433.                                                   dt, nonKeplerianAcceleration);
  434.             final Vector3D fixedA  = new Vector3D(-getMu() / (fixedR2 * fixedR), shiftedP,
  435.                                                   1, nonKeplerianAcceleration);

  436.             return new PVCoordinates(fixedP, fixedV, fixedA);

  437.         } else {
  438.             // don't include acceleration,
  439.             // so the shifted orbit is not considered to have derivatives
  440.             return new PVCoordinates(shiftedP, shiftedV);
  441.         }

  442.     }

  443.     /** Compute shifted position and velocity in hyperbolic case.
  444.      * @param dt time shift
  445.      * @return shifted position and velocity
  446.      */
  447.     private PVCoordinates shiftPVHyperbolic(final double dt) {

  448.         final PVCoordinates pv = getPVCoordinates();
  449.         final Vector3D pvP   = pv.getPosition();
  450.         final Vector3D pvV   = pv.getVelocity();
  451.         final Vector3D pvM   = pv.getMomentum();
  452.         final double r2      = pvP.getNormSq();
  453.         final double r       = FastMath.sqrt(r2);
  454.         final double rV2OnMu = r * pvV.getNormSq() / getMu();
  455.         final double a       = getA();
  456.         final double muA     = getMu() * a;
  457.         final double e       = FastMath.sqrt(1 - Vector3D.dotProduct(pvM, pvM) / muA);
  458.         final double sqrt    = FastMath.sqrt((e + 1) / (e - 1));

  459.         // compute mean anomaly
  460.         final double eSH     = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(-muA);
  461.         final double eCH     = rV2OnMu - 1;
  462.         final double H0      = FastMath.log((eCH + eSH) / (eCH - eSH)) / 2;
  463.         final double M0      = e * FastMath.sinh(H0) - H0;

  464.         // find canonical 2D frame with p pointing to perigee
  465.         final double v0      = 2 * FastMath.atan(sqrt * FastMath.tanh(H0 / 2));
  466.         final Vector3D p     = new Rotation(pvM, v0, RotationConvention.FRAME_TRANSFORM).applyTo(pvP).normalize();
  467.         final Vector3D q     = Vector3D.crossProduct(pvM, p).normalize();

  468.         // compute shifted eccentric anomaly
  469.         final double M1      = M0 + getKeplerianMeanMotion() * dt;
  470.         final double H1      = KeplerianAnomalyUtility.hyperbolicMeanToEccentric(e, M1);

  471.         // compute shifted in-plane Cartesian coordinates
  472.         final double cH     = FastMath.cosh(H1);
  473.         final double sH     = FastMath.sinh(H1);
  474.         final double sE2m1  = FastMath.sqrt((e - 1) * (e + 1));

  475.         // coordinates of position and velocity in the orbital plane
  476.         final double x      = a * (cH - e);
  477.         final double y      = -a * sE2m1 * sH;
  478.         final double factor = FastMath.sqrt(getMu() / -a) / (e * cH - 1);
  479.         final double xDot   = -factor * sH;
  480.         final double yDot   =  factor * sE2m1 * cH;

  481.         final Vector3D shiftedP = new Vector3D(x, p, y, q);
  482.         final Vector3D shiftedV = new Vector3D(xDot, p, yDot, q);
  483.         if (hasNonKeplerianAcceleration) {

  484.             // extract non-Keplerian part of the initial acceleration
  485.             final Vector3D nonKeplerianAcceleration = new Vector3D(1, getPVCoordinates().getAcceleration(),
  486.                                                                    getMu() / (r2 * r), pvP);

  487.             // add the quadratic motion due to the non-Keplerian acceleration to the Keplerian motion
  488.             final Vector3D fixedP   = new Vector3D(1, shiftedP,
  489.                                                    0.5 * dt * dt, nonKeplerianAcceleration);
  490.             final double   fixedR2 = fixedP.getNormSq();
  491.             final double   fixedR  = FastMath.sqrt(fixedR2);
  492.             final Vector3D fixedV  = new Vector3D(1, shiftedV,
  493.                                                   dt, nonKeplerianAcceleration);
  494.             final Vector3D fixedA  = new Vector3D(-getMu() / (fixedR2 * fixedR), shiftedP,
  495.                                                   1, nonKeplerianAcceleration);

  496.             return new PVCoordinates(fixedP, fixedV, fixedA);

  497.         } else {
  498.             // don't include acceleration,
  499.             // so the shifted orbit is not considered to have derivatives
  500.             return new PVCoordinates(shiftedP, shiftedV);
  501.         }

  502.     }

  503.     /** Create a 6x6 identity matrix.
  504.      * @return 6x6 identity matrix
  505.      */
  506.     private double[][] create6x6Identity() {
  507.         // this is the fastest way to set the 6x6 identity matrix
  508.         final double[][] identity = new double[6][6];
  509.         for (int i = 0; i < 6; i++) {
  510.             identity[i][i] = 1.0;
  511.         }
  512.         return identity;
  513.     }

  514.     @Override
  515.     protected double[][] computeJacobianMeanWrtCartesian() {
  516.         return create6x6Identity();
  517.     }

  518.     @Override
  519.     protected double[][] computeJacobianEccentricWrtCartesian() {
  520.         return create6x6Identity();
  521.     }

  522.     @Override
  523.     protected double[][] computeJacobianTrueWrtCartesian() {
  524.         return create6x6Identity();
  525.     }

  526.     /** {@inheritDoc} */
  527.     public void addKeplerContribution(final PositionAngle type, final double gm,
  528.                                       final double[] pDot) {

  529.         final PVCoordinates pv = getPVCoordinates();

  530.         // position derivative is velocity
  531.         final Vector3D velocity = pv.getVelocity();
  532.         pDot[0] += velocity.getX();
  533.         pDot[1] += velocity.getY();
  534.         pDot[2] += velocity.getZ();

  535.         // velocity derivative is Newtonian acceleration
  536.         final Vector3D position = pv.getPosition();
  537.         final double r2         = position.getNormSq();
  538.         final double coeff      = -gm / (r2 * FastMath.sqrt(r2));
  539.         pDot[3] += coeff * position.getX();
  540.         pDot[4] += coeff * position.getY();
  541.         pDot[5] += coeff * position.getZ();

  542.     }

  543.     /**  Returns a string representation of this Orbit object.
  544.      * @return a string representation of this object
  545.      */
  546.     public String toString() {
  547.         // use only the six defining elements, like the other Orbit.toString() methods
  548.         final String comma = ", ";
  549.         final PVCoordinates pv = getPVCoordinates();
  550.         final Vector3D p = pv.getPosition();
  551.         final Vector3D v = pv.getVelocity();
  552.         return "Cartesian parameters: {P(" +
  553.                 p.getX() + comma +
  554.                 p.getY() + comma +
  555.                 p.getZ() + "), V(" +
  556.                 v.getX() + comma +
  557.                 v.getY() + comma +
  558.                 v.getZ() + ")}";
  559.     }

  560.     /** Replace the instance with a data transfer object for serialization.
  561.      * <p>
  562.      * This intermediate class serializes all needed parameters,
  563.      * including position-velocity which are <em>not</em> serialized by parent
  564.      * {@link Orbit} class.
  565.      * </p>
  566.      * @return data transfer object that will be serialized
  567.      */
  568.     @DefaultDataContext
  569.     private Object writeReplace() {
  570.         return new DTO(this);
  571.     }

  572.     /** Internal class used only for serialization. */
  573.     @DefaultDataContext
  574.     private static class DTO implements Serializable {

  575.         /** Serializable UID. */
  576.         private static final long serialVersionUID = 20170414L;

  577.         /** Double values. */
  578.         private double[] d;

  579.         /** Frame in which are defined the orbital parameters. */
  580.         private final Frame frame;

  581.         /** Simple constructor.
  582.          * @param orbit instance to serialize
  583.          */
  584.         private DTO(final CartesianOrbit orbit) {

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

  586.             // decompose date
  587.             final AbsoluteDate j2000Epoch =
  588.                     DataContext.getDefault().getTimeScales().getJ2000Epoch();
  589.             final double epoch  = FastMath.floor(pv.getDate().durationFrom(j2000Epoch));
  590.             final double offset = pv.getDate().durationFrom(j2000Epoch.shiftedBy(epoch));

  591.             if (orbit.hasDerivatives()) {
  592.                 this.d = new double[] {
  593.                     epoch, offset, orbit.getMu(),
  594.                     pv.getPosition().getX(),     pv.getPosition().getY(),     pv.getPosition().getZ(),
  595.                     pv.getVelocity().getX(),     pv.getVelocity().getY(),     pv.getVelocity().getZ(),
  596.                     pv.getAcceleration().getX(), pv.getAcceleration().getY(), pv.getAcceleration().getZ()
  597.                 };
  598.             } else {
  599.                 this.d = new double[] {
  600.                     epoch, offset, orbit.getMu(),
  601.                     pv.getPosition().getX(),     pv.getPosition().getY(),     pv.getPosition().getZ(),
  602.                     pv.getVelocity().getX(),     pv.getVelocity().getY(),     pv.getVelocity().getZ()
  603.                 };
  604.             }

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

  606.         }

  607.         /** Replace the deserialized data transfer object with a {@link CartesianOrbit}.
  608.          * @return replacement {@link CartesianOrbit}
  609.          */
  610.         private Object readResolve() {
  611.             final AbsoluteDate j2000Epoch =
  612.                     DataContext.getDefault().getTimeScales().getJ2000Epoch();
  613.             if (d.length >= 12) {
  614.                 // we have derivatives
  615.                 return new CartesianOrbit(new TimeStampedPVCoordinates(j2000Epoch.shiftedBy(d[0]).shiftedBy(d[1]),
  616.                                                                        new Vector3D(d[3], d[ 4], d[ 5]),
  617.                                                                        new Vector3D(d[6], d[ 7], d[ 8]),
  618.                                                                        new Vector3D(d[9], d[10], d[11])),
  619.                                           frame, d[2]);
  620.             } else {
  621.                 // we don't have derivatives
  622.                 return new CartesianOrbit(new TimeStampedPVCoordinates(j2000Epoch.shiftedBy(d[0]).shiftedBy(d[1]),
  623.                                                                        new Vector3D(d[3], d[ 4], d[ 5]),
  624.                                                                        new Vector3D(d[6], d[ 7], d[ 8])),
  625.                                           frame, d[2]);
  626.             }
  627.         }

  628.     }

  629. }