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