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.
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 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] 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 , 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 , , , and 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.
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] , where , is the mean motion, and 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 in the case of the unperturbed motion. is the disturbing potential due to the disturbing bodies . 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 and velocity .
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 and its eccentricity , , where is the focal semidistance , and the minor semiaxis is defined as . 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 by means of where . defines a rotation around the -axis.
From (7) we have and from (10), it is easy to obtain . 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 , 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 , the secondary true anomaly and, for , the true anomaly . The anomaly is related to the mean anomaly by where the function is defined as
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 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 and . 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] where , and so The secondary true anomaly verifies 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
The mean anomaly is related to the eccentric anomaly through the Kepler equation . 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].
The parameter is related to by an so, can be developed as the value of verifies .
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 up to an arbitrary order in . As an example, up to fourth order in , we have
Replacing (33) in (29) finally we obtain
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 , and the estimated error in astronomical units for each iteration. The planetary elements are taken from [26].
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 , 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 So, 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] (Km, , , , , ), 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.
| (a) |
| (b) |
| (c) |
| (d) |
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 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.
A least square analytical approach to the optimal value of can be written as Figure 3 shows the value of , where the error reaches its minimum for each value of and its analytical approximation.
Figures 4, 5, 6 and 7 show the local integration errors, in position and velocity, for a satellite with Km and for the values of . 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 .
(a) erx
(b) ery
(c) ervx
(d) ervy
(a) erx
(b) ery
(c) ervx
(d) ervy
(a) erx
(b) ery
(c) ervx
(d) ervy
(a) erx
(b) ery
(c) ervx
(d) ervy
The solution of the equation is obtained solving, for each value of , (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 eccentricity . 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 days 4 h 43 m 11.5 s, orbital radius of the Moon Km, and mass 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 Km in position in an integration over 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 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.
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 , 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, 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 integration errors in a revolution depend on . This value increases with the eccentricity, and it is near for lower eccentricities and for higher ones .
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 of Km in 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This research has been partially supported by Grant P1.1B2012-47 from University Jaume I of Castellón.