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.
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: where is the radius vector of the satellite, is the gravitational constant, is the Earth mass, is the spaceflight constant, 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 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 and , where .
The second way involves the solution of the two-body problem: 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 be the third set of elements of Brouwer and Clemence [1] where is the major semiaxis, the eccentricity, the inclination of the orbit with respect to the plane, the argument of the ascending node, the argument of the perigee, and the mean anomaly at an epoch . In the elliptic case, let be the mean motion, with . The mean anomaly is connected to the mean motion through the formula , where is the epoch of the perigee passage. The quantities , , , , 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 anomaly . This anomaly is related to the mean anomaly through the Kepler equation The coordinates of the secondary in the orbital coordinate system are where is the center of mass of the Earth, the axis running to the perigee, the axis pointing to the angular momentum, and the axis making a direct trihedron. The radius vector of the satellite is given by The spatial coordinates are related to the orbital coordinates by means of where . defines a rotation around the -axis. The spatial speed of the satellite is given by and the orbital velocity 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) with the initial conditions , .
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 by means of 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 , where . This family includes the most common anomalies for appropriate values of and , the mean anomaly for , , the eccentric anomaly for , the true anomaly for , , and the Nacozy intermediate anomaly for and [4].
The interest in this issue can be seen in later papers, for instance, the ones by Ferrándiz et al. [6] who suggested a more complex family of transformations, where . Furthermore, Brumberg [7] suggested the use of the regularized length of arc as temporal variable that provides the function given by
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 Section 2 the general Sundman anomaly properties. We also obtain the differential equations of motion using an arbitrary general Sundman anomaly as temporal variable.
In Section 3 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
Let be a generalized Sundman transformation and let us define the generalized Sundman anomaly as the regularized value of given by where if , for all , and . The constant value of has been computed by López Ortí et al. [8] where is the hypergeometric function, , , and . For each we have so,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 obtained through the derivative of the Kepler equation . Replacing in (11), we obtain the equation and solving this implicit integral equation we get the value . 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 Figure 1, twenty points for an ellipse with and . Values of are .
| (a) |
| (b) |
| (c) |
| (d) |
To 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 Km, , , , , and . The period of the satellite is days and for the Earth the spaceflight constant is .
Table 1 shows the global error in position for one revolution, in units of Km, and the error in velocity in units of . 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 . In the same conditions, the Brumberg method produces the errors: Km and , for the regularized length of arc parameter . Compare, for example, with the results for in Table 1.
Figures 2, 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 for the values . These errors has been obtained by solving for each value of the equation (15). From the value of 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 . 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.
| (a) |
| (b) |
| (c) |
| (d) |
| (a) |
| (b) |
| (c) |
| (d) |
| (a) |
| (b) |
| (c) |
| (d) |
A least square analytical approach to the optimal value of can be written as Figure 5 shows the values of where the error reaches its minimum for each value of (dots) and its analytical approximation (continuous line).
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 mass is placed at the origin of coordinates. The second fix mass is placed at .(2)The initial osculating elements of the satellite are the same elements used in the unperturbed problem.(3)The values and have been chosen as and , where is the Earth mass and 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 using as temporal variable, it is necessary to determinate the value of , where depends on . For this purpose it is neccesary to solve the boundary problem: with the boundary conditions given by , , and . The value of the disturbing potential is given by
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 Km in position is minimum. Table 3 shows this required number of steps in an integration over 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 .
To evaluate the global error we have compared the results for the position using and . 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 (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 for the lowest eccentricities and for the highest.
To study the performance of the method for high eccentricities () 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 but the global error is better for the value .
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 a Km precision in position. The best results are obtained when we use for each step the optimal value for the osculating eccentricity.
Appendix
The Two-Fix-Center Problem
Let and be two punctual masses placed at the fix points of coordinates , , respectively, and let be the mass of a punctual body in the plane. This body is in motion in this plane under the action of the gravitational potential created by and . This is an integrable problem [12, 13].
To obtain the solution of this problem it is convenient to introduce the planar ellipsoidal coordinates defined as The metric element in this coordinate system is given by and consequently the kinetic energy can be written as
The potential energy at a point of coordinates can be written as [12] and operating as The Lagrangian function is Given that the motion is independent of the value, we take for simplicity. The Canonical momentum is defined as , and so and the kinetic energy is
The Hamiltonian function can be written as where and .
Hamiltonian does not depend explicitly on time , so it is a first integral . In this condition the Hamilton-Jacobi equation of the problem can be written as Multiplying the last equation by , we obtain This equation is separable, so it can be replaced by the system where and are constant and their values only depend on the initial conditions.
The action function is given by
The function provides a free Canonical transformation into defined by the equations The Hamiltonian function in the new variables can be written as and applying the Jacobi theorem [13] we have
The exact solution of the two-fix-center problem can be obtained from (A.15) and (A.14) after solving (A.13).
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-06I455.01/1 from Bancaja.