From 70ceeb4555389bc7d0c78cc7a6224f2caec126dd Mon Sep 17 00:00:00 2001 From: Andrew Goetz Date: Tue, 4 Oct 2022 10:38:11 -0400 Subject: [PATCH 1/2] Extracted Kepler anomaly conversions to separate classes. Deprecated old methods in CartesianOrbit and FieldCartesianOrbit for removal. --- .../org/orekit/orbits/CartesianOrbit.java | 58 +-- .../orekit/orbits/FieldCartesianOrbit.java | 63 +--- .../orbits/FieldKeplerianAnomalyUtility.java | 342 ++++++++++++++++++ .../orekit/orbits/FieldKeplerianOrbit.java | 236 +++--------- .../orbits/KeplerianAnomalyUtility.java | 306 ++++++++++++++++ .../org/orekit/orbits/KeplerianOrbit.java | 220 ++--------- .../FieldKeplerianAnomalyUtilityTest.java | 232 ++++++++++++ .../orbits/FieldKeplerianOrbitTest.java | 16 - .../orbits/KeplerianAnomalyUtilityTest.java | 177 +++++++++ .../org/orekit/orbits/KeplerianOrbitTest.java | 12 - 10 files changed, 1140 insertions(+), 522 deletions(-) create mode 100644 src/main/java/org/orekit/orbits/FieldKeplerianAnomalyUtility.java create mode 100644 src/main/java/org/orekit/orbits/KeplerianAnomalyUtility.java create mode 100644 src/test/java/org/orekit/orbits/FieldKeplerianAnomalyUtilityTest.java create mode 100644 src/test/java/org/orekit/orbits/KeplerianAnomalyUtilityTest.java diff --git a/src/main/java/org/orekit/orbits/CartesianOrbit.java b/src/main/java/org/orekit/orbits/CartesianOrbit.java index bdb2b0861..da35099a7 100644 --- a/src/main/java/org/orekit/orbits/CartesianOrbit.java +++ b/src/main/java/org/orekit/orbits/CartesianOrbit.java @@ -30,7 +30,6 @@ 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; @@ -527,7 +526,7 @@ public class CartesianOrbit extends Orbit { // compute shifted eccentric anomaly final double M1 = M0 + getKeplerianMeanMotion() * dt; - final double H1 = meanToHyperbolicEccentric(M1, e); + final double H1 = KeplerianAnomalyUtility.hyperbolicMeanToEccentric(e, M1); // compute shifted in-plane Cartesian coordinates final double cH = FastMath.cosh(H1); @@ -606,61 +605,6 @@ public class CartesianOrbit extends Orbit { } - /** Computes the hyperbolic eccentric anomaly from the mean anomaly. - *

- * The algorithm used here for solving hyperbolic Kepler equation is - * Danby's iterative method (3rd order) with Vallado's initial guess. - *

- * @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 */ diff --git a/src/main/java/org/orekit/orbits/FieldCartesianOrbit.java b/src/main/java/org/orekit/orbits/FieldCartesianOrbit.java index dc6fd76c4..bb86ece8d 100644 --- a/src/main/java/org/orekit/orbits/FieldCartesianOrbit.java +++ b/src/main/java/org/orekit/orbits/FieldCartesianOrbit.java @@ -20,8 +20,8 @@ package org.orekit.orbits; import java.util.Arrays; import java.util.stream.Stream; -import org.hipparchus.Field; import org.hipparchus.CalculusFieldElement; +import org.hipparchus.Field; import org.hipparchus.analysis.differentiation.FieldUnivariateDerivative2; import org.hipparchus.exception.LocalizedCoreFormats; import org.hipparchus.exception.MathIllegalStateException; @@ -32,7 +32,6 @@ import org.hipparchus.geometry.euclidean.threed.Vector3D; import org.hipparchus.util.FastMath; import org.hipparchus.util.FieldSinCos; import org.hipparchus.util.MathArrays; -import org.orekit.errors.OrekitMessages; import org.orekit.frames.Frame; import org.orekit.time.AbsoluteDate; import org.orekit.time.FieldAbsoluteDate; @@ -554,7 +553,7 @@ public class FieldCartesianOrbit> extends Fiel // compute shifted eccentric anomaly final T M1 = M0.add(getKeplerianMeanMotion().multiply(dt)); - final T H1 = meanToHyperbolicEccentric(M1, e); + final T H1 = FieldKeplerianAnomalyUtility.hyperbolicMeanToEccentric(e, M1); // compute shifted in-plane Cartesian coordinates final T cH = H1.cosh(); @@ -633,64 +632,6 @@ public class FieldCartesianOrbit> extends Fiel } - /** Computes the hyperbolic eccentric anomaly from the mean anomaly. - *

- * The algorithm used here for solving hyperbolic Kepler equation is - * Danby's iterative method (3rd order) with Vallado's initial guess. - *

- * @param M mean anomaly (rad) - * @param ecc eccentricity - * @return the hyperbolic eccentric anomaly - */ - private T meanToHyperbolicEccentric(final T M, final T ecc) { - - // Resolution of hyperbolic Kepler equation for Keplerian parameters - - // Field value of pi - final T pi = ecc.getPi(); - - // Initial guess - T H; - if (ecc.getReal() < 1.6) { - if (-pi.getReal() < M.getReal() && M.getReal() < 0. || M.getReal() > pi.getReal()) { - H = M.subtract(ecc); - } else { - H = M.add(ecc); - } - } else { - if (ecc.getReal() < 3.6 && FastMath.abs(M.getReal()) > pi.getReal()) { - H = M.subtract(ecc.copySign(M)); - } else { - H = M.divide(ecc.subtract(1.)); - } - } - - // Iterative computation - int iter = 0; - do { - final T f3 = ecc.multiply(H.cosh()); - final T f2 = ecc.multiply(H.sinh()); - final T f1 = f3.subtract(1.); - final T f0 = f2.subtract(H).subtract(M); - final T f12 = f1.multiply(2); - final T d = f0.divide(f12); - final T fdf = f1.subtract(d.multiply(f2)); - final T ds = f0.divide(fdf); - - final T shift = f0.divide(fdf.add(ds.multiply(ds).multiply(f3.divide(6.)))); - - H = H.subtract(shift); - - if (FastMath.abs(shift.getReal()) <= 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 */ diff --git a/src/main/java/org/orekit/orbits/FieldKeplerianAnomalyUtility.java b/src/main/java/org/orekit/orbits/FieldKeplerianAnomalyUtility.java new file mode 100644 index 000000000..a4afef804 --- /dev/null +++ b/src/main/java/org/orekit/orbits/FieldKeplerianAnomalyUtility.java @@ -0,0 +1,342 @@ +/* Copyright 2002-2022 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 org.hipparchus.CalculusFieldElement; +import org.hipparchus.exception.MathIllegalArgumentException; +import org.hipparchus.util.FastMath; +import org.hipparchus.util.FieldSinCos; +import org.hipparchus.util.MathUtils; +import org.hipparchus.util.Precision; +import org.orekit.errors.OrekitMessages; + +/** + * Utility methods for converting between different Keplerian anomalies. + * @author Luc Maisonobe + * @author Andrea Antolino + * @author Andrew Goetz + */ +public class FieldKeplerianAnomalyUtility { + + /** First coefficient to compute elliptic Kepler equation solver starter. */ + private static final double A; + + /** Second coefficient to compute elliptic Kepler equation solver starter. */ + private static final double B; + + static { + final double k1 = 3 * FastMath.PI + 2; + final double k2 = FastMath.PI - 1; + final double k3 = 6 * FastMath.PI - 1; + A = 3 * k2 * k2 / k1; + B = k3 * k3 / (6 * k1); + } + + private FieldKeplerianAnomalyUtility() { + } + + /** + * Computes the elliptic true anomaly from the elliptic mean anomaly. + * @param field type + * @param e eccentricity such that 0 ≤ e < 1 + * @param M elliptic mean anomaly (rad) + * @return elliptic true anomaly (rad) + */ + public static > T ellipticMeanToTrue(final T e, final T M) { + final T E = ellipticMeanToEccentric(e, M); + final T v = ellipticEccentricToTrue(e, E); + return v; + } + + /** + * Computes the elliptic mean anomaly from the elliptic true anomaly. + * @param field type + * @param e eccentricity such that 0 ≤ e < 1 + * @param v elliptic true anomaly (rad) + * @return elliptic mean anomaly (rad) + */ + public static > T ellipticTrueToMean(final T e, final T v) { + final T E = ellipticTrueToEccentric(e, v); + final T M = ellipticEccentricToMean(e, E); + return M; + } + + /** + * Computes the elliptic true anomaly from the elliptic eccentric anomaly. + * @param field type + * @param e eccentricity such that 0 ≤ e < 1 + * @param E elliptic eccentric anomaly (rad) + * @return elliptic true anomaly (rad) + */ + public static > T ellipticEccentricToTrue(final T e, final T E) { + final T beta = e.divide(e.multiply(e).negate().add(1).sqrt().add(1)); + final FieldSinCos scE = FastMath.sinCos(E); + return E.add(beta.multiply(scE.sin()).divide(beta.multiply(scE.cos()).subtract(1).negate()).atan().multiply(2)); + } + + /** + * Computes the elliptic eccentric anomaly from the elliptic true anomaly. + * @param field type + * @param e eccentricity such that 0 ≤ e < 1 + * @param v elliptic true anomaly (rad) + * @return elliptic eccentric anomaly (rad) + */ + public static > T ellipticTrueToEccentric(final T e, final T v) { + final T beta = e.divide(e.multiply(e).negate().add(1).sqrt().add(1)); + final FieldSinCos scv = FastMath.sinCos(v); + return v.subtract((beta.multiply(scv.sin()).divide(beta.multiply(scv.cos()).add(1))).atan().multiply(2)); + } + + /** + * Computes the elliptic eccentric anomaly from the elliptic mean anomaly. + *

+ * The algorithm used here for solving hyperbolic Kepler equation is from Odell, + * A.W., Gooding, R.H. "Procedures for solving Kepler's equation." Celestial + * Mechanics 38, 307–334 (1986). https://doi.org/10.1007/BF01238923 + *

+ * @param field type + * @param e eccentricity such that 0 ≤ e < 1 + * @param M elliptic mean anomaly (rad) + * @return elliptic eccentric anomaly (rad) + */ + public static > T ellipticMeanToEccentric(final T e, final T M) { + // reduce M to [-PI PI) interval + final T reducedM = MathUtils.normalizeAngle(M, M.getField().getZero()); + + // compute start value according to A. W. Odell and R. H. Gooding S12 starter + T E; + if (reducedM.abs().getReal() < 1.0 / 6.0) { + if (FastMath.abs(reducedM.getReal()) < Precision.SAFE_MIN) { + // this is an Orekit change to the S12 starter, mainly used when T is some kind + // of derivative structure. If reducedM is 0.0, the derivative of cbrt is + // infinite which induces NaN appearing later in the computation. As in this + // case E and M are almost equal, we initialize E with reducedM + E = reducedM; + } else { + // this is the standard S12 starter + E = reducedM.add(e.multiply((reducedM.multiply(6).cbrt()).subtract(reducedM))); + } + } else { + final T pi = e.getPi(); + if (reducedM.getReal() < 0) { + final T w = reducedM.add(pi); + E = reducedM.add(e.multiply(w.multiply(A).divide(w.negate().add(B)).subtract(pi).subtract(reducedM))); + } else { + final T w = reducedM.negate().add(pi); + E = reducedM + .add(e.multiply(w.multiply(A).divide(w.negate().add(B)).negate().subtract(reducedM).add(pi))); + } + } + + final T e1 = e.negate().add(1); + final boolean noCancellationRisk = (e1.getReal() + E.getReal() * E.getReal() / 6) >= 0.1; + + // perform two iterations, each consisting of one Halley step and one + // Newton-Raphson step + for (int j = 0; j < 2; ++j) { + + final T f; + T fd; + final FieldSinCos scE = FastMath.sinCos(E); + final T fdd = e.multiply(scE.sin()); + final T fddd = e.multiply(scE.cos()); + + if (noCancellationRisk) { + + f = (E.subtract(fdd)).subtract(reducedM); + fd = fddd.negate().add(1); + } else { + + f = eMeSinE(e, E).subtract(reducedM); + final T s = E.multiply(0.5).sin(); + fd = e1.add(e.multiply(s).multiply(s).multiply(2)); + } + final T dee = f.multiply(fd).divide(f.multiply(fdd).multiply(0.5).subtract(fd.multiply(fd))); + + // update eccentric anomaly, using expressions that limit underflow problems + final T w = fd.add(dee.multiply(fdd.add(dee.multiply(fddd.divide(3)))).multiply(0.5)); + fd = fd.add(dee.multiply(fdd.add(dee.multiply(fddd).multiply(0.5)))); + E = E.subtract(f.subtract(dee.multiply(fd.subtract(w))).divide(fd)); + + } + + // expand the result back to original range + E = E.add(M).subtract(reducedM); + return E; + } + + /** + * Accurate computation of E - e sin(E). + *

+ * This method is used when E is close to 0 and e close to 1, i.e. near the + * perigee of almost parabolic orbits + *

+ * @param field type + * @param e eccentricity + * @param E eccentric anomaly (rad) + * @return E - e sin(E) + */ + private static > T eMeSinE(final T e, final T E) { + T x = (e.negate().add(1)).multiply(E.sin()); + final T mE2 = E.negate().multiply(E); + T term = E; + double d = 0; + // the inequality test below IS intentional and should NOT be replaced by a + // check with a small tolerance + for (T x0 = E.getField().getZero().add(Double.NaN); !Double.valueOf(x.getReal()) + .equals(Double.valueOf(x0.getReal()));) { + d += 2; + term = term.multiply(mE2.divide(d * (d + 1))); + x0 = x; + x = x.subtract(term); + } + return x; + } + + /** + * Computes the elliptic mean anomaly from the elliptic eccentric anomaly. + * @param field type + * @param e eccentricity such that 0 ≤ e < 1 + * @param E elliptic eccentric anomaly (rad) + * @return elliptic mean anomaly (rad) + */ + public static > T ellipticEccentricToMean(final T e, final T E) { + return E.subtract(e.multiply(E.sin())); + } + + /** + * Computes the hyperbolic true anomaly from the hyperbolic mean anomaly. + * @param field type + * @param e eccentricity > 1 + * @param M hyperbolic mean anomaly + * @return hyperbolic true anomaly (rad) + */ + public static > T hyperbolicMeanToTrue(final T e, final T M) { + final T H = hyperbolicMeanToEccentric(e, M); + final T v = hyperbolicEccentricToTrue(e, H); + return v; + } + + /** + * Computes the hyperbolic mean anomaly from the hyperbolic true anomaly. + * @param field type + * @param e eccentricity > 1 + * @param v hyperbolic true anomaly (rad) + * @return hyperbolic mean anomaly + */ + public static > T hyperbolicTrueToMean(final T e, final T v) { + final T H = hyperbolicTrueToEccentric(e, v); + final T M = hyperbolicEccentricToMean(e, H); + return M; + } + + /** + * Computes the hyperbolic true anomaly from the hyperbolic eccentric anomaly. + * @param field type + * @param e eccentricity > 1 + * @param H hyperbolic eccentric anomaly + * @return hyperbolic true anomaly (rad) + */ + public static > T hyperbolicEccentricToTrue(final T e, final T H) { + final T s = e.add(1).divide(e.subtract(1)).sqrt(); + final T tanH = H.multiply(0.5).tanh(); + return s.multiply(tanH).atan().multiply(2); + } + + /** + * Computes the hyperbolic eccentric anomaly from the hyperbolic true anomaly. + * @param field type + * @param e eccentricity > 1 + * @param v hyperbolic true anomaly (rad) + * @return hyperbolic eccentric anomaly + */ + public static > T hyperbolicTrueToEccentric(final T e, final T v) { + final FieldSinCos scv = FastMath.sinCos(v); + final T sinhH = e.multiply(e).subtract(1).sqrt().multiply(scv.sin()).divide(e.multiply(scv.cos()).add(1)); + return sinhH.asinh(); + } + + /** + * Computes the hyperbolic eccentric anomaly from the hyperbolic mean anomaly. + *

+ * The algorithm used here for solving hyperbolic Kepler equation is Danby's + * iterative method (3rd order) with Vallado's initial guess. + *

+ * @param field type + * @param e eccentricity > 1 + * @param M hyperbolic mean anomaly + * @return hyperbolic eccentric anomaly + */ + public static > T hyperbolicMeanToEccentric(final T e, final T M) { + // Resolution of hyperbolic Kepler equation for Keplerian parameters + + // Field value of pi + final T pi = e.getPi(); + + // Initial guess + T H; + if (e.getReal() < 1.6) { + if (-pi.getReal() < M.getReal() && M.getReal() < 0. || M.getReal() > pi.getReal()) { + H = M.subtract(e); + } else { + H = M.add(e); + } + } else { + if (e.getReal() < 3.6 && M.abs().getReal() > pi.getReal()) { + H = M.subtract(e.copySign(M)); + } else { + H = M.divide(e.subtract(1)); + } + } + + // Iterative computation + int iter = 0; + do { + final T f3 = e.multiply(H.cosh()); + final T f2 = e.multiply(H.sinh()); + final T f1 = f3.subtract(1); + final T f0 = f2.subtract(H).subtract(M); + final T f12 = f1.multiply(2); + final T d = f0.divide(f12); + final T fdf = f1.subtract(d.multiply(f2)); + final T ds = f0.divide(fdf); + + final T shift = f0.divide(fdf.add(ds.multiply(ds.multiply(f3.divide(6))))); + + H = H.subtract(shift); + + if ((shift.abs().getReal()) <= 1.0e-12) { + return H; + } + + } while (++iter < 50); + + throw new MathIllegalArgumentException(OrekitMessages.UNABLE_TO_COMPUTE_HYPERBOLIC_ECCENTRIC_ANOMALY, iter); + } + + /** + * Computes the hyperbolic mean anomaly from the hyperbolic eccentric anomaly. + * @param field type + * @param e eccentricity > 1 + * @param H hyperbolic eccentric anomaly + * @return hyperbolic mean anomaly + */ + public static > T hyperbolicEccentricToMean(final T e, final T H) { + return e.multiply(H.sinh()).subtract(H); + } + +} diff --git a/src/main/java/org/orekit/orbits/FieldKeplerianOrbit.java b/src/main/java/org/orekit/orbits/FieldKeplerianOrbit.java index 9e2166f04..d881451e7 100644 --- a/src/main/java/org/orekit/orbits/FieldKeplerianOrbit.java +++ b/src/main/java/org/orekit/orbits/FieldKeplerianOrbit.java @@ -24,13 +24,11 @@ import java.util.stream.Stream; import org.hipparchus.CalculusFieldElement; import org.hipparchus.analysis.differentiation.FieldUnivariateDerivative1; import org.hipparchus.analysis.interpolation.FieldHermiteInterpolator; -import org.hipparchus.exception.MathIllegalArgumentException; import org.hipparchus.geometry.euclidean.threed.FieldVector3D; import org.hipparchus.util.FastMath; import org.hipparchus.util.FieldSinCos; import org.hipparchus.util.MathArrays; import org.hipparchus.util.MathUtils; -import org.hipparchus.util.Precision; import org.orekit.errors.OrekitException; import org.orekit.errors.OrekitIllegalArgumentException; import org.orekit.errors.OrekitInternalError; @@ -86,20 +84,6 @@ public class FieldKeplerianOrbit> extends Fiel /** Name of the eccentricity parameter. */ private static final String ECCENTRICITY = "eccentricity"; - /** First coefficient to compute Kepler equation solver starter. */ - private static final double A; - - /** Second coefficient to compute Kepler equation solver starter. */ - private static final double B; - - static { - final double k1 = 3 * FastMath.PI + 2; - final double k2 = FastMath.PI - 1; - final double k3 = 6 * FastMath.PI - 1; - A = 3 * k2 * k2 / k1; - B = k3 * k3 / (6 * k1); - } - /** Semi-major axis (m). */ private final T a; @@ -238,13 +222,13 @@ public class FieldKeplerianOrbit> extends Fiel switch (type) { case MEAN : vUD = (a.getReal() < 0) ? - hyperbolicEccentricToTrue(meanToHyperbolicEccentric(anomalyUD, eUD), eUD) : - ellipticEccentricToTrue(meanToEllipticEccentric(anomalyUD, eUD), eUD); + FieldKeplerianAnomalyUtility.hyperbolicMeanToTrue(eUD, anomalyUD) : + FieldKeplerianAnomalyUtility.ellipticMeanToTrue(eUD, anomalyUD); break; case ECCENTRIC : vUD = (a.getReal() < 0) ? - hyperbolicEccentricToTrue(anomalyUD, eUD) : - ellipticEccentricToTrue(anomalyUD, eUD); + FieldKeplerianAnomalyUtility.hyperbolicEccentricToTrue(eUD, anomalyUD) : + FieldKeplerianAnomalyUtility.ellipticEccentricToTrue(eUD, anomalyUD); break; case TRUE : vUD = anomalyUD; @@ -259,14 +243,14 @@ public class FieldKeplerianOrbit> extends Fiel case MEAN : this.v = (a.getReal() < 0) ? - hyperbolicEccentricToTrue(meanToHyperbolicEccentric(anomaly, e), e) : - ellipticEccentricToTrue(meanToEllipticEccentric(anomaly, e), e); + FieldKeplerianAnomalyUtility.hyperbolicMeanToTrue(e, anomaly) : + FieldKeplerianAnomalyUtility.ellipticMeanToTrue(e, anomaly); break; case ECCENTRIC : this.v = (a.getReal() < 0) ? - hyperbolicEccentricToTrue(anomaly, e) : - ellipticEccentricToTrue(anomaly, e); + FieldKeplerianAnomalyUtility.hyperbolicEccentricToTrue(e, anomaly) : + FieldKeplerianAnomalyUtility.ellipticEccentricToTrue(e, anomaly); break; case TRUE : @@ -367,13 +351,13 @@ public class FieldKeplerianOrbit> extends Fiel final T eSE = FieldVector3D.dotProduct(pvP, pvV).divide(muA.sqrt()); final T eCE = rV2OnMu.subtract(1); e = (eSE.multiply(eSE).add(eCE.multiply(eCE))).sqrt(); - v = ellipticEccentricToTrue(eSE.atan2(eCE), e); + v = FieldKeplerianAnomalyUtility.ellipticEccentricToTrue(e, eSE.atan2(eCE)); } else { // hyperbolic orbit final T eSH = FieldVector3D.dotProduct(pvP, pvV).divide(muA.negate().sqrt()); final T eCH = rV2OnMu.subtract(1); e = (m2.negate().divide(muA).add(1)).sqrt(); - v = hyperbolicEccentricToTrue((eCH.add(eSH)).divide(eCH.subtract(eSH)).log().divide(2), e); + v = FieldKeplerianAnomalyUtility.hyperbolicEccentricToTrue(e, (eCH.add(eSH)).divide(eCH.subtract(eSH)).log().divide(2)); } // Checking eccentricity range @@ -411,8 +395,8 @@ public class FieldKeplerianOrbit> extends Fiel final FieldUnivariateDerivative1 eUD = new FieldUnivariateDerivative1<>(e, eDot); final FieldUnivariateDerivative1 MUD = new FieldUnivariateDerivative1<>(getMeanAnomaly(), MDot); final FieldUnivariateDerivative1 vUD = (a.getReal() < 0) ? - FieldKeplerianOrbit.hyperbolicEccentricToTrue(FieldKeplerianOrbit.meanToHyperbolicEccentric(MUD, eUD), eUD) : - FieldKeplerianOrbit.ellipticEccentricToTrue(FieldKeplerianOrbit.meanToEllipticEccentric(MUD, eUD), eUD); + FieldKeplerianAnomalyUtility.hyperbolicMeanToTrue(eUD, MUD) : + FieldKeplerianAnomalyUtility.ellipticMeanToTrue(eUD, MUD); vDot = vUD.getDerivative(1); } else { @@ -547,7 +531,9 @@ public class FieldKeplerianOrbit> extends Fiel * @return eccentric anomaly (rad) */ public T getEccentricAnomaly() { - return (a.getReal() < 0) ? trueToHyperbolicEccentric(v, e) : trueToEllipticEccentric(v, e); + return (a.getReal() < 0) ? + FieldKeplerianAnomalyUtility.hyperbolicTrueToEccentric(e, v) : + FieldKeplerianAnomalyUtility.ellipticTrueToEccentric(e, v); } /** Get the eccentric anomaly derivative. @@ -565,8 +551,8 @@ public class FieldKeplerianOrbit> extends Fiel final FieldUnivariateDerivative1 eUD = new FieldUnivariateDerivative1<>(e, eDot); final FieldUnivariateDerivative1 vUD = new FieldUnivariateDerivative1<>(v, vDot); final FieldUnivariateDerivative1 EUD = (a.getReal() < 0) ? - trueToHyperbolicEccentric(vUD, eUD) : - trueToEllipticEccentric(vUD, eUD); + FieldKeplerianAnomalyUtility.hyperbolicTrueToEccentric(eUD, vUD) : + FieldKeplerianAnomalyUtility.ellipticTrueToEccentric(eUD, vUD); return EUD.getDerivative(1); } @@ -576,8 +562,8 @@ public class FieldKeplerianOrbit> extends Fiel */ public T getMeanAnomaly() { return (a.getReal() < 0) ? - hyperbolicEccentricToMean(trueToHyperbolicEccentric(v, e), e) : - ellipticEccentricToMean(trueToEllipticEccentric(v, e), e); + FieldKeplerianAnomalyUtility.hyperbolicTrueToMean(e, v) : + FieldKeplerianAnomalyUtility.ellipticTrueToMean(e, v); } /** Get the mean anomaly derivative. @@ -595,8 +581,8 @@ public class FieldKeplerianOrbit> extends Fiel final FieldUnivariateDerivative1 eUD = new FieldUnivariateDerivative1<>(e, eDot); final FieldUnivariateDerivative1 vUD = new FieldUnivariateDerivative1<>(v, vDot); final FieldUnivariateDerivative1 MUD = (a.getReal() < 0) ? - hyperbolicEccentricToMean(trueToHyperbolicEccentric(vUD, eUD), eUD) : - ellipticEccentricToMean(trueToEllipticEccentric(vUD, eUD), eUD); + FieldKeplerianAnomalyUtility.hyperbolicTrueToMean(eUD, vUD) : + FieldKeplerianAnomalyUtility.ellipticTrueToMean(eUD, vUD); return MUD.getDerivative(1); } @@ -635,11 +621,11 @@ public class FieldKeplerianOrbit> extends Fiel * @param e eccentricity * @param type of the field elements * @return v the true anomaly + * @deprecated As of 11.3, replaced by + * {@link FieldKeplerianAnomalyUtility#ellipticEccentricToMean(CalculusFieldElement, CalculusFieldElement)}. */ public static > T ellipticEccentricToTrue(final T E, final T e) { - final T beta = e.divide(e.multiply(e).negate().add(1).sqrt().add(1)); - final FieldSinCos scE = FastMath.sinCos(E); - return E.add(beta.multiply(scE.sin()).divide(beta.multiply(scE.cos()).subtract(1).negate()).atan().multiply(2)); + return FieldKeplerianAnomalyUtility.ellipticEccentricToMean(e, E); } /** Computes the elliptic eccentric anomaly from the true anomaly. @@ -647,11 +633,11 @@ public class FieldKeplerianOrbit> extends Fiel * @param e eccentricity * @param type of the field elements * @return E the elliptic eccentric anomaly + * @deprecated As of 11.3, replaced by + * {@link FieldKeplerianAnomalyUtility#ellipticTrueToEccentric(CalculusFieldElement, CalculusFieldElement)}. */ public static > T trueToEllipticEccentric(final T v, final T e) { - final T beta = e.divide(e.multiply(e).negate().add(1).sqrt().add(1)); - final FieldSinCos scv = FastMath.sinCos(v); - return v.subtract((beta.multiply(scv.sin()).divide(beta.multiply(scv.cos()).add(1))).atan().multiply(2)); + return FieldKeplerianAnomalyUtility.ellipticTrueToEccentric(e, v); } /** Computes the true anomaly from the hyperbolic eccentric anomaly. @@ -659,11 +645,11 @@ public class FieldKeplerianOrbit> extends Fiel * @param e eccentricity * @param type of the field elements * @return v the true anomaly + * @deprecated As of 11.3, replaced by + * {@link FieldKeplerianAnomalyUtility#hyperbolicEccentricToTrue(CalculusFieldElement, CalculusFieldElement)}. */ public static > T hyperbolicEccentricToTrue(final T H, final T e) { - final T s = e.add(1).divide(e.subtract(1)).sqrt(); - final T tanH = H.multiply(0.5).tanh(); - return s.multiply(tanH).atan().multiply(2); + return FieldKeplerianAnomalyUtility.hyperbolicEccentricToTrue(e, H); } /** Computes the hyperbolic eccentric anomaly from the true anomaly. @@ -671,11 +657,11 @@ public class FieldKeplerianOrbit> extends Fiel * @param e eccentricity * @param type of the field elements * @return H the hyperbolic eccentric anomaly + * @deprecated As of 11.3, replaced by + * {@link FieldKeplerianAnomalyUtility#hyperbolicTrueToEccentric(CalculusFieldElement, CalculusFieldElement)}. */ public static > T trueToHyperbolicEccentric(final T v, final T e) { - final FieldSinCos scv = FastMath.sinCos(v); - final T sinhH = e.multiply(e).subtract(1).sqrt().multiply(scv.sin()).divide(e.multiply(scv.cos()).add(1)); - return sinhH.asinh(); + return FieldKeplerianAnomalyUtility.hyperbolicTrueToEccentric(e, v); } /** Computes the mean anomaly from the hyperbolic eccentric anomaly. @@ -683,169 +669,35 @@ public class FieldKeplerianOrbit> extends Fiel * @param e eccentricity * @param type of the field elements * @return M the mean anomaly + * @deprecated As of 11.3, replaced by + * {@link FieldKeplerianAnomalyUtility#hyperbolicEccentricToMean(CalculusFieldElement, CalculusFieldElement)}. */ public static > T hyperbolicEccentricToMean(final T H, final T e) { - return e.multiply(H.sinh()).subtract(H); + return FieldKeplerianAnomalyUtility.hyperbolicEccentricToMean(e, H); } /** Computes the elliptic eccentric anomaly from the mean anomaly. - *

- * The algorithm used here for solving Kepler equation has been published - * in: "Procedures for solving Kepler's Equation", A. W. Odell and - * R. H. Gooding, Celestial Mechanics 38 (1986) 307-334 - *

* @param M mean anomaly (rad) * @param e eccentricity * @param type of the field elements * @return E the eccentric anomaly + * @deprecated As of 11.3, replaced by + * {@link FieldKeplerianAnomalyUtility#ellipticMeanToEccentric(CalculusFieldElement, CalculusFieldElement)}. */ public static > T meanToEllipticEccentric(final T M, final T e) { - // reduce M to [-PI PI) interval - final T reducedM = MathUtils.normalizeAngle(M, M.getField().getZero()); - - // compute start value according to A. W. Odell and R. H. Gooding S12 starter - T E; - if (reducedM.abs().getReal() < 1.0 / 6.0) { - if (FastMath.abs(reducedM.getReal()) < Precision.SAFE_MIN) { - // this is an Orekit change to the S12 starter, mainly used when T is some kind of derivative structure. - // If reducedM is 0.0, the derivative of cbrt is infinite which induces NaN appearing later in - // the computation. As in this case E and M are almost equal, we initialize E with reducedM - E = reducedM; - } else { - // this is the standard S12 starter - E = reducedM.add(e.multiply( (reducedM.multiply(6).cbrt()).subtract(reducedM))); - } - } else { - final T pi = e.getPi(); - if (reducedM.getReal() < 0) { - final T w = reducedM.add(pi); - E = reducedM.add(e.multiply(w.multiply(A).divide(w.negate().add(B)).subtract(pi).subtract(reducedM))); - } else { - final T w = reducedM.negate().add(pi); - E = reducedM.add(e.multiply(w.multiply(A).divide(w.negate().add(B)).negate().subtract(reducedM).add(pi))); - } - } - - final T e1 = e.negate().add(1); - final boolean noCancellationRisk = (e1.getReal() + E.getReal() * E.getReal() / 6) >= 0.1; - - // perform two iterations, each consisting of one Halley step and one Newton-Raphson step - for (int j = 0; j < 2; ++j) { - - final T f; - T fd; - final FieldSinCos scE = FastMath.sinCos(E); - final T fdd = e.multiply(scE.sin()); - final T fddd = e.multiply(scE.cos()); - - if (noCancellationRisk) { - - f = (E.subtract(fdd)).subtract(reducedM); - fd = fddd.negate().add(1); - } else { - - - f = eMeSinE(E, e).subtract(reducedM); - final T s = E.multiply(0.5).sin(); - fd = e1.add(e.multiply(s).multiply(s).multiply(2)); - } - final T dee = f.multiply(fd).divide(f.multiply(fdd).multiply(0.5).subtract(fd.multiply(fd))); - - // update eccentric anomaly, using expressions that limit underflow problems - final T w = fd.add(dee.multiply(fdd.add(dee.multiply(fddd.divide(3)))).multiply(0.5)); - fd = fd.add(dee.multiply(fdd.add(dee.multiply(fddd).multiply(0.5)))); - E = E.subtract(f.subtract(dee.multiply(fd.subtract(w))).divide(fd)); - - } - - // expand the result back to original range - E = E.add(M).subtract(reducedM); - return E; - } - - /** Accurate computation of E - e sin(E). - *

- * This method is used when E is close to 0 and e close to 1, - * i.e. near the perigee of almost parabolic orbits - *

- * @param E eccentric anomaly - * @param e eccentricity - * @param Type of the field elements - * @return E - e sin(E) - */ - private static > T eMeSinE(final T E, final T e) { - - T x = (e.negate().add(1)).multiply(E.sin()); - final T mE2 = E.negate().multiply(E); - T term = E; - double d = 0; - // the inequality test below IS intentional and should NOT be replaced by a check with a small tolerance - for (T x0 = E.getField().getZero().add(Double.NaN); !Double.valueOf(x.getReal()).equals(Double.valueOf(x0.getReal()));) { - d += 2; - term = term.multiply(mE2.divide(d * (d + 1))); - x0 = x; - x = x.subtract(term); - } - return x; + return FieldKeplerianAnomalyUtility.ellipticMeanToEccentric(e, M); } /** Computes the hyperbolic eccentric anomaly from the mean anomaly. - *

- * The algorithm used here for solving hyperbolic Kepler equation is - * Danby's iterative method (3rd order) with Vallado's initial guess. - *

* @param M mean anomaly (rad) * @param e eccentricity * @param Type of the field elements * @return H the hyperbolic eccentric anomaly + * @deprecated As of 11.3, replaced by + * {@link FieldKeplerianAnomalyUtility#hyperbolicMeanToEccentric(CalculusFieldElement, CalculusFieldElement)}. */ public static > T meanToHyperbolicEccentric(final T M, final T e) { - - // Resolution of hyperbolic Kepler equation for Keplerian parameters - - // Field value of pi - final T pi = e.getPi(); - - // Initial guess - T H; - if (e.getReal() < 1.6) { - if (-pi.getReal() < M.getReal() && M.getReal() < 0. || M.getReal() > pi.getReal()) { - H = M.subtract(e); - } else { - H = M.add(e); - } - } else { - if (e.getReal() < 3.6 && M.abs().getReal() > pi.getReal()) { - H = M.subtract(e.copySign(M)); - } else { - H = M.divide(e.subtract(1)); - } - } - - // Iterative computation - int iter = 0; - do { - final T f3 = e.multiply(H.cosh()); - final T f2 = e.multiply(H.sinh()); - final T f1 = f3.subtract(1); - final T f0 = f2.subtract(H).subtract(M); - final T f12 = f1.multiply(2); - final T d = f0.divide(f12); - final T fdf = f1.subtract(d.multiply(f2)); - final T ds = f0.divide(fdf); - - final T shift = f0.divide(fdf.add(ds.multiply(ds.multiply(f3.divide(6))))); - - H = H.subtract(shift); - - if ((shift.abs().getReal()) <= 1.0e-12) { - return H; - } - - } while (++iter < 50); - - throw new MathIllegalArgumentException(OrekitMessages.UNABLE_TO_COMPUTE_HYPERBOLIC_ECCENTRIC_ANOMALY, - iter); + return FieldKeplerianAnomalyUtility.hyperbolicMeanToEccentric(e, M); } /** Computes the mean anomaly from the elliptic eccentric anomaly. @@ -853,9 +705,11 @@ public class FieldKeplerianOrbit> extends Fiel * @param e eccentricity * @param type of the field elements * @return M the mean anomaly + * @deprecated As of 11.3, replaced by + * {@link FieldKeplerianAnomalyUtility#ellipticEccentricToMean(CalculusFieldElement, CalculusFieldElement)}. */ public static > T ellipticEccentricToMean(final T E, final T e) { - return E.subtract(e.multiply(E.sin())); + return FieldKeplerianAnomalyUtility.ellipticEccentricToMean(e, E); } /** {@inheritDoc} */ diff --git a/src/main/java/org/orekit/orbits/KeplerianAnomalyUtility.java b/src/main/java/org/orekit/orbits/KeplerianAnomalyUtility.java new file mode 100644 index 000000000..5120cfdc7 --- /dev/null +++ b/src/main/java/org/orekit/orbits/KeplerianAnomalyUtility.java @@ -0,0 +1,306 @@ +/* Copyright 2002-2022 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 org.hipparchus.exception.MathIllegalStateException; +import org.hipparchus.util.FastMath; +import org.hipparchus.util.MathUtils; +import org.hipparchus.util.SinCos; +import org.orekit.errors.OrekitMessages; + +/** + * Utility methods for converting between different Keplerian anomalies. + * @author Luc Maisonobe + * @author Andrew Goetz + */ +public final class KeplerianAnomalyUtility { + + /** First coefficient to compute elliptic Kepler equation solver starter. */ + private static final double A; + + /** Second coefficient to compute elliptic Kepler equation solver starter. */ + private static final double B; + + static { + final double k1 = 3 * FastMath.PI + 2; + final double k2 = FastMath.PI - 1; + final double k3 = 6 * FastMath.PI - 1; + A = 3 * k2 * k2 / k1; + B = k3 * k3 / (6 * k1); + } + + private KeplerianAnomalyUtility() { + } + + /** + * Computes the elliptic true anomaly from the elliptic mean anomaly. + * @param e eccentricity such that 0 ≤ e < 1 + * @param M elliptic mean anomaly (rad) + * @return elliptic true anomaly (rad) + */ + public static double ellipticMeanToTrue(final double e, final double M) { + final double E = ellipticMeanToEccentric(e, M); + final double v = ellipticEccentricToTrue(e, E); + return v; + } + + /** + * Computes the elliptic mean anomaly from the elliptic true anomaly. + * @param e eccentricity such that 0 ≤ e < 1 + * @param v elliptic true anomaly (rad) + * @return elliptic mean anomaly (rad) + */ + public static double ellipticTrueToMean(final double e, final double v) { + final double E = ellipticTrueToEccentric(e, v); + final double M = ellipticEccentricToMean(e, E); + return M; + } + + /** + * Computes the elliptic true anomaly from the elliptic eccentric anomaly. + * @param e eccentricity such that 0 ≤ e < 1 + * @param E elliptic eccentric anomaly (rad) + * @return elliptic true anomaly (rad) + */ + public static double ellipticEccentricToTrue(final double e, final double E) { + final double beta = e / (1 + FastMath.sqrt((1 - e) * (1 + e))); + final SinCos scE = FastMath.sinCos(E); + return E + 2 * FastMath.atan(beta * scE.sin() / (1 - beta * scE.cos())); + } + + /** + * Computes the elliptic eccentric anomaly from the elliptic true anomaly. + * @param e eccentricity such that 0 ≤ e < 1 + * @param v elliptic true anomaly (rad) + * @return elliptic eccentric anomaly (rad) + */ + public static double ellipticTrueToEccentric(final double e, final double v) { + final double beta = e / (1 + FastMath.sqrt(1 - e * e)); + final SinCos scv = FastMath.sinCos(v); + return v - 2 * FastMath.atan(beta * scv.sin() / (1 + beta * scv.cos())); + } + + /** + * Computes the elliptic eccentric anomaly from the elliptic mean anomaly. + *

+ * The algorithm used here for solving hyperbolic Kepler equation is from Odell, + * A.W., Gooding, R.H. "Procedures for solving Kepler's equation." Celestial + * Mechanics 38, 307–334 (1986). https://doi.org/10.1007/BF01238923 + *

+ * @param e eccentricity such that 0 ≤ e < 1 + * @param M elliptic mean anomaly (rad) + * @return elliptic eccentric anomaly (rad) + */ + public static double ellipticMeanToEccentric(final double e, final double M) { + // reduce M to [-PI PI) interval + final double reducedM = MathUtils.normalizeAngle(M, 0.0); + + // compute start value according to A. W. Odell and R. H. Gooding S12 starter + double E; + if (FastMath.abs(reducedM) < 1.0 / 6.0) { + E = reducedM + e * (FastMath.cbrt(6 * reducedM) - reducedM); + } else { + if (reducedM < 0) { + final double w = FastMath.PI + reducedM; + E = reducedM + e * (A * w / (B - w) - FastMath.PI - reducedM); + } else { + final double w = FastMath.PI - reducedM; + E = reducedM + e * (FastMath.PI - A * w / (B - w) - reducedM); + } + } + + final double e1 = 1 - e; + final boolean noCancellationRisk = (e1 + E * E / 6) >= 0.1; + + // perform two iterations, each consisting of one Halley step and one + // Newton-Raphson step + for (int j = 0; j < 2; ++j) { + final double f; + double fd; + final SinCos sc = FastMath.sinCos(E); + final double fdd = e * sc.sin(); + final double fddd = e * sc.cos(); + if (noCancellationRisk) { + f = (E - fdd) - reducedM; + fd = 1 - fddd; + } else { + f = eMeSinE(e, E) - reducedM; + final double s = FastMath.sin(0.5 * E); + fd = e1 + 2 * e * s * s; + } + final double dee = f * fd / (0.5 * f * fdd - fd * fd); + + // update eccentric anomaly, using expressions that limit underflow problems + final double w = fd + 0.5 * dee * (fdd + dee * fddd / 3); + fd += dee * (fdd + 0.5 * dee * fddd); + E -= (f - dee * (fd - w)) / fd; + + } + + // expand the result back to original range + E += M - reducedM; + + return E; + } + + /** + * Accurate computation of E - e sin(E). + *

+ * This method is used when E is close to 0 and e close to 1, i.e. near the + * perigee of almost parabolic orbits + *

+ * @param e eccentricity + * @param E eccentric anomaly (rad) + * @return E - e sin(E) + */ + private static double eMeSinE(final double e, final double E) { + double x = (1 - e) * FastMath.sin(E); + final double mE2 = -E * E; + double term = E; + double d = 0; + // the inequality test below IS intentional and should NOT be replaced by a + // check with a small tolerance + for (double x0 = Double.NaN; !Double.valueOf(x).equals(Double.valueOf(x0));) { + d += 2; + term *= mE2 / (d * (d + 1)); + x0 = x; + x = x - term; + } + return x; + } + + /** + * Computes the elliptic mean anomaly from the elliptic eccentric anomaly. + * @param e eccentricity such that 0 ≤ e < 1 + * @param E elliptic eccentric anomaly (rad) + * @return elliptic mean anomaly (rad) + */ + public static double ellipticEccentricToMean(final double e, final double E) { + return E - e * FastMath.sin(E); + } + + /** + * Computes the hyperbolic true anomaly from the hyperbolic mean anomaly. + * @param e eccentricity > 1 + * @param M hyperbolic mean anomaly + * @return hyperbolic true anomaly (rad) + */ + public static double hyperbolicMeanToTrue(final double e, final double M) { + final double H = hyperbolicMeanToEccentric(e, M); + final double v = hyperbolicEccentricToTrue(e, H); + return v; + } + + /** + * Computes the hyperbolic mean anomaly from the hyperbolic true anomaly. + * @param e eccentricity > 1 + * @param v hyperbolic true anomaly (rad) + * @return hyperbolic mean anomaly + */ + public static double hyperbolicTrueToMean(final double e, final double v) { + final double H = hyperbolicTrueToEccentric(e, v); + final double M = hyperbolicEccentricToMean(e, H); + return M; + } + + /** + * Computes the hyperbolic true anomaly from the hyperbolic eccentric anomaly. + * @param e eccentricity > 1 + * @param H hyperbolic eccentric anomaly + * @return hyperbolic true anomaly (rad) + */ + public static double hyperbolicEccentricToTrue(final double e, final double H) { + return 2 * FastMath.atan(FastMath.sqrt((e + 1) / (e - 1)) * FastMath.tanh(0.5 * H)); + } + + /** + * Computes the hyperbolic eccentric anomaly from the hyperbolic true anomaly. + * @param e eccentricity > 1 + * @param v hyperbolic true anomaly (rad) + * @return hyperbolic eccentric anomaly + */ + public static double hyperbolicTrueToEccentric(final double e, final double v) { + final SinCos scv = FastMath.sinCos(v); + final double sinhH = FastMath.sqrt(e * e - 1) * scv.sin() / (1 + e * scv.cos()); + return FastMath.asinh(sinhH); + } + + /** + * Computes the hyperbolic eccentric anomaly from the hyperbolic mean anomaly. + *

+ * The algorithm used here for solving hyperbolic Kepler equation is Danby's + * iterative method (3rd order) with Vallado's initial guess. + *

+ * @param e eccentricity > 1 + * @param M hyperbolic mean anomaly + * @return hyperbolic eccentric anomaly + */ + public static double hyperbolicMeanToEccentric(final double e, final double M) { + // Resolution of hyperbolic Kepler equation for Keplerian parameters + + // Initial guess + double H; + if (e < 1.6) { + if (-FastMath.PI < M && M < 0. || M > FastMath.PI) { + H = M - e; + } else { + H = M + e; + } + } else { + if (e < 3.6 && FastMath.abs(M) > FastMath.PI) { + H = M - FastMath.copySign(e, M); + } else { + H = M / (e - 1.); + } + } + + // Iterative computation + int iter = 0; + do { + final double f3 = e * FastMath.cosh(H); + final double f2 = e * 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); + } + + /** + * Computes the hyperbolic mean anomaly from the hyperbolic eccentric anomaly. + * @param e eccentricity > 1 + * @param H hyperbolic eccentric anomaly + * @return hyperbolic mean anomaly + */ + public static double hyperbolicEccentricToMean(final double e, final double H) { + return e * FastMath.sinh(H) - H; + } + +} diff --git a/src/main/java/org/orekit/orbits/KeplerianOrbit.java b/src/main/java/org/orekit/orbits/KeplerianOrbit.java index 5ad29fa60..00f585de2 100644 --- a/src/main/java/org/orekit/orbits/KeplerianOrbit.java +++ b/src/main/java/org/orekit/orbits/KeplerianOrbit.java @@ -23,7 +23,6 @@ import java.util.stream.Stream; import org.hipparchus.analysis.differentiation.UnivariateDerivative1; import org.hipparchus.analysis.interpolation.HermiteInterpolator; -import org.hipparchus.exception.MathIllegalStateException; import org.hipparchus.geometry.euclidean.threed.Vector3D; import org.hipparchus.util.FastMath; import org.hipparchus.util.MathUtils; @@ -87,20 +86,6 @@ public class KeplerianOrbit extends Orbit { /** Name of the eccentricity parameter. */ private static final String ECCENTRICITY = "eccentricity"; - /** First coefficient to compute Kepler equation solver starter. */ - private static final double A; - - /** Second coefficient to compute Kepler equation solver starter. */ - private static final double B; - - static { - final double k1 = 3 * FastMath.PI + 2; - final double k2 = FastMath.PI - 1; - final double k3 = 6 * FastMath.PI - 1; - A = 3 * k2 * k2 / k1; - B = k3 * k3 / (6 * k1); - } - /** Semi-major axis (m). */ private final double a; @@ -223,13 +208,13 @@ public class KeplerianOrbit extends Orbit { switch (type) { case MEAN : vUD = (a < 0) ? - FieldKeplerianOrbit.hyperbolicEccentricToTrue(FieldKeplerianOrbit.meanToHyperbolicEccentric(anomalyUD, eUD), eUD) : - FieldKeplerianOrbit.ellipticEccentricToTrue(FieldKeplerianOrbit.meanToEllipticEccentric(anomalyUD, eUD), eUD); + FieldKeplerianAnomalyUtility.hyperbolicMeanToTrue(eUD, anomalyUD) : + FieldKeplerianAnomalyUtility.ellipticMeanToTrue(eUD, anomalyUD); break; case ECCENTRIC : vUD = (a < 0) ? - FieldKeplerianOrbit.hyperbolicEccentricToTrue(anomalyUD, eUD) : - FieldKeplerianOrbit.ellipticEccentricToTrue(anomalyUD, eUD); + FieldKeplerianAnomalyUtility.hyperbolicEccentricToTrue(eUD, anomalyUD) : + FieldKeplerianAnomalyUtility.ellipticEccentricToTrue(eUD, anomalyUD); break; case TRUE : vUD = anomalyUD; @@ -243,13 +228,13 @@ public class KeplerianOrbit extends Orbit { switch (type) { case MEAN : this.v = (a < 0) ? - hyperbolicEccentricToTrue(meanToHyperbolicEccentric(anomaly, e), e) : - ellipticEccentricToTrue(meanToEllipticEccentric(anomaly, e), e); + KeplerianAnomalyUtility.hyperbolicMeanToTrue(e, anomaly) : + KeplerianAnomalyUtility.ellipticMeanToTrue(e, anomaly); break; case ECCENTRIC : this.v = (a < 0) ? - hyperbolicEccentricToTrue(anomaly, e) : - ellipticEccentricToTrue(anomaly, e); + KeplerianAnomalyUtility.hyperbolicEccentricToTrue(e, anomaly) : + KeplerianAnomalyUtility.ellipticEccentricToTrue(e, anomaly); break; case TRUE : this.v = anomaly; @@ -339,13 +324,13 @@ public class KeplerianOrbit extends Orbit { final double eSE = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(muA); final double eCE = rV2OnMu - 1; e = FastMath.sqrt(eSE * eSE + eCE * eCE); - v = ellipticEccentricToTrue(FastMath.atan2(eSE, eCE), e); + v = KeplerianAnomalyUtility.ellipticEccentricToTrue(e, FastMath.atan2(eSE, eCE)); } else { // hyperbolic orbit final double eSH = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(-muA); final double eCH = rV2OnMu - 1; e = FastMath.sqrt(1 - m2 / muA); - v = hyperbolicEccentricToTrue(FastMath.log((eCH + eSH) / (eCH - eSH)) / 2, e); + v = KeplerianAnomalyUtility.hyperbolicEccentricToTrue(e, FastMath.log((eCH + eSH) / (eCH - eSH)) / 2); } // Checking eccentricity range @@ -383,8 +368,8 @@ public class KeplerianOrbit extends Orbit { final UnivariateDerivative1 eUD = new UnivariateDerivative1(e, eDot); final UnivariateDerivative1 MUD = new UnivariateDerivative1(getMeanAnomaly(), MDot); final UnivariateDerivative1 vUD = (a < 0) ? - FieldKeplerianOrbit.hyperbolicEccentricToTrue(FieldKeplerianOrbit.meanToHyperbolicEccentric(MUD, eUD), eUD) : - FieldKeplerianOrbit.ellipticEccentricToTrue(FieldKeplerianOrbit.meanToEllipticEccentric(MUD, eUD), eUD); + FieldKeplerianAnomalyUtility.hyperbolicMeanToTrue(eUD, MUD) : + FieldKeplerianAnomalyUtility.ellipticMeanToTrue(eUD, MUD); vDot = vUD.getDerivative(1); } else { @@ -518,7 +503,9 @@ public class KeplerianOrbit extends Orbit { * @return eccentric anomaly (rad) */ public double getEccentricAnomaly() { - return (a < 0) ? trueToHyperbolicEccentric(v, e) : trueToEllipticEccentric(v, e); + return (a < 0) ? + KeplerianAnomalyUtility.hyperbolicTrueToEccentric(e, v) : + KeplerianAnomalyUtility.ellipticTrueToEccentric(e, v); } /** Get the eccentric anomaly derivative. @@ -529,8 +516,8 @@ public class KeplerianOrbit extends Orbit { final UnivariateDerivative1 eUD = new UnivariateDerivative1(e, eDot); final UnivariateDerivative1 vUD = new UnivariateDerivative1(v, vDot); final UnivariateDerivative1 EUD = (a < 0) ? - FieldKeplerianOrbit.trueToHyperbolicEccentric(vUD, eUD) : - FieldKeplerianOrbit.trueToEllipticEccentric(vUD, eUD); + FieldKeplerianAnomalyUtility.hyperbolicTrueToEccentric(eUD, vUD) : + FieldKeplerianAnomalyUtility.ellipticTrueToEccentric(eUD, vUD); return EUD.getDerivative(1); } @@ -539,8 +526,8 @@ public class KeplerianOrbit extends Orbit { */ public double getMeanAnomaly() { return (a < 0) ? - hyperbolicEccentricToMean(trueToHyperbolicEccentric(v, e), e) : - ellipticEccentricToMean(trueToEllipticEccentric(v, e), e); + KeplerianAnomalyUtility.hyperbolicTrueToMean(e, v) : + KeplerianAnomalyUtility.ellipticTrueToMean(e, v); } /** Get the mean anomaly derivative. @@ -551,8 +538,8 @@ public class KeplerianOrbit extends Orbit { final UnivariateDerivative1 eUD = new UnivariateDerivative1(e, eDot); final UnivariateDerivative1 vUD = new UnivariateDerivative1(v, vDot); final UnivariateDerivative1 MUD = (a < 0) ? - FieldKeplerianOrbit.hyperbolicEccentricToMean(FieldKeplerianOrbit.trueToHyperbolicEccentric(vUD, eUD), eUD) : - FieldKeplerianOrbit.ellipticEccentricToMean(FieldKeplerianOrbit.trueToEllipticEccentric(vUD, eUD), eUD); + FieldKeplerianAnomalyUtility.hyperbolicTrueToMean(eUD, vUD) : + FieldKeplerianAnomalyUtility.ellipticTrueToMean(eUD, vUD); return MUD.getDerivative(1); } @@ -581,207 +568,70 @@ public class KeplerianOrbit extends Orbit { * @param E eccentric anomaly (rad) * @param e eccentricity * @return v the true anomaly - * @since 9.0 + * @deprecated As of 11.3, replaced by {@link KeplerianAnomalyUtility#ellipticEccentricToTrue(double, double)}. */ public static double ellipticEccentricToTrue(final double E, final double e) { - final double beta = e / (1 + FastMath.sqrt((1 - e) * (1 + e))); - final SinCos scE = FastMath.sinCos(E); - return E + 2 * FastMath.atan(beta * scE.sin() / (1 - beta * scE.cos())); + return KeplerianAnomalyUtility.ellipticEccentricToTrue(e, E); } /** Computes the elliptic eccentric anomaly from the true anomaly. * @param v true anomaly (rad) * @param e eccentricity * @return E the elliptic eccentric anomaly - * @since 9.0 + * @deprecated As of 11.3, replaced by {@link KeplerianAnomalyUtility#ellipticTrueToEccentric(double, double)}. */ public static double trueToEllipticEccentric(final double v, final double e) { - final double beta = e / (1 + FastMath.sqrt(1 - e * e)); - final SinCos scv = FastMath.sinCos(v); - return v - 2 * FastMath.atan(beta * scv.sin() / (1 + beta * scv.cos())); + return KeplerianAnomalyUtility.ellipticTrueToEccentric(e, v); } /** Computes the true anomaly from the hyperbolic eccentric anomaly. * @param H hyperbolic eccentric anomaly (rad) * @param e eccentricity * @return v the true anomaly + * @deprecated As of 11.3, replaced by {@link KeplerianAnomalyUtility#hyperbolicEccentricToTrue(double, double)}. */ public static double hyperbolicEccentricToTrue(final double H, final double e) { - return 2 * FastMath.atan(FastMath.sqrt((e + 1) / (e - 1)) * FastMath.tanh(H / 2)); + return KeplerianAnomalyUtility.hyperbolicEccentricToTrue(e, H); } /** Computes the hyperbolic eccentric anomaly from the true anomaly. * @param v true anomaly (rad) * @param e eccentricity * @return H the hyperbolic eccentric anomaly - * @since 9.0 + * @deprecated As of 11.3, replaced by {@link KeplerianAnomalyUtility#hyperbolicTrueToEccentric(double, double)}. */ public static double trueToHyperbolicEccentric(final double v, final double e) { - final SinCos scv = FastMath.sinCos(v); - final double sinhH = FastMath.sqrt(e * e - 1) * scv.sin() / (1 + e * scv.cos()); - return FastMath.asinh(sinhH); + return KeplerianAnomalyUtility.hyperbolicTrueToEccentric(e, v); } /** Computes the elliptic eccentric anomaly from the mean anomaly. - *

- * The algorithm used here for solving Kepler equation has been published - * in: "Procedures for solving Kepler's Equation", A. W. Odell and - * R. H. Gooding, Celestial Mechanics 38 (1986) 307-334 - *

* @param M mean anomaly (rad) * @param e eccentricity * @return E the eccentric anomaly + * @deprecated As of 11.3, replaced by {@link KeplerianAnomalyUtility#ellipticMeanToEccentric(double, double)}. */ public static double meanToEllipticEccentric(final double M, final double e) { - - // reduce M to [-PI PI) interval - final double reducedM = MathUtils.normalizeAngle(M, 0.0); - - // compute start value according to A. W. Odell and R. H. Gooding S12 starter - double E; - if (FastMath.abs(reducedM) < 1.0 / 6.0) { - E = reducedM + e * (FastMath.cbrt(6 * reducedM) - reducedM); - } else { - if (reducedM < 0) { - final double w = FastMath.PI + reducedM; - E = reducedM + e * (A * w / (B - w) - FastMath.PI - reducedM); - } else { - final double w = FastMath.PI - reducedM; - E = reducedM + e * (FastMath.PI - A * w / (B - w) - reducedM); - } - } - - final double e1 = 1 - e; - final boolean noCancellationRisk = (e1 + E * E / 6) >= 0.1; - - // perform two iterations, each consisting of one Halley step and one Newton-Raphson step - for (int j = 0; j < 2; ++j) { - final double f; - double fd; - final SinCos sc = FastMath.sinCos(E); - final double fdd = e * sc.sin(); - final double fddd = e * sc.cos(); - if (noCancellationRisk) { - f = (E - fdd) - reducedM; - fd = 1 - fddd; - } else { - f = eMeSinE(E, e) - reducedM; - final double s = FastMath.sin(0.5 * E); - fd = e1 + 2 * e * s * s; - } - final double dee = f * fd / (0.5 * f * fdd - fd * fd); - - // update eccentric anomaly, using expressions that limit underflow problems - final double w = fd + 0.5 * dee * (fdd + dee * fddd / 3); - fd += dee * (fdd + 0.5 * dee * fddd); - E -= (f - dee * (fd - w)) / fd; - - } - - // expand the result back to original range - E += M - reducedM; - - return E; - - } - - /** Accurate computation of E - e sin(E). - *

- * This method is used when E is close to 0 and e close to 1, - * i.e. near the perigee of almost parabolic orbits - *

- * @param E eccentric anomaly - * @param e eccentricity - * @return E - e sin(E) - */ - private static double eMeSinE(final double E, final double e) { - double x = (1 - e) * FastMath.sin(E); - final double mE2 = -E * E; - double term = E; - double d = 0; - // the inequality test below IS intentional and should NOT be replaced by a check with a small tolerance - for (double x0 = Double.NaN; !Double.valueOf(x).equals(Double.valueOf(x0));) { - d += 2; - term *= mE2 / (d * (d + 1)); - x0 = x; - x = x - term; - } - return x; - } - - /** Computes the hyperbolic eccentric anomaly from the mean anomaly. - *

- * The algorithm used here for solving hyperbolic Kepler equation is - * Danby's iterative method (3rd order) with Vallado's initial guess. - *

- * @param M mean anomaly (rad) - * @param ecc eccentricity - * @return H the hyperbolic eccentric anomaly - */ - public static 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); + return KeplerianAnomalyUtility.ellipticEccentricToMean(e, M); } /** Computes the mean anomaly from the elliptic eccentric anomaly. * @param E eccentric anomaly (rad) * @param e eccentricity * @return M the mean anomaly - * @since 9.0 + * @deprecated As of 11.3, replaced by {@link KeplerianAnomalyUtility#ellipticEccentricToMean(double, double)}. */ public static double ellipticEccentricToMean(final double E, final double e) { - return E - e * FastMath.sin(E); + return KeplerianAnomalyUtility.ellipticEccentricToMean(e, E); } /** Computes the mean anomaly from the hyperbolic eccentric anomaly. * @param H hyperbolic eccentric anomaly (rad) * @param e eccentricity * @return M the mean anomaly - * @since 9.0 + * @deprecated As of 11.3, replaced by {@link KeplerianAnomalyUtility#hyperbolicEccentricToMean(double, double)}. */ public static double hyperbolicEccentricToMean(final double H, final double e) { - return e * FastMath.sinh(H) - H; + return KeplerianAnomalyUtility.hyperbolicEccentricToMean(e, H); } /** {@inheritDoc} */ diff --git a/src/test/java/org/orekit/orbits/FieldKeplerianAnomalyUtilityTest.java b/src/test/java/org/orekit/orbits/FieldKeplerianAnomalyUtilityTest.java new file mode 100644 index 000000000..cd33eae87 --- /dev/null +++ b/src/test/java/org/orekit/orbits/FieldKeplerianAnomalyUtilityTest.java @@ -0,0 +1,232 @@ +/* Copyright 2002-2022 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 org.hipparchus.CalculusFieldElement; +import org.hipparchus.Field; +import org.hipparchus.util.Decimal64Field; +import org.junit.jupiter.api.Assertions; +import org.junit.jupiter.api.Test; + +public class FieldKeplerianAnomalyUtilityTest { + + @Test + public void testEllipticMeanToTrue() { + doTestEllipticMeanToTrue(Decimal64Field.getInstance()); + } + + @Test + public void testEllipticTrueToMean() { + doTestEllipticTrueToMean(Decimal64Field.getInstance()); + } + + @Test + public void testEllipticEccentricToTrue() { + doTestEllipticEccentricToTrue(Decimal64Field.getInstance()); + } + + @Test + public void testEllipticTrueToEccentric() { + doTestEllipticTrueToEccentric(Decimal64Field.getInstance()); + } + + @Test + public void testEllipticMeanToEccentric() { + doTestEllipticMeanToEccentric(Decimal64Field.getInstance()); + } + + @Test + public void testEllipticEccentricToMean() { + doTestEllipticEccentricToMean(Decimal64Field.getInstance()); + } + + @Test + public void testHyperbolicMeanToTrue() { + doTestHyperbolicMeanToTrue(Decimal64Field.getInstance()); + } + + @Test + public void testHyperbolicTrueToMean() { + doTestHyperbolicTrueToMean(Decimal64Field.getInstance()); + } + + @Test + public void testHyperbolicEccentricToTrue() { + doTestHyperbolicEccentricToTrue(Decimal64Field.getInstance()); + } + + @Test + public void testHyperbolicTrueToEccentric() { + doTestHyperbolicTrueToEccentric(Decimal64Field.getInstance()); + } + + @Test + public void testHyperbolicMeanToEccentric() { + doTestHyperbolicMeanToEccentric(Decimal64Field.getInstance()); + } + + @Test + public void testHyperbolicEccentricToMean() { + doTestHyperbolicEccentricToMean(Decimal64Field.getInstance()); + } + + @Test + public void testIssue544() { + doTestIssue544(Decimal64Field.getInstance()); + } + + private > void doTestEllipticMeanToTrue(final Field field) { + final T e = field.getZero().add(0.231); + final T M = field.getZero().add(2.045); + final T v = FieldKeplerianAnomalyUtility.ellipticMeanToTrue(e, M); + Assertions.assertEquals(2.4004986679372027, v.getReal(), 1e-14); + } + + private > void doTestEllipticTrueToMean(final Field field) { + final T e = field.getZero().add(0.487); + final T v = field.getZero().add(1.386); + final T M = FieldKeplerianAnomalyUtility.ellipticTrueToMean(e, v); + Assertions.assertEquals(0.5238159114936672, M.getReal(), 1e-14); + } + + private > void doTestEllipticEccentricToTrue(final Field field) { + final T e = field.getZero().add(0.687); + final T E = field.getZero().add(4.639); + final T v = FieldKeplerianAnomalyUtility.ellipticEccentricToTrue(e, E); + Assertions.assertEquals(3.903008140176819, v.getReal(), 1e-14); + } + + private > void doTestEllipticTrueToEccentric(final Field field) { + final T e = field.getZero().add(0.527); + final T v = field.getZero().add(0.768); + final T E = FieldKeplerianAnomalyUtility.ellipticTrueToEccentric(e, v); + Assertions.assertEquals(0.44240462411915754, E.getReal(), 1e-14); + } + + private > void doTestEllipticMeanToEccentric(final Field field) { + final T e1 = field.getZero().add(0.726); + final T M1 = field.getZero().add(0.); + final T E1 = FieldKeplerianAnomalyUtility.ellipticMeanToEccentric(e1, M1); + Assertions.assertEquals(0.0, E1.getReal(), 1e-14); + + final T e2 = field.getZero().add(0.065); + final T M2 = field.getZero().add(4.586); + final T E2 = FieldKeplerianAnomalyUtility.ellipticMeanToEccentric(e2, M2); + Assertions.assertEquals(4.522172385101093, E2.getReal(), 1e-14); + + final T e3 = field.getZero().add(0.403); + final T M3 = field.getZero().add(0.121); + final T E3 = FieldKeplerianAnomalyUtility.ellipticMeanToEccentric(e3, M3); + Assertions.assertEquals(0.20175794699115656, E3.getReal(), 1e-14); + + final T e4 = field.getZero().add(0.999); + final T M4 = field.getZero().add(0.028); + final T E4 = FieldKeplerianAnomalyUtility.ellipticMeanToEccentric(e4, M4); + Assertions.assertEquals(0.5511071508829587, E4.getReal(), 1e-14); + } + + private > void doTestEllipticEccentricToMean(final Field field) { + final T e = field.getZero().add(0.192); + final T E = field.getZero().add(2.052); + final T M = FieldKeplerianAnomalyUtility.ellipticEccentricToMean(e, E); + Assertions.assertEquals(1.881803817764882, M.getReal(), 1e-14); + } + + private > void doTestHyperbolicMeanToTrue(final Field field) { + final T e = field.getZero().add(1.027); + final T M = field.getZero().add(1.293); + final T v = FieldKeplerianAnomalyUtility.hyperbolicMeanToTrue(e, M); + Assertions.assertEquals(2.8254185280004855, v.getReal(), 1e-14); + } + + private > void doTestHyperbolicTrueToMean(final Field field) { + final T e = field.getZero().add(1.161); + final T v = field.getZero().add(-2.469); + final T M = FieldKeplerianAnomalyUtility.hyperbolicTrueToMean(e, v); + Assertions.assertEquals(-2.5499244818919915, M.getReal(), 1e-14); + } + + private > void doTestHyperbolicEccentricToTrue(final Field field) { + final T e = field.getZero().add(2.161); + final T E = field.getZero().add(-1.204); + final T v = FieldKeplerianAnomalyUtility.hyperbolicEccentricToTrue(e, E); + Assertions.assertEquals(-1.4528528149658333, v.getReal(), 1e-14); + } + + private > void doTestHyperbolicTrueToEccentric(final Field field) { + final T e = field.getZero().add(1.595); + final T v = field.getZero().add(0.298); + final T E = FieldKeplerianAnomalyUtility.hyperbolicTrueToEccentric(e, v); + Assertions.assertEquals(0.1440079208139455, E.getReal(), 1e-14); + } + + private > void doTestHyperbolicMeanToEccentric(final Field field) { + final T e1 = field.getZero().add(1.201); + final T M1 = field.getZero(); + final T E1 = FieldKeplerianAnomalyUtility.hyperbolicMeanToEccentric(e1, M1); + Assertions.assertEquals(0.0, E1.getReal(), 1e-14); + + final T e2 = field.getZero().add(1.127); + final T M2 = field.getZero().add(-3.624); + final T E2 = FieldKeplerianAnomalyUtility.hyperbolicMeanToEccentric(e2, M2); + Assertions.assertEquals(-2.3736718687722265, E2.getReal(), 1e-14); + + final T e3 = field.getZero().add(1.338); + final T M3 = field.getZero().add(-0.290); + final T E3 = FieldKeplerianAnomalyUtility.hyperbolicMeanToEccentric(e3, M3); + Assertions.assertEquals(-0.6621795141831807, E3.getReal(), 1e-14); + + final T e4 = field.getZero().add(1.044); + final T M4 = field.getZero().add(3.996); + final T E4 = FieldKeplerianAnomalyUtility.hyperbolicMeanToEccentric(e4, M4); + Assertions.assertEquals(2.532614977388778, E4.getReal(), 1e-14); + + final T e5 = field.getZero().add(2.052); + final T M5 = field.getZero().add(4.329); + final T E5 = FieldKeplerianAnomalyUtility.hyperbolicMeanToEccentric(e5, M5); + Assertions.assertEquals(1.816886788278918, E5.getReal(), 1e-14); + + final T e6 = field.getZero().add(2.963); + final T M6 = field.getZero().add(-1.642); + final T E6 = FieldKeplerianAnomalyUtility.hyperbolicMeanToEccentric(e6, M6); + Assertions.assertEquals(-0.7341946491456494, E6.getReal(), 1e-14); + + final T e7 = field.getZero().add(4.117); + final T M7 = field.getZero().add(-0.286); + final T E7 = FieldKeplerianAnomalyUtility.hyperbolicMeanToEccentric(e7, M7); + Assertions.assertEquals(-0.09158570899196887, E7.getReal(), 1e-14); + } + + private > void doTestHyperbolicEccentricToMean(final Field field) { + final T e = field.getZero().add(1.801); + final T E = field.getZero().add(3.287); + final T M = FieldKeplerianAnomalyUtility.hyperbolicEccentricToMean(e, E); + Assertions.assertEquals(20.77894350750361, M.getReal(), 1e-14); + } + + private > void doTestIssue544(final Field field) { + // Initial parameters + // In order to test the issue, we voluntarily set the anomaly at T.NaN. + T e = field.getZero().add(0.7311); + T anomaly = field.getZero().add(Double.NaN); + // Computes the elliptic eccentric anomaly + T E = FieldKeplerianAnomalyUtility.ellipticMeanToEccentric(e, anomaly); + // Verify that an infinite loop did not occur + Assertions.assertTrue(Double.isNaN(E.getReal())); + } + +} diff --git a/src/test/java/org/orekit/orbits/FieldKeplerianOrbitTest.java b/src/test/java/org/orekit/orbits/FieldKeplerianOrbitTest.java index 857b203da..ca79f1b4e 100644 --- a/src/test/java/org/orekit/orbits/FieldKeplerianOrbitTest.java +++ b/src/test/java/org/orekit/orbits/FieldKeplerianOrbitTest.java @@ -256,11 +256,6 @@ public class FieldKeplerianOrbitTest { doTestCopyNonKeplerianAcceleration(Decimal64Field.getInstance()); } - @Test - public void testIssue544() { - doTestIssue544(Decimal64Field.getInstance()); - } - @Test public void testIssue674() { doTestIssue674(Decimal64Field.getInstance()); @@ -1898,17 +1893,6 @@ public class FieldKeplerianOrbitTest { } - private > void doTestIssue544(Field field) { - // Initial parameters - // In order to test the issue, we volontary set the anomaly at Double.NaN. - T e= field.getZero().add(0.7311); - T anomaly= field.getZero().add(Double.NaN); - // Computes the elliptic eccentric anomaly - T E = FieldKeplerianOrbit.meanToEllipticEccentric(anomaly, e); - // Verify that an infinite loop did not occur - Assertions.assertTrue(Double.isNaN(E.getReal())); - } - private > void doTestIssue674(Field field) { try { final T zero = field.getZero(); diff --git a/src/test/java/org/orekit/orbits/KeplerianAnomalyUtilityTest.java b/src/test/java/org/orekit/orbits/KeplerianAnomalyUtilityTest.java new file mode 100644 index 000000000..e925dec77 --- /dev/null +++ b/src/test/java/org/orekit/orbits/KeplerianAnomalyUtilityTest.java @@ -0,0 +1,177 @@ +/* Copyright 2002-2022 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 org.junit.jupiter.api.Assertions; +import org.junit.jupiter.api.Test; + +public class KeplerianAnomalyUtilityTest { + + @Test + public void testEllipticMeanToTrue() { + final double e = 0.231; + final double M = 2.045; + final double v = KeplerianAnomalyUtility.ellipticMeanToTrue(e, M); + Assertions.assertEquals(2.4004986679372027, v, 1e-14); + } + + @Test + public void testEllipticTrueToMean() { + final double e = 0.487; + final double v = 1.386; + final double M = KeplerianAnomalyUtility.ellipticTrueToMean(e, v); + Assertions.assertEquals(0.5238159114936672, M, 1e-14); + } + + @Test + public void testEllipticEccentricToTrue() { + final double e = 0.687; + final double E = 4.639; + final double v = KeplerianAnomalyUtility.ellipticEccentricToTrue(e, E); + Assertions.assertEquals(3.903008140176819, v, 1e-14); + } + + @Test + public void testEllipticTrueToEccentric() { + final double e = 0.527; + final double v = 0.768; + final double E = KeplerianAnomalyUtility.ellipticTrueToEccentric(e, v); + Assertions.assertEquals(0.44240462411915754, E, 1e-14); + } + + @Test + public void testEllipticMeanToEccentric() { + final double e1 = 0.726; + final double M1 = 0.; + final double E1 = KeplerianAnomalyUtility.ellipticMeanToEccentric(e1, M1); + Assertions.assertEquals(0.0, E1, 1e-14); + + final double e2 = 0.065; + final double M2 = 4.586; + final double E2 = KeplerianAnomalyUtility.ellipticMeanToEccentric(e2, M2); + Assertions.assertEquals(4.522172385101093, E2, 1e-14); + + final double e3 = 0.403; + final double M3 = 0.121; + final double E3 = KeplerianAnomalyUtility.ellipticMeanToEccentric(e3, M3); + Assertions.assertEquals(0.20175794699115656, E3, 1e-14); + + final double e4 = 0.999; + final double M4 = 0.028; + final double E4 = KeplerianAnomalyUtility.ellipticMeanToEccentric(e4, M4); + Assertions.assertEquals(0.5511071508829587, E4, 1e-14); + } + + @Test + public void testEllipticEccentricToMean() { + final double e = 0.192; + final double E = 2.052; + final double M = KeplerianAnomalyUtility.ellipticEccentricToMean(e, E); + Assertions.assertEquals(1.881803817764882, M, 1e-14); + } + + @Test + public void testHyperbolicMeanToTrue() { + final double e = 1.027; + final double M = 1.293; + final double v = KeplerianAnomalyUtility.hyperbolicMeanToTrue(e, M); + Assertions.assertEquals(2.8254185280004855, v, 1e-14); + } + + @Test + public void testHyperbolicTrueToMean() { + final double e = 1.161; + final double v = -2.469; + final double M = KeplerianAnomalyUtility.hyperbolicTrueToMean(e, v); + Assertions.assertEquals(-2.5499244818919915, M, 1e-14); + } + + @Test + public void testHyperbolicEccentricToTrue() { + final double e = 2.161; + final double E = -1.204; + final double v = KeplerianAnomalyUtility.hyperbolicEccentricToTrue(e, E); + Assertions.assertEquals(-1.4528528149658333, v, 1e-14); + } + + @Test + public void testHyperbolicTrueToEccentric() { + final double e = 1.595; + final double v = 0.298; + final double E = KeplerianAnomalyUtility.hyperbolicTrueToEccentric(e, v); + Assertions.assertEquals(0.1440079208139455, E, 1e-14); + } + + @Test + public void testHyperbolicMeanToEccentric() { + final double e1 = 1.201; + final double M1 = 0.0; + final double E1 = KeplerianAnomalyUtility.hyperbolicMeanToEccentric(e1, M1); + Assertions.assertEquals(0.0, E1, 1e-14); + + final double e2 = 1.127; + final double M2 = -3.624; + final double E2 = KeplerianAnomalyUtility.hyperbolicMeanToEccentric(e2, M2); + Assertions.assertEquals(-2.3736718687722265, E2, 1e-14); + + final double e3 = 1.338; + final double M3 = -0.290; + final double E3 = KeplerianAnomalyUtility.hyperbolicMeanToEccentric(e3, M3); + Assertions.assertEquals(-0.6621795141831807, E3, 1e-14); + + final double e4 = 1.044; + final double M4 = 3.996; + final double E4 = KeplerianAnomalyUtility.hyperbolicMeanToEccentric(e4, M4); + Assertions.assertEquals(2.532614977388778, E4, 1e-14); + + final double e5 = 2.052; + final double M5 = 4.329; + final double E5 = KeplerianAnomalyUtility.hyperbolicMeanToEccentric(e5, M5); + Assertions.assertEquals(1.816886788278918, E5, 1e-14); + + final double e6 = 2.963; + final double M6 = -1.642; + final double E6 = KeplerianAnomalyUtility.hyperbolicMeanToEccentric(e6, M6); + Assertions.assertEquals(-0.7341946491456494, E6, 1e-14); + + final double e7 = 4.117; + final double M7 = -0.286; + final double E7 = KeplerianAnomalyUtility.hyperbolicMeanToEccentric(e7, M7); + Assertions.assertEquals(-0.09158570899196887, E7, 1e-14); + } + + @Test + public void testHyperbolicEccentricToMean() { + final double e = 1.801; + final double E = 3.287; + final double M = KeplerianAnomalyUtility.hyperbolicEccentricToMean(e, E); + Assertions.assertEquals(20.77894350750361, M, 1e-14); + } + + @Test + public void testIssue544() { + // Initial parameters + // In order to test the issue, we voluntarily set the anomaly at Double.NaN. + double e = 0.7311; + double anomaly = Double.NaN; + // Computes the elliptic eccentric anomaly + double E = KeplerianAnomalyUtility.ellipticMeanToEccentric(e, anomaly); + // Verify that an infinite loop did not occur + Assertions.assertTrue(Double.isNaN(E)); + } + +} diff --git a/src/test/java/org/orekit/orbits/KeplerianOrbitTest.java b/src/test/java/org/orekit/orbits/KeplerianOrbitTest.java index 2f88a9ae4..68d5457be 100644 --- a/src/test/java/org/orekit/orbits/KeplerianOrbitTest.java +++ b/src/test/java/org/orekit/orbits/KeplerianOrbitTest.java @@ -1587,18 +1587,6 @@ public class KeplerianOrbitTest { } - @Test - public void testIssue544() { - // Initial parameters - // In order to test the issue, we volontary set the anomaly at Double.NaN. - double e = 0.7311; - double anomaly = Double.NaN; - // Computes the elliptic eccentric anomaly - double E = KeplerianOrbit.meanToEllipticEccentric(anomaly, e); - // Verify that an infinite loop did not occur - Assertions.assertTrue(Double.isNaN(E)); - } - @Test public void testIssue674() { try { -- GitLab From 2faf13ff0cb377e2bbbfde41870f4141d1918611 Mon Sep 17 00:00:00 2001 From: Andrew Goetz Date: Tue, 4 Oct 2022 10:38:26 -0400 Subject: [PATCH 2/2] Incorporated Gooding and Odell algorithm for solving the hyperbolic Kepler equation. --- src/changes/changes.xml | 6 + .../orbits/FieldKeplerianAnomalyUtility.java | 119 +++++++++++------- .../orbits/KeplerianAnomalyUtility.java | 103 +++++++++------ .../FieldKeplerianAnomalyUtilityTest.java | 6 + .../orbits/FieldKeplerianOrbitTest.java | 2 +- .../orbits/KeplerianAnomalyUtilityTest.java | 6 + .../org/orekit/orbits/KeplerianOrbitTest.java | 2 +- 7 files changed, 162 insertions(+), 82 deletions(-) diff --git a/src/changes/changes.xml b/src/changes/changes.xml index e6566df7f..3691ab9a0 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -21,6 +21,12 @@ + + Moved Keplerian anomaly conversion methods to KeplerianAnomalyUtility + and FieldKeplerianAnomalyUtility, deprecating the methods in + KeplerianOrbit and FieldKeplerianOrbit. Incorporated Gooding and Odell + algorithm for solving the hyperbolic Kepler equation. + Added Unscented Semi-analytical Kalman Estimator. diff --git a/src/main/java/org/orekit/orbits/FieldKeplerianAnomalyUtility.java b/src/main/java/org/orekit/orbits/FieldKeplerianAnomalyUtility.java index a4afef804..f189d270a 100644 --- a/src/main/java/org/orekit/orbits/FieldKeplerianAnomalyUtility.java +++ b/src/main/java/org/orekit/orbits/FieldKeplerianAnomalyUtility.java @@ -17,7 +17,8 @@ package org.orekit.orbits; import org.hipparchus.CalculusFieldElement; -import org.hipparchus.exception.MathIllegalArgumentException; +import org.hipparchus.Field; +import org.hipparchus.exception.MathIllegalStateException; import org.hipparchus.util.FastMath; import org.hipparchus.util.FieldSinCos; import org.hipparchus.util.MathUtils; @@ -273,8 +274,10 @@ public class FieldKeplerianAnomalyUtility { /** * Computes the hyperbolic eccentric anomaly from the hyperbolic mean anomaly. *

- * The algorithm used here for solving hyperbolic Kepler equation is Danby's - * iterative method (3rd order) with Vallado's initial guess. + * The algorithm used here for solving hyperbolic Kepler equation is from + * Gooding, R.H., Odell, A.W. "The hyperbolic Kepler equation (and the elliptic + * equation revisited)." Celestial Mechanics 44, 267–282 (1988). + * https://doi.org/10.1007/BF01235540 *

* @param field type * @param e eccentricity > 1 @@ -282,50 +285,82 @@ public class FieldKeplerianAnomalyUtility { * @return hyperbolic eccentric anomaly */ public static > T hyperbolicMeanToEccentric(final T e, final T M) { - // Resolution of hyperbolic Kepler equation for Keplerian parameters - - // Field value of pi - final T pi = e.getPi(); - - // Initial guess - T H; - if (e.getReal() < 1.6) { - if (-pi.getReal() < M.getReal() && M.getReal() < 0. || M.getReal() > pi.getReal()) { - H = M.subtract(e); + final Field field = e.getField(); + final T zero = field.getZero(); + final T one = field.getOne(); + final T two = zero.add(2.0); + final T three = zero.add(3.0); + final T half = zero.add(0.5); + final T onePointFive = zero.add(1.5); + final T fourThirds = zero.add(4.0).divide(zero.add(3.0)); + + // Solve L = S - g * asinh(S) for S = sinh(H). + final T L = M.divide(e); + final T g = e.reciprocal(); + final T g1 = one.subtract(g); + + // Starter based on Lagrange's theorem. + T S = L; + if (L.isZero()) { + return M.getField().getZero(); + } + final T cl = L.multiply(L).add(one).sqrt(); + final T al = L.asinh(); + final T w = g.multiply(g).multiply(al).divide(cl.multiply(cl).multiply(cl)); + S = one.subtract(g.divide(cl)); + S = L.add(g.multiply(al).divide(S.multiply(S).multiply(S) + .add(w.multiply(L).multiply(onePointFive.subtract(fourThirds.multiply(g)))).cbrt())); + + // Two iterations (at most) of Halley-then-Newton process. + for (int i = 0; i < 2; ++i) { + final T s0 = S.multiply(S); + final T s1 = s0.add(one); + final T s2 = s1.sqrt(); + final T s3 = s1.multiply(s2); + final T fdd = g.multiply(S).divide(s3); + final T fddd = g.multiply(one.subtract(two.multiply(s0))).divide(s1.multiply(s3)); + + T f; + T fd; + if (s0.divide(zero.add(6.0)).add(g1).getReal() >= 0.5) { + f = S.subtract(g.multiply(S.asinh())).subtract(L); + fd = one.subtract(g.divide(s2)); } else { - H = M.add(e); + // Accurate computation of S - (1 - g1) * asinh(S) + // when (g1, S) is close to (0, 0). + final T t = S.divide(one.add(one.add(S.multiply(S)).sqrt())); + final T tsq = t.multiply(t); + T x = S.multiply(g1.add(g.multiply(tsq))); + T term = two.multiply(g).multiply(t); + T twoI1 = one; + T x0; + int j = 0; + do { + if (++j == 1000000) { + // This isn't expected to happen, but it protects against an infinite loop. + throw new MathIllegalStateException( + OrekitMessages.UNABLE_TO_COMPUTE_HYPERBOLIC_ECCENTRIC_ANOMALY, j); + } + twoI1 = twoI1.add(2.0); + term = term.multiply(tsq); + x0 = x; + x = x.subtract(term.divide(twoI1)); + } while (x.getReal() != x0.getReal()); + f = x.subtract(L); + fd = s0.divide(s2.add(one)).add(g1).divide(s2); } - } else { - if (e.getReal() < 3.6 && M.abs().getReal() > pi.getReal()) { - H = M.subtract(e.copySign(M)); - } else { - H = M.divide(e.subtract(1)); + final T ds = f.multiply(fd).divide(half.multiply(f).multiply(fdd).subtract(fd.multiply(fd))); + final T stemp = S.add(ds); + if (S.getReal() == stemp.getReal()) { + break; } + f = f.add(ds.multiply(fd.add(half.multiply(ds.multiply(fdd.add(ds.divide(three).multiply(fddd))))))); + fd = fd.add(ds.multiply(fdd.add(half.multiply(ds).multiply(fddd)))); + S = stemp.subtract(f.divide(fd)); } - // Iterative computation - int iter = 0; - do { - final T f3 = e.multiply(H.cosh()); - final T f2 = e.multiply(H.sinh()); - final T f1 = f3.subtract(1); - final T f0 = f2.subtract(H).subtract(M); - final T f12 = f1.multiply(2); - final T d = f0.divide(f12); - final T fdf = f1.subtract(d.multiply(f2)); - final T ds = f0.divide(fdf); - - final T shift = f0.divide(fdf.add(ds.multiply(ds.multiply(f3.divide(6))))); - - H = H.subtract(shift); - - if ((shift.abs().getReal()) <= 1.0e-12) { - return H; - } - - } while (++iter < 50); - - throw new MathIllegalArgumentException(OrekitMessages.UNABLE_TO_COMPUTE_HYPERBOLIC_ECCENTRIC_ANOMALY, iter); + final T H = S.asinh(); + return H; } /** diff --git a/src/main/java/org/orekit/orbits/KeplerianAnomalyUtility.java b/src/main/java/org/orekit/orbits/KeplerianAnomalyUtility.java index 5120cfdc7..16ccec0bb 100644 --- a/src/main/java/org/orekit/orbits/KeplerianAnomalyUtility.java +++ b/src/main/java/org/orekit/orbits/KeplerianAnomalyUtility.java @@ -242,55 +242,82 @@ public final class KeplerianAnomalyUtility { /** * Computes the hyperbolic eccentric anomaly from the hyperbolic mean anomaly. *

- * The algorithm used here for solving hyperbolic Kepler equation is Danby's - * iterative method (3rd order) with Vallado's initial guess. + * The algorithm used here for solving hyperbolic Kepler equation is from + * Gooding, R.H., Odell, A.W. "The hyperbolic Kepler equation (and the elliptic + * equation revisited)." Celestial Mechanics 44, 267–282 (1988). + * https://doi.org/10.1007/BF01235540 *

* @param e eccentricity > 1 * @param M hyperbolic mean anomaly * @return hyperbolic eccentric anomaly */ public static double hyperbolicMeanToEccentric(final double e, final double M) { - // Resolution of hyperbolic Kepler equation for Keplerian parameters - - // Initial guess - double H; - if (e < 1.6) { - if (-FastMath.PI < M && M < 0. || M > FastMath.PI) { - H = M - e; + // Solve L = S - g * asinh(S) for S = sinh(H). + final double L = M / e; + final double g = 1.0 / e; + final double g1 = 1.0 - g; + + // Starter based on Lagrange's theorem. + double S = L; + if (L == 0.0) { + return 0.0; + } + final double cl = FastMath.sqrt(1.0 + L * L); + final double al = FastMath.asinh(L); + final double w = g * g * al / (cl * cl * cl); + S = 1.0 - g / cl; + S = L + g * al / FastMath.cbrt(S * S * S + w * L * (1.5 - (4.0 / 3.0) * g)); + + // Two iterations (at most) of Halley-then-Newton process. + for (int i = 0; i < 2; ++i) { + final double s0 = S * S; + final double s1 = s0 + 1.0; + final double s2 = FastMath.sqrt(s1); + final double s3 = s1 * s2; + final double fdd = g * S / s3; + final double fddd = g * (1.0 - 2.0 * s0) / (s1 * s3); + + double f; + double fd; + if (s0 / 6.0 + g1 >= 0.5) { + f = (S - g * FastMath.asinh(S)) - L; + fd = 1.0 - g / s2; } else { - H = M + e; + // Accurate computation of S - (1 - g1) * asinh(S) + // when (g1, S) is close to (0, 0). + final double t = S / (1.0 + FastMath.sqrt(1.0 + S * S)); + final double tsq = t * t; + double x = S * (g1 + g * tsq); + double term = 2.0 * g * t; + double twoI1 = 1.0; + double x0; + int j = 0; + do { + if (++j == 1000000) { + // This isn't expected to happen, but it protects against an infinite loop. + throw new MathIllegalStateException( + OrekitMessages.UNABLE_TO_COMPUTE_HYPERBOLIC_ECCENTRIC_ANOMALY, j); + } + twoI1 += 2.0; + term *= tsq; + x0 = x; + x -= term / twoI1; + } while (x != x0); + f = x - L; + fd = (s0 / (s2 + 1.0) + g1) / s2; } - } else { - if (e < 3.6 && FastMath.abs(M) > FastMath.PI) { - H = M - FastMath.copySign(e, M); - } else { - H = M / (e - 1.); + final double ds = f * fd / (0.5 * f * fdd - fd * fd); + final double stemp = S + ds; + if (S == stemp) { + break; } + f += ds * (fd + 0.5 * ds * (fdd + ds / 3.0 * fddd)); + fd += ds * (fdd + 0.5 * ds * fddd); + S = stemp - f / fd; } - // Iterative computation - int iter = 0; - do { - final double f3 = e * FastMath.cosh(H); - final double f2 = e * 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); + final double H = FastMath.asinh(S); + return H; } /** diff --git a/src/test/java/org/orekit/orbits/FieldKeplerianAnomalyUtilityTest.java b/src/test/java/org/orekit/orbits/FieldKeplerianAnomalyUtilityTest.java index cd33eae87..5a4b517d7 100644 --- a/src/test/java/org/orekit/orbits/FieldKeplerianAnomalyUtilityTest.java +++ b/src/test/java/org/orekit/orbits/FieldKeplerianAnomalyUtilityTest.java @@ -209,6 +209,12 @@ public class FieldKeplerianAnomalyUtilityTest { final T M7 = field.getZero().add(-0.286); final T E7 = FieldKeplerianAnomalyUtility.hyperbolicMeanToEccentric(e7, M7); Assertions.assertEquals(-0.09158570899196887, E7.getReal(), 1e-14); + + // Issue 951. + final T e8 = field.getZero().add(1.251844925917281); + final T M8 = field.getZero().add(54.70111712786907); + final T E8 = FieldKeplerianAnomalyUtility.hyperbolicMeanToEccentric(e8, M8); + Assertions.assertEquals(4.550432282228856, E8.getReal(), 1e-14); } private > void doTestHyperbolicEccentricToMean(final Field field) { diff --git a/src/test/java/org/orekit/orbits/FieldKeplerianOrbitTest.java b/src/test/java/org/orekit/orbits/FieldKeplerianOrbitTest.java index ca79f1b4e..ead639700 100644 --- a/src/test/java/org/orekit/orbits/FieldKeplerianOrbitTest.java +++ b/src/test/java/org/orekit/orbits/FieldKeplerianOrbitTest.java @@ -1739,7 +1739,7 @@ public class FieldKeplerianOrbitTest { MatcherAssert.assertThat(rebuilt.getPerigeeArgumentDot().getReal(), relativelyCloseTo(orbit.getPerigeeArgumentDot().getReal(), 1)); MatcherAssert.assertThat(rebuilt.getRightAscensionOfAscendingNodeDot().getReal(), relativelyCloseTo(orbit.getRightAscensionOfAscendingNodeDot().getReal(), 1)); for (PositionAngle type2 : PositionAngle.values()) { - MatcherAssert.assertThat(rebuilt.getAnomaly(type2).getReal(), relativelyCloseTo(orbit.getAnomaly(type2).getReal(), 2)); + MatcherAssert.assertThat(rebuilt.getAnomaly(type2).getReal(), relativelyCloseTo(orbit.getAnomaly(type2).getReal(), 3)); MatcherAssert.assertThat(rebuilt.getAnomalyDot(type2).getReal(), relativelyCloseTo(orbit.getAnomalyDot(type2).getReal(), 4)); } } diff --git a/src/test/java/org/orekit/orbits/KeplerianAnomalyUtilityTest.java b/src/test/java/org/orekit/orbits/KeplerianAnomalyUtilityTest.java index e925dec77..387b5524f 100644 --- a/src/test/java/org/orekit/orbits/KeplerianAnomalyUtilityTest.java +++ b/src/test/java/org/orekit/orbits/KeplerianAnomalyUtilityTest.java @@ -152,6 +152,12 @@ public class KeplerianAnomalyUtilityTest { final double M7 = -0.286; final double E7 = KeplerianAnomalyUtility.hyperbolicMeanToEccentric(e7, M7); Assertions.assertEquals(-0.09158570899196887, E7, 1e-14); + + // Issue 951. + final double e8 = 1.251844925917281; + final double M8 = 54.70111712786907; + final double E8 = KeplerianAnomalyUtility.hyperbolicMeanToEccentric(e8, M8); + Assertions.assertEquals(4.550432282228856, E8, 1e-14); } @Test diff --git a/src/test/java/org/orekit/orbits/KeplerianOrbitTest.java b/src/test/java/org/orekit/orbits/KeplerianOrbitTest.java index 68d5457be..c666faff4 100644 --- a/src/test/java/org/orekit/orbits/KeplerianOrbitTest.java +++ b/src/test/java/org/orekit/orbits/KeplerianOrbitTest.java @@ -1380,7 +1380,7 @@ public class KeplerianOrbitTest { MatcherAssert.assertThat(rebuilt.getPerigeeArgumentDot(), relativelyCloseTo(orbit.getPerigeeArgumentDot(), 1)); MatcherAssert.assertThat(rebuilt.getRightAscensionOfAscendingNodeDot(), relativelyCloseTo(orbit.getRightAscensionOfAscendingNodeDot(), 1)); for (PositionAngle type2 : PositionAngle.values()) { - MatcherAssert.assertThat(rebuilt.getAnomaly(type2), relativelyCloseTo(orbit.getAnomaly(type2), 2)); + MatcherAssert.assertThat(rebuilt.getAnomaly(type2), relativelyCloseTo(orbit.getAnomaly(type2), 3)); MatcherAssert.assertThat(rebuilt.getAnomalyDot(type2), relativelyCloseTo(orbit.getAnomalyDot(type2), 4)); } } -- GitLab