# A Study about the Integration of the Elliptical Orbital Motion Based on a Special One-Parametric Family of Anomalies

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

---

## Abstract

This paper aimed to address the study of a new family of anomalies, called natural anomalies, defined as a one-parameter convex linear combination of the true and secondary anomalies, measured from the primary and the secondary focus of the ellipse, and its use in the study of analytical and numerical solutions of perturbed two-body problem.  We take two approaches: first, the study of the analytical development of the basic quantities of the two-body problem to be used in the analytical theories of the planetary motion and second, the study of the minimization of the errors in the numerical integration by an appropriate choice of parameters in our family for each value of the eccentricity. The use of an appropriate value of the parameter can improve the length of the developments in the analytical theories and reduce the errors in the case of the numerical integration.


---

## Body

## 1. Introduction

One of the main problems in celestial mechanics is the study of the motion in the solar system. In this regard, it is interesting to study the planetary theories and the motion of an artificial satellite around the Earth.To construct a planetary theory, we can follow two main ways: one, the use of a numerical integrator, and two, the use of analytical methods to integrate the problem [1–4]. This paper is focused on planetary theories in the second sense. The analytical and semianalytical planetary theories involve the management of very long Poisson series developments [5] using, as temporal variables, the mean longitudes of the planets. In this sense, the study of a change of the temporal variable in order to obtain developments that are more compact than the ones obtained if the mean longitudes (or mean anomalies) are used [6–9] can be interesting. To construct an analytical or semianalytical planetary theory, it is necessary to obtain a set of developments for the two-body problem, and from that, we can evaluate the inverse of the distance between the two planets and so to develop the second member of the differential equations of motion. Obtaining precise developments of the inverse of the distance is one of most important problems to construct the analytical theories.To study the motion of an artificial satellite, it is more convenient to use a numerical integration method. An excellent overview on numerical integrators can be seen in [10] containing a theory of symplectic and symmetric integrator, including Runge-Kutta (composition, splitting,
…) and multistep and specially designed integrators, also their construction and practical merits are discussed. However, the history of modern numerical integrators begins long ago. Adams introduced a multistep method to study the perturbed motion of the couple Jupiter-Saturn, Kutta introduced at the beginning of the 20th century the well known family of Runge-Kutta integrators, Bulirsch and Stoer developed extrapolation methods that were used by the Institute of Applied Astronomy of St. Petersburg to obtain the minor planet ephemerides, and so forth.Numerical integrators are appropriate to study the planetary theories and the satellite motion. In the first case, it is very important for the errors of the methods to be bounded for a very large time span and for this purpose, symplectic integrators can be adequate. One of the most important problems in the study of a satellite motion is the case of very high eccentricities. The main problem with using the mean anomaly in this case is the non uniformity on the point distribution on the ellipse due to the second Kepler law; this, for a constant step size, implies a minor concentration of points in the perigee and so larger errors in this region. In this case, it can be convenient to use a technique known as analytical regularization of step size. This technique involves a change of the temporal variable in order to have,for the new variable, a point distribution containing a major density of points in the orbital regions, where the velocity is higher. In this context, Sundman [11] introduces a new temporal variable τ related to the time t through d
t
=
C
r
d
τ. Nacozy [12], Janin and Bond [13, 14], and Velez and Hilinski [15] extend this technique defining a new one-parameter family of transformations α called generalized Sundman transformations:
(1)
d
t
=
Q
(
r
,
α
)
d
τ
α
,
where Q
(
r
,
α
)
=
C
α
r
α is the called partition function. A more complicated transformation was introduced by Ferrándiz et al. [16] Q
(
r
)
=
r
2
/
3
(
a
0
+
a
1
r
)
-
1
/
2 and Brumberg [17] who proposed the use of the regularized length of arc:
(2)
Q
(
r
)
=
1
2
G
M
r
1
-
r
/
2
a
,
where G
M is the spaceflight constant of the Earth and a the major semiaxe of the osculating elliptic orbit of the satellite.In this paper, we follow this technique. In particular, we propose a new one-parametric family of anomalies, called natural anomalies, defined by a convex linear combinationΨ
α
=
α
f
+
(
1
-
α
)
f
′, where f is the true anomaly and f
′ the polar angle referred to the secondary focus called secondary true anomaly.The rest of this article is organized as follows.In Section2, we introduce the general background. This section contains the equations of the perturbed motion in two ways: first, to study the analytical planetary theories and second, to use appropriate numerical methods.In Section3, the properties of natural family of anomalies are also described. These properties include the connection between the quantities r and r
′ with f
′, the connection between the natural anomaly Ψ
α and the eccentric anomaly E, and the study of the partition function denoted by Q
α
(
r
) for this case.In Section4, the analytical developments according to the natural anomaly are studied. This study contains the developments of E, sin
E, cos
⁡
E, and 1
/
r with respect to the natural anomaly. So, an appropriate use of them is enough to obtain the development of the mean anomaly according to the natural anomaly, that is the Kepler equation. In this section, we show the length of the analytical development of the inverse of the distance for the couple Jupiter-Saturn using several values of the α.In Section5, we obtain the differential equations of motion using an arbitrary anomaly from the natural family. A set of numerical examples about the two-body problem are analyzed. We carry out a perturbed problem in order to analyze the robustness of the method.In Section6, the main conclusions and remarks of this paper are exposed.


## 2. Basic Equations

The analytical methods are based on the solution of the two-body problem (Sun planet) through a set of orbital elements, for example, the third set of Brouwer and Clemence [18] (
a
,
e
,
i
,
Ω
,
ω
,
M
), where M
=
M
0
+
n
(
t
-
t
0
), n is the mean motion, and t
0 is the initial epoch being (
a
,
e
,
i
,
Ω
,
ω
,
n
) constants in the unperturbed two-body problem. This solution can be used as a first approximation of the perturbed problem and we can use the Lagrange method of variation of constants to replace the first elements by the osculating ones given by the Lagrange planetary equations [19]:
(3)
d
a
d
t
=
2
n
a
∂
R
∂
σ
,
d
e
d
t
=
-
1
-
e
2
n
a
2
e
∂
R
∂
ω
+
1
-
e
2
n
a
2
e
∂
R
∂
σ
,
d
i
d
t
=
-
1
n
a
2
1
-
e
2
sin
i
∂
R
∂
Ω
+
ctg
i
n
a
2
1
-
e
2
∂
R
∂
ω
,
d
Ω
d
t
=
1
n
a
2
1
-
e
2
sin
i
∂
R
∂
i
,
d
ω
d
t
=
1
-
e
2
n
a
2
e
∂
R
∂
e
-
cos
⁡
i
n
a
2
1
-
e
2
sin
i
∂
R
∂
i
,
d
σ
d
t
=
-
2
n
a
∂
R
∂
a
-
1
-
e
2
n
a
2
e
∂
R
∂
e
,
where σ is a new variable defined by the equation:
(4)
M
=
σ
+
∫
t
0
t
n
d
t
and it coincides with M
0 in the case of the unperturbed motion. R is the disturbing potential R
=
∑
k
=
1
N
R
i due to the disturbing bodies i
=
1
,
…
,
N. It is defined as [19]
(5)
R
=
∑
k
=
1
N
G
m
k
[
(
1
Δ
k
)
-
x
·
x
k
+
y
·
y
k
+
z
·
z
k
r
k
3
]
,
where r
→
=
(
x
,
y
,
z
) and r
→
k
=
(
x
k
,
y
k
,
z
k
) are the heliocentric vector position of the secondary body and the kth disturbing body, respectively, Δ
k is the distance between the secondary body and the disturbing body, and m
k is the mass of the disturbing body.In order to integrate the Lagrange planetary equations through analytical methods, it is necessary to develop the second member of the Lagrange planetary equations as truncated Fourier series, which is a classical problem in celestial mechanics [4, 5, 18, 20, 21]. The analytical methods are appropriated to study planetary motion because the eccentricities of the planetary orbits are small. Despite this, analytical methods provide very long series solution and it is convenient to obtain more compact developments using as temporal variable an appropriate anomaly [6, 7].To obtain the developments with respect to an anomalyΨ
i, it is necessary to obtain for each planet i the developments of the coordinates and the inverse of the radius in Fourier series of Ψ
i. Then, the integration of the Lagrange planetary equations with respect to the Ψ
i anomalies requires to obtain solution of the corresponding Kepler equation M
i
=
M
i
(
Ψ
i
), for each planet.To use numerical integration methods it is more appropriate to consider the equation of motion in the form of the second Newton law. To study the efficiency of the numerical integration methods using an anomalyΨ as temporal variable, we select the problem of the motion of an artificial satellite around the Earth. The relative motion of the secondary with respect to the Earth is defined by the second order differential equations
(6)
d
2
r
→
d
t
2
=
-
G
M
r
→
r
3
-
∇
→
U
-
F
→
,
where r
→ is the radius vector of the satellite, U the potential from which the perturbative conservative forces are derived, and F
→ includes the nonconservative forces. To integrate system (6), it is necessary to know the initial values of the radius vector r
→
0 and velocity v
→
0.In order to uniformize the truncation errors when a numerical integrator is used, there are three main ways: the use of a very small stepsize, the use of a adaptive stepsize method, and the use of a change in the temporal variable to get an appropriate distribution of the points in the orbit so that the points are mostly concentrated in the regions where the speed and curvature are maxima. In this paper, we follow the third way, as previously stated.


## 3. Natural Anomalies

In this section, a new family of anomalies depending on a parameter is defined. Let us represent in Figure1 the elliptic orbit corresponding to the motion of the two-body problem. This ellipse is defined by its major semiaxis a
=
O
Q
¯ and its eccentricity e
=
c
/
a, 0
≤
e
<
1, where c is the focal semidistance c
=
F
F
′
¯
/
2, and the minor semiaxis b is defined as b
=
a
1
-
e
2. Let O be the center of the ellipse, let F be the primary focus, let F
′ be the secondary focus (also called equality point), let Q be the pericenter, and let P be the position of the secondary in the orbit. Let us define the coordinates (
ξ
,
η
) referred to the primary focus and let r and r
′ be the distance between the secondary P and the primary focus and the secondary focus, F and F
′, respectively. The angle f is called the true anomaly and for the angle f
′, we propose the name “secondary true anomaly”.Figure 1
Elliptic motion.The quantitiesr and r
′ satisfy
(7)
r
+
r
′
=
2
a
,
r
=
a
2
(
1
-
e
2
)
1
+
e
cos
⁡
f
,
r
′
=
a
2
(
1
-
e
2
)
1
-
e
cos
⁡
f
′
.
On the other hand, the coordinates of the secondary with respect to the primary, in the orbital plane, are
(8)
ξ
=
r
sin
f
=
r
′
sin
f
′
=
a
(
e
-
cos
⁡
E
)
,
η
=
r
cos
⁡
f
=
r
′
cos
⁡
f
′
-
2
a
e
=
a
1
-
e
2
sin
E
.The spatial orbital coordinatesr
→
=
(
x
,
y
,
x
)
t are related to the orbital coordinates (
ξ
,
η
,
0
)
t by means of
(9)
r
→
=
A
r
→
orb
,
where A
=
R
1
(
-
Ω
)
R
3
(
-
i
)
R
1
(
-
ω
). R
i defines a rotation around the i-axis.From (7) we have
(10)
d
r
+
d
r
′
=
0
,
d
r
=
-
a
2
(
1
-
e
2
)
sin
f
d
f
(
1
+
e
cos
⁡
f
)
2
,
d
f
′
=
a
2
(
1
-
e
2
)
sin
f
′
d
f
′
(
1
-
e
cos
⁡
f
′
)
2
and from (10), it is easy to obtain r
2
sin
f
d
f
=
r
′
2
sin
f
′
d
f
′. Taking into account (8), we obtain
(11)
r
d
f
=
r
′
d
f
′
.The radiir and r
′ are related to the eccentric anomaly E through
(12)
r
=
a
(
1
-
e
cos
⁡
E
)
,
r
′
=
a
(
1
+
e
cos
⁡
E
)
.
The anomalies f and f
′ are related to the eccentric anomaly E by
(13)
tan
f
2
=
1
+
e
1
-
e
tan
E
2
,
tan
f
′
2
=
1
-
e
1
+
e
tan
E
2
.
To relate the anomaly f
′ to the mean anomaly M we use the integral of areas d
f
=
(
a
2
1
-
e
2
/
r
2
)
d
M, and after replacing in (11), we get
(14)
d
f
′
=
a
2
1
-
e
2
r
r
′
d
M
=
1
-
e
2
1
-
e
2
cos
⁡
2
E
d
M
.
To compare f
′ and M up to second order in e we have
(15)
f
′
=
M
+
e
2
4
sin
2
M
+
O
(
e
3
)
and so, for small values of eccentricity, the anomaly f
′ is near M.Let us define the natural family of anomaliesΨ
α as
(16)
Ψ
α
=
α
f
+
(
1
-
α
)
f
′
,
0
≤
α
≤
1
this family includes, for α
=
0, the secondary true anomaly f
′ and, for α
=
1, the true anomaly f. The anomaly Ψ
α is related to the mean anomaly by
(17)
d
M
=
Q
α
(
r
)
d
Ψ
α
,
where the function Q
α
(
r
) is defined as
(18)
Q
α
(
r
)
=
1
a
2
1
-
e
2
r
[
α
r
+
1
-
α
2
a
-
r
]
-
1
.


## 4. Analytical Developments

To integrate the planetary Lagrange equations, it is necessary to develop the second member of the Lagrange planetary equations as Fourier series with respect to the selected anomalies for each couple of planets. For this task, it is necessary to obtain for each planet the developments according to the selected anomaly of the orbital coordinatesξ, and η, the inverse of the radius 1
/
r and the mean anomaly M.For this purpose, first of all, we obtain the development of an arbitrary anomalyΨ
(
α
) in the natural anomalies family. In the future, we will denote Ψ
(
α
) as Ψ, if it does not lead to confusion. The developments of the orbital coordinates ξ, η (8) are completely determined by the ones of sin
E and cos
⁡
E. To obtain these developments, it is necessary first to expand Ψ with respect to E [22].To relate the natural anomaly to the eccentric anomaly, it is convenient to obtain the development of the natural anomaly as Fourier series of eccentric anomaly. To this aim, it is convenient to use the classic formula [4]
(19)
tan
f
2
=
1
+
e
1
-
e
tan
E
2
=
1
+
q
1
-
q
tan
E
2
,
where q
=
tan
(
ϕ
/
2
), e
=
sin
ϕ and so
(20)
f
=
E
+
∑
k
=
1
∞
2
k
q
k
sin
k
E
.
The secondary true anomaly verifies
(21)
tan
f
′
2
=
1
-
e
1
+
e
tan
E
2
=
1
-
q
1
+
q
tan
E
2
and consequently
(22)
f
′
=
E
+
∑
k
=
1
∞
(
-
1
)
k
2
k
q
k
sin
k
E
.
Replacing (20) and (22) in (16) we obtain
(23)
Ψ
=
E
+
∑
k
=
1
∞
[
1
+
(
-
1
+
(
-
1
)
k
)
q
k
α
]
2
k
sin
k
E
.From this development, we can obtain up to an arbitrary orderk in the developments, with respect to Ψ, of E and an arbitrary analytical function f
(
E
) using the Deprit inversion algorithm [23]. As an example, up to fourth order in q, we have
(24)
E
=
Ψ
+
{
2
(
1
-
2
α
)
q
+
(
8
α
3
-
12
α
2
+
4
α
)
q
3
}
sin
Ψ
+
{
(
8
α
2
-
8
α
+
1
)
q
2
(
-
128
3
α
4
+
256
α
3
3
-
160
α
2
3
+
32
α
3
)
2
(
256
α
3
3
)
2
+
(
-
128
3
α
4
+
256
α
3
3
-
160
α
2
3
+
32
α
3
)
q
4
}
×
sin
2
Ψ
+
{
-
24
α
3
+
36
α
2
-
40
α
3
+
2
3
}
q
3
sin
3
Ψ
+
{
256
α
4
3
-
512
α
3
3
+
320
α
2
3
-
64
α
3
+
1
2
}
q
4
×
sin
4
Ψ
(25)
sin
E
=
{
+
(
4
α
4
3
-
8
α
3
3
+
8
α
2
3
-
4
α
3
)
1
+
(
-
2
α
2
+
2
α
-
1
)
q
2
+
(
4
α
4
3
-
8
α
3
3
+
8
α
2
3
-
4
α
3
)
q
4
}
sin
Ψ
+
{
(
32
α
3
3
-
16
α
2
+
22
α
3
-
q
3
+
)
(
1
-
2
α
)
q
+
(
32
α
3
3
-
16
α
2
+
22
α
3
-
1
)
q
3
}
sin
2
Ψ
+
{
(
6
α
2
-
6
α
+
1
)
q
2
+
(
-
54
α
4
+
108
α
3
-
72
α
2
+
18
α
-
1
)
q
4
}
×
sin
3
Ψ
+
{
-
64
3
α
3
+
32
α
2
-
38
α
3
+
1
}
q
3
sin
4
Ψ
+
{
250
α
4
3
-
500
α
3
3
+
320
α
2
3
-
70
α
3
+
1
}
q
4
×
sin
5
Ψ
(26)
cos
⁡
E
=
(
1
-
2
α
)
q
+
{
+
(
20
α
4
3
-
40
α
3
3
+
28
α
2
3
-
8
3
)
1
+
(
-
1
+
6
α
-
6
α
2
)
q
2
+
(
20
α
4
3
-
40
α
3
3
+
28
α
2
3
-
8
3
)
q
4
}
cos
⁡
Ψ
+
{
q
+
(
64
α
3
3
-
32
α
2
+
38
α
3
-
2
α
q
-
1
)
q
3
}
×
cos
⁡
2
Ψ
+
{
(
1
-
6
α
+
6
α
2
)
q
+
(
-
90
α
4
+
180
α
3
-
116
α
2
+
26
α
-
1
)
q
4
}
×
cos
⁡
3
Ψ
+
{
-
64
3
α
3
+
32
α
2
-
38
α
3
+
1
}
q
3
cos
⁡
(
4
Ψ
)
+
{
250
α
4
3
-
500
α
3
3
+
320
α
2
3
-
70
α
3
+
1
}
q
4
×
cos
⁡
(
5
Ψ
)
.The mean anomalyM is related to the eccentric anomaly through the Kepler equation M
=
E
-
e
sin
E. Replacing (24) and (25) in the Kepler equation we obtain
(27)
M
=
Ψ
+
{
-
4
α
q
+
(
8
α
3
-
8
α
2
+
4
)
q
3
}
sin
Ψ
+
{
+
(
-
128
3
α
4
+
64
α
3
-
64
α
2
3
8
α
+
4
)
(
8
α
2
-
4
α
-
1
)
q
+
(
-
128
3
α
4
+
64
α
3
-
64
α
2
3
8
α
+
4
)
q
4
}
sin
2
Ψ
+
{
-
24
α
3
+
24
α
2
-
4
α
3
-
4
3
}
q
3
sin
3
Ψ
+
{
256
α
4
3
-
128
α
3
+
128
α
2
3
+
4
α
-
3
2
}
q
4
sin
4
Ψ
.This is the Kepler equation for the natural anomalyΨ and it is necessary to integrate the second member of planetary equation of Lagrange [8, 24].The parameterq is related to e by q
=
(
1
-
1
-
e
2
)
/
e an so, q can be developed as
(28)
q
=
∑
k
=
1
∞
(
2
k
-
1
)
!
!
2
k
k
!
e
2
k
-
1
,
the value of q verifies 0
≤
q
=
e
/
(
1
+
1
-
e
2
)
≤
e.To developa
/
r with respect to Ψ, it is more convenient to use (7) and then we have
(29)
a
r
=
1
+
e
cos
⁡
f
1
-
e
2
.It is easy to see from (19) and (21) that
(30)
tan
f
′
2
=
1
-
e
1
+
e
tan
f
2
and so
(31)
f
′
=
f
+
∑
k
=
1
∞
(
-
1
)
k
2
k
e
k
sin
k
f
.Consequently,(32)
Ψ
=
f
+
∑
k
=
1
∞
(
-
1
)
k
(
1
-
α
)
2
k
q
k
sin
k
f
.Applying Deprit algorithm [23], we obtain the development of cos
⁡
f up to an arbitrary order in e. As an example, up to fourth order in e, we have
(33)
cos
⁡
f
=
(
α
-
1
)
e
+
{
1
+
(
5
α
2
-
3
α
2
2
+
1
)
e
2
+
(
5
α
4
12
-
7
α
3
6
+
13
α
2
12
-
α
3
)
e
4
}
cos
⁡
Ψ
+
{
(
1
-
α
)
e
+
(
8
α
3
3
-
6
α
2
+
13
α
3
-
1
)
e
3
}
cos
⁡
2
Ψ
+
{
(
3
α
2
2
-
5
α
2
+
)
e
2
+
(
-
45
8
α
4
+
63
α
3
4
-
127
α
2
8
+
27
α
4
-
1
)
e
4
}
×
cos
⁡
3
Ψ
+
{
-
8
3
α
3
+
6
α
2
-
13
α
3
+
1
}
e
3
cos
⁡
4
Ψ
+
{
125
α
4
24
-
175
α
3
12
+
355
α
2
24
-
77
α
12
+
1
}
e
4
cos
⁡
5
Ψ
.Replacing (33) in (29) finally we obtain
(34)
a
r
=
1
+
α
e
4
+
α
e
2
+
{
e
+
(
-
3
2
α
2
+
5
α
2
)
e
3
}
cos
⁡
Ψ
+
{
(
1
-
α
)
e
2
+
(
8
α
3
3
-
6
α
2
+
10
α
3
)
e
4
}
cos
⁡
2
Ψ
+
{
3
α
2
2
-
5
α
2
+
1
}
e
3
cos
⁡
3
Ψ
+
{
-
8
3
α
3
+
6
α
2
-
13
α
3
+
1
}
e
4
cos
⁡
4
Ψ
.The convergence of these developments is poor when the eccentricity is near to one. Fortunately, in the case of the planetary motion, the eccentricity values are small and the convergence rate is appropriate. On the other hand, for intermediate values of eccentricity, the length of the series using an appropriate anomalyΨ
α is lower than if the mean anomaly is used as temporal variable.The main problem to develop the second member of the Lagrange planetary equations is to develop the inverse of the distance. This development can be obtained using the Kovalevsky algorithm [25]. Table 1 shows, for the couple Jupiter-Saturn, the performance of the algorithm. This table contains the length of the expansion of the inverse of the distance using the mean anomaly M and the natural anomaly for α
=
0
,
0.25
,
0.50
,
0.75
,
1.0, and the estimated error in astronomical units for each iteration. The planetary elements are taken from [26].Table 1
Length of the1
/
Δ using M
,
Ψ
0
,
Ψ
0.25
,
Ψ
0.5
,
Ψ
0.75,  and Ψ
1.0.
Number of iteration
M
α
=
0.0
α
=
0.25
α
=
0.50
α
=
0.75
α
=
1.0
err
0
2197
1597
1438
1298
1095
419
1
2518
1882
1717
1570
1379
1172
1.01
·
10
-
2
2
2728
2144
1988
1855
1728
1707
1.36
·
10
-
3
3
2841
2294
2182
2102
2028
2056
3.67
·
10
-
5
4
2613
2052
1988
1957
1924
1912
3.50
·
10
-
8
5
2627
2112
2036
2007
1940
1952
6.27
·
10
-
14


## 5. Numerical Examples

In general, the perturbative forces are small, for this reason, it is convenient to test the methods by applying them to the well known two-body problem, referred to the orbital coordinate system(
x
,
y
,
0
), to select an appropriate new temporal variable in order to minimize the distribution of the truncation errors in the orbit. Let us define a generic family Ψ
α of anomalies depending on a parameter α as d
t
=
Q
α
(
r
)
d
Ψ
α. For each α we have
(35)
d
d
t
=
n
d
d
M
=
n
d
d
Ψ
α
d
Ψ
α
d
M
=
n
Q
α
(
r
)
d
d
Ψ
α
.
So,
(36)
d
x
d
Ψ
α
=
Q
α
(
r
)
n
v
x
,
d
v
x
d
Ψ
α
=
-
Q
α
(
r
)
n
[
G
M
x
r
3
+
∂
V
∂
x
-
F
x
]
,
d
y
d
Ψ
α
=
Q
α
(
r
)
n
v
y
,
d
v
y
d
Ψ
α
=
-
Q
α
(
r
)
n
[
G
M
y
r
3
+
∂
V
∂
y
-
F
y
]
.
In order to test the performance of this method, we use a fictitious artificial satellite with the same elements than HEOS II used by Brumberg [17] (a
=
118363.47Km, e
=
0.942572319, i
=
2
8
o
.
16096, Ω
=
18
5
o
.
07554, ω
=
27
0
o
.
07151,   M
0
=
0
o), except for its eccentricity, that it is changed to study the optimum value of α depending on the value of the eccentricity e. In Figure 2, we show a sample of twenty points for Ψ
α with homogeneous distribution over the orbit.Points distribution fore
=
0.7, α
=
0.0
,
0.5
,
0.75
,
1.0.
(a)
α
=
0.0
(b)
α
=
0.5
(c)
α
=
0.75
(d)
α
=
1.0Table2 shows the values of the different α, where the minimum of the errors for this fictitious satellite with the same (
a
,
i
,
Ω
,
ω
,
M
) elements than HEOSII and different values of eccentricity (
e
=
0.0
,
0.05
,
…
,
0.95
) is reached. This table has been carried out using a classic Runge Kutta integrator of fourth order with 1000 uniform steps. In this table, we can see that the value of α, where the errors in position reach their minimum, depends on the eccentricity.Table 2
Optimalα for e
|
Δ
r
→
| in 10
-
5 Km |
Δ
v
→
| in 10
-
8 Km/s.
e
α
|
Δ
r
→
|
|
Δ
v
→
|
e
α
|
Δ
r
→
|
|
Δ
v
→
|
0.00
0.7284
3.730
·
10
-
2
2.893
·
10
-
4
0.50
0.9127
2.357
·
10
-
2
8.684
·
10
-
2
0.05
0.7581
3.555
·
10
-
2
3.426
·
10
-
3
0.55
0.9225
2.111
·
10
-
2
9.858
·
10
-
2
0.10
0.7839
3.451
·
10
-
2
8.409
·
10
-
3
0.60
0.9310
1.852
·
10
-
2
1.091
·
10
-
1
0.15
0.8066
3.374
·
10
-
2
1.483
·
10
-
2
0.65
0.9381
1.643
·
10
-
2
1.170
·
10
-
1
0.20
0.8268
3.303
·
10
-
2
2.253
·
10
-
2
0.70
0.9435
1.479
·
10
-
2
1.199
·
10
-
1
0.25
0.8450
3.229
·
10
-
2
3.143
·
10
-
2
0.75
0.9463
1.353
·
10
-
2
1.094
·
10
-
1
0.30
0.8613
3.116
·
10
-
2
4.118
·
10
-
2
0.80
0.9461
1.264
·
10
-
2
6.894
·
10
-
2
0.35
0.8761
2.980
·
10
-
2
5.179
·
10
-
2
0.85
0.9454
6.821
·
10
-
3
4.313
·
10
-
2
0.40
0.8896
2.810
·
10
-
2
6.319
·
10
-
2
0.90
0.9512
6.064
·
10
-
2
3.600
·
10
-
1
0.45
0.9017
2.599
·
10
-
2
7.483
·
10
-
2
0.95
0.9661
1.121
·
10
0
2.185
·
10
0A least square analytical approach to the optimal value ofα can be written as
(37)
α
=
2.10826
e
5
-
4.77809
e
4
+
3.92513
e
3
-
1.71554
e
2
+
0.71606
e
+
0.72724
.
Figure 3 shows the value of α, where the error |
Δ
r
→
| reaches its minimum for each value of e
∈
[
0,0.95
] and its analytical approximation.Figure 3
Optimal value ofα for each value of e
∈
[
0,0.95
].Figures4, 5, 6 and 7 show the local integration errors, in position and velocity, for a satellite with a
=
118363.47Km and e
=
0.7 for the values of α
=
0
,
0.5,0.75
,
1.0. These errors have been obtained solving for each value of Ψ
α the equation Ψ
α
=
Ψ
α
(
E
). From the value of E, we compute, using the two-body equations, the exact solution for the position and velocity vectors. These values are the initial conditions for the numerical integrator. The local truncation errors are exactly evaluated by comparing the values obtained through integration in the next step with the corresponding ones evaluated solving the equation Ψ
α
+
h
=
Ψ
α
(
E
).Local integration errors distributione
=
0.7, α
=
0.0.
(a)
erx
(b)
ery
(c)
ervx
(d)
ervyLocal integration errors distributione
=
0.7, α
=
0.5.
(a)
erx
(b)
ery
(c)
ervx
(d)
ervyLocal integration errors distributione
=
0.7, α
=
0.75.
(a)
erx
(b)
ery
(c)
ervx
(d)
ervyDistribution of the local integration errorse
=
0.7, α
=
1.0.
(a)
erx
(b)
ery
(c)
ervx
(d)
ervyThe solution of the equationΨ
α
=
Ψ
α
(
E
) is obtained solving, for each value of Ψ
α
=
i
*
h, (16) and (13).To test the robustness of the method, we will study a perturbed case included in the planar restricted three-body problem, this problem includes the Earth, the Moon, and an artificial satellite with the same semiaxe that Heos II and eccentricity0.95. The problem uses the following approach to the motion of the satellite perturbed by the Moon.(1)
The satellite is in the orbital plane of the Moon.
(2)
The Moon motion is approached through a circular motion around the Earth with a sidereal period of27 days 4 h 43 m 11.5 s, orbital radius of the Moon R
L
=
384.400Km, and mass 0,01231 in units of Earth mass.
(3)
The couple Earth-Moon is not perturbed by the satellite motion.Table3 shows the number of steps necessary to provide an accuracy of 1
0
-
4Km in position in an integration over 100 days using as integrator a classic Runge-Kutta of fourth order and a Runge-Kutta of order eight [27]. To evaluate the global error, we have compared the position results using N and 1.1
N steps. In this table, we can see that the integration can be improved by the use of an appropriate choice of α in the natural family of anomalies given by α
(
e
) (37). It is important to emphasize that all the anomalies of this family improve the values obtained by the use of the mean anomaly.Table 3
Number of steps to get a precision of|
Δ
r
→
|
<
1
0
-
4 Km.
N
α
=
0
α
=
0.25
α
=
0.5
α
=
0.75
α
=
1.0
α
(
e
)
RK4
7084079
3292627
401504
251229
190193
158367
146393
RK8
305033
94362
14303
10833
9025
6840
6630


## 6. Concluding Remarks

In this paper, a new one-parametric family of anomalies, called natural anomalies, has been defined.This family of anomalies is adequate to be used in the construction of the analytical theories of planetary motion. In this sense, we have described a method to obtain the analytical development as Fourier series of the natural anomaly up an arbitrary order in eccentricity of the most common quantities (25), (26), and (29) of the two-body problem, and from these developments, we can obtain the development of the second member of the Lagrange planetary equation.To test the efficiency of the use of the variablesΨ
α as independent variables in the analytical developments, we show the performance of the Kovalevsky algorithm to obtain the inverse of the distance for the couple Jupiter-Saturn, obtaining more compact developments for values of α near to values given by (37). To integrate the Lagrange planetary equations by analytical methods using a generalized anomaly as temporal variable, it is necessary to use the generalized Kepler equation (27).The natural anomalies conform to a parametric family of anomalies that includes the true anomaly. Forα
=
0, we have the secondary true anomaly, its value being near to the mean anomaly for small values of the eccentricity and, in the case of an eccentricity that is not small, the spatial points distribution over the orbit is the most appropriate, if the mean anomaly is used.In the symmetrical case,Ψ
1
/
2 defines a point distribution over the orbit with a major concentration in the perigee than if the eccentric anomaly E is used.The two-body problem has been used to test the performance of the numerical integrators. In this study, the local and global errors depend on the value ofα. The global integration errors in a revolution depend on α. This value increases with the eccentricity, and it is near 0.72 for lower eccentricities and 0.9661 for higher ones (
e
=
0.95
).To test the robustness of the method, a simplified problem, included in the restricted three-body problem, has been analyzed. In this problem, we study the number of steps that is necessary to get a precision of1
·
10
-
4Km in |
Δ
r
→
| using a classical method of Runge-Kutta of fourth order and a Runge-Kutta of eighth order. The number of steps is minimum when we take for each step the value of α that minimizes the error in the osculating two-body problem (37). The result is similar using a classic RK4 method and using a RK8 method, obviously with a minor number of steps in the last case.



---

*Source: 162060-2014-02-24.xml*