# A Note on the Use of Generalized Sundman Anomalies in the Numerical Integration of the Elliptical Orbital Motion

Authors: José Antonio López Ortí; Francisco José Marco Castillo; María José Martínez Usó
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

The orbital motion around a central body is an interesting problem
that involves the theory of artificial satellites and the planetary theories
in the solar system. Nevertheless some difficult situations appear while
studying this apparently simple problem, depending on each particular
case. The real problem consists of searching the perturbed solution from
a basic two-body motion problem. In addition, the perturbed problem
must be solved using a numerical method and its efficiency depends on
the selected coordinate system and the corresponding time. In fact, local
and global errors are not necessarily homogeneously distributed over the
orbit. In other words, there is a strong relationship between the spatial
distribution of the selected points and the temporal independent variable.
This is particularly dramatic in specially difficult cases.
This issue leads us to consider different anomalies as temporal variables, searching for both precision and efficiency. Therefore, we are interested in the study of techniques to integrate the orbital motion equations
using different anomalies as temporal variables which are functions of one
or more parameters. The final aim of this paper is the minimization of the integration errors using an appropriate choice of the parameter depending on the eccentricity
value in the family of the generalized Sundman anomalies.


---

## Body

## 1. Introduction

One of the most important problems in spatial mechanics is the study of the motion of an artificial satellite around the Earth. The relative motion of the satellite with respect to a coordinate system fixed on the Earth can be modelized through the second order differential equation:(1)
d
2
r
→
d
t
2
=
-
μ
r
→
r
3
-
∇
→
U
+
F
→
,
where r
→ is the radius vector of the satellite, G is the gravitational constant, m is the Earth mass, μ
=
G
m is the spaceflight constant, U is the potential induced by the conservative perturbation forces, from the Sun, Moon, and other planets gravitational forces and the perturbations due to the nonsphericity of the Earth, and F
→ represents the nonconservative disturbing forces as the atmospheric friction, the radiation pressure, the solar wind, and so forth.There are two ways to solve (1): first, a direct integration of the equations with the initial conditions r
→
(
0
)
=
r
→
0 and v
→
(
0
)
=
v
→
0, where v
→
=
r
→
˙.The second way involves the solution of the two-body problem:(2)
d
2
r
→
d
t
2
=
-
μ
r
→
r
3
,
r
→
(
0
)
=
r
→
0
,
v
→
(
0
)
=
v
→
0
obtaining the solution of (2) using the Keplerian elements and, then, getting the solution of (1) using the variation of constants method.The two-body problem is a classical integrable problem whose solution can be provided using the orbital elementsσ
→. Let (
a
,
e
,
i
,
Ω
,
ω
,
M
) be the third set of elements of Brouwer and Clemence [1] where a is the major semiaxis, e the eccentricity, i the inclination of the orbit with respect to the O
X
Y plane, Ω the argument of the ascending node, ω the argument of the perigee, and M the mean anomaly at an epoch t. In the elliptic case, let n
=
μ
/
a
3 be the mean motion, with μ
=
G
m. The mean anomaly is connected to the mean motion through the formula M
=
n
(
t
-
T
0
), where T
0 is the epoch of the perigee passage. The quantities a, e, i, Ω, ω are constants in the case of unperturbed motion.To solve exactly the two-body problem in a closed form it is convenient to use the eccentric anomalyE. This anomaly is related to the mean anomaly through the Kepler equation
(3)
E
-
e
sin
E
=
M
.
The coordinates of the secondary in the orbital coordinate system r
→
orb
=
(
ξ
,
η
,
ζ
)
t are
(4)
ξ
=
a
(
cos
⁡
E
-
e
)
,
η
=
a
1
-
e
2
sin
E
,
ζ
=
0
,
where O is the center of mass of the Earth, the axis O
ξ running to the perigee, the axis O
ζ pointing to the angular momentum, and the axis O
η making a direct trihedron. The radius vector of the satellite is given by
(5)
r
=
a
(
1
-
e
cos
⁡
E
)
.
The spatial coordinates r
→
=
(
x
,
y
,
x
)
t are related to the orbital coordinates by means of
(6)
r
→
=
A
r
→
orb
,
where A
=
R
1
(
-
Ω
)
R
3
(
-
i
)
R
1
(
-
ω
). R
i defines a rotation around the i-axis. The spatial speed of the satellite is given by
(7)
ξ
˙
=
a
n
sin
E
1
-
e
cos
⁡
E
,
η
˙
=
n
a
1
-
e
2
cos
⁡
E
1
-
e
cos
⁡
E
,
ζ
˙
=
0
,
and the orbital velocity
(8)
r
→
˙
=
A
r
→
˙
orb
.
To study the perturbed motion through analytical methods it is adequate to use the motion equation in Gauss planetary equations form [2]. In order to integrate the problem by means of numerical methods it is more convenient to use the system obtained from (1)
(9)
x
˙
=
v
x
,
v
˙
x
=
-
μ
x
r
3
-
∂
U
∂
x
+
F
x
,
y
˙
=
v
y
,
v
˙
y
=
-
μ
y
r
3
-
∂
U
∂
y
+
F
y
,
with the initial conditions r
→
(
0
)
=
r
→
0, r
→
˙
(
0
)
=
r
→
˙
0.Currently, different numerical integration methods are used to study the equations of motion. Our interest lies in being able to integrate the not disturbed problem, whose exact solution (in a formal and closed form, including series, integrals, and implicit equations) is known, by means of a method that allows reducing the errors using the maximum norm, both for the position and for the velocity. With this aim, we know that a bad distribution of the points in the periapsis and the apoapsis leads inexorably to truncation errors of such magnitude that make impossible the task proposed regarding the reduction of errors.To solve this problem Velez and Hilinski [3] suggested three pathways: the use of a very small fixed step-size; the use of a variable step-size method; and the use of a fixed step-size method together with a appropriate change of the temporal variable in order to uniformize the truncation errors over the orbit.In terms of the latter, already in 1922, Sundman [4] introduced a fictitious time τ related to the time t by means of d
t
=
C
r
d
τ to regularize the equations of the three bodies problem. Some time after, and generalizing this idea, Janin and Bond [5] proposed a new family of time transformations depending on a parameter α, called generalized Sundman transformations, given by d
t
=
Q
(
r
,
α
)
d
τ
α, where Q
(
r
,
α
)
=
C
α
r
α. This family includes the most common anomalies for appropriate values of α and C
α, the mean anomaly M for α
=
0, C
α
=
1
/
n, the eccentric anomaly E for α
=
1, the true anomaly V for α
=
2, C
=
(
1
/
n
a
2
)
(
1
/
1
-
e
2
), and the Nacozy intermediate anomaly for α
=
3
/
2 and C
1
/
2
=
(
μ
)
-
1 [4].The interest in this issue can be seen in later papers, for instance, the ones by Ferrándiz et al. [6] d
t
=
Q
(
r
)
d
τ who suggested a more complex family of transformations, where Q
(
r
)
=
r
2
/
3
(
a
0
+
a
1
r
)
-
1
/
2. Furthermore, Brumberg [7] suggested the use of the regularized length of arc s
* as temporal variable that provides the Q
(
r
) function given by
(10)
Q
(
r
)
=
1
2
μ
r
1
-
r
/
2
a
.As previously stated, it is convenient to increase the number of points around the regions with major curvature, perigee, and apogee, in order to have position and velocity under control.In this sense, Brumberg transformations can be considered as a first step, because they sparse the integration points as a uniform spatial distribution over the ellipse. It is worthy to remark that it is not a time uniform distribution. This technique improves the integration errors because it implies a greater concentration of points around the perigee than if the mean anomaly is used, but the problems at the apogee remain.On the other hand, the family of Sundman generalized anomaly contains the mean anomaly and the true anomaly as particular cases: in the first case the points are concentrated around the apogee but the number of points is minimum in the perigee where the curvature and the orbital speed are maxima. In the true anomaly case, the maximum of points are concentrated around the perigee and the minimum around the apogee.Making use of the properties of the Sundman generalized anomalies, the main idea of our proposed method consists of combining the good properties of the use of different anomalies at the perigee and the apogee. As we will see later, the generalized Sundman family of anomalies are a good choice.After this brief introduction to the general problem and its background, we describe in Section2 the general Sundman anomaly properties. We also obtain the differential equations of motion using an arbitrary general Sundman anomaly as temporal variable.In Section3 we consider two numerical examples. First, an example of a simple two-body problem by means of the integration of the old HEOS I satellite. Secondly, we carry out the resolution of a disturbed problem: the two-fixed-center problem (the background theory is included in the Appendix). Both proposed examples have closed solution (including an integral formulation solution for the second example).Finally, we expose the main conclusions and remarks about this paper.


## 2. Algorithms to Use the Generalized Sundman Anomalies as Temporal Variables

Letd
t
=
Q
(
r
,
α
)
d
τ be a generalized Sundman transformation and let us define the generalized Sundman anomaly Ψ
α as the regularized value of τ
α given by
(11)
d
M
=
K
α
(
e
)
r
α
d
Ψ
α
,
where Ψ
α
=
M if M
=
k
π, for all k
∈
Z, and Ψ
α
(
-
M
)
=
-
Ψ
α
(
M
). The constant value of K
α
(
e
) has been computed by López Ortí et al. [8]
(12)
K
α
(
e
)
=
a
-
α
{
(
1
-
e
)
p
F
(
1
2
,
-
p
,
1
;
u
)
+
(
1
+
e
)
1
-
α
F
(
1
2
,
-
p
,
1
;
v
)
}
,
where F
(
a
,
b
,
c
;
z
) is the hypergeometric function, p
=
1
-
α, u
=
2
e
/
(
e
-
1
), and v
=
2
e
/
(
1
+
e
). For each α we have
(13)
d
d
t
=
n
d
d
M
=
n
d
d
Ψ
α
d
Ψ
α
d
M
=
n
K
α
(
e
)
r
-
α
d
d
Ψ
α
,
so,
(14)
d
x
d
Ψ
α
=
K
α
(
e
)
n
r
α
v
x
,
d
v
x
d
Ψ
α
=
-
K
α
(
e
)
n
r
α
[
μ
x
r
3

+

∂
V
∂
x

-

F
x
]
,
d
y
d
Ψ
α
=
K
α
(
e
)
n
r
α
v
y
,
d
v
y
d
Ψ
α
=
-
K
α
(
e
)
n
r
α
[
μ
y
r
3

+

∂
V
∂
y

-

F
y
]
.
In the most common problems of satellite motion, the disturbing forces are small, so the optimal step-size can be obtained from the unperturbed elliptic motion. To compute this solution we proceed by integration of (14) using a Gragg [9] integrator with polynomial extrapolation and a classic Runge-Kutta (4) method.The global error for an arbitrary value ofΨ
α can be computed from the value of the exact solution of the two bodies problem for the value of eccentric anomaly E
(
Ψ
α
) obtained through the derivative of the Kepler equation d
M
=
(
1
-
e
cos
⁡
E
)
d
E
=
(
r
/
a
)
d
E. Replacing in (11), we obtain the equation
(15)
Ψ
α
=
a
-
α
K
α
(
e
)
∫
0
E
(
Ψ
α
)
‍
(
1
-
e
cos
⁡
E
)
1
-
α
d
E
,
and solving this implicit integral equation we get the value E
(
Ψ
α
). Replacing this value in (4), (7), (6), and (8) we get the exact solution of the problem and then the global error in position and velocity of the numerical method.


## 3. Numerical Examples

To show the spatial distributions of the points over the orbit we represent, see Figure1, twenty points Ψ
α
(
e
)
=
{
j
π
/
10
;
j
=
1
,
…
,
20
} for an ellipse with a
=
1 and e
=
0.7. Values of α are α
=
0,1
,
1.5,2.Points distribution fore
=
0.7 and α
=
0.0, 0.5, 0.75, 1.0.
(a)
α
=
0.0
(b)
α
=
0.5
(c)
α
=
0.75
(d)
α
=
1.0To test the efficiency of the previous transformations, we use as an example a highly eccentric satellite: the old satellite HEOS I used by Brumberg [7].The elements of this satellite [10] are a
=
118363.47 Km, e
=
0.942572319, i
=
2
8
°
.
16096, Ω
=
18
5
°
.
07554, ω
=
27
0
°
.
07151, and M
0
=
0
°. The period of the satellite is 4.69 days and for the Earth the spaceflight constant is μ
=
3.986005
K
m
3
s
-
2.Table1 shows the global error in position |
r
→
-
r
→
0
| for one revolution, in units of 1.1
0
-
5 Km, and the error in velocity |
v
→
-
v
→
0
| in units of 1.1
0
-
8
Km
s
-
1. This table has been constructed using a Gragg-Bulirsch-Stoer method with polynomial extrapolation with 1000 uniform steps and 10 evaluations of the function for each step. In this table, we can see that the minimum in position and velocity errors is reached near the value α
=
1.9. In the same conditions, the Brumberg method produces the errors: |
r
→
-
r
→
0
|
=
5.467
×
1
0
-
5 Km and |
v
→
-
v
→
0
|
=
0.012
×
1
0
-
8
Km
s
-
1, for the regularized length of arc parameter s
*. Compare, for example, with the results for α
=
1.9 in Table 1.Table 1
Errors in position and velocity for HEOS I in1
0
-
5
Km and 1
0
-
8
Km
/
s, respectively.
α
|
r
→
-
r
→
0
|
|
v
→
-
v
→
0
|
α
|
r
→
-
r
→
0
|
|
v
→
-
v
→
0
|
0.0
953553.59323
770880.33913
1.6
0.00086
0.00078
0.1
258032.45040
208605.19468
1.7
0.00037
0.00037
0.2
66193.91357
53516.13695
1.8
0.00016
0.00020
0.3
16522.81621
13359.02237
1.9
0.00005
0.00011
0.4
4074.21292
3294.34188
2.0
0.00009
0.00003
0.5
1002.70466
810.86164
2.1
0.00032
0.00030
0.6
248.26083
200.79429
2.2
0.00112
0.00127
0.7
62.24464
50.35530
2.3
0.00320
0.00399
0.8
15.89750
12.86526
2.4
0.00961
0.01251
0.9
4.15863
3.36714
2.5
0.03029
0.03987
1.0
1.12001
0.90756
2.6
0.10106
0.13113
1.1
0.31248
0.25354
2.7
0.36030
0.45054
1.2
0.08991
0.07312
2.8
1.37345
1.62725
1.3
0.02677
0.02188
2.9
5.59493
6.21512
1.4
0.00880
0.00725
3.0
24.23229
25.17338
1.5
0.00282
0.00238
3.1
110.62982
107.96944Figures2, 3, and 4 show the local truncation errors for a fictitious body with the elements of HEOS satellite, except that the eccentricity is taken as e
=
0.95 for the values α
=
1.0,1.5
,
2.0. These errors has been obtained by solving for each value of Ψ
α the equation Ψ
α
=
Ψ
α
(
E
) (15). From the value of E we compute the exact value of the position and velocity vectors. These values are the initial conditions for the numerical integrator. To determine the exact local truncation errors we compare the values obtained in the next step of Ψ
α with the ones evaluated solving the equation Ψ
α
+
h
=
Ψ
α
(
E
). Table 2 shows the minimum of the errors for this fictitious satellite taking several values of eccentricity. This table has been built up using a classic Runge-Kutta integrator of fourth order with 1000 uniform steps. In this table we can see that the errors in position reach the minimum at a value of α depending on the eccentricity.Table 2
Optimalα for e  |
Δ
r
→
| in 1
0
-
5 Km  |
Δ
v
→
|  and in 1
0
-
8
Km
/
s.
e
α
|
Δ
r
→
|
|
Δ
v
→
|
e
α
|
Δ
r
→
|
|
Δ
v
→
|
0.00
1.539
3.74
e
-
07
2.92
e
-
12
0.50
1.671
1.88
e
-
07
6.57
e
-
10
0.05
1.561
3.69
e
-
07
5.71
e
-
11
0.55
1.681
1.60
e
-
07
7.13
e
-
10
0.10
1.578
3.56
e
-
07
1.16
e
-
10
0.60
1.692
1.41
e
-
07
7.60
e
-
10
0.15
1.593
3.44
e
-
07
1.79
e
-
10
0.65
1.704
1.32
e
-
07
7.86
e
-
10
0.20
1.606
3.31
e
-
07
2.43
e
-
10
0.70
1.718
1.06
e
-
07
7.77
e
-
10
0.25
1.618
3.15
e
-
07
3.11
e
-
10
0.75
1.735
9.81
e
-
08
6.97
e
-
10
0.30
1.630
2.98
e
-
07
3.83
e
-
10
0.80
1.757
1.35
e
-
07
4.65
e
-
10
0.35
1.640
2.66
e
-
07
4.50
e
-
10
0.85
1.790
2.31
e
-
07
1.10
e
-
10
0.40
1.650
2.50
e
-
07
5.18
e
-
10
0.90
1.847
5.79
e
-
07
1.66
e
-
09
0.45
1.661
2.25
e
-
07
5.93
e
-
10
0.95
1.917
2.70
e
-
06
9.83
e
-
09Local integration errors distributionα
=
1.0.
(a)
e
r
x
·
1
0
5
(b)
e
r
y
·
1
0
5
(c)
e
r
v
x
·
1
0
8
(d)
e
r
v
y
·
1
0
8Local integration errors distributionα
=
1.5.
(a)
e
r
x
·
1
0
5
(b)
e
r
y
·
1
0
5
(c)
e
r
v
x
·
1
0
8
(d)
e
r
v
y
·
1
0
8Local integration errors distributionα
=
2.0.
(a)
e
r
x
·
1
0
5
(b)
e
r
y
·
1
0
5
(c)
e
r
v
x
·
1
0
8
(d)
e
r
v
y
·
1
0
8A least square analytical approach to the optimal value ofα can be written as
(16)
α
=
3.38992
e
5
-
6.49697
e
4
+
4.78192
e
3
-
1.73234
e
2
+
0.5381
e
+
1.53836
.
Figure 5 shows the values of α where the error |
Δ
r
→
| reaches its minimum for each value of e
∈
[
0,0.95
] (dots) and its analytical approximation (continuous line).Figure 5
Optimal value ofα for each value of e
∈
[
0,0.95
].To test the robustness of the method we use as perturbed problem a particular case of the two-fix centers with the following characteristics.(1)
The first fix massM is placed at the origin of coordinates. The second fix mass M
′ is placed at (
0,2
c
).
(2)
The initial osculating elements of the satellite are the same elements used in the unperturbed problem.
(3)
The valuesM
′ and c have been chosen as M
′
=
1
0
-
2
M and c
=
10
a, where M is the Earth mass and a is the initial osculating semiaxis. Notice that the exact resolution of the two-fix centers problem is quite complex. See more details in the Appendix. To solve this problem in a time span [
0
,
T
] using Ψ
α as temporal variable, it is necessary to determinate the value of Ψ
α
(
T
), where e depends on t. For this purpose it is neccesary to solve the boundary problem:
(17)
d
x
d
Ψ
α
=
K
α
(
e
)
n
r
α
v
x
,
d
v
x
d
Ψ
α
=
-
K
α
(
e
)
n
r
α
[
μ
x
r
3
+
∂
V
∂
x
]
,
d
y
d
Ψ
α
=
K
α
(
e
)
n
r
α
v
y
,
d
v
y
d
Ψ
α
=
-
K
α
(
e
)
n
r
α
[
μ
y
r
3
+
∂
V
∂
y
]
,
d
t
d
Ψ
α
=
K
α
(
e
)
n
r
α
,
Ψ
α
∈
[
0
,
Ψ
α
(
T
)
]
,
with the boundary conditions given by r
→
(
0
)
=
r
→
0, v
→
(
0
)
=
v
→
0, and t
(
Ψ
α
(
T
)
)
=
T. The value of the disturbing potential V is given by
(18)
V
=
-
G
M
′
x
2
+
(
y
-
2
c
)
2
.To solve this problem a shooting method has been used in order to obtain theα value for which the number of steps to assure an accuracy of 1
0
-
4 Km in position is minimum. Table 3 shows this required number of steps in an integration over 100 initial periods of revolution using as integrator a classic Runge-Kutta of fourth order and a Runge-Kutta of eighth order [11]. Notice that the value is reached for α
(
e
).Table 3
Number of steps to get a precision of|
Δ
r
→
|
<
1
0
-
4 Km.
α
=
0
α
=
0.5
α
=
1.0
α
=
1.5
α
=
2.0
α
(
e
)
rk4
2287352
1562132
969798
662225
601888
595921
rk8
73994
41764
21428
13302
16091
12091To evaluate the global error we have compared the results for the position usingN and 1.1
N. In Table 3 we can see that the integration is improved using an appropriate choice of α in the generalized Sundman family of anomalies given by α
(
e
) (16).


## 4. Concluding Remarks

The generalized family of Sundman anomalies is adequate to be used in the numerical integration of the perturbed two-body problem. To study the optimal value ofα in order to increase the performance of the numerical methods, we have carried out some numerical experiments on the unperturbed two-body problem.The local and global errors depend on the value of theα parameter of the selected anomaly. The optimal value of α in order to minimize the global integration errors in a revolution increases with the eccentricity. This value is near 1.5 for the lowest eccentricities and 1.9 for the highest.To study the performance of the method for high eccentricities (e
=
0.95) a numerical integration using 1000 steps for uniformly distributed Ψ
α has been carried out, taking the initial conditions of the two-body problem as initial conditions and evaluating the difference between the computed position and speed with the ones corresponding to the two-body problem. The minimum of the local error is reached near α
=
1.5 but the global error is better for the value 1.9.The use of an appropriate value ofα in the generalized Sundman family allows obtaining a more uniform spatial and temporal distribution for the points over the orbit than if the mean anomaly is used. The value of α when the minimum of global error is reached depends strongly on the value of eccentricity and both of them increasing at the same time.The robustness of the method has been tested using as perturbed problem a particular case of the two-fix center problem, which is a well known integrable problem whose solution has been used to evaluate the necessary steps to reach a1
0
-
4 Km precision in position. The best results are obtained when we use for each step the optimal α value for the osculating eccentricity.



---

*Source: 691926-2014-09-11.xml*