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

and Applied Analysis 3 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 equationM 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 d 2 ⃗ r dt 2 = −GM ⃗ r r 3 − ⃗ ∇U − ⃗ F, (6) 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 Figure 1 the elliptic orbit corresponding to the motion of the two-body problem. This ellipse is defined by its major semiaxis a = OQ and its eccentricity e = c/a, 0 ≤ e < 1, where c is the focal semidistance c = FF󸀠/2, and the minor semiaxis b is defined as b = a1 − e2. 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 andF, respectively. P


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][2][3][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][7][8][9] can be interesting.To construct an analytical or semianalytical planetary theory, it is necessary to obtain a set of developments for the twobody 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, 2 Abstract and Applied Analysis 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  through  = .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: where (, ) =     is the called partition function.A more complicated transformation was introduced by Ferrándiz et al. [16] () =  2/3 ( 0 +  1 ) −1/2 and Brumberg [17] who proposed the use of the regularized length of arc: where  is the spaceflight constant of the Earth and  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 Ψ  =  + (1 − )  , where  is the true anomaly and   the polar angle referred to the secondary focus called secondary true anomaly.
The rest of this article is organized as follows.
In Section 2, 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 Section 3, the properties of natural family of anomalies are also described.These properties include the connection between the quantities  and   with   , the connection between the natural anomaly Ψ  and the eccentric anomaly , and the study of the partition function denoted by   () for this case.
In Section 4, the analytical developments according to the natural anomaly are studied.This study contains the developments of , sin , cos , and 1/ 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 Section 5, 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 Section 6, the main conclusions and remarks of this paper are exposed.

Basic Equations
The analytical methods are based on the solution of the twobody problem (Sun planet) through a set of orbital elements, for example, the third set of Brouwer and Clemence [18] (, , , Ω, , ), where  =  0 + ( −  0 ),  is the mean motion, and  0 is the initial epoch being (, , , Ω, , ) 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]: where  is a new variable defined by the equation: and it coincides with  0 in the case of the unperturbed motion. is the disturbing potential  = ∑  =1   due to the disturbing bodies  = 1, . . ., .It is defined as [19] where ⃗  = (, , ) and ⃗   = (  ,   ,   ) are the heliocentric vector position of the secondary body and the th disturbing body, respectively, Δ  is the distance between the secondary body and the disturbing body, and   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 Ψ  , it is necessary to obtain for each planet  the developments of the coordinates and the inverse of the radius in Fourier series of Ψ  .Then, the integration of the Lagrange planetary equations with respect to the Ψ  anomalies requires to obtain solution of the corresponding Kepler equation   =   (Ψ  ), 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 where ⃗  is the radius vector of the satellite,  the potential from which the perturbative conservative forces are derived, and ⃗  includes the nonconservative forces.To integrate system (6), it is necessary to know the initial values of the radius vector ⃗  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.

Natural Anomalies
In this section, a new family of anomalies depending on a parameter is defined.Let us represent in Figure 1 the elliptic orbit corresponding to the motion of the two-body problem.This ellipse is defined by its major semiaxis  =  and its eccentricity  = /, 0 ≤  < 1, where  is the focal semidistance  =   /2, and the minor semiaxis  is defined as  =  √ 1 −  2 .Let  be the center of the ellipse, let  be the primary focus, let   be the secondary focus (also called equality point), let  be the pericenter, and let  be the position of the secondary in the orbit.Let us define the coordinates (, ) referred to the primary focus and let  and   be the distance between the secondary  and the primary focus and the secondary focus,  and   , respectively.
The angle  is called the true anomaly and for the angle   , we propose the name "secondary true anomaly".
The quantities  and   satisfy On the other hand, the coordinates of the secondary with respect to the primary, in the orbital plane, are The spatial orbital coordinates ⃗  = (, , )  are related to the orbital coordinates (, , 0)  by means of where  =  1 (−Ω) 3 (−) 1 (−).  defines a rotation around the -axis.From (7) we have and from (10), it is easy to obtain  2 sin  =  2 sin     .Taking into account (8), we obtain The radii  and   are related to the eccentric anomaly  through The anomalies  and   are related to the eccentric anomaly  by To relate the anomaly   to the mean anomaly  we use the integral of areas  = ( 2 √ 1 −  2 / 2 ), and after replacing in (11), we get To compare   and  up to second order in  we have and so, for small values of eccentricity, the anomaly   is near .
Let us define the natural family of anomalies Ψ  as this family includes, for  = 0, the secondary true anomaly   and, for  = 1, the true anomaly .The anomaly Ψ  is related to the mean anomaly by where the function   () is defined as

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/ and the mean anomaly .
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  and cos .To obtain these developments, it is necessary first to expand Ψ with respect to  [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] tan where  = tan(/2),  = sin  and so The secondary true anomaly verifies tan and consequently Replacing ( 20) and ( 22) in ( 16) we obtain From this development, we can obtain up to an arbitrary order  in the developments, with respect to Ψ, of  and an arbitrary analytical function () using the Deprit inversion algorithm [23].As an example, up to fourth order in , we have Abstract and Applied Analysis 5 The mean anomaly  is related to the eccentric anomaly through the Kepler equation  =  −  sin .Replacing (24) and (25) in the Kepler equation we obtain 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].
To develop / with respect to Ψ, it is more convenient to use (7) and then we have It is easy to see from ( 19) and ( 21) that and so Consequently, Applying Deprit algorithm [23], we obtain the development of cos  up to an arbitrary order in .As an example, up to fourth order in , we have cos  = ( − 1) (34)  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  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].

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 Abstract and Applied Analysis 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] ( = 118363.47Km,  = 0.942572319,  = 28 o .16096,Ω = 185 o .07554, = 270 o .07151, 0 = 0 o ), except for its eccentricity, that it is changed to study the optimum value of  depending on the value of the eccentricity .In Figure 2, we show a sample of twenty points for Ψ  with homogeneous distribution over the orbit.
Table 2 shows the values of the different , where the minimum of the errors for this fictitious satellite with the same (, , Ω, , ) elements than HEOSII and different values of eccentricity ( = 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.
(37) Figure 3 shows the value of , where the error |Δ ⃗ | reaches its minimum for each value of  ∈ [0, 0.95] and its analytical approximation.
Figures 4, 5, 6 and 7 show the local integration errors, in position and velocity, for a satellite with  = 118363.47Km and  = 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 Ψ  = Ψ  ().From the value of , 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 Ψ  + ℎ = Ψ  ().
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 eccentricity 0.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 of 27 days 4 h 43 m 11.5 s, orbital radius of the Moon   = 384.400Km, and mass 0, 01231 in units of Earth mass.
(3) The couple Earth-Moon is not perturbed by the satellite motion.
Table 3 shows the number of steps necessary to provide an accuracy of 10 −4 Km 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  and 1.1  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 () (37).It is important to emphasize that all the anomalies of this family improve the values obtained by the use of the mean anomaly.

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  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
coordinate system (, , 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  =   ()Ψ  .For each  we have