CartesianOrbit.java

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

  19. import org.hipparchus.analysis.differentiation.UnivariateDerivative2;
  20. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  21. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  22. import org.hipparchus.linear.MatrixUtils;
  23. import org.hipparchus.util.FastMath;
  24. import org.orekit.frames.Frame;
  25. import org.orekit.frames.KinematicTransform;
  26. import org.orekit.time.AbsoluteDate;
  27. import org.orekit.time.TimeOffset;
  28. import org.orekit.utils.FieldPVCoordinates;
  29. import org.orekit.utils.PVCoordinates;
  30. import org.orekit.utils.TimeStampedPVCoordinates;


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

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

  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 {

  70.     /** 6x6 identity matrix. */
  71.     private static final double[][] SIX_BY_SIX_IDENTITY = MatrixUtils.createRealIdentityMatrix(6).getData();

  73.     /** Indicator for non-Keplerian derivatives. */
  74.     private final transient boolean hasNonKeplerianAcceleration;

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

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

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

  122.     /** Constructor from any kind of orbital parameters.
  123.      * @param op orbital parameters to copy
  124.      */
  125.     public CartesianOrbit(final Orbit op) {
  126.         super(op.getPVCoordinates(), op.getFrame(), op.getMu());
  127.         hasNonKeplerianAcceleration = op.hasNonKeplerianAcceleration();
  128.         if (op instanceof EquinoctialOrbit) {
  129.             equinoctial = (EquinoctialOrbit) op;
  130.         } else if (op instanceof CartesianOrbit) {
  131.             equinoctial = ((CartesianOrbit) op).equinoctial;
  132.         } else {
  133.             equinoctial = null;
  134.         }
  135.     }

  137.     /** {@inheritDoc} */
  138.     public OrbitType getType() {
  139.         return OrbitType.CARTESIAN;
  140.     }

  142.     /** {@inheritDoc} */
  143.     @Override
  144.     protected Vector3D nonKeplerianAcceleration() {
  145.         final double norm = getPosition().getNorm();
  146.         return getPVCoordinates().getAcceleration().add(new Vector3D(getMu() / (norm * norm * norm), getPosition()));
  147.     }

  149.     /** Lazy evaluation of equinoctial parameters. */
  150.     private void initEquinoctial() {
  151.         if (equinoctial == null) {
  152.             if (hasNonKeplerianAcceleration()) {
  153.                 // getPVCoordinates includes accelerations that will be interpreted as derivatives
  154.                 equinoctial = new EquinoctialOrbit(getPVCoordinates(), getFrame(), getDate(), getMu());
  155.             } else {
  156.                 // get rid of Keplerian acceleration so we don't assume
  157.                 // we have derivatives when in fact we don't have them
  158.                 equinoctial = new EquinoctialOrbit(new PVCoordinates(getPosition(),
  159.                                                                      getPVCoordinates().getVelocity()),
  160.                                                    getFrame(), getDate(), getMu());
  161.             }
  162.         }
  163.     }

  165.     /** Get the position/velocity with derivatives.
  166.      * @return position/velocity with derivatives
  167.      * @since 10.2
  168.      */
  169.     private FieldPVCoordinates<UnivariateDerivative2> getPVDerivatives() {
  170.         // PVA coordinates
  171.         final PVCoordinates pva = getPVCoordinates();
  172.         final Vector3D      p   = pva.getPosition();
  173.         final Vector3D      v   = pva.getVelocity();
  174.         final Vector3D      a   = pva.getAcceleration();
  175.         // Field coordinates
  176.         final FieldVector3D<UnivariateDerivative2> pG = new FieldVector3D<>(new UnivariateDerivative2(p.getX(), v.getX(), a.getX()),
  177.                                                                new UnivariateDerivative2(p.getY(), v.getY(), a.getY()),
  178.                                                                new UnivariateDerivative2(p.getZ(), v.getZ(), a.getZ()));
  179.         final FieldVector3D<UnivariateDerivative2> vG = new FieldVector3D<>(new UnivariateDerivative2(v.getX(), a.getX(), 0.0),
  180.                                                                new UnivariateDerivative2(v.getY(), a.getY(), 0.0),
  181.                                                                new UnivariateDerivative2(v.getZ(), a.getZ(), 0.0));
  182.         return new FieldPVCoordinates<>(pG, vG);
  183.     }

  185.     /** {@inheritDoc} */
  186.     public double getA() {
  187.         final double r  = getPosition().getNorm();
  188.         final double V2 = getPVCoordinates().getVelocity().getNormSq();
  189.         return r / (2 - r * V2 / getMu());
  190.     }

  192.     /** {@inheritDoc} */
  193.     public double getADot() {
  194.         if (hasNonKeplerianAcceleration) {
  195.             final FieldPVCoordinates<UnivariateDerivative2> pv = getPVDerivatives();
  196.             final UnivariateDerivative2 r  = pv.getPosition().getNorm();
  197.             final UnivariateDerivative2 V2 = pv.getVelocity().getNormSq();
  198.             final UnivariateDerivative2 a  = r.divide(r.multiply(V2).divide(getMu()).subtract(2).negate());
  199.             return a.getDerivative(1);
  200.         } else {
  201.             return 0.;
  202.         }
  203.     }

  205.     /** {@inheritDoc} */
  206.     public double getE() {
  207.         final double muA = getMu() * getA();
  208.         if (isElliptical()) {
  209.             // elliptic or circular orbit
  210.             final Vector3D pvP   = getPosition();
  211.             final Vector3D pvV   = getPVCoordinates().getVelocity();
  212.             final double rV2OnMu = pvP.getNorm() * pvV.getNormSq() / getMu();
  213.             final double eSE     = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(muA);
  214.             final double eCE     = rV2OnMu - 1;
  215.             return FastMath.sqrt(eCE * eCE + eSE * eSE);
  216.         } else {
  217.             // hyperbolic orbit
  218.             final Vector3D pvM = getPVCoordinates().getMomentum();
  219.             return FastMath.sqrt(1 - pvM.getNormSq() / muA);
  220.         }
  221.     }

  223.     /** {@inheritDoc} */
  224.     public double getEDot() {
  225.         if (hasNonKeplerianAcceleration) {
  226.             final FieldPVCoordinates<UnivariateDerivative2> pv = getPVDerivatives();
  227.             final FieldVector3D<UnivariateDerivative2> pvP   = pv.getPosition();
  228.             final FieldVector3D<UnivariateDerivative2> pvV   = pv.getVelocity();
  229.             final UnivariateDerivative2 r       = pvP.getNorm();
  230.             final UnivariateDerivative2 V2      = pvV.getNormSq();
  231.             final UnivariateDerivative2 rV2OnMu = r.multiply(V2).divide(getMu());
  232.             final UnivariateDerivative2 a       = r.divide(rV2OnMu.negate().add(2));
  233.             final UnivariateDerivative2 eSE     = FieldVector3D.dotProduct(pvP, pvV).divide(a.multiply(getMu()).sqrt());
  234.             final UnivariateDerivative2 eCE     = rV2OnMu.subtract(1);
  235.             final UnivariateDerivative2 e       = eCE.multiply(eCE).add(eSE.multiply(eSE)).sqrt();
  236.             return e.getDerivative(1);
  237.         } else {
  238.             return 0.;
  239.         }
  240.     }

  242.     /** {@inheritDoc} */
  243.     public double getI() {
  244.         return Vector3D.angle(Vector3D.PLUS_K, getPVCoordinates().getMomentum());
  245.     }

  247.     /** {@inheritDoc} */
  248.     public double getIDot() {
  249.         if (hasNonKeplerianAcceleration) {
  250.             final FieldPVCoordinates<UnivariateDerivative2> pv = getPVDerivatives();
  251.             final FieldVector3D<UnivariateDerivative2> momentum =
  252.                             FieldVector3D.crossProduct(pv.getPosition(), pv.getVelocity());
  253.             final UnivariateDerivative2 i = FieldVector3D.angle(Vector3D.PLUS_K, momentum);
  254.             return i.getDerivative(1);
  255.         } else {
  256.             return 0.;
  257.         }
  258.     }

  260.     /** {@inheritDoc} */
  261.     public double getEquinoctialEx() {
  262.         initEquinoctial();
  263.         return equinoctial.getEquinoctialEx();
  264.     }

  266.     /** {@inheritDoc} */
  267.     public double getEquinoctialExDot() {
  268.         initEquinoctial();
  269.         return equinoctial.getEquinoctialExDot();
  270.     }

  272.     /** {@inheritDoc} */
  273.     public double getEquinoctialEy() {
  274.         initEquinoctial();
  275.         return equinoctial.getEquinoctialEy();
  276.     }

  278.     /** {@inheritDoc} */
  279.     public double getEquinoctialEyDot() {
  280.         initEquinoctial();
  281.         return equinoctial.getEquinoctialEyDot();
  282.     }

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

  294.     /** {@inheritDoc} */
  295.     public double getHxDot() {
  296.         if (hasNonKeplerianAcceleration) {
  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 hx = w.getY().negate().divide(w.getZ().add(1));
  308.             return hx.getDerivative(1);
  309.         } else {
  310.             return 0.;
  311.         }
  312.     }

  314.     /** {@inheritDoc} */
  315.     public double getHy() {
  316.         final Vector3D w = getPVCoordinates().getMomentum().normalize();
  317.         // Check for equatorial retrograde orbit
  318.         if ((w.getX() * w.getX() + w.getY() * w.getY()) == 0 && w.getZ() < 0) {
  319.             return Double.NaN;
  320.         }
  321.         return  w.getX() / (1 + w.getZ());
  322.     }

  324.     /** {@inheritDoc} */
  325.     public double getHyDot() {
  326.         if (hasNonKeplerianAcceleration) {
  327.             final FieldPVCoordinates<UnivariateDerivative2> pv = getPVDerivatives();
  328.             final FieldVector3D<UnivariateDerivative2> w =
  329.                             FieldVector3D.crossProduct(pv.getPosition(), pv.getVelocity()).normalize();
  330.             // Check for equatorial retrograde orbit
  331.             final double x = w.getX().getValue();
  332.             final double y = w.getY().getValue();
  333.             final double z = w.getZ().getValue();
  334.             if ((x * x + y * y) == 0 && z < 0) {
  335.                 return Double.NaN;
  336.             }
  337.             final UnivariateDerivative2 hy = w.getX().divide(w.getZ().add(1));
  338.             return hy.getDerivative(1);
  339.         } else {
  340.             return 0.;
  341.         }
  342.     }

  344.     /** {@inheritDoc} */
  345.     public double getLv() {
  346.         initEquinoctial();
  347.         return equinoctial.getLv();
  348.     }

  350.     /** {@inheritDoc} */
  351.     public double getLvDot() {
  352.         initEquinoctial();
  353.         return equinoctial.getLvDot();
  354.     }

  356.     /** {@inheritDoc} */
  357.     public double getLE() {
  358.         initEquinoctial();
  359.         return equinoctial.getLE();
  360.     }

  362.     /** {@inheritDoc} */
  363.     public double getLEDot() {
  364.         initEquinoctial();
  365.         return equinoctial.getLEDot();
  366.     }

  368.     /** {@inheritDoc} */
  369.     public double getLM() {
  370.         initEquinoctial();
  371.         return equinoctial.getLM();
  372.     }

  374.     /** {@inheritDoc} */
  375.     public double getLMDot() {
  376.         initEquinoctial();
  377.         return equinoctial.getLMDot();
  378.     }

  380.     /** {@inheritDoc} */
  381.     @Override
  382.     public boolean hasNonKeplerianAcceleration() {
  383.         return hasNonKeplerianAcceleration;
  384.     }

  386.     /** {@inheritDoc} */
  387.     protected Vector3D initPosition() {
  388.         // nothing to do here, as the canonical elements are already the Cartesian ones
  389.         return getPVCoordinates().getPosition();
  390.     }

  392.     /** {@inheritDoc} */
  393.     protected TimeStampedPVCoordinates initPVCoordinates() {
  394.         // nothing to do here, as the canonical elements are already the Cartesian ones
  395.         return getPVCoordinates();
  396.     }

  398.     /** {@inheritDoc} */
  399.     @Override
  400.     public CartesianOrbit inFrame(final Frame inertialFrame) {
  401.         if (hasNonKeplerianAcceleration()) {
  402.             return new CartesianOrbit(getPVCoordinates(inertialFrame), inertialFrame, getMu());
  403.         } else {
  404.             final KinematicTransform transform = getFrame().getKinematicTransformTo(inertialFrame, getDate());
  405.             return new CartesianOrbit(transform.transformOnlyPV(getPVCoordinates()), inertialFrame, getDate(), getMu());
  406.         }
  407.     }

  409.     /** {@inheritDoc} */
  410.     public CartesianOrbit shiftedBy(final double dt) {
  411.         final PVCoordinates shiftedPV = shiftPV(dt);
  412.         return new CartesianOrbit(shiftedPV, getFrame(), getDate().shiftedBy(dt), getMu());
  413.     }

  415.     /** {@inheritDoc} */
  416.     public CartesianOrbit shiftedBy(final TimeOffset dt) {
  417.         final PVCoordinates shiftedPV = shiftPV(dt.toDouble());
  418.         return new CartesianOrbit(shiftedPV, getFrame(), getDate().shiftedBy(dt), getMu());
  419.     }

  421.     /** Compute shifted position and velocity.
  422.      * @param dt time shift
  423.      * @return shifted position and velocity
  424.      */
  425.     private PVCoordinates shiftPV(final double dt) {

  427.         final Vector3D pvP = getPosition();
  428.         final PVCoordinates shiftedPV = KeplerianMotionCartesianUtility.predictPositionVelocity(dt, pvP,
  429.             getPVCoordinates().getVelocity(), getMu());

  431.         if (dt != 0. && hasNonKeplerianAcceleration) {

  433.             // extract non-Keplerian part of the initial acceleration
  434.             final double r2 = pvP.getNormSq();
  435.             final double r = FastMath.sqrt(r2);
  436.             final Vector3D nonKeplerianAcceleration = new Vector3D(1, getPVCoordinates().getAcceleration(),
  437.                                                                    getMu() / (r2 * r), pvP);

  439.             // add the quadratic motion due to the non-Keplerian acceleration to the Keplerian motion
  440.             final Vector3D shiftedP = shiftedPV.getPosition();
  441.             final Vector3D shiftedV = shiftedPV.getVelocity();
  442.             final Vector3D fixedP   = new Vector3D(1, shiftedP,
  443.                                                    0.5 * dt * dt, nonKeplerianAcceleration);
  444.             final double   fixedR2 = fixedP.getNormSq();
  445.             final double   fixedR  = FastMath.sqrt(fixedR2);
  446.             final Vector3D fixedV  = new Vector3D(1, shiftedV,
  447.                                                   dt, nonKeplerianAcceleration);
  448.             final Vector3D fixedA  = new Vector3D(-getMu() / (fixedR2 * fixedR), shiftedP,
  449.                                                   1, nonKeplerianAcceleration);

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

  453.         } else {
  454.             // don't include acceleration,
  455.             // so the shifted orbit is not considered to have derivatives
  456.             return shiftedPV;
  457.         }

  459.     }

  461.     @Override
  462.     protected double[][] computeJacobianMeanWrtCartesian() {
  463.         return SIX_BY_SIX_IDENTITY;
  464.     }

  466.     @Override
  467.     protected double[][] computeJacobianEccentricWrtCartesian() {
  468.         return SIX_BY_SIX_IDENTITY;
  469.     }

  471.     @Override
  472.     protected double[][] computeJacobianTrueWrtCartesian() {
  473.         return SIX_BY_SIX_IDENTITY;
  474.     }

  476.     /** {@inheritDoc} */
  477.     public void addKeplerContribution(final PositionAngleType type, final double gm,
  478.                                       final double[] pDot) {

  480.         final PVCoordinates pv = getPVCoordinates();

  482.         // position derivative is velocity
  483.         final Vector3D velocity = pv.getVelocity();
  484.         pDot[0] += velocity.getX();
  485.         pDot[1] += velocity.getY();
  486.         pDot[2] += velocity.getZ();

  488.         // velocity derivative is Newtonian acceleration
  489.         final Vector3D position = pv.getPosition();
  490.         final double r2         = position.getNormSq();
  491.         final double coeff      = -gm / (r2 * FastMath.sqrt(r2));
  492.         pDot[3] += coeff * position.getX();
  493.         pDot[4] += coeff * position.getY();
  494.         pDot[5] += coeff * position.getZ();

  496.     }

  498.     /**  Returns a string representation of this Orbit object.
  499.      * @return a string representation of this object
  500.      */
  501.     public String toString() {
  502.         // use only the six defining elements, like the other Orbit.toString() methods
  503.         final String comma = ", ";
  504.         final PVCoordinates pv = getPVCoordinates();
  505.         final Vector3D p = pv.getPosition();
  506.         final Vector3D v = pv.getVelocity();
  507.         return "Cartesian parameters: {P(" +
  508.                 p.getX() + comma +
  509.                 p.getY() + comma +
  510.                 p.getZ() + "), V(" +
  511.                 v.getX() + comma +
  512.                 v.getY() + comma +
  513.                 v.getZ() + ")}";
  514.     }

  516. }
Created with JaCoCo 0.8.12.202403310830