Mathematical Methods Applied to the Celestial Mechanics of Arti๏ฌcial Satellites
View this Special IssueResearch Article  Open Access
Martin Lara, Juan F. San Juan, Luis M. Lรณpez, "Semianalytic Integration of HighAltitude Orbits under Lunisolar Effects", Mathematical Problems in Engineering, vol. 2012, Article ID 659396, 17 pages, 2012. https://doi.org/10.1155/2012/659396
Semianalytic Integration of HighAltitude Orbits under Lunisolar Effects
Abstract
The longterm effect of lunisolar perturbations on highaltitude orbits is studied after a double averaging procedure that removes both the mean anomaly of the satellite and that of the moon. Lunisolar effects acting on highaltitude orbits are comparable in magnitude to the Earthโs oblateness perturbation. Hence, their accurate modeling does not allow for the usual truncation of the expansion of the thirdbody disturbing function up to the second degree. Using canonical perturbation theory, the averaging is carried out up to the order where secondorder terms in the Earth oblateness coefficient are apparent. This truncation order forces to take into account up to the fifth degree in the expansion of the lunar disturbing function. The small values of the moonโs orbital eccentricity and inclination with respect to the ecliptic allow for some simplification. Nevertheless, as far as the averaging is carried out in closed form of the satelliteโs orbit eccentricity, it is not restricted to loweccentricity orbits.
1. Introduction
In an increasingly saturated space about the Earth, aerospace engineers confront the mathematical problem of accurately predicting the position of Earthโs artificial satellites. This is required not only for the correct operation of satellites but also for preserving the integrity of space assets. Thus, operational satellites are threatened by the not so remote possibility of a collision with a defunct satellite [1] but most probably by the impact with other uncontrolled manmade space objects as spent rocket stages or collision fragmentsโall of them commonly called space debris.
Precise predictions require the integration of complete force models including both gravitational and nongravitational effects, like a highdegree and order geopotential, ephemeridesbased lunisolar perturbations, drag, solar radiation pressure including eclipses, and so forth (see [2, 3], for instance). The most accurate integration is expected from numerical methods, although precision ephemeris can also be obtained by means of semianalytical integration [4]. In fact, both approaches, numerical and semianalytical, do not need to enter a competition. Thus, for instance, while semianalytical methods may be efficient in keeping a running catalog of hundreds of thousands of space objects within a reasonable accuracy, if the probability of impact with an operational satellite is detected to surpass a certain threshold, then a more accurate numerical integration will help control engineers in deciding whether a collision avoidance maneuver is required or, on the contrary, the integrity of the satellite is not in risk [5].
In a semianalytical approach the highest frequencies of the motion, which normally have small amplitudes, are filtered analytically via averaging procedures. This filtering allows the numerical integration of the averaged system to proceed with very long step sizes. Then, the shortperiod terms are recovered analytically, if desired, at any step of the numerical integration [6โ8].
Averaging techniques are also used for exploring questions affecting stability, such as those derived from tesseral resonances or thirdbody perturbations, in a reducedphase space [9]. In this respect, much attention has been recently paid to the longterm evolution of classical GNSS constellations, either for operational or disposal orbits [10].
While the noncentralities of the Earth gravitational potential play a key role in the motion of low altitude satellites, thirdbody perturbations have also a decisive influence in the longterm evolution of medium and highaltitude Earth orbits. The thirdbody disturbing function is commonly given by a series expansion in Legendre polynomials. Often, the series is truncated to the first term in the expansion [11], but this early truncation is not always accurate enough to represent the real dynamics [4, 12, 13]. Nevertheless, recursions in the literature allow to extend the Legendre polynomial expansion to any desired order either in classical or nonsingular elements [14โ16].
Because of the physical characteristics of the Earth gravitational potential, where the secondorder zonal coefficient () clearly dominates over all other harmonics, secondorder effects of may be important and must be taken into account when the effect of higherorder harmonics is studied. Correspondingly, the truncation in the expansion of thirdbody perturbations must include terms of magnitude comparable to squared. Because the thirdbody disturbing function is expanded in the ratio semimajor axis of the satelliteโs orbit to semimajor axis of the thirdbodyโs orbit, the degree at which the expansion must be truncated depends on the altitude of the satellites to be studied.
In this paper we study the effect of lunisolar perturbations on highaltitude orbits about a noncentral Earth, which is assumed to be oblate although without equatorial symmetry. More specifically, we are interested in the semianalytical integration of satellites on altitudes of classical global navigation satellite system (GNSS) constellations such as GPS, Glonass, or Galileo. Note, however, that the assumption of an axisymmetric geopotential prevents to tackle the tesseral resonance problem that commonly suffer this kind of orbits. With respect to previous research [17], where we approached the general case of thirdbody perturbations via double averaging, we release here the common simplification of assuming the thirdbody in circular orbit. Also, we focus on the case of Earthโs artificial satellites dealing with more real lunisolar perturbations. We do that because recent results [18] seem to contradict the claim that taking the thirdbody in circular orbit does not produce any noticeable degradation in the longterm propagation of real Earth orbits [19]. Besides, for the actual values of the orbits of the sun and moon, neglecting terms in the eccentricity is not consistent with a higherorder expansion of the lunisolar disturbing function. Hence, both the eccentricity of the orbits of the sun and moon are taken into account in the present work. For the latter, the moon, the orbit is assumed to remain with constant inclination with respect to the ecliptic plane, over which the longitude of the ascending node moves with linear motion. The argument of the perigee of the moon is also assumed to evolve linearly, while we take the apparent orbit of the sun to be purely Keplerian.
We use canonical perturbation theory by Lie transforms [20โ22]. The order of each term of the disturbing function is determined by a virtual small parameter that is taken proportional to the ratio of the satelliteโs orbit semimajor axis to the moonโs orbit semimajor axis. Then, we check that the magnitude of second order terms due to the Earth oblateness is comparable to that of the fifth degree in the expansion of the moonโs thirdbody potential, and hence we truncate the moonโs thirdbody potential up to the fifth degree. On the contrary, the usual truncation up to the second degree is enough for modeling the sun disturbing function assumed that we limit the theory to the order of squared terms. The Hamiltonian model also takes into account the asymmetries in the Earth potential caused by the term, whose influence in the dynamics of the lower eccentricity orbits is clearly noted.
The initial Hamiltonian is of three degrees of freedom and time dependent. Note, however, that the time appears in different scales related to the mean motion of the moon, that of the sun, and the rate of variation of the moonโs argument of the perigee and longitude of the ascending node. In our averaging procedure, we eliminate the mean anomaly of the satellite and that of the moon to obtain a twodegreesoffreedom Hamiltonian, which still depends on time albeit only through very slowly varying quantities. As far as we do not find resonances between the mean motions of the satellite and the moon, the averaging can be done in two steps: the mean anomaly of the satellite is removed first by means of a classical Delaunay normalization [23]; then, a similar transformation is used to eliminate the mean anomaly of the moon. The splitting of the averaging has the advantage of simplifying computations. Specifically, the generating function of each transformation is obtained from the solution of simple quadratures.
The averaging is carried out not only in closed form of the satelliteโs orbit eccentricity, thus making the simplified Hamiltonian useful for studying the longterm evolution of any orbit, but also in closed form of the eccentricity of the moonโs orbit. Nevertheless, because of the small values of the eccentricity of the orbits of both perturbing bodies as well as the small inclination of the orbit of the moon with respect to the ecliptic, further simplifications can be done by neglecting terms of order higher than the truncation order of the theory.
Numerical experiments using the doubly averaged Hamiltonian demonstrate that the inclusion in the model of the orbital eccentricities of the sun and moon does not cause qualitative differences with respect to the circular orbit approximation. Besides, the recovery of the shortperiod effects by means of the analytical transformation equations of the averaging provides a quite reasonable precision in the longterm integrations.
2. Model
We study the motion of an Earth artificial satellite of negligible mass whose Keplerian motion is slightly distorted by the noncentralities of the geopotential due to the Earthโs oblateness and latitudinal asymmetry, as represented by the harmonic coefficient and , respectively, and under the pointmass attraction of the sun and moon. Thus, the motion of the satellite is derived from the potential where is the Earthโs gravitational parameter, its equatorial radius, is the zonal harmonic coefficient of degree , is distance from the origin, and is the satelliteโs coordinate in the direction of the symmetry axis of the Earth. In the masspoint approximation, the thirdbody disturbing potential , , is where is the thirdbodyโs gravitational parameter, and are the radius vector of the satellite and of the thirdbody, respectively, of corresponding modulus and . If the disturbing body is far away from the perturbed body, (2.2) can be expanded in power series of the ratio : where , is the mean motion of the thirdbody in its orbit of semimajor axis , are Legendre polynomials, and is the angle encompassed by and . The absence of the Keplerian term as a summand in (2.3) has no effect in the restricted problem approximation, in which the mass of the satellite has negligible effects on the thirdbodyโs orbit.
In our model we assume an Earthcentered frame with the origin defined by the intersection of the Earthโs equator and the ecliptic; we both of which consider fixed planes. Whereas we take the sunโs apparent orbit about the Earth to be purely Keplerian, we assume that the moon moves with Keplerian motion on a precessing ellipse that evolves with constant rate of the argument of perigee. Besides, the moonโs orbital plane is assumed to have constant inclination with respect to the ecliptic; yet it is also assumed to be precessing over the ecliptic with linear variation of the longitude of the ascending node.
Calling , , , , , and to semimajor axis, eccentricity, inclination, argument of the perigee, longitude of the ascending node, and mean anomaly, respectively, we use the approximate values for the sun, where notes the obliquity of the ecliptic, we assume a fictitious epoch in which , and take /sidereal year. For the moon, and take , , with We recall that and are referred to the ecliptic. Besides, for the sun and for the moon.
Remark that this simple model is not, of course, valid for computing precise ephemeris of the moon, who may be clearly out of the predicted position because of the amplitude of periodic oscillations in and , and also in and , comper [24]. Nevertheless, it will be enough for our purposes of investigating the qualitative features that lunisolar perturbations introduce in the longterm behavior of Earthโs artificial satellites.
3. Perturbation Approach
By using canonical variables we can study the problem in Hamiltonian form. Besides, in order to apply perturbation theory, we arrange the Hamiltonian as a power series in a small parameter:
We are interested in highaltitude orbits in the range of 20โ000 to 35โ000โkm, the socalled upper mean earth orbit (MEO) region. In order to take account of all relevant lunisolar perturbations, we take a virtual small parameter proportional to the ratio semimajor axis of the satellite semimajor axis of the moon. For the altitudes of interest, this ratio is . Then, we find that the Earthโs oblateness coefficient is a fourthorder quantity, and the effect appears at the eighth order. With respect to the moon the consecutive terms of the Legendre polynomials expansion of the moon disturbing function appear at consecutive orders of the Hamiltonian, starting by the fifth order. We include up to the fifth degree, whose effect may be comparable to secondorder effects of . For the sun, it is enough to take the first term in the Legendre polynomials expansion. Then, the zeroorder term is the Keplerian Besides, for , and where with where is the true anomaly and is the usual eccentricity function. The equations of motion are obtained from the Hamilton equations:
Recall that the Cartesian coordinates of the satellite are expressed in orbital elements by means of simple rotations. Thus, where is the argument of latitude and and are the usual rotation matrices about the  and axes, respectively.
Since the Hamiltonian must be expressed in canonical variables, we assume hereafter that the orbital elements of the satellite are always expressed as function of the Delaunay variables (, , , , , ), given by the known relations , , , , , and . Nevertheless, we use the orbital elements notation because of its immediate physical insight.
Delaunayโs variables are singular for zero eccentricity and/or zero inclination. In spite of that, the Lie transforms theory is naturally computed in Delaunay variables. Once the generating function of the averaging is computed, it can be applied to any function of the Delaunayโs variables. Specifically nonsingular variables are assumed to be functions of Delaunay variables to obtain the transformation equations of the averaging in nonsingular variables [25].
We note that the sun and moon coordinates are assumed to be known functions of time. Therefore, (3.1) is a timedependent Hamiltonian of three degrees of freedom. In our perturbation approach, we avoid dealing with time by moving to an โextendedโโ phase space (see, for instance, [26]). Since the timedependent variables evolve in quite different time scales, we find convenient to introduce four new pairs of canonical conjugate variables: for the mean anomaly of the sun and its conjugate momentum the same for the moon, for the moon argument of perigee and its conjugate momentum, and for the moon longitude of the ascending node and its conjugate momentum. The specific values of the introduced canonical momenta are irrelevant for our purposes, and we find convenient to make the ordering
Once we have set the Hamiltonian order, we will apply perturbation theory by the Lie transforms in order to filter the shortperiod effects in the potential equation (2.1). More precisely, we base our computations on Depritโs algorithm, which is specifically designed for automatic computation by machine [21, 27].
3.1. The Delaunay Normalization
We compute the canonical transformation that removes the mean anomaly of the satellite from the Hamiltonian up to the eighth order. Note that the mean anomaly does not appear explicitly but through its dependence on the true anomaly. This fact introduces some subtleties in the averaging procedure, but the usual differential relations between true, mean, and eccentric anomalies allow to carry out the normalization on closed form of the eccentricity.
After normalization, we get the new Hamiltonian where the terms are expressed in the prime variables. We find for , and where and stand for sine and cosine of the inclination, respectively; where is an integer division. The eccentricity coefficients , the inclination ones , and the thirdbody direction coefficients and are given in Tables 1, 2, and 3.



Note that except for squared terms the new Hamiltonian is obtained by the simple average over the mean anomaly of all the terms in (3.1). Besides, we introduced Kozaiโs arbitrary constant [28] in the solution of the fourthorder generating function to keep the prime variables as close as possible to the average value of corresponding osculating ones.
The semimajor axis of the satellite remains constant after averaging, as well as its canonical partner, the Delaunay action . Then, the time evolution of the mean anomaly decouples from the twodegreesoffreedom, timedependent system The numerical integration of this system can be done with longer step sizes than the original one because of the filtering of short periodic effects via averaging. At each step of the numerical integration, the osculating elements can be recovered analytically using the transformation equations computed also by the Lie transforms method.
3.2. Elimination of the Mean Anomaly of the Moon
A new canonical transformation removes the mean anomaly from the Hamiltonian. The algorithm starts now by setting the new Hamiltonian where . Then, we express the direction of the moon in terms of the moonโs orbital elements by the composition of rotations: where . The components of the moon direction vector are then replaced in coefficients , of Table 1. Finally, application of the Lie triangle provides the new Hamiltonian where the terms are expressed in the double prime variables, as well as the generating function of the transformation.
The new averaging is similar to the preceding Delaunay normalization in the sense that we base on the differential relations between the true and mean anomalies to perform the averaging. As before, up to the computed order in the perturbation theory, there is no coupling between the different Hamiltonian terms, which, therefore, reduce to their averaging over the mean anomaly of the moon . Thus, for , and where, from (3.11), with
To get some further simplification, we note that ~~. Thus, in our eighthorder theory, we neglect higherorder terms factored by such that . Under these simplifying assumptions, the coefficients and are presented in Tables 4 and 5, respectively.


As in the preceding averaging, adequate arbitrary constants have been introduced in the computation of the generating function to guarantee its average to zero.
After the double averaging, the Hamiltonian only depends on longperiod terms related to the sunโs apparent motion and very longperiod terms related to the precession of the nodes and recession of the line of apsides of the moonโs orbit. The numerical integration of corresponding Hamilton equations is, now, very much faster and efficient.
4. Numerical Experiments
For the numerical tests we used a higherorder RungeKutta method. Specifically, the numerical integration was performed with the Dormand and Prince implementation of the RungeKutta method coded in FORTRAN 77 by Hairer et al. [29].
To illustrate the significance of recovering the shortperiod effects up to higher orders, we first show in Figure 1 a sequence of the errors obtained in the semimajor axis after one month semianalytical propagation. For this example, we take the initial osculating elements , , deg., , deg., and ; besides we assumed that both the mean anomalies of the moon and the sun are zero at the origin of time. Corresponding elements in the single and doubleaveraged phase space depend on the truncation order of the perturbation theory and are presented in the top and bottom parts of Table 6, respectively.

In reference to Figure 1, the top plot shows a direct comparison between the numerical integration of the Hamilton equations of the original, nonaveraged problemโwhose disturbing potential is given in (2.1)โand that of the longterm Hamiltonian after the double averaging up to the eight order of the small parameter. Then, from top to bottom, we show the errors obtained at each step of the integration when recovering shortperiod terms by computing the transformation equations of the averaging up to the fourth, fifth, sixth, seventh, and eighth orders. For the latter, the amplitude of the periodic errors is reduced to several centimeter; note that, in addition to shortperiod errors related to the orbital period of the satellite, we can appreciate a twoweek modulation related to the moonโs mean motion. Remark that the amplitude of the periodic errors roughly divides by ten with each order of the transformation equations, which is consistent with the assumed magnitude of the virtual small parameter .
The shortperiod effects can be ignored in the study of the longterm orbital behavior, where the simple propagation of the doubleaveraged equations appears very fast and efficient. Sample propagations are shown in figures below, which show the important effects that have the initial right ascension of the ascending node and argument of perigee of the satelliteโs orbit in the long term and specifically in the eccentricity and inclination evolution of the satelliteโs orbit. Thus, Figure 2 shows the notably different evolution of the satelliteโs eccentricity and inclination for different initial nodes; for the other initial conditions we have taken the same as in the preceding shortterm propagations but now assuming directly that they are mean elements. In fact, the more relevant parameter is the difference between the node of the satelliteโs orbit and that of the moonโs orbit, as illustrated in Figure 3 where the initial longitude of the ascending node of the moonโs orbit over the ecliptic has been taken as instead of deg.
(a)
(b)
(a)
(b)
Figure 4 shows that the effect of the initial argument of the perigee of the satelliteโs orbit is also important in the long term, eccentricity evolution. The effect is almost negligible in the orbital inclination in the longterm and it is not presented.
Finally, we must mention that further tests demonstrated that there are no important qualitative differences in the long term when using lowerorder truncations of the theory, resulting in faster numerical integration of the mean elements. So the fifthorder truncation should be the preferred longterm Hamiltonian. Besides, we also checked for a variety of orbits that making does not either introduce qualitative differences in the longterm, in agreement with [19].
5. Conclusions
Modeling lunisolar perturbations on highaltitude Earth orbits requires to retain high degrees in the Legendre polynomials expansion of the thirdbody disturbing function of the moon. In consequence the eccentricity of the thirdbodyโs orbit cannot be neglected in the case of either the moon or the sun.
The longterm behavior of highaltitude Earth orbits is approached in a semianalytical way via averaging procedures, in which we take advantage of the different scales in which appears the time to do the averaging in the extended phase space. In addition, up to secondorder terms in the Earthโs oblateness coefficient, the averaging has been computed in closed form of the eccentricity, and, therefore, the semianalytical integration can be applied to any orbit.
Sample numerical propagations of test cases show that the more relevant parameter on the longterm behavior is the difference between the right ascension of the ascending node of the satelliteโs orbit and the longitude of the ascending node of the moonโs orbit over the ecliptic, having an apparent effect manifested by the almostsecular growing of the orbit eccentricity and also by verylongperiod oscillations of the inclination with an amplitude of several degrees. Also the initial argument of the perigee of the satelliteโs orbit has notable effects in the satelliteโs orbit, but only in what respects to the longterm evolution of the eccentricity.
Acknowledgments
Part of this research has been supported by the Government of Spain (Projects AYA 200911896, AYA 201018796, and Grant Gobierno de La Rioja Fomenta 2010/16).
References
 T. S. Kelso, โAnalysis of the Iridum 33—Cosmos 2251 Colission,โ Paper AAS 09368, American Astronautical Society, 2009. View at: Google Scholar
 O. Montenbruck and E. Gill, Satellite Orbits: Models, Methods, and Applications, Springer, Berlin, Germany, 2000. View at: Zentralblatt MATH
 D. A. Vallado, Fundamentals of Astrodynamics and Applications, Microcosm Press and Kluwer Academic, El Segundo, Calif, USA, 2007. View at: Zentralblatt MATH
 P. J. Cefola, V. Yurasov, Z. Folcik, E. Phelps, R. J. Proulx, and A. Nazarenko, โComparison of the DSST and the USM SemiAnalytical Orbit Propagators,โ Paper AAS 03236, American Astronautical Society, 2003. View at: Google Scholar
 E. J. Fletcher, โStatus and progress in the space surveillance and tracking segment of ESA’s space situational awareness programme,โ in Proceedings of the AMOS Conference, Maui, Hawaii, USA, September 2010. View at: Publisher Site  Google Scholar
 W. D. McClain, A Recursively Formulated FirstOrder Semianalytic Artificial Satellite Theory Based on the Generalized Method of Averaging : The Generalized Method of Averaging Applied to the Artificial Satellite Problem, vol. 1, Computer Sciences Corporation CSC/TR77/6010, 1977.
 S. Valk, A. Lemaître, and F. Deleflie, โSemianalytical theory of mean orbital motion for geosynchronous space debris under gravitational influence,โ Advances in Space Research, vol. 43, no. 7, pp. 1070โ1082, 2009. View at: Publisher Site  Google Scholar
 J. F. SanJuan, L. M. López, and R. López, โHigherorder analytical attitude propagation of an oblate rigid body under gravitygradient torque,โ Mathematical Problems in Engineering, vol. 2012, Article ID 123138, 15 pages, 2012. View at: Publisher Site  Google Scholar
 O. Winter, D. Mourão, C. de Melo, E. Macau, J. Ferreira, and J. Carvalho, โControlling the eccentricity of polar lunar orbits with lowthrust propulsion,โ Mathematical Problems in Engineering, vol. 2009, Article ID 159287, 10 pages, 2009. View at: Google Scholar  Zentralblatt MATH
 D. Merguizo Sanchez, T. Yokoyama, P. I. De Oliveira Brasil, and R. R. Cordeiro, โSome initial conditions for disposed satellites of the systems GPS and Galileo constellations,โ Mathematical Problems in Engineering, vol. 2009, Article ID 510759, 22 pages, 2009. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 R. A. Broucke, โLongterm thirdbody effects via double averaging,โ Journal of Guidance, Control, and Dynamics, vol. 26, no. 1, pp. 27โ32, 2003. View at: Publisher Site  Google Scholar
 G. Metris and P. Exertier, โSemianalytical theory of the mean orbital motion,โ Astronomy and Astrophysics, vol. 294, pp. 278โ286, 1995. View at: Google Scholar
 A. F. B. A. Prado, โThirdbody perturbation in orbits around natural satellites,โ Journal of Guidance, Control, and Dynamics, vol. 26, no. 1, pp. 33โ40, 2003. View at: Publisher Site  Google Scholar
 W. M. Kaula, โDevelopment of the lunar and solar disturbing functions for a close satellite,โ Astronomical Journal, vol. 67, pp. 300โ303, 1962. View at: Publisher Site  Google Scholar
 G. Giacaglia, โLunar perturbations of artificial satellites of the earth,โ Celestial Mechanics, vol. 9, pp. 239โ267, 1974. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 S. K. Collins and P. J. Cefola, โDouble Averaged Third Body Model for Prediction of SuperSynchronous Orbits over Long Time Spans,โ Paper AAS 79135, American Astronautical Society, 1979. View at: Google Scholar
 M. Lara, J. F. SanJuan, L. M. Lopez, and P. Cefola, โOn thirdbody perturbations on highaltitude orbits,โ submitted to. Celestial Mechanics and Dynamical Astronomy. View at: Publisher Site  Google Scholar
 R. C. Domingos, R. Vilhena de Moraes, and A. F. B. D. A. Prado, โThirdbody perturbation in the case of elliptic orbits for the disturbing body,โ Mathematical Problems in Engineering, vol. 2008, Article ID 763654, 14 pages, 2008. View at: Google Scholar  Zentralblatt MATH
 C. C. Chao and R. A. Gick, โLongterm evolution of navigation satellite orbits: GPS/GLONASS/GALILEO,โ Advances in Space Research, vol. 34, no. 5, pp. 1221โ1226, 2004. View at: Publisher Site  Google Scholar
 G.I. Hori, โTheory of general perturbation with unspecified canonical variable,โ Publications of the Astronomical Society of Japan, vol. 18, no. 4, pp. 287โ296, 1966. View at: Google Scholar
 A. Deprit, โCanonical transformations depending on a small parameter,โ Celestial Mechanics, vol. 1, pp. 12โ30, 1969. View at: Publisher Site  Google Scholar
 J. A. Campbell and W. H. Jefferys, โEquivalence of the perturbation theories of Hori and Deprit,โ Celestial Mechanics, vol. 2, no. 4, pp. 467โ473, 1970. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 A. Deprit, โDelaunay normalisations,โ Celestial Mechanics, vol. 26, no. 1, pp. 9โ21, 1982. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 J. M. A. Danby, Fundamentals of Celestial Mechanics, WillmannBell, Richmond, VA, USA, 1992.
 A. Deprit and A. Rom, โThe main problem of artificial satellite theory for small and moderate eccentricities,โ Celestial Mechanics, vol. 2, no. 2, pp. 166โ206, 1970. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 D. Brouwer and G. M. Clemence, Methods of Celestial Mechanics, Academic Press, New York, NY, USA, 1961.
 J. F. SanJuan, L. M. López, and R. López, โMathATESAT: a symbolicnumeric environment in astrodynamics and celestial mechanics,โ Lecture Notes in Computer Science, vol. 6783, part 2, pp. 436โ449, 2011. View at: Publisher Site  Google Scholar
 Y. Kozai, โSecondorder solution of artificial satellite theory without air drag,โ Astronomical Journal, vol. 67, no. 7, pp. 446โ461, 1962. View at: Publisher Site  Google Scholar
 E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations. I. NonStiff Problems, vol. 8, Springer, Berlin, Germany, 2nd edition, 1993.
Copyright
Copyright © 2012 Martin Lara et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.