- About this Journal ·
- Abstracting and Indexing ·
- Advance Access ·
- Aims and Scope ·
- Annual Issues ·
- Article Processing Charges ·
- Articles in Press ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents
Mathematical Problems in Engineering
Volume 2012 (2012), Article ID 659396, 17 pages
Semianalytic Integration of High-Altitude Orbits under Lunisolar Effects
1C/Columnas de Hércules 1, ES-11100 San Fernando, Spain
2Departamento de Matemáticas y Computación, Universidad de La Rioja, 26004 Logroño, Spain
3Departamento de Ingeniería Mecánica, Universidad de La Rioja, 26004 Logroño, Spain
Received 6 November 2011; Accepted 20 January 2012
Academic Editor: Josep Masdemont
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.
The long-term effect of lunisolar perturbations on high-altitude 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 high-altitude 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 third-body disturbing function up to the second degree. Using canonical perturbation theory, the averaging is carried out up to the order where second-order 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 low-eccentricity orbits.
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  but most probably by the impact with other uncontrolled man-made 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 high-degree and order geopotential, ephemerides-based 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 . In fact, both approaches, numerical and semi-analytical, do not need to enter a competition. Thus, for instance, while semi-analytical 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 .
In a semi-analytical 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 short-period 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 third-body perturbations, in a reduced-phase space . In this respect, much attention has been recently paid to the long-term evolution of classical GNSS constellations, either for operational or disposal orbits .
While the noncentralities of the Earth gravitational potential play a key role in the motion of low altitude satellites, third-body perturbations have also a decisive influence in the long-term evolution of medium- and high-altitude Earth orbits. The third-body disturbing function is commonly given by a series expansion in Legendre polynomials. Often, the series is truncated to the first term in the expansion , 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 second-order zonal coefficient () clearly dominates over all other harmonics, second-order effects of may be important and must be taken into account when the effect of higher-order harmonics is studied. Correspondingly, the truncation in the expansion of third-body perturbations must include terms of magnitude comparable to -squared. Because the third-body disturbing function is expanded in the ratio semi-major axis of the satellite’s orbit to semi-major axis of the third-body’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 high-altitude orbits about a noncentral Earth, which is assumed to be oblate although without equatorial symmetry. More specifically, we are interested in the semi-analytical 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 , where we approached the general case of third-body perturbations via double averaging, we release here the common simplification of assuming the third-body 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  seem to contradict the claim that taking the third-body in circular orbit does not produce any noticeable degradation in the long-term propagation of real Earth orbits . Besides, for the actual values of the orbits of the sun and moon, neglecting terms in the eccentricity is not consistent with a higher-order 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 semi-major 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 third-body potential, and hence we truncate the moon’s third-body 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 two-degrees-of-freedom 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 ; 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 long-term 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 short-period effects by means of the analytical transformation equations of the averaging provides a quite reasonable precision in the long-term integrations.
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 point-mass 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 mass-point approximation, the third-body disturbing potential , , is where is the third-body’s gravitational parameter, and are the radius vector of the satellite and of the third-body, 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 third-body in its orbit of semi-major 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 third-body’s orbit.
In our model we assume an Earth-centered 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 . Nevertheless, it will be enough for our purposes of investigating the qualitative features that lunisolar perturbations introduce in the long-term 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 high-altitude orbits in the range of 20 000 to 35 000 km, the so-called 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 fourth-order 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 second-order effects of . For the sun, it is enough to take the first term in the Legendre polynomials expansion. Then, the zero-order 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 .
We note that the sun and moon coordinates are assumed to be known functions of time. Therefore, (3.1) is a time-dependent 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, ). Since the time-dependent 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 short-period 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 third-body 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  in the solution of the fourth-order 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 two-degrees-of-freedom, time-dependent 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 eighth-order theory, we neglect higher-order 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 long-period terms related to the sun’s apparent motion and very long-period 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 higher-order Runge-Kutta method. Specifically, the numerical integration was performed with the Dormand and Prince implementation of the Runge-Kutta method coded in FORTRAN 77 by Hairer et al. .
To illustrate the significance of recovering the short-period effects up to higher orders, we first show in Figure 1 a sequence of the errors obtained in the semimajor axis after one month semi-analytical 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 double-averaged 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 long-term 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 short-period 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 short-period errors related to the orbital period of the satellite, we can appreciate a two-week 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 short-period effects can be ignored in the study of the long-term orbital behavior, where the simple propagation of the double-averaged 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 short-term 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.
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 long-term 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 lower-order truncations of the theory, resulting in faster numerical integration of the mean elements. So the fifth-order truncation should be the preferred long-term Hamiltonian. Besides, we also checked for a variety of orbits that making does not either introduce qualitative differences in the long-term, in agreement with .
Modeling lunisolar perturbations on high-altitude Earth orbits requires to retain high degrees in the Legendre polynomials expansion of the third-body disturbing function of the moon. In consequence the eccentricity of the third-body’s orbit cannot be neglected in the case of either the moon or the sun.
The long-term behavior of high-altitude 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 second-order 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 long-term 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 almost-secular growing of the orbit eccentricity and also by very-long-period 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 long-term evolution of the eccentricity.
Part of this research has been supported by the Government of Spain (Projects AYA 2009-11896, AYA 2010-18796, and Grant Gobierno de La Rioja Fomenta 2010/16).
- T. S. Kelso, “Analysis of the Iridum 33—Cosmos 2251 Colission,” Paper AAS 09-368, American Astronautical Society, 2009.
- O. Montenbruck and E. Gill, Satellite Orbits: Models, Methods, and Applications, Springer, Berlin, Germany, 2000.
- D. A. Vallado, Fundamentals of Astrodynamics and Applications, Microcosm Press and Kluwer Academic, El Segundo, Calif, USA, 2007.
- P. J. Cefola, V. Yurasov, Z. Folcik, E. Phelps, R. J. Proulx, and A. Nazarenko, “Comparison of the DSST and the USM Semi-Analytical Orbit Propagators,” Paper AAS 03-236, American Astronautical Society, 2003.
- 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.
- W. D. McClain, A Recursively Formulated First-Order 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/TR-77/6010, 1977.
- S. Valk, A. Lemaître, and F. Deleflie, “Semi-analytical theory of mean orbital motion for geosynchronous space debris under gravitational influence,” Advances in Space Research, vol. 43, no. 7, pp. 1070–1082, 2009.
- J. F. San-Juan, L. M. López, and R. López, “Higher-order analytical attitude propagation of an oblate rigid body under gravity-gradient torque,” Mathematical Problems in Engineering, vol. 2012, Article ID 123138, 15 pages, 2012.
- O. Winter, D. Mourão, C. de Melo, E. Macau, J. Ferreira, and J. Carvalho, “Controlling the eccentricity of polar lunar orbits with low-thrust propulsion,” Mathematical Problems in Engineering, vol. 2009, Article ID 159287, 10 pages, 2009.
- 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.
- R. A. Broucke, “Long-term third-body effects via double averaging,” Journal of Guidance, Control, and Dynamics, vol. 26, no. 1, pp. 27–32, 2003.
- G. Metris and P. Exertier, “Semi-analytical theory of the mean orbital motion,” Astronomy and Astrophysics, vol. 294, pp. 278–286, 1995.
- A. F. B. A. Prado, “Third-body perturbation in orbits around natural satellites,” Journal of Guidance, Control, and Dynamics, vol. 26, no. 1, pp. 33–40, 2003.
- W. M. Kaula, “Development of the lunar and solar disturbing functions for a close satellite,” Astronomical Journal, vol. 67, pp. 300–303, 1962.
- G. Giacaglia, “Lunar perturbations of artificial satellites of the earth,” Celestial Mechanics, vol. 9, pp. 239–267, 1974.
- S. K. Collins and P. J. Cefola, “Double Averaged Third Body Model for Prediction of Super-Synchronous Orbits over Long Time Spans,” Paper AAS 79-135, American Astronautical Society, 1979.
- M. Lara, J. F. San-Juan, L. M. Lopez, and P. Cefola, “On third-body perturbations on high-altitude orbits,” submitted to. Celestial Mechanics and Dynamical Astronomy.
- R. C. Domingos, R. Vilhena de Moraes, and A. F. B. D. A. Prado, “Third-body perturbation in the case of elliptic orbits for the disturbing body,” Mathematical Problems in Engineering, vol. 2008, Article ID 763654, 14 pages, 2008.
- C. C. Chao and R. A. Gick, “Long-term evolution of navigation satellite orbits: GPS/GLONASS/GALILEO,” Advances in Space Research, vol. 34, no. 5, pp. 1221–1226, 2004.
- 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.
- A. Deprit, “Canonical transformations depending on a small parameter,” Celestial Mechanics, vol. 1, pp. 12–30, 1969.
- 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.
- A. Deprit, “Delaunay normalisations,” Celestial Mechanics, vol. 26, no. 1, pp. 9–21, 1982.
- J. M. A. Danby, Fundamentals of Celestial Mechanics, Willmann-Bell, 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.
- D. Brouwer and G. M. Clemence, Methods of Celestial Mechanics, Academic Press, New York, NY, USA, 1961.
- J. F. San-Juan, L. M. López, and R. López, “MathATESAT: a symbolic-numeric environment in astrodynamics and celestial mechanics,” Lecture Notes in Computer Science, vol. 6783, part 2, pp. 436–449, 2011.
- Y. Kozai, “Second-order solution of artificial satellite theory without air drag,” Astronomical Journal, vol. 67, no. 7, pp. 446–461, 1962.
- E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations. I. Non-Stiff Problems, vol. 8, Springer, Berlin, Germany, 2nd edition, 1993.