From 9ed90178a8661b6fc0f5e65053cf34dec2e9aa7d Mon Sep 17 00:00:00 2001 From: ktyl Date: Wed, 3 Apr 2024 19:30:13 +0100 Subject: [PATCH 1/8] Add eccentric anomaly determination by binary search. --- include/astro/orbitalElementConversions.hpp | 73 +++++++++++++++++++++ 1 file changed, 73 insertions(+) diff --git a/include/astro/orbitalElementConversions.hpp b/include/astro/orbitalElementConversions.hpp index 56e86df..45d6c94 100644 --- a/include/astro/orbitalElementConversions.hpp +++ b/include/astro/orbitalElementConversions.hpp @@ -848,6 +848,78 @@ Real convertEllipticalMeanAnomalyToEccentricAnomaly( return eccentricAnomaly; } +// TODO: where the hell to put this ? +template int sign(T val) { + return (T(0) < val) - (val < T(0)); +} + +// TODO: nice and cool documentation ! +// what is the problem with the newton's method approach? +// +template +Real convertEllipticalMeanAnomalyToEccentricAnomalyBS( + const Real eccentricity, + const Real meanAnomaly) +{ + assert(eccentricity >= 0.0 && eccentricity < (1.0 - 1.0e-11)); + + const Real pi = 3.14159265358979323846; + + // Set mean anomaly to domain between 0 and 2pi. + Real meanAnomalyShifted = std::fmod(meanAnomaly, 2.0 * pi); + if (meanAnomalyShifted <0.0) + { + meanAnomalyShifted += 2.0 * pi; + } + + Real M = meanAnomalyShifted; + Real E = eccentricity; + // estimate for eccentric offset E + Real E0; + Real D; + + const int iterations = 33; + + // 110 F = SGN(M) : M = ABS(M)/(2*PI) + Real F = sign(M); + M = fabs(M) / (2.0 * pi); + // 120 M = (M-INT(M))*2*PI*F + M = (M - (int)M) * 2.0 * pi * F; + // 130 IF M < 0 THEN M = M+2*PI + if (M < 0) + { + M = M + 2.0 * pi; + } + + // 140 F = 1 + F = 1.0; + // 150 IF M > PI THEN F = -1 + // 160 IF M > PI THEN M = 2*PI - M + if (M > pi) + { + F = -1.0; + M = 2.0 * pi - M; + } + + // 170 E0 = PI/2 : D = PI/4 + E0 = pi * 0.5; + D = pi * 0.25; + // 180 FOR J = 1 TO 33 + for (int J = 0; J < iterations; J++) + { + // 190 M1 = E0 - E*SIN(E0) + Real M1 = E0 - E * sin(E0); + // 200 E0 = E0 + D*SGN(M-M1) : D = D/2 + E0 = E0 + D * sign(M - M1); + D *= 0.5; + } + + // 220 E0 = E0*F + E0 *= F; + + return E0; +} + } // namespace astro /*! @@ -856,4 +928,5 @@ Real convertEllipticalMeanAnomalyToEccentricAnomaly( * 2007. * Musegaas, P. Optimization of Space Trajectories Including Multiple Gravity Assists and Deep * Space Maneuvers. MSc thesis, Delft University of Technology, 2013. + * TODO: fill out reference! from Astronomical Algorithms p. 206 */ From eecc8714266f5bf4b2fefbb8a4529590f224d881 Mon Sep 17 00:00:00 2001 From: ktyl Date: Sat, 6 Apr 2024 11:44:11 +0100 Subject: [PATCH 2/8] Update reference. --- include/astro/orbitalElementConversions.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/include/astro/orbitalElementConversions.hpp b/include/astro/orbitalElementConversions.hpp index 45d6c94..50b281e 100644 --- a/include/astro/orbitalElementConversions.hpp +++ b/include/astro/orbitalElementConversions.hpp @@ -928,5 +928,9 @@ Real convertEllipticalMeanAnomalyToEccentricAnomalyBS( * 2007. * Musegaas, P. Optimization of Space Trajectories Including Multiple Gravity Assists and Deep * Space Maneuvers. MSc thesis, Delft University of Technology, 2013. - * TODO: fill out reference! from Astronomical Algorithms p. 206 + * TODO: in the text Meeus actually references Roger W. Sinnott's Sky & Telescope which was published + * a few years earlier. I have not yet been able to locate a copy to confirm the algorithm exists + * there as it does in AA, but this reference should likely be updated to the original developer + * of the algorithm rather than Meeus' redistribution of it. + * Meeus, J. Astronomical Algorithms. Second Edition, Willmann-Bell, Inc. 1991 */ From ddecdaee7dbb4b36c4d156215b6e48ee45a7a027 Mon Sep 17 00:00:00 2001 From: ktyl Date: Sun, 7 Apr 2024 19:47:30 +0100 Subject: [PATCH 3/8] Add function documentation. --- include/astro/orbitalElementConversions.hpp | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/include/astro/orbitalElementConversions.hpp b/include/astro/orbitalElementConversions.hpp index 50b281e..19a9e0e 100644 --- a/include/astro/orbitalElementConversions.hpp +++ b/include/astro/orbitalElementConversions.hpp @@ -853,15 +853,27 @@ template int sign(T val) { return (T(0) < val) - (val < T(0)); } -// TODO: nice and cool documentation ! -// what is the problem with the newton's method approach? -// +//! Convert elliptical mean anomaly to eccentric anomaly using a binary search. +/*! + * Converts mean anomaly to eccentric anomaly for elliptical orbits, + * for all eccentricities >= 0.0 and < 1.0. + * + * If the eccentricity falls outside the valid range, then a runtime + * exception is thrown. + * + * Also, note that the mean anomaly is automatically transformed to fit within the 0 to 2.0pi range. + * + * @tparam Real Real number type + * @param eccentricity Eccentricity [-] + * @param meanAnomaly Mean anomaly [rad] + * @return Eccentric anomaly [rad] + */ template Real convertEllipticalMeanAnomalyToEccentricAnomalyBS( const Real eccentricity, const Real meanAnomaly) { - assert(eccentricity >= 0.0 && eccentricity < (1.0 - 1.0e-11)); + assert(eccentricity >= 0.0 && eccentricity < 1.0); const Real pi = 3.14159265358979323846; From 0cec785587644859ef326c0fd9e190816fce7fd3 Mon Sep 17 00:00:00 2001 From: ktyl Date: Sun, 7 Apr 2024 19:50:58 +0100 Subject: [PATCH 4/8] Flatten utility function. --- include/astro/orbitalElementConversions.hpp | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/include/astro/orbitalElementConversions.hpp b/include/astro/orbitalElementConversions.hpp index 19a9e0e..4b65758 100644 --- a/include/astro/orbitalElementConversions.hpp +++ b/include/astro/orbitalElementConversions.hpp @@ -848,11 +848,6 @@ Real convertEllipticalMeanAnomalyToEccentricAnomaly( return eccentricAnomaly; } -// TODO: where the hell to put this ? -template int sign(T val) { - return (T(0) < val) - (val < T(0)); -} - //! Convert elliptical mean anomaly to eccentric anomaly using a binary search. /*! * Converts mean anomaly to eccentric anomaly for elliptical orbits, @@ -892,8 +887,8 @@ Real convertEllipticalMeanAnomalyToEccentricAnomalyBS( const int iterations = 33; - // 110 F = SGN(M) : M = ABS(M)/(2*PI) - Real F = sign(M); + // TODO: sign(M). This is used twice in the algorithm, and could be extracted to a separate function. + Real F = (Real(0) < M) - (M < Real(0)); M = fabs(M) / (2.0 * pi); // 120 M = (M-INT(M))*2*PI*F M = (M - (int)M) * 2.0 * pi * F; @@ -922,11 +917,11 @@ Real convertEllipticalMeanAnomalyToEccentricAnomalyBS( // 190 M1 = E0 - E*SIN(E0) Real M1 = E0 - E * sin(E0); // 200 E0 = E0 + D*SGN(M-M1) : D = D/2 - E0 = E0 + D * sign(M - M1); + Real val = M - M1; + E0 = E0 + D * ((Real(0) < val) - (val < Real(0))); D *= 0.5; } - // 220 E0 = E0*F E0 *= F; return E0; From 41b9ba0116658c9ee2485f5a1ec0387fc64fad00 Mon Sep 17 00:00:00 2001 From: ktyl Date: Sun, 7 Apr 2024 19:53:34 +0100 Subject: [PATCH 5/8] Remove redundant comments. --- include/astro/orbitalElementConversions.hpp | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/include/astro/orbitalElementConversions.hpp b/include/astro/orbitalElementConversions.hpp index 4b65758..44af878 100644 --- a/include/astro/orbitalElementConversions.hpp +++ b/include/astro/orbitalElementConversions.hpp @@ -881,7 +881,6 @@ Real convertEllipticalMeanAnomalyToEccentricAnomalyBS( Real M = meanAnomalyShifted; Real E = eccentricity; - // estimate for eccentric offset E Real E0; Real D; @@ -890,33 +889,24 @@ Real convertEllipticalMeanAnomalyToEccentricAnomalyBS( // TODO: sign(M). This is used twice in the algorithm, and could be extracted to a separate function. Real F = (Real(0) < M) - (M < Real(0)); M = fabs(M) / (2.0 * pi); - // 120 M = (M-INT(M))*2*PI*F M = (M - (int)M) * 2.0 * pi * F; - // 130 IF M < 0 THEN M = M+2*PI if (M < 0) { M = M + 2.0 * pi; } - // 140 F = 1 F = 1.0; - // 150 IF M > PI THEN F = -1 - // 160 IF M > PI THEN M = 2*PI - M if (M > pi) { F = -1.0; M = 2.0 * pi - M; } - // 170 E0 = PI/2 : D = PI/4 E0 = pi * 0.5; D = pi * 0.25; - // 180 FOR J = 1 TO 33 for (int J = 0; J < iterations; J++) { - // 190 M1 = E0 - E*SIN(E0) Real M1 = E0 - E * sin(E0); - // 200 E0 = E0 + D*SGN(M-M1) : D = D/2 Real val = M - M1; E0 = E0 + D * ((Real(0) < val) - (val < Real(0))); D *= 0.5; From 791fa3dc7de20260c2e5af09a25eb756d385182b Mon Sep 17 00:00:00 2001 From: ktyl Date: Sun, 7 Apr 2024 19:54:04 +0100 Subject: [PATCH 6/8] Add TODO. --- include/astro/orbitalElementConversions.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/astro/orbitalElementConversions.hpp b/include/astro/orbitalElementConversions.hpp index 44af878..68ee84f 100644 --- a/include/astro/orbitalElementConversions.hpp +++ b/include/astro/orbitalElementConversions.hpp @@ -884,6 +884,8 @@ Real convertEllipticalMeanAnomalyToEccentricAnomalyBS( Real E0; Real D; + // TODO: 33 iterations gets us to the limit of accuracy for a 10bit computer. We live in the + // future now, how many iterations do we need to get 64 bits of accuracy? const int iterations = 33; // TODO: sign(M). This is used twice in the algorithm, and could be extracted to a separate function. From 87c00e88f79e460c17bbd6208cc573ba5e0af5e0 Mon Sep 17 00:00:00 2001 From: ktyl Date: Mon, 8 Apr 2024 23:40:53 +0100 Subject: [PATCH 7/8] Configure BS iterations with function parameter. --- include/astro/orbitalElementConversions.hpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/include/astro/orbitalElementConversions.hpp b/include/astro/orbitalElementConversions.hpp index 68ee84f..2dc068b 100644 --- a/include/astro/orbitalElementConversions.hpp +++ b/include/astro/orbitalElementConversions.hpp @@ -856,17 +856,22 @@ Real convertEllipticalMeanAnomalyToEccentricAnomaly( * If the eccentricity falls outside the valid range, then a runtime * exception is thrown. * + * A binary search with 100 iterations should yield about 32 bits of precision. + * * Also, note that the mean anomaly is automatically transformed to fit within the 0 to 2.0pi range. * * @tparam Real Real number type + * @tparam Integer Integer type * @param eccentricity Eccentricity [-] * @param meanAnomaly Mean anomaly [rad] + * @param iterations Binary search iterations [-] * @return Eccentric anomaly [rad] */ -template +template Real convertEllipticalMeanAnomalyToEccentricAnomalyBS( - const Real eccentricity, - const Real meanAnomaly) + const Real eccentricity, + const Real meanAnomaly, + const Integer iterations = 100) { assert(eccentricity >= 0.0 && eccentricity < 1.0); @@ -884,10 +889,6 @@ Real convertEllipticalMeanAnomalyToEccentricAnomalyBS( Real E0; Real D; - // TODO: 33 iterations gets us to the limit of accuracy for a 10bit computer. We live in the - // future now, how many iterations do we need to get 64 bits of accuracy? - const int iterations = 33; - // TODO: sign(M). This is used twice in the algorithm, and could be extracted to a separate function. Real F = (Real(0) < M) - (M < Real(0)); M = fabs(M) / (2.0 * pi); From e12e008be28438ea30e49b06eb66452d83ad984e Mon Sep 17 00:00:00 2001 From: ktyl Date: Mon, 8 Apr 2024 23:51:09 +0100 Subject: [PATCH 8/8] Remove referencing TODO. The reference referred to in the TODO was: Sinnott, R. Sky and Telescope. Vol. 70, AAS Sky Publishing LLC. 1985. But my original citation includes the stated algorithm in full, so I think it is fine. --- include/astro/orbitalElementConversions.hpp | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/include/astro/orbitalElementConversions.hpp b/include/astro/orbitalElementConversions.hpp index 2dc068b..8b8d381 100644 --- a/include/astro/orbitalElementConversions.hpp +++ b/include/astro/orbitalElementConversions.hpp @@ -928,9 +928,5 @@ Real convertEllipticalMeanAnomalyToEccentricAnomalyBS( * 2007. * Musegaas, P. Optimization of Space Trajectories Including Multiple Gravity Assists and Deep * Space Maneuvers. MSc thesis, Delft University of Technology, 2013. - * TODO: in the text Meeus actually references Roger W. Sinnott's Sky & Telescope which was published - * a few years earlier. I have not yet been able to locate a copy to confirm the algorithm exists - * there as it does in AA, but this reference should likely be updated to the original developer - * of the algorithm rather than Meeus' redistribution of it. - * Meeus, J. Astronomical Algorithms. Second Edition, Willmann-Bell, Inc. 1991 + * Meeus, J. Astronomical Algorithms. Second Edition, Willmann-Bell, Inc. 1991. */