diff --git a/src/changes/changes.xml b/src/changes/changes.xml index e6566df7f787ed8481f833664f394d89c52dd786..3691ab9a026ed39872ccb9fd5f8b91018dd63636 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -21,6 +21,12 @@
- * 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 dc6fd76c4e96437af8aeb5035782d434a4a9ad90..bb86ece8dd0b297c07f52346f92769d1cf6025e4 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- * 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 0000000000000000000000000000000000000000..f189d270a3c587788dece543d0d5344758b6a3b3 --- /dev/null +++ b/src/main/java/org/orekit/orbits/FieldKeplerianAnomalyUtility.java @@ -0,0 +1,377 @@ +/* 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.exception.MathIllegalStateException; +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+ * 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+ * 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+ * 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- * 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- * 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- * 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+ * 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 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) { + // 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 { + // 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; + } + 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; + } + + final double H = FastMath.asinh(S); + return H; + } + + /** + * 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 5ad29fa60d83f36b1aa1f27349671453c4407ddc..00f585de20923f2d1d238f3560395ce2f45740c0 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 0000000000000000000000000000000000000000..5a4b517d7804a855844ea4914a8e62e1bd593adf --- /dev/null +++ b/src/test/java/org/orekit/orbits/FieldKeplerianAnomalyUtilityTest.java @@ -0,0 +1,238 @@ +/* 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