1 /* Copyright 2002-2026 CS GROUP
2 * Licensed to CS GROUP (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.orbits;
18
19 import org.hipparchus.exception.MathIllegalStateException;
20 import org.hipparchus.util.FastMath;
21 import org.hipparchus.util.MathUtils;
22 import org.hipparchus.util.SinCos;
23 import org.orekit.errors.OrekitException;
24 import org.orekit.errors.OrekitMessages;
25
26 /**
27 * Utility methods for converting between different Keplerian anomalies.
28 * @author Luc Maisonobe
29 * @author Andrew Goetz
30 * @author Alberto Fossa'
31 */
32 public final class KeplerianAnomalyUtility {
33
34 /**
35 * Maximum value of the hyperbolic sine in the initialization of the
36 * iterative computation of the hyperbolic eccentric anomaly difference.
37 */
38 static final double MAXIMUM_INITIAL_HYPERBOLIC_SINE = 1e3;
39
40 /**
41 * Convergence threshold in the iterative computation of the
42 * elliptic and hyperbolic eccentric anomalies differences.
43 */
44 static final double CONVERGENCE_THRESHOLD = 1e-10;
45
46 /** Tolerance for stopping criterion in iterative conversion from mean to eccentric angle. */
47 static final double RELATIVE_CONVERGENCE_THRESHOLD = 1.0e-13;
48
49 /**
50 * Maximum number of iterations in the iterative computation of
51 * the elliptic and hyperbolic eccentric anomalies differences.
52 */
53 static final int MAXIMUM_ITERATIONS = 50;
54
55 /** First coefficient to compute elliptic Kepler equation solver starter. */
56 private static final double A;
57
58 /** Second coefficient to compute elliptic Kepler equation solver starter. */
59 private static final double B;
60
61 static {
62 final double k1 = 3 * FastMath.PI + 2;
63 final double k2 = FastMath.PI - 1;
64 final double k3 = 6 * FastMath.PI - 1;
65 A = 3 * k2 * k2 / k1;
66 B = k3 * k3 / (6 * k1);
67 }
68
69 /** Private constructor for utility class. */
70 private KeplerianAnomalyUtility() {
71 // Nothing to do
72 }
73
74 /**
75 * Computes the elliptic true anomaly from the elliptic mean anomaly.
76 * @param e eccentricity such that 0 ≤ e < 1
77 * @param M elliptic mean anomaly (rad)
78 * @return elliptic true anomaly (rad)
79 */
80 public static double ellipticMeanToTrue(final double e, final double M) {
81 final double E = ellipticMeanToEccentric(e, M);
82 final double v = ellipticEccentricToTrue(e, E);
83 return v;
84 }
85
86 /**
87 * Computes the elliptic mean anomaly from the elliptic true anomaly.
88 * @param e eccentricity such that 0 ≤ e < 1
89 * @param v elliptic true anomaly (rad)
90 * @return elliptic mean anomaly (rad)
91 */
92 public static double ellipticTrueToMean(final double e, final double v) {
93 final double E = ellipticTrueToEccentric(e, v);
94 final double M = ellipticEccentricToMean(e, E);
95 return M;
96 }
97
98 /**
99 * Computes the elliptic true anomaly from the elliptic eccentric anomaly.
100 * @param e eccentricity such that 0 ≤ e < 1
101 * @param E elliptic eccentric anomaly (rad)
102 * @return elliptic true anomaly (rad)
103 */
104 public static double ellipticEccentricToTrue(final double e, final double E) {
105 final double beta = e / (1 + FastMath.sqrt((1 - e) * (1 + e)));
106 final SinCos scE = FastMath.sinCos(E);
107 return E + 2 * FastMath.atan(beta * scE.sin() / (1 - beta * scE.cos()));
108 }
109
110 /**
111 * Computes the elliptic eccentric anomaly from the elliptic true anomaly.
112 * @param e eccentricity such that 0 ≤ e < 1
113 * @param v elliptic true anomaly (rad)
114 * @return elliptic eccentric anomaly (rad)
115 */
116 public static double ellipticTrueToEccentric(final double e, final double v) {
117 final double beta = e / (1 + FastMath.sqrt(1 - e * e));
118 final SinCos scv = FastMath.sinCos(v);
119 return v - 2 * FastMath.atan(beta * scv.sin() / (1 + beta * scv.cos()));
120 }
121
122 /**
123 * Computes the elliptic eccentric anomaly from the elliptic mean anomaly.
124 * <p>
125 * The algorithm used here for solving hyperbolic Kepler equation is from Odell,
126 * A.W., Gooding, R.H. "Procedures for solving Kepler's equation." Celestial
127 * Mechanics 38, 307–334 (1986). https://doi.org/10.1007/BF01238923
128 * </p>
129 * @param e eccentricity such that 0 ≤ e < 1
130 * @param M elliptic mean anomaly (rad)
131 * @return elliptic eccentric anomaly (rad)
132 */
133 public static double ellipticMeanToEccentric(final double e, final double M) {
134 // reduce M to [-PI PI) interval
135 final double reducedM = MathUtils.normalizeAngle(M, 0.0);
136
137 // compute start value according to A. W. Odell and R. H. Gooding S12 starter
138 double E;
139 if (FastMath.abs(reducedM) < 1.0 / 6.0) {
140 E = reducedM + e * (FastMath.cbrt(6 * reducedM) - reducedM);
141 } else {
142 if (reducedM < 0) {
143 final double w = FastMath.PI + reducedM;
144 E = reducedM + e * (A * w / (B - w) - FastMath.PI - reducedM);
145 } else {
146 final double w = FastMath.PI - reducedM;
147 E = reducedM + e * (FastMath.PI - A * w / (B - w) - reducedM);
148 }
149 }
150
151 final double e1 = 1 - e;
152 final boolean noCancellationRisk = (e1 + E * E / 6) >= 0.1;
153
154 // perform two iterations, each consisting of one Halley step and one
155 // Newton-Raphson step
156 for (int j = 0; j < 2; ++j) {
157 final double f;
158 double fd;
159 final SinCos sc = FastMath.sinCos(E);
160 final double fdd = e * sc.sin();
161 final double fddd = e * sc.cos();
162 if (noCancellationRisk) {
163 f = (E - fdd) - reducedM;
164 fd = 1 - fddd;
165 } else {
166 f = eMeSinE(e, E) - reducedM;
167 final double s = FastMath.sin(0.5 * E);
168 fd = e1 + 2 * e * s * s;
169 }
170 final double dee = f * fd / (0.5 * f * fdd - fd * fd);
171
172 // update eccentric anomaly, using expressions that limit underflow problems
173 final double w = fd + 0.5 * dee * (fdd + dee * fddd / 3);
174 fd += dee * (fdd + 0.5 * dee * fddd);
175 E -= (f - dee * (fd - w)) / fd;
176
177 }
178
179 // expand the result back to original range
180 E += M - reducedM;
181
182 return E;
183 }
184
185 /**
186 * Accurate computation of E - e sin(E).
187 * <p>
188 * This method is used when E is close to 0 and e close to 1, i.e. near the
189 * perigee of almost parabolic orbits
190 * </p>
191 * @param e eccentricity
192 * @param E eccentric anomaly (rad)
193 * @return E - e sin(E)
194 */
195 private static double eMeSinE(final double e, final double E) {
196 double x = (1 - e) * FastMath.sin(E);
197 final double mE2 = -E * E;
198 double term = E;
199 double d = 0;
200 // the inequality test below IS intentional and should NOT be replaced by a
201 // check with a small tolerance
202 for (double x0 = Double.NaN; !Double.valueOf(x).equals(x0);) {
203 d += 2;
204 term *= mE2 / (d * (d + 1));
205 x0 = x;
206 x = x - term;
207 }
208 return x;
209 }
210
211 /**
212 * Computes the elliptic mean anomaly from the elliptic eccentric anomaly.
213 * @param e eccentricity such that 0 ≤ e < 1
214 * @param E elliptic eccentric anomaly (rad)
215 * @return elliptic mean anomaly (rad)
216 */
217 public static double ellipticEccentricToMean(final double e, final double E) {
218 return E - e * FastMath.sin(E);
219 }
220
221 /**
222 * Computes the hyperbolic true anomaly from the hyperbolic mean anomaly.
223 * @param e eccentricity > 1
224 * @param M hyperbolic mean anomaly
225 * @return hyperbolic true anomaly (rad)
226 */
227 public static double hyperbolicMeanToTrue(final double e, final double M) {
228 final double H = hyperbolicMeanToEccentric(e, M);
229 final double v = hyperbolicEccentricToTrue(e, H);
230 return v;
231 }
232
233 /**
234 * Computes the hyperbolic mean anomaly from the hyperbolic true anomaly.
235 * @param e eccentricity > 1
236 * @param v hyperbolic true anomaly (rad)
237 * @return hyperbolic mean anomaly
238 */
239 public static double hyperbolicTrueToMean(final double e, final double v) {
240 final double H = hyperbolicTrueToEccentric(e, v);
241 final double M = hyperbolicEccentricToMean(e, H);
242 return M;
243 }
244
245 /**
246 * Computes the hyperbolic true anomaly from the hyperbolic eccentric anomaly.
247 * @param e eccentricity > 1
248 * @param H hyperbolic eccentric anomaly
249 * @return hyperbolic true anomaly (rad)
250 */
251 public static double hyperbolicEccentricToTrue(final double e, final double H) {
252 return 2 * FastMath.atan(FastMath.sqrt((e + 1) / (e - 1)) * FastMath.tanh(0.5 * H));
253 }
254
255 /**
256 * Computes the hyperbolic eccentric anomaly from the hyperbolic true anomaly.
257 * @param e eccentricity > 1
258 * @param v hyperbolic true anomaly (rad)
259 * @return hyperbolic eccentric anomaly
260 */
261 public static double hyperbolicTrueToEccentric(final double e, final double v) {
262 final SinCos scv = FastMath.sinCos(v);
263 final double sinhH = FastMath.sqrt(e * e - 1) * scv.sin() / (1 + e * scv.cos());
264 return FastMath.asinh(sinhH);
265 }
266
267 /**
268 * Computes the hyperbolic eccentric anomaly from the hyperbolic mean anomaly.
269 * <p>
270 * The algorithm used here for solving hyperbolic Kepler equation is from
271 * Gooding, R.H., Odell, A.W. "The hyperbolic Kepler equation (and the elliptic
272 * equation revisited)." Celestial Mechanics 44, 267–282 (1988).
273 * https://doi.org/10.1007/BF01235540
274 * </p>
275 * @param e eccentricity > 1
276 * @param M hyperbolic mean anomaly
277 * @return hyperbolic eccentric anomaly
278 */
279 public static double hyperbolicMeanToEccentric(final double e, final double M) {
280 // Solve L = S - g * asinh(S) for S = sinh(H).
281 final double L = M / e;
282 final double g = 1.0 / e;
283 final double g1 = 1.0 - g;
284
285 // Starter based on Lagrange's theorem.
286 double S = L;
287 if (L == 0.0) {
288 return 0.0;
289 }
290 final double cl = FastMath.sqrt(1.0 + L * L);
291 final double al = FastMath.asinh(L);
292 final double w = g * g * al / (cl * cl * cl);
293 S = 1.0 - g / cl;
294 S = L + g * al / FastMath.cbrt(S * S * S + w * L * (1.5 - (4.0 / 3.0) * g));
295
296 // Two iterations (at most) of Halley-then-Newton process.
297 for (int i = 0; i < 2; ++i) {
298 final double s0 = S * S;
299 final double s1 = s0 + 1.0;
300 final double s2 = FastMath.sqrt(s1);
301 final double s3 = s1 * s2;
302 final double fdd = g * S / s3;
303 final double fddd = g * (1.0 - 2.0 * s0) / (s1 * s3);
304
305 double f;
306 double fd;
307 if (s0 / 6.0 + g1 >= 0.5) {
308 f = (S - g * FastMath.asinh(S)) - L;
309 fd = 1.0 - g / s2;
310 } else {
311 // Accurate computation of S - (1 - g1) * asinh(S)
312 // when (g1, S) is close to (0, 0).
313 final double t = S / (1.0 + FastMath.sqrt(1.0 + S * S));
314 final double tsq = t * t;
315 double x = S * (g1 + g * tsq);
316 double term = 2.0 * g * t;
317 double twoI1 = 1.0;
318 double x0;
319 int j = 0;
320 do {
321 if (++j == 1000000) {
322 // This isn't expected to happen, but it protects against an infinite loop.
323 throw new MathIllegalStateException(
324 OrekitMessages.UNABLE_TO_COMPUTE_HYPERBOLIC_ECCENTRIC_ANOMALY, j);
325 }
326 twoI1 += 2.0;
327 term *= tsq;
328 x0 = x;
329 x -= term / twoI1;
330 } while (x != x0);
331 f = x - L;
332 fd = (s0 / (s2 + 1.0) + g1) / s2;
333 }
334 final double ds = f * fd / (0.5 * f * fdd - fd * fd);
335 final double stemp = S + ds;
336 if (S == stemp) {
337 break;
338 }
339 f += ds * (fd + 0.5 * ds * (fdd + ds / 3.0 * fddd));
340 fd += ds * (fdd + 0.5 * ds * fddd);
341 S = stemp - f / fd;
342 }
343
344 final double H = FastMath.asinh(S);
345 return H;
346 }
347
348 /**
349 * Computes the hyperbolic mean anomaly from the hyperbolic eccentric anomaly.
350 * @param e eccentricity > 1
351 * @param H hyperbolic eccentric anomaly
352 * @return hyperbolic mean anomaly
353 */
354 public static double hyperbolicEccentricToMean(final double e, final double H) {
355 return e * FastMath.sinh(H) - H;
356 }
357
358 /**
359 * Convert anomaly.
360 * @param oldType old position angle type
361 * @param anomaly old value for anomaly
362 * @param e eccentricity
363 * @param newType new position angle type
364 * @return converted anomaly
365 * @since 12.2
366 */
367 public static double convertAnomaly(final PositionAngleType oldType, final double anomaly, final double e,
368 final PositionAngleType newType) {
369 return switch (newType) {
370 case MEAN -> getMeanAnomaly(e, anomaly, oldType);
371 case TRUE -> getTrueAnomaly(e, anomaly, oldType);
372 case ECCENTRIC -> getEccentricAnomaly(e, anomaly, oldType);
373 };
374 }
375
376 /** Get the true anomaly.
377 * @param e eccentricity
378 * @param anomaly anomaly (rad)
379 * @param positionAngleType anomaly type
380 * @return true anomaly (rad)
381 * @since 14.0
382 */
383 private static double getTrueAnomaly(final double e, final double anomaly,
384 final PositionAngleType positionAngleType) {
385 return switch (positionAngleType) {
386 case MEAN -> (e > 1.) ? hyperbolicMeanToTrue(e, anomaly) : ellipticMeanToTrue(e, anomaly);
387 case TRUE -> anomaly;
388 case ECCENTRIC -> (e > 1.) ? hyperbolicEccentricToTrue(e, anomaly) : ellipticEccentricToTrue(e, anomaly);
389 };
390 }
391
392 /** Get the eccentric anomaly.
393 * @param e eccentricity
394 * @param anomaly anomaly (rad)
395 * @param positionAngleType anomaly type
396 * @return eccentric anomaly (rad)
397 * @since 14.0
398 */
399 private static double getEccentricAnomaly(final double e, final double anomaly,
400 final PositionAngleType positionAngleType) {
401 return switch (positionAngleType) {
402 case MEAN -> (e > 1.) ? hyperbolicMeanToEccentric(e, anomaly) : ellipticMeanToEccentric(e, anomaly);
403 case TRUE -> (e > 1.) ? hyperbolicTrueToEccentric(e, anomaly) : ellipticTrueToEccentric(e, anomaly);
404 case ECCENTRIC -> anomaly;
405 };
406 }
407
408 /** Get the mean anomaly.
409 * @param e eccentricity
410 * @param anomaly anomaly (rad)
411 * @param positionAngleType anomaly type
412 * @return mean anomaly (rad)
413 * @since 14.0
414 */
415 private static double getMeanAnomaly(final double e, final double anomaly,
416 final PositionAngleType positionAngleType) {
417 return switch (positionAngleType) {
418 case MEAN -> anomaly;
419 case TRUE -> (e > 1.) ? hyperbolicTrueToMean(e, anomaly) : ellipticTrueToMean(e, anomaly);
420 case ECCENTRIC -> (e > 1.) ? hyperbolicEccentricToMean(e, anomaly) : ellipticEccentricToMean(e, anomaly);
421 };
422 }
423
424 /**
425 * Computes the difference in elliptic eccentric anomaly given that in elliptic mean anomaly.
426 * <p>This function solves Eq. 4.43 in Battin (1999).</p>
427 * @param s0 intermediate quantity {@code s0}
428 * @param c0 intermediate quantity {@code c0}
429 * @param deltaM difference in elliptic mean anomaly
430 * @return difference in elliptic eccentric anomaly
431 */
432 public static double ellipticMeanToEccentricDifference(final double s0, final double c0, final double deltaM) {
433
434 double deltaE = deltaM;
435 double delta;
436 boolean hasConverged;
437 int iter = 0;
438
439 do {
440 final double sE = FastMath.sin(deltaE);
441 final double cE = FastMath.cos(deltaE);
442
443 final double f = deltaE + s0 * (1.0 - cE) - c0 * sE - deltaM;
444 final double df = 1.0 + s0 * sE - c0 * cE;
445
446 delta = f / df;
447 deltaE -= delta;
448 hasConverged = FastMath.abs(delta) <= CONVERGENCE_THRESHOLD;
449 } while (++iter < MAXIMUM_ITERATIONS && !hasConverged);
450
451 if (!hasConverged) {
452 throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_ELLIPTIC_ECCENTRIC_ANOMALY, iter);
453 }
454
455 return deltaE;
456 }
457
458 /**
459 * Computes the difference in hyperbolic eccentric anomaly given that in hyperbolic mean anomaly.
460 * <p>This function solves Eq. 4.64 in Battin (1999).</p>
461 * @param s0 intermediate quantity {@code s0}
462 * @param c0 intermediate quantity {@code c0}
463 * @param deltaN difference in hyperbolic mean anomaly
464 * @return difference in hyperbolic eccentric anomaly
465 * @since 14.0
466 */
467 public static double hyperbolicMeanToEccentricDifference(final double s0, final double c0, final double deltaN) {
468
469 double shH = FastMath.signum(deltaN) * FastMath.min(FastMath.abs(FastMath.sinh(deltaN)), MAXIMUM_INITIAL_HYPERBOLIC_SINE);
470 double delta;
471 boolean hasConverged;
472 int iter = 0;
473
474 do {
475 final double deltaH = FastMath.asinh(shH);
476 final double chH = FastMath.sqrt(1.0 + shH * shH);
477
478 final double f = deltaH + s0 * (1.0 - chH) - c0 * shH + deltaN;
479 final double df = 1.0 / chH - (s0 * shH) / chH - c0;
480
481 delta = f / df;
482 shH -= delta;
483 hasConverged = FastMath.abs(delta) <= CONVERGENCE_THRESHOLD;
484 } while (++iter < MAXIMUM_ITERATIONS && !hasConverged);
485
486 if (!hasConverged) {
487 throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_HYPERBOLIC_ECCENTRIC_ANOMALY, iter);
488 }
489
490 return FastMath.asinh(shH);
491 }
492 }