CartesianOrbit.java
/* Copyright 2002-2020 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.orbits;
import java.io.Serializable;
import java.util.stream.Stream;
import org.hipparchus.analysis.differentiation.UnivariateDerivative2;
import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.exception.MathIllegalStateException;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.RotationConvention;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.SinCos;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.data.DataContext;
import org.orekit.errors.OrekitMessages;
import org.orekit.frames.Frame;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.CartesianDerivativesFilter;
import org.orekit.utils.FieldPVCoordinates;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.TimeStampedPVCoordinates;
/** This class holds Cartesian orbital parameters.
* <p>
* The parameters used internally are the Cartesian coordinates:
* <ul>
* <li>x</li>
* <li>y</li>
* <li>z</li>
* <li>xDot</li>
* <li>yDot</li>
* <li>zDot</li>
* </ul>
* contained in {@link PVCoordinates}.
*
* <p>
* Note that the implementation of this class delegates all non-Cartesian related
* computations ({@link #getA()}, {@link #getEquinoctialEx()}, ...) to an underlying
* instance of the {@link EquinoctialOrbit} class. This implies that using this class
* only for analytical computations which are always based on non-Cartesian
* parameters is perfectly possible but somewhat sub-optimal.
* </p>
* <p>
* The instance <code>CartesianOrbit</code> is guaranteed to be immutable.
* </p>
* @see Orbit
* @see KeplerianOrbit
* @see CircularOrbit
* @see EquinoctialOrbit
* @author Luc Maisonobe
* @author Guylaine Prat
* @author Fabien Maussion
* @author Véronique Pommier-Maurussane
*/
public class CartesianOrbit extends Orbit {
/** Serializable UID. */
private static final long serialVersionUID = 20170414L;
/** Indicator for non-Keplerian derivatives. */
private final transient boolean hasNonKeplerianAcceleration;
/** Underlying equinoctial orbit to which high-level methods are delegated. */
private transient EquinoctialOrbit equinoctial;
/** Constructor from Cartesian parameters.
*
* <p> The acceleration provided in {@code pvCoordinates} is accessible using
* {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
* use {@code mu} and the position to compute the acceleration, including
* {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
*
* @param pvaCoordinates the position, velocity and acceleration of the satellite.
* @param frame the frame in which the {@link PVCoordinates} are defined
* (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
* @param mu central attraction coefficient (m³/s²)
* @exception IllegalArgumentException if frame is not a {@link
* Frame#isPseudoInertial pseudo-inertial frame}
*/
public CartesianOrbit(final TimeStampedPVCoordinates pvaCoordinates,
final Frame frame, final double mu)
throws IllegalArgumentException {
super(pvaCoordinates, frame, mu);
hasNonKeplerianAcceleration = hasNonKeplerianAcceleration(pvaCoordinates, mu);
equinoctial = null;
}
/** Constructor from Cartesian parameters.
*
* <p> The acceleration provided in {@code pvCoordinates} is accessible using
* {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
* use {@code mu} and the position to compute the acceleration, including
* {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
*
* @param pvaCoordinates the position and velocity of the satellite.
* @param frame the frame in which the {@link PVCoordinates} are defined
* (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
* @param date date of the orbital parameters
* @param mu central attraction coefficient (m³/s²)
* @exception IllegalArgumentException if frame is not a {@link
* Frame#isPseudoInertial pseudo-inertial frame}
*/
public CartesianOrbit(final PVCoordinates pvaCoordinates, final Frame frame,
final AbsoluteDate date, final double mu)
throws IllegalArgumentException {
this(new TimeStampedPVCoordinates(date, pvaCoordinates), frame, mu);
}
/** Constructor from any kind of orbital parameters.
* @param op orbital parameters to copy
*/
public CartesianOrbit(final Orbit op) {
super(op.getPVCoordinates(), op.getFrame(), op.getMu());
hasNonKeplerianAcceleration = op.hasDerivatives();
if (op instanceof EquinoctialOrbit) {
equinoctial = (EquinoctialOrbit) op;
} else if (op instanceof CartesianOrbit) {
equinoctial = ((CartesianOrbit) op).equinoctial;
} else {
equinoctial = null;
}
}
/** {@inheritDoc} */
public OrbitType getType() {
return OrbitType.CARTESIAN;
}
/** Lazy evaluation of equinoctial parameters. */
private void initEquinoctial() {
if (equinoctial == null) {
if (hasDerivatives()) {
// getPVCoordinates includes accelerations that will be interpreted as derivatives
equinoctial = new EquinoctialOrbit(getPVCoordinates(), getFrame(), getDate(), getMu());
} else {
// get rid of Keplerian acceleration so we don't assume
// we have derivatives when in fact we don't have them
equinoctial = new EquinoctialOrbit(new PVCoordinates(getPVCoordinates().getPosition(),
getPVCoordinates().getVelocity()),
getFrame(), getDate(), getMu());
}
}
}
/** Get the position/velocity with derivatives.
* @return position/velocity with derivatives
* @since 10.2
*/
private FieldPVCoordinates<UnivariateDerivative2> getPVDerivatives() {
// PVA coordinates
final PVCoordinates pva = getPVCoordinates();
final Vector3D p = pva.getPosition();
final Vector3D v = pva.getVelocity();
final Vector3D a = pva.getAcceleration();
// Field coordinates
final FieldVector3D<UnivariateDerivative2> pG = new FieldVector3D<>(new UnivariateDerivative2(p.getX(), v.getX(), a.getX()),
new UnivariateDerivative2(p.getY(), v.getY(), a.getY()),
new UnivariateDerivative2(p.getZ(), v.getZ(), a.getZ()));
final FieldVector3D<UnivariateDerivative2> vG = new FieldVector3D<>(new UnivariateDerivative2(v.getX(), a.getX(), 0.0),
new UnivariateDerivative2(v.getY(), a.getY(), 0.0),
new UnivariateDerivative2(v.getZ(), a.getZ(), 0.0));
return new FieldPVCoordinates<>(pG, vG);
}
/** {@inheritDoc} */
public double getA() {
final double r = getPVCoordinates().getPosition().getNorm();
final double V2 = getPVCoordinates().getVelocity().getNormSq();
return r / (2 - r * V2 / getMu());
}
/** {@inheritDoc} */
public double getADot() {
if (hasDerivatives()) {
final FieldPVCoordinates<UnivariateDerivative2> pv = getPVDerivatives();
final UnivariateDerivative2 r = pv.getPosition().getNorm();
final UnivariateDerivative2 V2 = pv.getVelocity().getNormSq();
final UnivariateDerivative2 a = r.divide(r.multiply(V2).divide(getMu()).subtract(2).negate());
return a.getDerivative(1);
} else {
return Double.NaN;
}
}
/** {@inheritDoc} */
public double getE() {
final double muA = getMu() * getA();
if (muA > 0) {
// elliptic or circular orbit
final Vector3D pvP = getPVCoordinates().getPosition();
final Vector3D pvV = getPVCoordinates().getVelocity();
final double rV2OnMu = pvP.getNorm() * pvV.getNormSq() / getMu();
final double eSE = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(muA);
final double eCE = rV2OnMu - 1;
return FastMath.sqrt(eCE * eCE + eSE * eSE);
} else {
// hyperbolic orbit
final Vector3D pvM = getPVCoordinates().getMomentum();
return FastMath.sqrt(1 - pvM.getNormSq() / muA);
}
}
/** {@inheritDoc} */
public double getEDot() {
if (hasDerivatives()) {
final FieldPVCoordinates<UnivariateDerivative2> pv = getPVDerivatives();
final FieldVector3D<UnivariateDerivative2> pvP = pv.getPosition();
final FieldVector3D<UnivariateDerivative2> pvV = pv.getVelocity();
final UnivariateDerivative2 r = pvP.getNorm();
final UnivariateDerivative2 V2 = pvV.getNormSq();
final UnivariateDerivative2 rV2OnMu = r.multiply(V2).divide(getMu());
final UnivariateDerivative2 a = r.divide(rV2OnMu.negate().add(2));
final UnivariateDerivative2 eSE = FieldVector3D.dotProduct(pvP, pvV).divide(a.multiply(getMu()).sqrt());
final UnivariateDerivative2 eCE = rV2OnMu.subtract(1);
final UnivariateDerivative2 e = eCE.multiply(eCE).add(eSE.multiply(eSE)).sqrt();
return e.getDerivative(1);
} else {
return Double.NaN;
}
}
/** {@inheritDoc} */
public double getI() {
return Vector3D.angle(Vector3D.PLUS_K, getPVCoordinates().getMomentum());
}
/** {@inheritDoc} */
public double getIDot() {
if (hasDerivatives()) {
final FieldPVCoordinates<UnivariateDerivative2> pv = getPVDerivatives();
final FieldVector3D<UnivariateDerivative2> momentum =
FieldVector3D.crossProduct(pv.getPosition(), pv.getVelocity());
final UnivariateDerivative2 i = FieldVector3D.angle(Vector3D.PLUS_K, momentum);
return i.getDerivative(1);
} else {
return Double.NaN;
}
}
/** {@inheritDoc} */
public double getEquinoctialEx() {
initEquinoctial();
return equinoctial.getEquinoctialEx();
}
/** {@inheritDoc} */
public double getEquinoctialExDot() {
initEquinoctial();
return equinoctial.getEquinoctialExDot();
}
/** {@inheritDoc} */
public double getEquinoctialEy() {
initEquinoctial();
return equinoctial.getEquinoctialEy();
}
/** {@inheritDoc} */
public double getEquinoctialEyDot() {
initEquinoctial();
return equinoctial.getEquinoctialEyDot();
}
/** {@inheritDoc} */
public double getHx() {
final Vector3D w = getPVCoordinates().getMomentum().normalize();
// Check for equatorial retrograde orbit
if (((w.getX() * w.getX() + w.getY() * w.getY()) == 0) && w.getZ() < 0) {
return Double.NaN;
}
return -w.getY() / (1 + w.getZ());
}
/** {@inheritDoc} */
public double getHxDot() {
if (hasDerivatives()) {
final FieldPVCoordinates<UnivariateDerivative2> pv = getPVDerivatives();
final FieldVector3D<UnivariateDerivative2> w =
FieldVector3D.crossProduct(pv.getPosition(), pv.getVelocity()).normalize();
// Check for equatorial retrograde orbit
final double x = w.getX().getValue();
final double y = w.getY().getValue();
final double z = w.getZ().getValue();
if (((x * x + y * y) == 0) && z < 0) {
return Double.NaN;
}
final UnivariateDerivative2 hx = w.getY().negate().divide(w.getZ().add(1));
return hx.getDerivative(1);
} else {
return Double.NaN;
}
}
/** {@inheritDoc} */
public double getHy() {
final Vector3D w = getPVCoordinates().getMomentum().normalize();
// Check for equatorial retrograde orbit
if (((w.getX() * w.getX() + w.getY() * w.getY()) == 0) && w.getZ() < 0) {
return Double.NaN;
}
return w.getX() / (1 + w.getZ());
}
/** {@inheritDoc} */
public double getHyDot() {
if (hasDerivatives()) {
final FieldPVCoordinates<UnivariateDerivative2> pv = getPVDerivatives();
final FieldVector3D<UnivariateDerivative2> w =
FieldVector3D.crossProduct(pv.getPosition(), pv.getVelocity()).normalize();
// Check for equatorial retrograde orbit
final double x = w.getX().getValue();
final double y = w.getY().getValue();
final double z = w.getZ().getValue();
if (((x * x + y * y) == 0) && z < 0) {
return Double.NaN;
}
final UnivariateDerivative2 hy = w.getX().divide(w.getZ().add(1));
return hy.getDerivative(1);
} else {
return Double.NaN;
}
}
/** {@inheritDoc} */
public double getLv() {
initEquinoctial();
return equinoctial.getLv();
}
/** {@inheritDoc} */
public double getLvDot() {
initEquinoctial();
return equinoctial.getLvDot();
}
/** {@inheritDoc} */
public double getLE() {
initEquinoctial();
return equinoctial.getLE();
}
/** {@inheritDoc} */
public double getLEDot() {
initEquinoctial();
return equinoctial.getLEDot();
}
/** {@inheritDoc} */
public double getLM() {
initEquinoctial();
return equinoctial.getLM();
}
/** {@inheritDoc} */
public double getLMDot() {
initEquinoctial();
return equinoctial.getLMDot();
}
/** {@inheritDoc} */
public boolean hasDerivatives() {
return hasNonKeplerianAcceleration;
}
/** {@inheritDoc} */
protected TimeStampedPVCoordinates initPVCoordinates() {
// nothing to do here, as the canonical elements are already the Cartesian ones
return getPVCoordinates();
}
/** {@inheritDoc} */
public CartesianOrbit shiftedBy(final double dt) {
final PVCoordinates shiftedPV = (getA() < 0) ? shiftPVHyperbolic(dt) : shiftPVElliptic(dt);
return new CartesianOrbit(shiftedPV, getFrame(), getDate().shiftedBy(dt), getMu());
}
/** {@inheritDoc}
* <p>
* The interpolated instance is created by polynomial Hermite interpolation
* ensuring velocity remains the exact derivative of position.
* </p>
* <p>
* As this implementation of interpolation is polynomial, it should be used only
* with small samples (about 10-20 points) in order to avoid <a
* href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
* and numerical problems (including NaN appearing).
* </p>
* <p>
* If orbit interpolation on large samples is needed, using the {@link
* org.orekit.propagation.analytical.Ephemeris} class is a better way than using this
* low-level interpolation. The Ephemeris class automatically handles selection of
* a neighboring sub-sample with a predefined number of point from a large global sample
* in a thread-safe way.
* </p>
*/
public CartesianOrbit interpolate(final AbsoluteDate date, final Stream<Orbit> sample) {
final TimeStampedPVCoordinates interpolated =
TimeStampedPVCoordinates.interpolate(date, CartesianDerivativesFilter.USE_PVA,
sample.map(orbit -> orbit.getPVCoordinates()));
return new CartesianOrbit(interpolated, getFrame(), date, getMu());
}
/** Compute shifted position and velocity in elliptic case.
* @param dt time shift
* @return shifted position and velocity
*/
private PVCoordinates shiftPVElliptic(final double dt) {
// preliminary computation
final Vector3D pvP = getPVCoordinates().getPosition();
final Vector3D pvV = getPVCoordinates().getVelocity();
final double r2 = pvP.getNormSq();
final double r = FastMath.sqrt(r2);
final double rV2OnMu = r * pvV.getNormSq() / getMu();
final double a = getA();
final double eSE = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(getMu() * a);
final double eCE = rV2OnMu - 1;
final double e2 = eCE * eCE + eSE * eSE;
// we can use any arbitrary reference 2D frame in the orbital plane
// in order to simplify some equations below, we use the current position as the u axis
final Vector3D u = pvP.normalize();
final Vector3D v = Vector3D.crossProduct(getPVCoordinates().getMomentum(), u).normalize();
// the following equations rely on the specific choice of u explained above,
// some coefficients that vanish to 0 in this case have already been removed here
final double ex = (eCE - e2) * a / r;
final double ey = -FastMath.sqrt(1 - e2) * eSE * a / r;
final double beta = 1 / (1 + FastMath.sqrt(1 - e2));
final double thetaE0 = FastMath.atan2(ey + eSE * beta * ex, r / a + ex - eSE * beta * ey);
final SinCos scThetaE0 = FastMath.sinCos(thetaE0);
final double thetaM0 = thetaE0 - ex * scThetaE0.sin() + ey * scThetaE0.cos();
// compute in-plane shifted eccentric argument
final double thetaM1 = thetaM0 + getKeplerianMeanMotion() * dt;
final double thetaE1 = meanToEccentric(thetaM1, ex, ey);
final SinCos scTE = FastMath.sinCos(thetaE1);
final double cTE = scTE.cos();
final double sTE = scTE.sin();
// compute shifted in-plane Cartesian coordinates
final double exey = ex * ey;
final double exCeyS = ex * cTE + ey * sTE;
final double x = a * ((1 - beta * ey * ey) * cTE + beta * exey * sTE - ex);
final double y = a * ((1 - beta * ex * ex) * sTE + beta * exey * cTE - ey);
final double factor = FastMath.sqrt(getMu() / a) / (1 - exCeyS);
final double xDot = factor * (-sTE + beta * ey * exCeyS);
final double yDot = factor * ( cTE - beta * ex * exCeyS);
final Vector3D shiftedP = new Vector3D(x, u, y, v);
final Vector3D shiftedV = new Vector3D(xDot, u, yDot, v);
if (hasNonKeplerianAcceleration) {
// extract non-Keplerian part of the initial acceleration
final Vector3D nonKeplerianAcceleration = new Vector3D(1, getPVCoordinates().getAcceleration(),
getMu() / (r2 * r), pvP);
// add the quadratic motion due to the non-Keplerian acceleration to the Keplerian motion
final Vector3D fixedP = new Vector3D(1, shiftedP,
0.5 * dt * dt, nonKeplerianAcceleration);
final double fixedR2 = fixedP.getNormSq();
final double fixedR = FastMath.sqrt(fixedR2);
final Vector3D fixedV = new Vector3D(1, shiftedV,
dt, nonKeplerianAcceleration);
final Vector3D fixedA = new Vector3D(-getMu() / (fixedR2 * fixedR), shiftedP,
1, nonKeplerianAcceleration);
return new PVCoordinates(fixedP, fixedV, fixedA);
} else {
// don't include acceleration,
// so the shifted orbit is not considered to have derivatives
return new PVCoordinates(shiftedP, shiftedV);
}
}
/** Compute shifted position and velocity in hyperbolic case.
* @param dt time shift
* @return shifted position and velocity
*/
private PVCoordinates shiftPVHyperbolic(final double dt) {
final PVCoordinates pv = getPVCoordinates();
final Vector3D pvP = pv.getPosition();
final Vector3D pvV = pv.getVelocity();
final Vector3D pvM = pv.getMomentum();
final double r2 = pvP.getNormSq();
final double r = FastMath.sqrt(r2);
final double rV2OnMu = r * pvV.getNormSq() / getMu();
final double a = getA();
final double muA = getMu() * a;
final double e = FastMath.sqrt(1 - Vector3D.dotProduct(pvM, pvM) / muA);
final double sqrt = FastMath.sqrt((e + 1) / (e - 1));
// compute mean anomaly
final double eSH = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(-muA);
final double eCH = rV2OnMu - 1;
final double H0 = FastMath.log((eCH + eSH) / (eCH - eSH)) / 2;
final double M0 = e * FastMath.sinh(H0) - H0;
// find canonical 2D frame with p pointing to perigee
final double v0 = 2 * FastMath.atan(sqrt * FastMath.tanh(H0 / 2));
final Vector3D p = new Rotation(pvM, v0, RotationConvention.FRAME_TRANSFORM).applyTo(pvP).normalize();
final Vector3D q = Vector3D.crossProduct(pvM, p).normalize();
// compute shifted eccentric anomaly
final double M1 = M0 + getKeplerianMeanMotion() * dt;
final double H1 = meanToHyperbolicEccentric(M1, e);
// compute shifted in-plane Cartesian coordinates
final double cH = FastMath.cosh(H1);
final double sH = FastMath.sinh(H1);
final double sE2m1 = FastMath.sqrt((e - 1) * (e + 1));
// coordinates of position and velocity in the orbital plane
final double x = a * (cH - e);
final double y = -a * sE2m1 * sH;
final double factor = FastMath.sqrt(getMu() / -a) / (e * cH - 1);
final double xDot = -factor * sH;
final double yDot = factor * sE2m1 * cH;
final Vector3D shiftedP = new Vector3D(x, p, y, q);
final Vector3D shiftedV = new Vector3D(xDot, p, yDot, q);
if (hasNonKeplerianAcceleration) {
// extract non-Keplerian part of the initial acceleration
final Vector3D nonKeplerianAcceleration = new Vector3D(1, getPVCoordinates().getAcceleration(),
getMu() / (r2 * r), pvP);
// add the quadratic motion due to the non-Keplerian acceleration to the Keplerian motion
final Vector3D fixedP = new Vector3D(1, shiftedP,
0.5 * dt * dt, nonKeplerianAcceleration);
final double fixedR2 = fixedP.getNormSq();
final double fixedR = FastMath.sqrt(fixedR2);
final Vector3D fixedV = new Vector3D(1, shiftedV,
dt, nonKeplerianAcceleration);
final Vector3D fixedA = new Vector3D(-getMu() / (fixedR2 * fixedR), shiftedP,
1, nonKeplerianAcceleration);
return new PVCoordinates(fixedP, fixedV, fixedA);
} else {
// don't include acceleration,
// so the shifted orbit is not considered to have derivatives
return new PVCoordinates(shiftedP, shiftedV);
}
}
/** Computes the eccentric in-plane argument from the mean in-plane argument.
* @param thetaM = mean in-plane argument (rad)
* @param ex first component of eccentricity vector
* @param ey second component of eccentricity vector
* @return the eccentric in-plane argument.
*/
private double meanToEccentric(final double thetaM, final double ex, final double ey) {
// Generalization of Kepler equation to in-plane parameters
// with thetaE = eta + E and
// thetaM = eta + M = thetaE - ex.sin(thetaE) + ey.cos(thetaE)
// and eta being counted from an arbitrary reference in the orbital plane
double thetaE = thetaM;
double thetaEMthetaM = 0.0;
int iter = 0;
do {
final SinCos scThetaE = FastMath.sinCos(thetaE);
final double f2 = ex * scThetaE.sin() - ey * scThetaE.cos();
final double f1 = 1.0 - ex * scThetaE.cos() - ey * scThetaE.sin();
final double f0 = thetaEMthetaM - f2;
final double f12 = 2.0 * f1;
final double shift = f0 * f12 / (f1 * f12 - f0 * f2);
thetaEMthetaM -= shift;
thetaE = thetaM + thetaEMthetaM;
if (FastMath.abs(shift) <= 1.0e-12) {
return thetaE;
}
} while (++iter < 50);
throw new MathIllegalStateException(LocalizedCoreFormats.CONVERGENCE_FAILED);
}
/** Computes the hyperbolic eccentric anomaly from the mean anomaly.
* <p>
* The algorithm used here for solving hyperbolic Kepler equation is
* Danby's iterative method (3rd order) with Vallado's initial guess.
* </p>
* @param M mean anomaly (rad)
* @param ecc eccentricity
* @return the hyperbolic eccentric anomaly
*/
private double meanToHyperbolicEccentric(final double M, final double ecc) {
// Resolution of hyperbolic Kepler equation for Keplerian parameters
// Initial guess
double H;
if (ecc < 1.6) {
if ((-FastMath.PI < M && M < 0.) || M > FastMath.PI) {
H = M - ecc;
} else {
H = M + ecc;
}
} else {
if (ecc < 3.6 && FastMath.abs(M) > FastMath.PI) {
H = M - FastMath.copySign(ecc, M);
} else {
H = M / (ecc - 1.);
}
}
// Iterative computation
int iter = 0;
do {
final double f3 = ecc * FastMath.cosh(H);
final double f2 = ecc * FastMath.sinh(H);
final double f1 = f3 - 1.;
final double f0 = f2 - H - M;
final double f12 = 2. * f1;
final double d = f0 / f12;
final double fdf = f1 - d * f2;
final double ds = f0 / fdf;
final double shift = f0 / (fdf + ds * ds * f3 / 6.);
H -= shift;
if (FastMath.abs(shift) <= 1.0e-12) {
return H;
}
} while (++iter < 50);
throw new MathIllegalStateException(OrekitMessages.UNABLE_TO_COMPUTE_HYPERBOLIC_ECCENTRIC_ANOMALY,
iter);
}
/** Create a 6x6 identity matrix.
* @return 6x6 identity matrix
*/
private double[][] create6x6Identity() {
// this is the fastest way to set the 6x6 identity matrix
final double[][] identity = new double[6][6];
for (int i = 0; i < 6; i++) {
identity[i][i] = 1.0;
}
return identity;
}
@Override
protected double[][] computeJacobianMeanWrtCartesian() {
return create6x6Identity();
}
@Override
protected double[][] computeJacobianEccentricWrtCartesian() {
return create6x6Identity();
}
@Override
protected double[][] computeJacobianTrueWrtCartesian() {
return create6x6Identity();
}
/** {@inheritDoc} */
public void addKeplerContribution(final PositionAngle type, final double gm,
final double[] pDot) {
final PVCoordinates pv = getPVCoordinates();
// position derivative is velocity
final Vector3D velocity = pv.getVelocity();
pDot[0] += velocity.getX();
pDot[1] += velocity.getY();
pDot[2] += velocity.getZ();
// velocity derivative is Newtonian acceleration
final Vector3D position = pv.getPosition();
final double r2 = position.getNormSq();
final double coeff = -gm / (r2 * FastMath.sqrt(r2));
pDot[3] += coeff * position.getX();
pDot[4] += coeff * position.getY();
pDot[5] += coeff * position.getZ();
}
/** Returns a string representation of this Orbit object.
* @return a string representation of this object
*/
public String toString() {
// use only the six defining elements, like the other Orbit.toString() methods
final String comma = ", ";
final PVCoordinates pv = getPVCoordinates();
final Vector3D p = pv.getPosition();
final Vector3D v = pv.getVelocity();
return "Cartesian parameters: {P(" +
p.getX() + comma +
p.getY() + comma +
p.getZ() + "), V(" +
v.getX() + comma +
v.getY() + comma +
v.getZ() + ")}";
}
/** Replace the instance with a data transfer object for serialization.
* <p>
* This intermediate class serializes all needed parameters,
* including position-velocity which are <em>not</em> serialized by parent
* {@link Orbit} class.
* </p>
* @return data transfer object that will be serialized
*/
@DefaultDataContext
private Object writeReplace() {
return new DTO(this);
}
/** Internal class used only for serialization. */
@DefaultDataContext
private static class DTO implements Serializable {
/** Serializable UID. */
private static final long serialVersionUID = 20170414L;
/** Double values. */
private double[] d;
/** Frame in which are defined the orbital parameters. */
private final Frame frame;
/** Simple constructor.
* @param orbit instance to serialize
*/
private DTO(final CartesianOrbit orbit) {
final TimeStampedPVCoordinates pv = orbit.getPVCoordinates();
// decompose date
final AbsoluteDate j2000Epoch =
DataContext.getDefault().getTimeScales().getJ2000Epoch();
final double epoch = FastMath.floor(pv.getDate().durationFrom(j2000Epoch));
final double offset = pv.getDate().durationFrom(j2000Epoch.shiftedBy(epoch));
if (orbit.hasDerivatives()) {
this.d = new double[] {
epoch, offset, orbit.getMu(),
pv.getPosition().getX(), pv.getPosition().getY(), pv.getPosition().getZ(),
pv.getVelocity().getX(), pv.getVelocity().getY(), pv.getVelocity().getZ(),
pv.getAcceleration().getX(), pv.getAcceleration().getY(), pv.getAcceleration().getZ()
};
} else {
this.d = new double[] {
epoch, offset, orbit.getMu(),
pv.getPosition().getX(), pv.getPosition().getY(), pv.getPosition().getZ(),
pv.getVelocity().getX(), pv.getVelocity().getY(), pv.getVelocity().getZ()
};
}
this.frame = orbit.getFrame();
}
/** Replace the deserialized data transfer object with a {@link CartesianOrbit}.
* @return replacement {@link CartesianOrbit}
*/
private Object readResolve() {
final AbsoluteDate j2000Epoch =
DataContext.getDefault().getTimeScales().getJ2000Epoch();
if (d.length >= 12) {
// we have derivatives
return new CartesianOrbit(new TimeStampedPVCoordinates(j2000Epoch.shiftedBy(d[0]).shiftedBy(d[1]),
new Vector3D(d[3], d[ 4], d[ 5]),
new Vector3D(d[6], d[ 7], d[ 8]),
new Vector3D(d[9], d[10], d[11])),
frame, d[2]);
} else {
// we don't have derivatives
return new CartesianOrbit(new TimeStampedPVCoordinates(j2000Epoch.shiftedBy(d[0]).shiftedBy(d[1]),
new Vector3D(d[3], d[ 4], d[ 5]),
new Vector3D(d[6], d[ 7], d[ 8])),
frame, d[2]);
}
}
}
}