diff --git a/docs/html/astroAll_8hpp.html b/docs/html/astroAll_8hpp.html index c0de964..3f2cb53 100644 --- a/docs/html/astroAll_8hpp.html +++ b/docs/html/astroAll_8hpp.html @@ -105,42 +105,40 @@
- + - - - - - + + + + + - - - + + + - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + - - - + + +
diff --git a/docs/html/astroAll_8hpp__incl.map b/docs/html/astroAll_8hpp__incl.map index b853c30..97adb0d 100644 --- a/docs/html/astroAll_8hpp__incl.map +++ b/docs/html/astroAll_8hpp__incl.map @@ -1,38 +1,36 @@ - + - - - - - + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/html/astroAll_8hpp__incl.md5 b/docs/html/astroAll_8hpp__incl.md5 index fd1f73c..182acdf 100644 --- a/docs/html/astroAll_8hpp__incl.md5 +++ b/docs/html/astroAll_8hpp__incl.md5 @@ -1 +1 @@ -a2b28307db35fbd06780c502b423f098 \ No newline at end of file +635858882283ac34f0213486c71b7173 \ No newline at end of file diff --git a/docs/html/astroAll_8hpp__incl.png b/docs/html/astroAll_8hpp__incl.png index 8c026b2..11ca694 100644 Binary files a/docs/html/astroAll_8hpp__incl.png and b/docs/html/astroAll_8hpp__incl.png differ diff --git a/docs/html/namespaceastro.html b/docs/html/namespaceastro.html index fd0f857..f9639b8 100644 --- a/docs/html/namespaceastro.html +++ b/docs/html/namespaceastro.html @@ -502,9 +502,9 @@

Definition at line 82 of file orbitalElementConversions.hpp.

+

Definition at line 79 of file orbitalElementConversions.hpp.

-

References trueLongitudeOfPeriapsis, xPositionIndex, xVelocityIndex, yPositionIndex, yVelocityIndex, zPositionIndex, and zVelocityIndex.

+

References argumentOfLatitude, trueLongitude, trueLongitudeOfPeriapsis, xPositionIndex, xVelocityIndex, yPositionIndex, yVelocityIndex, zPositionIndex, and zVelocityIndex.

@@ -574,7 +574,7 @@

Definition at line 345 of file orbitalElementConversions.hpp.

+

Definition at line 316 of file orbitalElementConversions.hpp.

References argumentOfPeriapsisIndex, eccentricityIndex, inclinationIndex, longitudeOfAscendingNodeIndex, semiMajorAxisIndex, trueAnomalyIndex, xPositionIndex, xVelocityIndex, yPositionIndex, yVelocityIndex, zPositionIndex, and zVelocityIndex.

@@ -625,9 +625,9 @@

Returns
Elliptical eccentric anomaly [rad]
-

Definition at line 449 of file orbitalElementConversions.hpp.

+

Definition at line 420 of file orbitalElementConversions.hpp.

-

Referenced by convertTrueAnomalyToEccentricAnomaly().

+

Referenced by convertTrueAnomalyToEccentricAnomaly().

@@ -676,9 +676,9 @@

Returns
Hyperbolic eccentric anomaly [rad]
-

Definition at line 482 of file orbitalElementConversions.hpp.

+

Definition at line 453 of file orbitalElementConversions.hpp.

-

Referenced by convertTrueAnomalyToEccentricAnomaly().

+

Referenced by convertTrueAnomalyToEccentricAnomaly().

@@ -728,9 +728,9 @@

Returns
Eccentric anomaly [rad]
-

Definition at line 524 of file orbitalElementConversions.hpp.

+

Definition at line 495 of file orbitalElementConversions.hpp.

-

References convertTrueAnomalyToEllipticalEccentricAnomaly(), and convertTrueAnomalyToHyperbolicEccentricAnomaly().

+

References convertTrueAnomalyToEllipticalEccentricAnomaly(), and convertTrueAnomalyToHyperbolicEccentricAnomaly().

@@ -779,9 +779,9 @@

Returns
Mean anomaly [rad]
-

Definition at line 563 of file orbitalElementConversions.hpp.

+

Definition at line 534 of file orbitalElementConversions.hpp.

-

Referenced by convertEccentricAnomalyToMeanAnomaly().

+

Referenced by convertEccentricAnomalyToMeanAnomaly().

@@ -830,9 +830,9 @@

Returns
Mean anomaly [rad]
-

Definition at line 587 of file orbitalElementConversions.hpp.

+

Definition at line 558 of file orbitalElementConversions.hpp.

-

Referenced by convertEccentricAnomalyToMeanAnomaly().

+

Referenced by convertEccentricAnomalyToMeanAnomaly().

@@ -882,9 +882,9 @@

Returns
Mean anomaly [rad]
-

Definition at line 615 of file orbitalElementConversions.hpp.

+

Definition at line 586 of file orbitalElementConversions.hpp.

-

References convertEllipticalEccentricAnomalyToMeanAnomaly(), and convertHyperbolicEccentricAnomalyToMeanAnomaly().

+

References convertEllipticalEccentricAnomalyToMeanAnomaly(), and convertHyperbolicEccentricAnomalyToMeanAnomaly().

@@ -931,9 +931,9 @@

Returns
True anomaly [rad]
-

Definition at line 651 of file orbitalElementConversions.hpp.

+

Definition at line 622 of file orbitalElementConversions.hpp.

-

Referenced by convertEccentricAnomalyToTrueAnomaly().

+

Referenced by convertEccentricAnomalyToTrueAnomaly().

@@ -980,9 +980,9 @@

Returns
True anomaly [rad]
-

Definition at line 676 of file orbitalElementConversions.hpp.

+

Definition at line 647 of file orbitalElementConversions.hpp.

-

Referenced by convertEccentricAnomalyToTrueAnomaly().

+

Referenced by convertEccentricAnomalyToTrueAnomaly().

@@ -1032,9 +1032,9 @@

Returns
True anomaly [rad]
-

Definition at line 711 of file orbitalElementConversions.hpp.

+

Definition at line 682 of file orbitalElementConversions.hpp.

-

References convertEllipticalEccentricAnomalyToTrueAnomaly(), and convertHyperbolicEccentricAnomalyToTrueAnomaly().

+

References convertEllipticalEccentricAnomalyToTrueAnomaly(), and convertHyperbolicEccentricAnomalyToTrueAnomaly().

@@ -1097,9 +1097,9 @@

Returns
Kepler equation value for given elliptical orbit [rad]
-

Definition at line 757 of file orbitalElementConversions.hpp.

+

Definition at line 728 of file orbitalElementConversions.hpp.

-

Referenced by convertEllipticalMeanAnomalyToEccentricAnomaly().

+

Referenced by convertEllipticalMeanAnomalyToEccentricAnomaly().

@@ -1153,9 +1153,9 @@

Returns
First-derivative of Kepler's function for elliptical orbits [rad]
-

Definition at line 783 of file orbitalElementConversions.hpp.

+

Definition at line 754 of file orbitalElementConversions.hpp.

-

Referenced by convertEllipticalMeanAnomalyToEccentricAnomaly().

+

Referenced by convertEllipticalMeanAnomalyToEccentricAnomaly().

@@ -1221,9 +1221,9 @@

Returns
Eccentric anomaly [rad]
-

Definition at line 813 of file orbitalElementConversions.hpp.

+

Definition at line 784 of file orbitalElementConversions.hpp.

-

References computeEllipticalKeplerFunction(), and computeFirstDerivativeEllipticalKeplerFunction().

+

References computeEllipticalKeplerFunction(), and computeFirstDerivativeEllipticalKeplerFunction().

diff --git a/docs/html/orbitalElementConversions_8hpp.html b/docs/html/orbitalElementConversions_8hpp.html index 52ec6bb..73bf55a 100644 --- a/docs/html/orbitalElementConversions_8hpp.html +++ b/docs/html/orbitalElementConversions_8hpp.html @@ -103,29 +103,26 @@ #include <stdexcept>
#include <vector>
#include "astro/stateVectorIndices.hpp"
-#include <iostream>
Include dependency graph for orbitalElementConversions.hpp:
- + - + - + - + - + - + - + - - - +
diff --git a/docs/html/orbitalElementConversions_8hpp__incl.map b/docs/html/orbitalElementConversions_8hpp__incl.map index e8525bf..353def8 100644 --- a/docs/html/orbitalElementConversions_8hpp__incl.map +++ b/docs/html/orbitalElementConversions_8hpp__incl.map @@ -1,19 +1,17 @@ - + - + - + - + - + - + - + - - - + diff --git a/docs/html/orbitalElementConversions_8hpp__incl.md5 b/docs/html/orbitalElementConversions_8hpp__incl.md5 index a5e8a28..554a8f7 100644 --- a/docs/html/orbitalElementConversions_8hpp__incl.md5 +++ b/docs/html/orbitalElementConversions_8hpp__incl.md5 @@ -1 +1 @@ -309d890323d94d11feb0af9ef6329cc6 \ No newline at end of file +53a07e996ba6c52db59723bba567e423 \ No newline at end of file diff --git a/docs/html/orbitalElementConversions_8hpp__incl.png b/docs/html/orbitalElementConversions_8hpp__incl.png index 253eac1..d4fcd7f 100644 Binary files a/docs/html/orbitalElementConversions_8hpp__incl.png and b/docs/html/orbitalElementConversions_8hpp__incl.png differ diff --git a/docs/html/orbitalElementConversions_8hpp_source.html b/docs/html/orbitalElementConversions_8hpp_source.html index f0f6dcb..d675210 100644 --- a/docs/html/orbitalElementConversions_8hpp_source.html +++ b/docs/html/orbitalElementConversions_8hpp_source.html @@ -113,627 +113,600 @@
18
20
-
21// (debug)
-
22#include <iostream>
+
21namespace astro
+
22{
23
-
24namespace astro
-
25{
-
26
-
28
-
81template <typename Real, typename Vector6>
- -
83 const Vector6& cartesianElements,
-
84 const Real gravitationalParameter,
-
85 const Real tolerance = 10.0 * std::numeric_limits<Real>::epsilon())
-
86{
-
87 assert(cartesianElements.size() == 6);
-
88 assert(gravitationalParameter > 0.0);
+
25
+
78template <typename Real, typename Vector6>
+ +
80 const Vector6& cartesianElements,
+
81 const Real gravitationalParameter,
+
82 const Real tolerance = 10.0 * std::numeric_limits<Real>::epsilon())
+
83{
+
84 assert(cartesianElements.size() == 6);
+
85 assert(gravitationalParameter > 0.0);
+
86
+
87 const Real pi = 3.14159265358979323846;
+
88 const Real gravitationalParameterInverse = 1.0 / gravitationalParameter;
89
-
90 const Real pi = 3.14159265358979323846;
-
91 const Real gravitationalParameterInverse = 1.0 / gravitationalParameter;
+
90 typedef std::vector<Real> Vector;
+
91 Vector6 keplerianElements = cartesianElements;
92
-
93 typedef std::vector<Real> Vector;
-
94 Vector6 keplerianElements = cartesianElements;
-
95
-
96 // Cartesian position
-
97 const Vector position = Vector({cartesianElements[xPositionIndex],
-
98 cartesianElements[yPositionIndex],
-
99 cartesianElements[zPositionIndex]});
-
100 const Real positionNormSquared = position[0] * position[0]
-
101 + position[1] * position[1]
-
102 + position[2] * position[2];
-
103 const Real positionNorm = std::sqrt(positionNormSquared);
-
104 const Real positionNormInverse = 1.0 / positionNorm;
-
105
-
106std::cout << "(debug) r: " << positionNorm << std::endl;
-
107
-
108 // Cartesian velocity
-
109 const Vector velocity = Vector({cartesianElements[xVelocityIndex],
-
110 cartesianElements[yVelocityIndex],
-
111 cartesianElements[zVelocityIndex]});
-
112 const Real velocityNormSquared = velocity[0] * velocity[0]
-
113 + velocity[1] * velocity[1]
-
114 + velocity[2] * velocity[2];
-
115 const Real velocityNorm = std::sqrt(velocityNormSquared);
-
116
-
117std::cout << "(debug) v: " << std::sqrt(velocityNormSquared) << std::endl;
-
118
-
119 // Angular momentum
-
120 const Vector angularMomentum
-
121 = Vector({position[1] * velocity[2] - position[2] * velocity[1],
-
122 position[2] * velocity[0] - position[0] * velocity[2],
-
123 position[0] * velocity[1] - position[1] * velocity[0]});
-
124 const Real angularMomentumNormSquared = angularMomentum[0] * angularMomentum[0]
-
125 + angularMomentum[1] * angularMomentum[1]
-
126 + angularMomentum[2] * angularMomentum[2];
-
127 const Real angularMomentumNorm = std::sqrt(angularMomentumNormSquared);
-
128 const Vector angularMomentumUnitVector
-
129 = Vector({angularMomentum[0] / angularMomentumNorm,
-
130 angularMomentum[1] / angularMomentumNorm,
-
131 angularMomentum[2] / angularMomentumNorm});
-
132std::cout << "(debug) h: " << angularMomentum[0] << ","
-
133 << angularMomentum[1] << ","
-
134 << angularMomentum[2] << std::endl;
+
93 // Cartesian position
+
94 const Vector position = Vector({cartesianElements[xPositionIndex],
+
95 cartesianElements[yPositionIndex],
+
96 cartesianElements[zPositionIndex]});
+
97 const Real positionNormSquared = position[0] * position[0]
+
98 + position[1] * position[1]
+
99 + position[2] * position[2];
+
100 const Real positionNorm = std::sqrt(positionNormSquared);
+
101 const Real positionNormInverse = 1.0 / positionNorm;
+
102
+
103 // Cartesian velocity
+
104 const Vector velocity = Vector({cartesianElements[xVelocityIndex],
+
105 cartesianElements[yVelocityIndex],
+
106 cartesianElements[zVelocityIndex]});
+
107 const Real velocityNormSquared = velocity[0] * velocity[0]
+
108 + velocity[1] * velocity[1]
+
109 + velocity[2] * velocity[2];
+
110 const Real velocityNorm = std::sqrt(velocityNormSquared);
+
111
+
112 // Angular momentum
+
113 const Vector angularMomentum
+
114 = Vector({position[1] * velocity[2] - position[2] * velocity[1],
+
115 position[2] * velocity[0] - position[0] * velocity[2],
+
116 position[0] * velocity[1] - position[1] * velocity[0]});
+
117 const Real angularMomentumNormSquared = angularMomentum[0] * angularMomentum[0]
+
118 + angularMomentum[1] * angularMomentum[1]
+
119 + angularMomentum[2] * angularMomentum[2];
+
120 const Real angularMomentumNorm = std::sqrt(angularMomentumNormSquared);
+
121 const Vector angularMomentumUnitVector
+
122 = Vector({angularMomentum[0] / angularMomentumNorm,
+
123 angularMomentum[1] / angularMomentumNorm,
+
124 angularMomentum[2] / angularMomentumNorm});
+
125
+
126 // Eccentricity
+
127 const Real eccentricityVectorFirsTermMultiplier
+
128 = velocityNormSquared * gravitationalParameterInverse - positionNormInverse;
+
129
+
130 const Real positionDotVelocity = position[0] * velocity[0]
+
131 + position[1] * velocity[1]
+
132 + position[2] * velocity[2];
+
133
+
134 const Real positionDotVelocityScaled = positionDotVelocity / gravitationalParameter;
135
-
136std::cout << "(debug) hx: " << angularMomentum[0] << std::endl;
-
137std::cout << "(debug) hy: " << angularMomentum[1] << std::endl;
-
138std::cout << "(debug) hz: " << angularMomentum[2] << std::endl;
-
139std::cout << "(debug) h: " << angularMomentumNorm << std::endl;std::cout << "(debug) hx: " << angularMomentum[0] << std::endl;
-
140std::cout << "(debug) hy: " << angularMomentum[1] << std::endl;
-
141std::cout << "(debug) hz: " << angularMomentum[2] << std::endl;
-
142std::cout << "(debug) h: " << angularMomentumNorm << std::endl;
+
136 const Vector eccentricityVector
+
137 = Vector({eccentricityVectorFirsTermMultiplier * position[0]
+
138 - positionDotVelocity * gravitationalParameterInverse * velocity[0],
+
139 eccentricityVectorFirsTermMultiplier * position[1]
+
140 - positionDotVelocity * gravitationalParameterInverse * velocity[1],
+
141 eccentricityVectorFirsTermMultiplier * position[2]
+
142 - positionDotVelocity * gravitationalParameterInverse * velocity[2]});
143
-
144 // Eccentricity
-
145 const Real eccentricityVectorFirsTermMultiplier
-
146 = velocityNormSquared * gravitationalParameterInverse - positionNormInverse;
+
144 const Real eccentricityNormSquared = eccentricityVector[0] * eccentricityVector[0]
+
145 + eccentricityVector[1] * eccentricityVector[1]
+
146 + eccentricityVector[2] * eccentricityVector[2];
147
-
148 const Real positionDotVelocity = position[0] * velocity[0]
-
149 + position[1] * velocity[1]
-
150 + position[2] * velocity[2];
-
151
-
152 const Real positionDotVelocityScaled = positionDotVelocity / gravitationalParameter;
-
153
-
154 const Vector eccentricityVector
-
155 = Vector({eccentricityVectorFirsTermMultiplier * position[0]
-
156 - positionDotVelocity * gravitationalParameterInverse * velocity[0],
-
157 eccentricityVectorFirsTermMultiplier * position[1]
-
158 - positionDotVelocity * gravitationalParameterInverse * velocity[1],
-
159 eccentricityVectorFirsTermMultiplier * position[2]
-
160 - positionDotVelocity * gravitationalParameterInverse * velocity[2]});
-
161
-
162 const Real eccentricityNormSquared = eccentricityVector[0] * eccentricityVector[0]
-
163 + eccentricityVector[1] * eccentricityVector[1]
-
164 + eccentricityVector[2] * eccentricityVector[2];
-
165
-
166 const Real eccentricity = std::sqrt(eccentricityNormSquared);
-
167 keplerianElements[1] = eccentricity;
-
168
-
169std::cout << "(debug) ex: " << eccentricityVector[0] << std::endl;
-
170std::cout << "(debug) ey: " << eccentricityVector[1] << std::endl;
-
171std::cout << "(debug) ez: " << eccentricityVector[2] << std::endl;
-
172std::cout << "(debug) e: " << keplerianElements[1] << std::endl;
-
173
-
174 // Specifc total energy using the vis-viva equation
-
175 const Real specificTotalEnergy
-
176 = 0.5 * velocityNormSquared - gravitationalParameter / positionNorm;
-
177
-
178std::cout << "(debug) Esp: " << specificTotalEnergy << std::endl;
+
148 const Real eccentricity = std::sqrt(eccentricityNormSquared);
+
149 keplerianElements[1] = eccentricity;
+
150
+
151 // Specifc total energy using the vis-viva equation
+
152 const Real specificTotalEnergy
+
153 = 0.5 * velocityNormSquared - gravitationalParameter / positionNorm;
+
154
+
155 // Semi-major axis
+
156 Real semiMajorAxis = std::numeric_limits<Real>::quiet_NaN();
+
157 Real semiLatusRectum = std::numeric_limits<Real>::quiet_NaN();
+
158 if (std::fabs(eccentricity - 1.0) > tolerance)
+
159 {
+
160 semiMajorAxis = -0.5 * gravitationalParameter / specificTotalEnergy;
+
161 semiLatusRectum = semiMajorAxis * (1.0 - eccentricityNormSquared);
+
162 keplerianElements[0] = semiMajorAxis;
+
163 }
+
164 else
+
165 {
+
166 semiMajorAxis = std::numeric_limits<Real>::infinity();
+
167 semiLatusRectum = angularMomentumNormSquared * gravitationalParameterInverse;
+
168 keplerianElements[0] = semiLatusRectum;
+
169 }
+
170 assert(semiMajorAxis > 0.0);
+
171
+
172 // Inclination
+
173 const Real inclination = std::acos(angularMomentumUnitVector[2]);
+
174 keplerianElements[2] = inclination;
+
175
+
176 // Ascending node vector
+
177 const Vector ascendingNodeVector
+
178 = Vector({-angularMomentum[1], angularMomentum[0], 0.0});
179
-
180 // Semi-major axis
-
181 Real semiMajorAxis = std::numeric_limits<Real>::quiet_NaN();
-
182 Real semiLatusRectum = std::numeric_limits<Real>::quiet_NaN();
-
183 if (std::fabs(eccentricity - 1.0) > tolerance)
-
184 {
-
185 semiMajorAxis = -0.5 * gravitationalParameter / specificTotalEnergy;
-
186 semiLatusRectum = semiMajorAxis * (1.0 - eccentricityNormSquared);
-
187 keplerianElements[0] = semiMajorAxis;
-
188 }
-
189 else
-
190 {
-
191 semiMajorAxis = std::numeric_limits<Real>::infinity();
-
192 semiLatusRectum = angularMomentumNormSquared * gravitationalParameterInverse;
-
193 keplerianElements[0] = semiLatusRectum;
-
194 }
-
195 assert(semiMajorAxis > 0.0);
-
196std::cout << "(debug) a: " << semiMajorAxis << std::endl;
-
197
-
198std::cout << "(debug) a: " << keplerianElements[0] << std::endl;
-
199std::cout << "(debug) p: " << semiLatusRectum << std::endl;
-
200
-
201 // Inclination
-
202 const Real inclination = std::acos(angularMomentumUnitVector[2]);
-
203 keplerianElements[2] = inclination;
-
204
-
205std::cout << "(debug) i: " << keplerianElements[2] / pi * 180.0 << std::endl;
-
206
-
207 // Ascending node vector
-
208 const Vector ascendingNodeVector
-
209 = Vector({-angularMomentum[1], angularMomentum[0], 0.0});
-
210std::cout << "(debug) asc. node vec.: " << ascendingNodeVector[0] << ", "
-
211 << ascendingNodeVector[1] << ", "
-
212 << ascendingNodeVector[2] << std::endl;
-
213
-
214 const Real ascendingNodeVectorNormSquared
-
215 = ascendingNodeVector[0] * ascendingNodeVector[0]
-
216 + ascendingNodeVector[1] * ascendingNodeVector[1]
-
217 + ascendingNodeVector[2] * ascendingNodeVector[2];
-
218
-
219 const Real ascendingNodeVectorNorm = std::sqrt(ascendingNodeVectorNormSquared);
-
220
-
221 const Vector ascendingNodeUnitVector
-
222 = Vector({ascendingNodeVector[0] / ascendingNodeVectorNorm,
-
223 ascendingNodeVector[1] / ascendingNodeVectorNorm,
-
224 ascendingNodeVector[2] / ascendingNodeVectorNorm });
-
225
-
226std::cout << "(debug) nx: " << ascendingNodeVector[0] << std::endl;
-
227std::cout << "(debug) ny: " << ascendingNodeVector[1] << std::endl;
-
228std::cout << "(debug) nz: " << ascendingNodeVector[2] << std::endl;
-
229std::cout << "(debug) n: " << ascendingNodeVectorNorm << std::endl;
-
230
-
231 // Longitude of ascending node
-
232 Real longitudeOfAscendingNode = std::acos(ascendingNodeUnitVector[0]);
-
233
-
234 if (ascendingNodeVector[1] < 0.0)
-
235 {
-
236 longitudeOfAscendingNode = 2.0 * pi - longitudeOfAscendingNode;
-
237 }
-
238std::cout << "(debug) RAAN: " << longitudeOfAscendingNode << std::endl;
+
180 const Real ascendingNodeVectorNormSquared
+
181 = ascendingNodeVector[0] * ascendingNodeVector[0]
+
182 + ascendingNodeVector[1] * ascendingNodeVector[1]
+
183 + ascendingNodeVector[2] * ascendingNodeVector[2];
+
184
+
185 const Real ascendingNodeVectorNorm = std::sqrt(ascendingNodeVectorNormSquared);
+
186
+
187 const Vector ascendingNodeUnitVector
+
188 = Vector({ascendingNodeVector[0] / ascendingNodeVectorNorm,
+
189 ascendingNodeVector[1] / ascendingNodeVectorNorm,
+
190 ascendingNodeVector[2] / ascendingNodeVectorNorm });
+
191
+
192 // Longitude of ascending node
+
193 Real longitudeOfAscendingNode = std::acos(ascendingNodeUnitVector[0]);
+
194
+
195 if (ascendingNodeVector[1] < 0.0)
+
196 {
+
197 longitudeOfAscendingNode = 2.0 * pi - longitudeOfAscendingNode;
+
198 }
+
199
+
200 keplerianElements[4] = longitudeOfAscendingNode;
+
201
+
202 // Argument of periapsis
+
203 const Real ascendingNodeVectorDotEccentricityVector
+
204 = ascendingNodeVector[0] * eccentricityVector[0]
+
205 + ascendingNodeVector[1] * eccentricityVector[1]
+
206 + ascendingNodeVector[2] * eccentricityVector[2];
+
207
+
208 Real argumentOfPeriapsis = std::acos(ascendingNodeVectorDotEccentricityVector
+
209 / (ascendingNodeVectorNorm * eccentricity));
+
210
+
211 if (eccentricityVector[2] < 0.0)
+
212 {
+
213 argumentOfPeriapsis = 2.0 * pi - argumentOfPeriapsis;
+
214 }
+
215
+
216 keplerianElements[3] = argumentOfPeriapsis;
+
217
+
218 // True anomaly
+
219 const Real eccentricityVectorDotPosition = eccentricityVector[0] * position[0]
+
220 + eccentricityVector[1] * position[1]
+
221 + eccentricityVector[2] * position[2];
+
222
+
223 Real trueAnomaly = std::acos(eccentricityVectorDotPosition / (eccentricity * positionNorm));
+
224
+
225 if (positionDotVelocity < 0.0)
+
226 {
+
227 trueAnomaly = 2.0 * pi - trueAnomaly;
+
228 }
+
229
+
230 keplerianElements[5] = trueAnomaly;
+
231
+
232 // True longitude of periapsis
+
233 Real trueLongitudeOfPeriapsis = std::acos(eccentricityVector[0] / eccentricity);
+
234
+
235 if (eccentricityVector[1] < tolerance)
+
236 {
+ +
238 }
239
-
240 keplerianElements[4] = longitudeOfAscendingNode;
-
241
-
242std::cout << "(debug) raan: " << keplerianElements[4] / pi * 180.0 << std::endl;
-
243
-
244 // Argument of periapsis
-
245 const Real ascendingNodeVectorDotEccentricityVector
-
246 = ascendingNodeVector[0] * eccentricityVector[0]
-
247 + ascendingNodeVector[1] * eccentricityVector[1]
-
248 + ascendingNodeVector[2] * eccentricityVector[2];
-
249
-
250 Real argumentOfPeriapsis = std::acos(ascendingNodeVectorDotEccentricityVector
-
251 / (ascendingNodeVectorNorm * eccentricity));
-
252
-
253 if (eccentricityVector[2] < 0.0)
-
254 {
-
255 argumentOfPeriapsis = 2.0 * pi - argumentOfPeriapsis;
-
256 }
-
257
-
258 keplerianElements[3] = argumentOfPeriapsis;
+
240 // Special case: elliptical, equatorial
+
241 if(std::fabs(eccentricity) > tolerance && std::fabs(inclination) < tolerance)
+
242 {
+
243 longitudeOfAscendingNode = std::numeric_limits<Real>::quiet_NaN();
+
244 keplerianElements[4] = trueLongitudeOfPeriapsis;
+
245 }
+
246
+
247 // Argument of latitude
+
248 const Real ascendingNodeVectorDotPosition = ascendingNodeVector[0] * position[0]
+
249 + ascendingNodeVector[1] * position[1]
+
250 + ascendingNodeVector[2] * position[2];
+
251
+
252 Real argumentOfLatitude = std::acos(ascendingNodeVectorDotPosition
+
253 / (ascendingNodeVectorNorm * positionNorm));
+
254
+
255 if (position[2] < 0.0)
+
256 {
+ +
258 }
259
-
260std::cout << "(debug) omega: " << keplerianElements[3] / pi * 180.0 << std::endl;
-
261
-
262 // True anomaly
-
263 const Real eccentricityVectorDotPosition = eccentricityVector[0] * position[0]
-
264 + eccentricityVector[1] * position[1]
-
265 + eccentricityVector[2] * position[2];
+
260 // Special case: circular, inclined
+
261 if(std::fabs(eccentricity) < tolerance && std::fabs(inclination) > tolerance)
+
262 {
+
263 argumentOfPeriapsis = std::numeric_limits<Real>::quiet_NaN();
+
264 keplerianElements[3] = argumentOfLatitude;
+
265 }
266
-
267 Real trueAnomaly = std::acos(eccentricityVectorDotPosition / (eccentricity * positionNorm));
-
268
-
269 if (positionDotVelocity < 0.0)
+
267 // True longitude
+
268 Real trueLongitude = std::acos(position[0] / positionNorm);
+
269 if (position[1] < 0.0)
270 {
-
271 trueAnomaly = 2.0 * pi - trueAnomaly;
+
271 trueLongitude = 2.0 * pi - trueLongitude;
272 }
273
-
274 keplerianElements[5] = trueAnomaly;
-
275
-
276std::cout << "(debug) ta: " << keplerianElements[5] / pi * 180.0 << std::endl;
-
277
-
278 // Special cases
-
279 // Elliptical, equatorial
-
280 Real trueLongitudeOfPeriapsis = std::acos(eccentricityVector[0]/eccentricity);
-
281 if (std::fabs(eccentricityVector[1]) < tolerance)
-
282 {
- -
284 longitudeOfAscendingNode = std::numeric_limits<Real>::quiet_NaN();
-
285 keplerianElements[4] = trueLongitudeOfPeriapsis;
-
286 }
-
287
-
288 // // Circular, inclined
-
289 // const Real ascendingNodeVectorDotPosition = ascendingNodeVector[0] * position[0]
-
290 // + ascendingNodeVector[1] * position[1]
-
291 // + ascendingNodeVector[2] * position[2];
-
292 // Real argumentOfLatitude = std::acos(ascendingNodeVectorDotPosition
-
293 // / (ascendingNodeVectorNorm * positionNorm));
-
294 // if (position[2] < 0.0)
-
295 // {
-
296 // argumentOfLatitude = 2.0 * pi - argumentOfLatitude;
-
297 // longitudeOfAscendingNode = std::numeric_limits<Real>::quiet_NaN();
-
298 // keplerianElements[4] = argumentOfLatitude;
-
299 // }
-
300
-
301 // // Circular, equatorial
-
302 // Real trueLongitude = std::acos(position[0] / positionNorm);
-
303 // if (position[1] < 0.0)
-
304 // {
-
305 // trueLongitude = 2.0 * pi - trueLongitude;
-
306 // longitudeOfAscendingNode = std::numeric_limits<Real>::quiet_NaN();
-
307 // argumentOfPeriapsis = std::numeric_limits<Real>::quiet_NaN();
-
308 // keplerianElements[5] = trueLongitude;
-
309 // }
-
310
-
311 return keplerianElements;
-
312}
-
313
-
315
-
344template <typename Real, typename Vector6>
- -
346 const Vector6& keplerianElements, const Real gravitationalParameter,
-
347 const Real tolerance = 10.0 * std::numeric_limits<Real>::epsilon())
-
348{
-
349 Vector6 cartesianElements = keplerianElements;
-
350
-
351 const Real semiMajorAxis = keplerianElements[semiMajorAxisIndex];
-
352 const Real eccentricity = keplerianElements[eccentricityIndex];
-
353 const Real inclination = keplerianElements[inclinationIndex];
-
354 const Real argumentOfPeriapsis = keplerianElements[argumentOfPeriapsisIndex];
-
355 const Real longitudeOfAscendingNode = keplerianElements[longitudeOfAscendingNodeIndex];
-
356 const Real trueAnomaly = keplerianElements[trueAnomalyIndex];
-
357
-
358 // Pre-compute sines and cosines of angles for efficient computation.
-
359 const Real cosineOfInclination = std::cos(inclination);
-
360 const Real sineOfInclination = std::sin(inclination);
-
361 const Real cosineOfArgumentOfPeriapsis = std::cos(argumentOfPeriapsis);
-
362 const Real sineOfArgumentOfPeriapsis = std::sin(argumentOfPeriapsis);
-
363 const Real cosineOfLongitudeOfAscendingNode = std::cos(longitudeOfAscendingNode);
-
364 const Real sineOfLongitudeOfAscendingNode = std::sin(longitudeOfAscendingNode);
-
365 const Real cosineOfTrueAnomaly = std::cos(trueAnomaly);
-
366 const Real sineOfTrueAnomaly = std::sin(trueAnomaly);
-
367
-
368 // Compute semi-latus rectum in the case the orbit is not a parabola.
-
369 Real semiLatusRectum = 0.0;
-
370 if (std::fabs(eccentricity - 1.0) > tolerance)
-
371 {
-
372 semiLatusRectum = semiMajorAxis * (1.0 - eccentricity * eccentricity);
-
373 }
-
374
-
375 // Else set the semi-latus rectum as the first element in the vector of Keplerian elements.
-
376 else
-
377 {
-
378 semiLatusRectum = keplerianElements[0];
-
379 }
-
380
-
381 // Compute the magnitude of the orbital radius, measured from the focal point.
-
382 const Real radiusMagnitude = semiLatusRectum / (1.0 + eccentricity * cosineOfTrueAnomaly);
-
383
-
384 // Define position and velocity in the perifocal coordinate system.
-
385 const Real xPositionPerifocal = radiusMagnitude * cosineOfTrueAnomaly;
-
386 const Real yPositionPerifocal = radiusMagnitude * sineOfTrueAnomaly;
-
387 const Real xVelocityPerifocal
-
388 = -std::sqrt(gravitationalParameter / semiLatusRectum) * sineOfTrueAnomaly;
-
389 const Real yVelocityPerifocal
-
390 = std::sqrt(gravitationalParameter / semiLatusRectum)
-
391 * (eccentricity + cosineOfTrueAnomaly);
+
274 // // Special case: circular, equatorial
+
275 if (std::fabs(eccentricity) < tolerance && std::fabs(inclination) < tolerance)
+
276 {
+
277 argumentOfPeriapsis = std::numeric_limits<Real>::quiet_NaN();
+
278 longitudeOfAscendingNode = std::numeric_limits<Real>::quiet_NaN();
+
279 keplerianElements[5] = trueLongitude;
+
280 }
+
281
+
282 return keplerianElements;
+
283}
+
284
+
286
+
315template <typename Real, typename Vector6>
+ +
317 const Vector6& keplerianElements, const Real gravitationalParameter,
+
318 const Real tolerance = 10.0 * std::numeric_limits<Real>::epsilon())
+
319{
+
320 Vector6 cartesianElements = keplerianElements;
+
321
+
322 const Real semiMajorAxis = keplerianElements[semiMajorAxisIndex];
+
323 const Real eccentricity = keplerianElements[eccentricityIndex];
+
324 const Real inclination = keplerianElements[inclinationIndex];
+
325 const Real argumentOfPeriapsis = keplerianElements[argumentOfPeriapsisIndex];
+
326 const Real longitudeOfAscendingNode = keplerianElements[longitudeOfAscendingNodeIndex];
+
327 const Real trueAnomaly = keplerianElements[trueAnomalyIndex];
+
328
+
329 // Pre-compute sines and cosines of angles for efficient computation.
+
330 const Real cosineOfInclination = std::cos(inclination);
+
331 const Real sineOfInclination = std::sin(inclination);
+
332 const Real cosineOfArgumentOfPeriapsis = std::cos(argumentOfPeriapsis);
+
333 const Real sineOfArgumentOfPeriapsis = std::sin(argumentOfPeriapsis);
+
334 const Real cosineOfLongitudeOfAscendingNode = std::cos(longitudeOfAscendingNode);
+
335 const Real sineOfLongitudeOfAscendingNode = std::sin(longitudeOfAscendingNode);
+
336 const Real cosineOfTrueAnomaly = std::cos(trueAnomaly);
+
337 const Real sineOfTrueAnomaly = std::sin(trueAnomaly);
+
338
+
339 // Compute semi-latus rectum in the case the orbit is not a parabola.
+
340 Real semiLatusRectum = 0.0;
+
341 if (std::fabs(eccentricity - 1.0) > tolerance)
+
342 {
+
343 semiLatusRectum = semiMajorAxis * (1.0 - eccentricity * eccentricity);
+
344 }
+
345
+
346 // Else set the semi-latus rectum as the first element in the vector of Keplerian elements.
+
347 else
+
348 {
+
349 semiLatusRectum = keplerianElements[0];
+
350 }
+
351
+
352 // Compute the magnitude of the orbital radius, measured from the focal point.
+
353 const Real radiusMagnitude = semiLatusRectum / (1.0 + eccentricity * cosineOfTrueAnomaly);
+
354
+
355 // Define position and velocity in the perifocal coordinate system.
+
356 const Real xPositionPerifocal = radiusMagnitude * cosineOfTrueAnomaly;
+
357 const Real yPositionPerifocal = radiusMagnitude * sineOfTrueAnomaly;
+
358 const Real xVelocityPerifocal
+
359 = -std::sqrt(gravitationalParameter / semiLatusRectum) * sineOfTrueAnomaly;
+
360 const Real yVelocityPerifocal
+
361 = std::sqrt(gravitationalParameter / semiLatusRectum)
+
362 * (eccentricity + cosineOfTrueAnomaly);
+
363
+
364 // Compute scalar components of rotation matrix to rotate from periforcal to Earth-Centered
+
365 // Inertial (ECI) frame.
+
366 const Real rotationMatrixComponent11
+
367 = (cosineOfLongitudeOfAscendingNode * cosineOfArgumentOfPeriapsis
+
368 -sineOfLongitudeOfAscendingNode * sineOfArgumentOfPeriapsis * cosineOfInclination);
+
369 const Real rotationMatrixComponent12
+
370 = (-cosineOfLongitudeOfAscendingNode * sineOfArgumentOfPeriapsis
+
371 -sineOfLongitudeOfAscendingNode * cosineOfArgumentOfPeriapsis * cosineOfInclination);
+
372
+
373 const Real rotationMatrixComponent21
+
374 = (sineOfLongitudeOfAscendingNode * cosineOfArgumentOfPeriapsis
+
375 + cosineOfLongitudeOfAscendingNode * sineOfArgumentOfPeriapsis * cosineOfInclination);
+
376 const Real rotationMatrixComponent22
+
377 = (-sineOfLongitudeOfAscendingNode * sineOfArgumentOfPeriapsis
+
378 + cosineOfLongitudeOfAscendingNode * cosineOfArgumentOfPeriapsis * cosineOfInclination);
+
379
+
380 const Real rotationMatrixComponent31 = (sineOfArgumentOfPeriapsis * sineOfInclination);
+
381 const Real rotationMatrixComponent32 = (cosineOfArgumentOfPeriapsis * sineOfInclination);
+
382
+
383 // Compute Cartesian position and velocities.
+
384 cartesianElements[xPositionIndex] = rotationMatrixComponent11 * xPositionPerifocal
+
385 + rotationMatrixComponent12 * yPositionPerifocal;
+
386
+
387 cartesianElements[yPositionIndex] = rotationMatrixComponent21 * xPositionPerifocal
+
388 + rotationMatrixComponent22 * yPositionPerifocal;
+
389
+
390 cartesianElements[zPositionIndex] = rotationMatrixComponent31 * xPositionPerifocal
+
391 + rotationMatrixComponent32 * yPositionPerifocal;
392
-
393 // Compute scalar components of rotation matrix to rotate from periforcal to Earth-Centered
-
394 // Inertial (ECI) frame.
-
395 const Real rotationMatrixComponent11
-
396 = (cosineOfLongitudeOfAscendingNode * cosineOfArgumentOfPeriapsis
-
397 -sineOfLongitudeOfAscendingNode * sineOfArgumentOfPeriapsis * cosineOfInclination);
-
398 const Real rotationMatrixComponent12
-
399 = (-cosineOfLongitudeOfAscendingNode * sineOfArgumentOfPeriapsis
-
400 -sineOfLongitudeOfAscendingNode * cosineOfArgumentOfPeriapsis * cosineOfInclination);
+
393 cartesianElements[xVelocityIndex] = rotationMatrixComponent11 * xVelocityPerifocal
+
394 + rotationMatrixComponent12 * yVelocityPerifocal;
+
395
+
396 cartesianElements[yVelocityIndex] = rotationMatrixComponent21 * xVelocityPerifocal
+
397 + rotationMatrixComponent22 * yVelocityPerifocal;
+
398
+
399 cartesianElements[zVelocityIndex] = rotationMatrixComponent31 * xVelocityPerifocal
+
400 + rotationMatrixComponent32 * yVelocityPerifocal;
401
-
402 const Real rotationMatrixComponent21
-
403 = (sineOfLongitudeOfAscendingNode * cosineOfArgumentOfPeriapsis
-
404 + cosineOfLongitudeOfAscendingNode * sineOfArgumentOfPeriapsis * cosineOfInclination);
-
405 const Real rotationMatrixComponent22
-
406 = (-sineOfLongitudeOfAscendingNode * sineOfArgumentOfPeriapsis
-
407 + cosineOfLongitudeOfAscendingNode * cosineOfArgumentOfPeriapsis * cosineOfInclination);
-
408
-
409 const Real rotationMatrixComponent31 = (sineOfArgumentOfPeriapsis * sineOfInclination);
-
410 const Real rotationMatrixComponent32 = (cosineOfArgumentOfPeriapsis * sineOfInclination);
-
411
-
412 // Compute Cartesian position and velocities.
-
413 cartesianElements[xPositionIndex] = rotationMatrixComponent11 * xPositionPerifocal
-
414 + rotationMatrixComponent12 * yPositionPerifocal;
-
415
-
416 cartesianElements[yPositionIndex] = rotationMatrixComponent21 * xPositionPerifocal
-
417 + rotationMatrixComponent22 * yPositionPerifocal;
-
418
-
419 cartesianElements[zPositionIndex] = rotationMatrixComponent31 * xPositionPerifocal
-
420 + rotationMatrixComponent32 * yPositionPerifocal;
-
421
-
422 cartesianElements[xVelocityIndex] = rotationMatrixComponent11 * xVelocityPerifocal
-
423 + rotationMatrixComponent12 * yVelocityPerifocal;
+
402 return cartesianElements;
+
403}
+
404
+
406
+
419template <typename Real>
+ +
421 const Real eccentricity)
+
422{
+
423 assert(eccentricity >= 0.0 && eccentricity < 1.0);
424
-
425 cartesianElements[yVelocityIndex] = rotationMatrixComponent21 * xVelocityPerifocal
-
426 + rotationMatrixComponent22 * yVelocityPerifocal;
-
427
-
428 cartesianElements[zVelocityIndex] = rotationMatrixComponent31 * xVelocityPerifocal
-
429 + rotationMatrixComponent32 * yVelocityPerifocal;
-
430
-
431 return cartesianElements;
-
432}
-
433
-
435
-
448template <typename Real>
- -
450 const Real eccentricity)
-
451{
-
452 assert(eccentricity >= 0.0 && eccentricity < 1.0);
-
453
-
454 // Compute sine and cosine of eccentric anomaly.
-
455 const Real sineOfEccentricAnomaly
-
456 = std::sqrt(1.0 - std::pow(eccentricity, 2.0))
-
457 * std::sin(trueAnomaly) / (1.0 + eccentricity * std::cos(trueAnomaly));
-
458 const Real cosineOfEccentricAnomaly
-
459 = (eccentricity + std::cos(trueAnomaly))
-
460 / (1.0 + eccentricity * std::cos(trueAnomaly));
-
461
-
462 // Return elliptical eccentric anomaly.
-
463 return std::atan2(sineOfEccentricAnomaly, cosineOfEccentricAnomaly);
-
464}
+
425 // Compute sine and cosine of eccentric anomaly.
+
426 const Real sineOfEccentricAnomaly
+
427 = std::sqrt(1.0 - std::pow(eccentricity, 2.0))
+
428 * std::sin(trueAnomaly) / (1.0 + eccentricity * std::cos(trueAnomaly));
+
429 const Real cosineOfEccentricAnomaly
+
430 = (eccentricity + std::cos(trueAnomaly))
+
431 / (1.0 + eccentricity * std::cos(trueAnomaly));
+
432
+
433 // Return elliptical eccentric anomaly.
+
434 return std::atan2(sineOfEccentricAnomaly, cosineOfEccentricAnomaly);
+
435}
+
436
+
438
+
452template <typename Real>
+ +
454 const Real eccentricity)
+
455{
+
456 assert(eccentricity > 1.0);
+
457
+
458 // Compute hyperbolic sine and hyperbolic cosine of hyperbolic eccentric anomaly.
+
459 const Real hyperbolicSineOfHyperbolicEccentricAnomaly
+
460 = std::sqrt(std::pow(eccentricity, 2.0) - 1.0)
+
461 * std::sin(trueAnomaly) / (1.0 + std::cos(trueAnomaly));
+
462
+
463 const Real hyperbolicCosineOfHyperbolicEccentricAnomaly
+
464 = (std::cos(trueAnomaly) + eccentricity) / (1.0 + std::cos(trueAnomaly));
465
-
467
-
481template <typename Real>
- -
483 const Real eccentricity)
-
484{
-
485 assert(eccentricity > 1.0);
-
486
-
487 // Compute hyperbolic sine and hyperbolic cosine of hyperbolic eccentric anomaly.
-
488 const Real hyperbolicSineOfHyperbolicEccentricAnomaly
-
489 = std::sqrt(std::pow(eccentricity, 2.0) - 1.0)
-
490 * std::sin(trueAnomaly) / (1.0 + std::cos(trueAnomaly));
-
491
-
492 const Real hyperbolicCosineOfHyperbolicEccentricAnomaly
-
493 = (std::cos(trueAnomaly) + eccentricity) / (1.0 + std::cos(trueAnomaly));
-
494
-
495 // Return hyperbolic eccentric anomaly.
-
496 // The inverse hyperbolic tangent is computed here manually, since the atanh() function is not
-
497 // available in older C++ compilers.
-
498 const Real angle
-
499 = hyperbolicSineOfHyperbolicEccentricAnomaly
-
500 / hyperbolicCosineOfHyperbolicEccentricAnomaly;
-
501 return 0.5 * (std::log(1.0 + angle) - std::log(1.0 - angle));
-
502}
-
503
-
505
-
523template <typename Real>
-
524Real convertTrueAnomalyToEccentricAnomaly(const Real trueAnomaly, const Real eccentricity)
-
525{
-
526 assert(eccentricity >= 0.0
-
527 && (std::fabs(eccentricity - 1.0) > std::numeric_limits<Real>::epsilon()));
-
528
-
529 Real eccentricAnomaly = 0.0;
-
530
-
531 // Check if orbit is elliptical and compute eccentric anomaly.
-
532 if (eccentricity >= 0.0 && eccentricity < 1.0)
-
533 {
-
534 eccentricAnomaly
-
535 = convertTrueAnomalyToEllipticalEccentricAnomaly(trueAnomaly, eccentricity);
-
536 }
-
537
-
538 // Check if orbit is hyperbolic and compute eccentric anomaly.
-
539 else if (eccentricity > 1.0)
-
540 {
-
541 eccentricAnomaly
-
542 = convertTrueAnomalyToHyperbolicEccentricAnomaly(trueAnomaly, eccentricity);
-
543 }
-
544
-
545 return eccentricAnomaly;
-
546}
-
547
-
549
-
562template <typename Real>
-
563Real convertEllipticalEccentricAnomalyToMeanAnomaly(const Real ellipticalEccentricAnomaly,
-
564 const Real eccentricity)
-
565{
-
566 assert(eccentricity >= 0.0 && eccentricity < 1.0);
+
466 // Return hyperbolic eccentric anomaly.
+
467 // The inverse hyperbolic tangent is computed here manually, since the atanh() function is not
+
468 // available in older C++ compilers.
+
469 const Real angle
+
470 = hyperbolicSineOfHyperbolicEccentricAnomaly
+
471 / hyperbolicCosineOfHyperbolicEccentricAnomaly;
+
472 return 0.5 * (std::log(1.0 + angle) - std::log(1.0 - angle));
+
473}
+
474
+
476
+
494template <typename Real>
+
495Real convertTrueAnomalyToEccentricAnomaly(const Real trueAnomaly, const Real eccentricity)
+
496{
+
497 assert(eccentricity >= 0.0
+
498 && (std::fabs(eccentricity - 1.0) > std::numeric_limits<Real>::epsilon()));
+
499
+
500 Real eccentricAnomaly = 0.0;
+
501
+
502 // Check if orbit is elliptical and compute eccentric anomaly.
+
503 if (eccentricity >= 0.0 && eccentricity < 1.0)
+
504 {
+
505 eccentricAnomaly
+
506 = convertTrueAnomalyToEllipticalEccentricAnomaly(trueAnomaly, eccentricity);
+
507 }
+
508
+
509 // Check if orbit is hyperbolic and compute eccentric anomaly.
+
510 else if (eccentricity > 1.0)
+
511 {
+
512 eccentricAnomaly
+
513 = convertTrueAnomalyToHyperbolicEccentricAnomaly(trueAnomaly, eccentricity);
+
514 }
+
515
+
516 return eccentricAnomaly;
+
517}
+
518
+
520
+
533template <typename Real>
+
534Real convertEllipticalEccentricAnomalyToMeanAnomaly(const Real ellipticalEccentricAnomaly,
+
535 const Real eccentricity)
+
536{
+
537 assert(eccentricity >= 0.0 && eccentricity < 1.0);
+
538
+
539 return ellipticalEccentricAnomaly - eccentricity * std::sin(ellipticalEccentricAnomaly);
+
540}
+
541
+
543
+
557template <typename Real>
+ +
559 const Real hyperbolicEccentricAnomaly, const Real eccentricity)
+
560{
+
561 assert(eccentricity > 1.0);
+
562
+
563 return eccentricity * std::sinh(hyperbolicEccentricAnomaly) - hyperbolicEccentricAnomaly;
+
564}
+
565
567
-
568 return ellipticalEccentricAnomaly - eccentricity * std::sin(ellipticalEccentricAnomaly);
-
569}
-
570
-
572
-
586template <typename Real>
- -
588 const Real hyperbolicEccentricAnomaly, const Real eccentricity)
-
589{
-
590 assert(eccentricity > 1.0);
+
585template <typename Real>
+ +
587 const Real eccentricAnomaly, const Real eccentricity)
+
588{
+
589 assert(eccentricity >= 0.0
+
590 && std::fabs(eccentricity - 1.0) > std::numeric_limits<Real>::epsilon());
591
-
592 return eccentricity * std::sinh(hyperbolicEccentricAnomaly) - hyperbolicEccentricAnomaly;
-
593}
-
594
-
596
-
614template <typename Real>
- -
616 const Real eccentricAnomaly, const Real eccentricity)
-
617{
-
618 assert(eccentricity >= 0.0
-
619 && std::fabs(eccentricity - 1.0) > std::numeric_limits<Real>::epsilon());
-
620
-
621 Real meanAnomaly = 0.0;
-
622
-
623 // Check if orbit is elliptical and compute mean anomaly.
-
624 if (eccentricity >= 0.0 && eccentricity < 1.0)
-
625 {
-
626 meanAnomaly
-
627 = convertEllipticalEccentricAnomalyToMeanAnomaly(eccentricAnomaly, eccentricity);
-
628 }
-
629
-
630 // Check if orbit is hyperbolic and compute mean anomaly.
-
631 else if (eccentricity > 1.0)
-
632 {
-
633 meanAnomaly
-
634 = convertHyperbolicEccentricAnomalyToMeanAnomaly(eccentricAnomaly, eccentricity);
-
635 }
-
636
-
637 return meanAnomaly;
-
638}
-
639
-
641
-
650template <typename Real>
-
651Real convertEllipticalEccentricAnomalyToTrueAnomaly(const Real ellipticEccentricAnomaly,
-
652 const Real eccentricity)
-
653{
-
654 const Real sineOfTrueAnomaly = std::sqrt(1.0 - eccentricity * eccentricity)
-
655 * std::sin(ellipticEccentricAnomaly)
-
656 / (1.0 - eccentricity * std::cos(ellipticEccentricAnomaly));
-
657
-
658 const Real cosineOfTrueAnomaly
-
659 = (std::cos(ellipticEccentricAnomaly) - eccentricity)
-
660 / (1.0 - eccentricity * std::cos(ellipticEccentricAnomaly));
+
592 Real meanAnomaly = 0.0;
+
593
+
594 // Check if orbit is elliptical and compute mean anomaly.
+
595 if (eccentricity >= 0.0 && eccentricity < 1.0)
+
596 {
+
597 meanAnomaly
+
598 = convertEllipticalEccentricAnomalyToMeanAnomaly(eccentricAnomaly, eccentricity);
+
599 }
+
600
+
601 // Check if orbit is hyperbolic and compute mean anomaly.
+
602 else if (eccentricity > 1.0)
+
603 {
+
604 meanAnomaly
+
605 = convertHyperbolicEccentricAnomalyToMeanAnomaly(eccentricAnomaly, eccentricity);
+
606 }
+
607
+
608 return meanAnomaly;
+
609}
+
610
+
612
+
621template <typename Real>
+
622Real convertEllipticalEccentricAnomalyToTrueAnomaly(const Real ellipticEccentricAnomaly,
+
623 const Real eccentricity)
+
624{
+
625 const Real sineOfTrueAnomaly = std::sqrt(1.0 - eccentricity * eccentricity)
+
626 * std::sin(ellipticEccentricAnomaly)
+
627 / (1.0 - eccentricity * std::cos(ellipticEccentricAnomaly));
+
628
+
629 const Real cosineOfTrueAnomaly
+
630 = (std::cos(ellipticEccentricAnomaly) - eccentricity)
+
631 / (1.0 - eccentricity * std::cos(ellipticEccentricAnomaly));
+
632
+
633 return std::atan2(sineOfTrueAnomaly, cosineOfTrueAnomaly);
+
634}
+
635
+
637
+
646template <typename Real>
+
647Real convertHyperbolicEccentricAnomalyToTrueAnomaly(const Real hyperbolicEccentricAnomaly,
+
648 const Real eccentricity)
+
649{
+
650 const Real sineOfTrueAnomaly
+
651 = std::sqrt(eccentricity * eccentricity - 1.0)
+
652 * std::sinh(hyperbolicEccentricAnomaly)
+
653 / (eccentricity * std::cosh(hyperbolicEccentricAnomaly) - 1.0);
+
654
+
655 const Real cosineOfTrueAnomaly
+
656 = (eccentricity - std::cosh(hyperbolicEccentricAnomaly))
+
657 / (eccentricity * std::cosh(hyperbolicEccentricAnomaly) - 1.0);
+
658
+
659 return std::atan2(sineOfTrueAnomaly, cosineOfTrueAnomaly);
+
660}
661
-
662 return std::atan2(sineOfTrueAnomaly, cosineOfTrueAnomaly);
-
663}
-
664
-
666
-
675template <typename Real>
-
676Real convertHyperbolicEccentricAnomalyToTrueAnomaly(const Real hyperbolicEccentricAnomaly,
-
677 const Real eccentricity)
-
678{
-
679 const Real sineOfTrueAnomaly
-
680 = std::sqrt(eccentricity * eccentricity - 1.0)
-
681 * std::sinh(hyperbolicEccentricAnomaly)
-
682 / (eccentricity * std::cosh(hyperbolicEccentricAnomaly) - 1.0);
-
683
-
684 const Real cosineOfTrueAnomaly
-
685 = (eccentricity - std::cosh(hyperbolicEccentricAnomaly))
-
686 / (eccentricity * std::cosh(hyperbolicEccentricAnomaly) - 1.0);
-
687
-
688 return std::atan2(sineOfTrueAnomaly, cosineOfTrueAnomaly);
-
689}
-
690
-
692
-
710template <typename Real>
-
711Real convertEccentricAnomalyToTrueAnomaly(const Real eccentricAnomaly, const Real eccentricity)
-
712{
-
713 assert(eccentricity >= 0.0
-
714 && std::fabs(eccentricity - 1.0) > std::numeric_limits<Real>::epsilon());
-
715
-
716 Real trueAnomaly = 0.0;
-
717
-
718 // Check if orbit is elliptical and compute true anomaly.
-
719 if (eccentricity < 1.0)
-
720 {
-
721 trueAnomaly = convertEllipticalEccentricAnomalyToTrueAnomaly(eccentricAnomaly,
-
722 eccentricity);
-
723 }
-
724
-
725 else if (eccentricity > 1.0)
-
726 {
-
727 trueAnomaly = convertHyperbolicEccentricAnomalyToTrueAnomaly(eccentricAnomaly,
-
728 eccentricity);
-
729 }
-
730
-
731 return trueAnomaly;
-
732}
-
733
-
735
-
756template <typename Real>
-
757Real computeEllipticalKeplerFunction(const Real eccentricAnomaly,
-
758 const Real eccentricity,
-
759 const Real meanAnomaly)
-
760{
-
761 return eccentricAnomaly - eccentricity * std::sin(eccentricAnomaly) - meanAnomaly;
-
762}
-
763
-
765
-
782template <typename Real>
-
783Real computeFirstDerivativeEllipticalKeplerFunction(const Real eccentricAnomaly,
-
784 const Real eccentricity)
-
785{
-
786 return 1.0 - eccentricity * std::cos(eccentricAnomaly);
-
787}
-
788
-
790
-
812template <typename Real, typename Integer>
- -
814 const Real eccentricity,
-
815 const Real meanAnomaly,
-
816 const Real rootFindingTolerance = 1.0e-3 * std::numeric_limits<Real>::epsilon(),
-
817 const Integer maximumIterations = 100)
-
818{
-
819 assert(eccentricity >= 0.0 && eccentricity < (1.0 - 1.0e-11));
-
820
-
821 const Real pi = 3.14159265358979323846;
-
822
-
823 // Set mean anomaly to domain between 0 and 2pi.
-
824 Real meanAnomalyShifted = std::fmod(meanAnomaly, 2.0 * pi);
-
825 if (meanAnomalyShifted <0.0)
-
826 {
-
827 meanAnomalyShifted += 2.0 * pi;
-
828 }
-
829
-
830 Real eccentricAnomaly = std::numeric_limits<Real>::quiet_NaN();
-
831
-
832 // Set the initial guess for the eccentric anomaly.
-
833 // !!!!!!!!!!!!! IMPORTANT !!!!!!!!!!!!!
-
834 // If this scheme is changed, please run a very extensive test suite. The Newton-Raphson
-
835 // root-finding scheme tends to be chaotic for very specific combinations of mean anomaly
-
836 // and eccentricity. Various random tests of 100,000,000 samples were done to verify the
-
837 // default scheme for the initial guess. (Musegaas, 2012).
-
838 Real initialGuess = 0.0;
-
839 if (meanAnomalyShifted > pi)
-
840 {
-
841 initialGuess = meanAnomalyShifted - eccentricity;
-
842 }
-
843 else
-
844 {
-
845 initialGuess = meanAnomalyShifted + eccentricity;
-
846 }
-
847
-
848 // Execute Newton-Raphson root-finding algorithm.
-
849 eccentricAnomaly = initialGuess;
-
850 for (int i = 0; i < maximumIterations + 1; ++i)
-
851 {
-
852 if (i == maximumIterations)
-
853 {
-
854 throw std::runtime_error(
-
855 "ERROR: Maximum iterations for Newton-Raphson root-finding exceeded!");
-
856 }
-
857
-
858 const Real nextEccentricAnomaly
-
859 = eccentricAnomaly
-
860 - computeEllipticalKeplerFunction(eccentricAnomaly,
-
861 eccentricity,
-
862 meanAnomalyShifted)
- -
864 eccentricity);
-
865
-
866 const Real eccentricAnomalyDifference = eccentricAnomaly - nextEccentricAnomaly;
-
867
-
868 eccentricAnomaly = nextEccentricAnomaly;
-
869
-
870 if (eccentricAnomalyDifference < rootFindingTolerance)
-
871 {
-
872 break;
-
873 }
-
874 }
-
875
-
876 // Return eccentric anomaly.
-
877 return eccentricAnomaly;
-
878}
-
879
-
880} // namespace astro
-
881
+
663
+
681template <typename Real>
+
682Real convertEccentricAnomalyToTrueAnomaly(const Real eccentricAnomaly, const Real eccentricity)
+
683{
+
684 assert(eccentricity >= 0.0
+
685 && std::fabs(eccentricity - 1.0) > std::numeric_limits<Real>::epsilon());
+
686
+
687 Real trueAnomaly = 0.0;
+
688
+
689 // Check if orbit is elliptical and compute true anomaly.
+
690 if (eccentricity < 1.0)
+
691 {
+
692 trueAnomaly = convertEllipticalEccentricAnomalyToTrueAnomaly(eccentricAnomaly,
+
693 eccentricity);
+
694 }
+
695
+
696 else if (eccentricity > 1.0)
+
697 {
+
698 trueAnomaly = convertHyperbolicEccentricAnomalyToTrueAnomaly(eccentricAnomaly,
+
699 eccentricity);
+
700 }
+
701
+
702 return trueAnomaly;
+
703}
+
704
+
706
+
727template <typename Real>
+
728Real computeEllipticalKeplerFunction(const Real eccentricAnomaly,
+
729 const Real eccentricity,
+
730 const Real meanAnomaly)
+
731{
+
732 return eccentricAnomaly - eccentricity * std::sin(eccentricAnomaly) - meanAnomaly;
+
733}
+
734
+
736
+
753template <typename Real>
+
754Real computeFirstDerivativeEllipticalKeplerFunction(const Real eccentricAnomaly,
+
755 const Real eccentricity)
+
756{
+
757 return 1.0 - eccentricity * std::cos(eccentricAnomaly);
+
758}
+
759
+
761
+
783template <typename Real, typename Integer>
+ +
785 const Real eccentricity,
+
786 const Real meanAnomaly,
+
787 const Real rootFindingTolerance = 1.0e-3 * std::numeric_limits<Real>::epsilon(),
+
788 const Integer maximumIterations = 100)
+
789{
+
790 assert(eccentricity >= 0.0 && eccentricity < (1.0 - 1.0e-11));
+
791
+
792 const Real pi = 3.14159265358979323846;
+
793
+
794 // Set mean anomaly to domain between 0 and 2pi.
+
795 Real meanAnomalyShifted = std::fmod(meanAnomaly, 2.0 * pi);
+
796 if (meanAnomalyShifted <0.0)
+
797 {
+
798 meanAnomalyShifted += 2.0 * pi;
+
799 }
+
800
+
801 Real eccentricAnomaly = std::numeric_limits<Real>::quiet_NaN();
+
802
+
803 // Set the initial guess for the eccentric anomaly.
+
804 // !!!!!!!!!!!!! IMPORTANT !!!!!!!!!!!!!
+
805 // If this scheme is changed, please run a very extensive test suite. The Newton-Raphson
+
806 // root-finding scheme tends to be chaotic for very specific combinations of mean anomaly
+
807 // and eccentricity. Various random tests of 100,000,000 samples were done to verify the
+
808 // default scheme for the initial guess. (Musegaas, 2012).
+
809 Real initialGuess = 0.0;
+
810 if (meanAnomalyShifted > pi)
+
811 {
+
812 initialGuess = meanAnomalyShifted - eccentricity;
+
813 }
+
814 else
+
815 {
+
816 initialGuess = meanAnomalyShifted + eccentricity;
+
817 }
+
818
+
819 // Execute Newton-Raphson root-finding algorithm.
+
820 eccentricAnomaly = initialGuess;
+
821 for (int i = 0; i < maximumIterations + 1; ++i)
+
822 {
+
823 if (i == maximumIterations)
+
824 {
+
825 throw std::runtime_error(
+
826 "ERROR: Maximum iterations for Newton-Raphson root-finding exceeded!");
+
827 }
+
828
+
829 const Real nextEccentricAnomaly
+
830 = eccentricAnomaly
+
831 - computeEllipticalKeplerFunction(eccentricAnomaly,
+
832 eccentricity,
+
833 meanAnomalyShifted)
+ +
835 eccentricity);
+
836
+
837 const Real eccentricAnomalyDifference = eccentricAnomaly - nextEccentricAnomaly;
+
838
+
839 eccentricAnomaly = nextEccentricAnomaly;
+
840
+
841 if (eccentricAnomalyDifference < rootFindingTolerance)
+
842 {
+
843 break;
+
844 }
+
845 }
+
846
+
847 // Return eccentric anomaly.
+
848 return eccentricAnomaly;
+
849}
+
850
+
851} // namespace astro
+
852
-
Real convertEccentricAnomalyToTrueAnomaly(const Real eccentricAnomaly, const Real eccentricity)
Convert eccentric anomaly to true anomaly.
-
Vector6 convertKeplerianToCartesianElements(const Vector6 &keplerianElements, const Real gravitationalParameter, const Real tolerance=10.0 *std::numeric_limits< Real >::epsilon())
Convert Keplerian elements to Cartesian elements.
-
Real convertHyperbolicEccentricAnomalyToMeanAnomaly(const Real hyperbolicEccentricAnomaly, const Real eccentricity)
Convert hyperbolic eccentric anomaly to mean anomaly.
-
Real convertEllipticalEccentricAnomalyToTrueAnomaly(const Real ellipticEccentricAnomaly, const Real eccentricity)
Convert elliptical eccentric anomaly to true anomaly.
-
Real convertTrueAnomalyToEccentricAnomaly(const Real trueAnomaly, const Real eccentricity)
Convert true anomaly to eccentric anomaly.
-
Real convertHyperbolicEccentricAnomalyToTrueAnomaly(const Real hyperbolicEccentricAnomaly, const Real eccentricity)
Convert hyperbolic eccentric anomaly to true anomaly.
-
Real convertTrueAnomalyToEllipticalEccentricAnomaly(const Real trueAnomaly, const Real eccentricity)
Convert true anomaly to elliptical eccentric anomaly.
-
Real computeFirstDerivativeEllipticalKeplerFunction(const Real eccentricAnomaly, const Real eccentricity)
Compute 1st-derivative of Kepler function for elliptical orbits.
-
Real convertTrueAnomalyToHyperbolicEccentricAnomaly(const Real trueAnomaly, const Real eccentricity)
Convert true anomaly to hyperbolic eccentric anomaly.
-
Real convertEccentricAnomalyToMeanAnomaly(const Real eccentricAnomaly, const Real eccentricity)
Convert eccentric anomaly to mean anomaly.
-
Real computeEllipticalKeplerFunction(const Real eccentricAnomaly, const Real eccentricity, const Real meanAnomaly)
Compute Kepler function for elliptical orbits.
-
Vector6 convertCartesianToKeplerianElements(const Vector6 &cartesianElements, const Real gravitationalParameter, const Real tolerance=10.0 *std::numeric_limits< Real >::epsilon())
Convert Cartesian elements to Keplerian elements.
-
Real convertEllipticalEccentricAnomalyToMeanAnomaly(const Real ellipticalEccentricAnomaly, const Real eccentricity)
Convert elliptical eccentric anomaly to mean anomaly.
+
Real convertEccentricAnomalyToTrueAnomaly(const Real eccentricAnomaly, const Real eccentricity)
Convert eccentric anomaly to true anomaly.
+
Vector6 convertKeplerianToCartesianElements(const Vector6 &keplerianElements, const Real gravitationalParameter, const Real tolerance=10.0 *std::numeric_limits< Real >::epsilon())
Convert Keplerian elements to Cartesian elements.
+
Real convertHyperbolicEccentricAnomalyToMeanAnomaly(const Real hyperbolicEccentricAnomaly, const Real eccentricity)
Convert hyperbolic eccentric anomaly to mean anomaly.
+
Real convertEllipticalEccentricAnomalyToTrueAnomaly(const Real ellipticEccentricAnomaly, const Real eccentricity)
Convert elliptical eccentric anomaly to true anomaly.
+
Real convertTrueAnomalyToEccentricAnomaly(const Real trueAnomaly, const Real eccentricity)
Convert true anomaly to eccentric anomaly.
+
Real convertHyperbolicEccentricAnomalyToTrueAnomaly(const Real hyperbolicEccentricAnomaly, const Real eccentricity)
Convert hyperbolic eccentric anomaly to true anomaly.
+
Real convertTrueAnomalyToEllipticalEccentricAnomaly(const Real trueAnomaly, const Real eccentricity)
Convert true anomaly to elliptical eccentric anomaly.
+
Real computeFirstDerivativeEllipticalKeplerFunction(const Real eccentricAnomaly, const Real eccentricity)
Compute 1st-derivative of Kepler function for elliptical orbits.
+
Real convertTrueAnomalyToHyperbolicEccentricAnomaly(const Real trueAnomaly, const Real eccentricity)
Convert true anomaly to hyperbolic eccentric anomaly.
+
Real convertEccentricAnomalyToMeanAnomaly(const Real eccentricAnomaly, const Real eccentricity)
Convert eccentric anomaly to mean anomaly.
+
Real computeEllipticalKeplerFunction(const Real eccentricAnomaly, const Real eccentricity, const Real meanAnomaly)
Compute Kepler function for elliptical orbits.
+
Vector6 convertCartesianToKeplerianElements(const Vector6 &cartesianElements, const Real gravitationalParameter, const Real tolerance=10.0 *std::numeric_limits< Real >::epsilon())
Convert Cartesian elements to Keplerian elements.
+
Real convertEllipticalEccentricAnomalyToMeanAnomaly(const Real ellipticalEccentricAnomaly, const Real eccentricity)
Convert elliptical eccentric anomaly to mean anomaly.
+ +
@ longitudeOfAscendingNodeIndex
-
Real convertEllipticalMeanAnomalyToEccentricAnomaly(const Real eccentricity, const Real meanAnomaly, const Real rootFindingTolerance=1.0e-3 *std::numeric_limits< Real >::epsilon(), const Integer maximumIterations=100)
Convert elliptical mean anomaly to eccentric anomaly.
+
Real convertEllipticalMeanAnomalyToEccentricAnomaly(const Real eccentricity, const Real meanAnomaly, const Real rootFindingTolerance=1.0e-3 *std::numeric_limits< Real >::epsilon(), const Integer maximumIterations=100)
Convert elliptical mean anomaly to eccentric anomaly.