`Mathematical Problems in EngineeringVolume 2009 (2009), Article ID 740460, 24 pageshttp://dx.doi.org/10.1155/2009/740460`
Research Article

## Nonsphericity of the Moon and Near Sun-Synchronous Polar Lunar Orbits

2Division of Space Mechanics and Control—INPE, CEP 12227-010, São José dos Campos, SP, Brazil

Received 30 July 2009; Revised 28 October 2009; Accepted 3 December 2009

Copyright © 2009 Jean Paulo dos Santos Carvalho 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.

#### Abstract

Herein, we consider the problem of a lunar artificial satellite perturbed by the nonuniform distribution of mass of the Moon taking into account the oblateness () and the equatorial ellipticity (sectorial term ). Using Lie-Hori method up to the second order short-period terms of the Hamiltonian are eliminated. A study is done for the critical inclination in first and second order of the disturbing potential. Coupling terms due to the nonuniform distribution of mass of the Moon are analyzed. Numerical simulations are presented with the disturbing potential of first and second order is. It an approach for the behavior of the longitude of the ascending node of a near Sun-synchronous polar lunar orbit is presented.

#### 1. Introduction

In this paper, we consider the problem of a lunar artificial satellite of low altitude taking into account the oblateness (J2) and the equatorial ellipticity (sectorial term C22) of the Moon. The Lie-Hori [1] perturbation theory method up to the second order is applied to eliminate the short-period terms of the disturbing potential. The perturbation method up to the second order is applied to analyze coupling terms. In this work, the long-period term of the disturbing potential is analyzed. A formula is developed to compute the critical inclination when the perturbations due to the nonsphericity of the Moon as a function of the terms of the zonal and sectorial harmonics occur.

An approach is done for a special type of orbit, denominated Sun-synchronous orbit of Moon’s artificial satellites. The Sun-synchronous orbit is a particular case of an almost polar orbit. The satellite travels from the North Pole to the South Pole and vice versa, but its orbital plane is always fixed for an observer that is posted in the Sun. Thus the satellite always passes approximately on the same point of the surface of the Moon every day in the same hour. In such a way the satellite can transmit all the data collected for a lunar fixed antenna, during its orbits. An analysis of Sun-synchronous orbits considering the nonuniform distribution of mass of the Moon is done for the longitude of the ascending node with an approach based on Park and Junkins [2].

In [35] an analytical theory, is developed to study the orbital motion of lunar artificial satellites using the method of transformation of Lie [6, 7] as a perturbation method. The main perturbation is due to the nonspherical gravitational field of the Moon and the attraction of the Earth. The disturbing body is in circular orbit with the disturbing function developed in polynomial of Legendre up to the second order. In [811] an analytical theory, is developed with numerical applications taking into account the nonuniform distribution of mass of the Moon and the perturbation of the third body in elliptical orbit (Earth is considered). The disturbing function is expanded in Legendre associated functions up to the fourth order.

This paper is developed based on [2, 4], where the perturbation theory method of Lie-Hori up to the second order is used. Our contribution is characterized by (a) We developed of a new formula for the critical inclination of second order; (b) we fix g and h to assure the condition of frozen orbits; (c) we showed that the coefficients J2 and C22 affect the variation of the eccentricity strongly (it affects the eccentricity directly) in the second order contributing (especially the C22 term) to increasing the variation of the eccentricity mainly for small inclinations; (d) the coupled perturbations (nonuniform distribution of mass of the Moon (J2 and C22) and third-body (P2)) help to control the variation of the eccentricity for low-altitude polar orbits; (e) we presented a new formula to compute inclinations for Sun-synchronous orbits when it is taking into account the harmonic J2 and C22 in the first-order potential.

This paper has seven sections. In Section 2, the terms due to nonsphericity of the Moon are presented while Section 3 is devoted to the Hamiltonian of the system. In Section 4, an approach concerning the critical inclinations is used. In Section 5, an approach concerning the Sun-synchronous lunar orbits is used. Numeric applications are done in Section 6. Section 7 is devoted to the conclusions.

#### 2. Nonsphericity of the Moon

Besides the fact that the Moon is much less flattened than the Earth, it also causes perturbations in space vehicles. Table 1 presents orders of magnitude for some zonal and sectorial harmonics compared with the same parameters for the Earth. The term describes the equatorial bulge of the Moon, often referred to as the oblateness. The coefficient measures the ellipticity of the equator.

Table 1: Magnitude orders for and .

The space vehicle is a point of mass in a three-dimensional orbit with orbital elements: (semimajor axis), (eccentricity), (inclination), (argument of periapsis), (longitude of the ascending node), and (mean motion) given by the third Kepler’s law , where is the mass of the Moon.

Then, we will present the Hamiltonian formalism using Delaunay canonical variables [7] defined as , , mean anomaly, argument of periapsis, and longitude of the ascending node.

The force function, the negative of the total energy as used in physics, is given by

here, V is the negative of the potential energy, and T is the kinetic energy. The force function can be put as [12]:

or

the function , comprising all terms of V except the central term, is known as the disturbing potential. The term due to the unperturbed potential is given by

Considering the lunar equatorial plane as the reference plane, the disturbing potential can be written in the form [13] of

where is the Lunar gravitational constant, R0 is the equatorial radius of the Moon ( km), Pn represent the Legendre polynomial, Pnm represent the associated Legendre polynomial, the angle is the latitude of the orbit with respect to the equator of the Moon, the angle is the longitude measured from the direction of the longest axis of the Moon, where , since is the longitude reckoned from any fixed direction, and is the longitude of the Moon’s longest meridian from the same fixed direction. However, will contain the time explicitly (see [4, 13], for a detailed discussion). Using spherical trigonometry, we have where f is the true anomaly.

The following assumptions have been made [4, 14] and will be used in this work: (a) the motion of the Moon is uniform (librations are neglected); (b) the lunar equator lies in the ecliptic (we neglect the inclination of about 1.5° of the lunar equator to the ecliptic, and the inclination of the lunar orbit to the ecliptic of about 5°); (c) the longitude of the lunar longest meridian is equal to the mean longitude of the Earth (librations are neglected); (d) the mean longitude of the Earth, , is equal to . See De Saedeleer [4] for a detailed discussion.

Since the variables and appear only as a combination of , where , with being the lunar mean motion, the degree of freedom can be reduced by choosing as a new variable . A new term must then be added to the Hamiltonian in order to get . The Hamiltonian is still time-dependent through . Since the longest meridian is always pointing toward the Earth, it is possible to choose a rotating system whose -axis passes through this meridian.

Regarding the Earth’s potential, the dominant coefficient is J2. The rest are of higher-order terms [15]. In contrast to the Earth, the first harmonics of the Lunar potential are all almost of the same order (see Table 1). This fact complicates the choice of the harmonic where the potential can be truncated and this makes its choice a little arbitrary. The influence of the Earth and of the nonsphericity of the Moon on the stability of lunar satellites was also pointed out by [16] but sectorial harmonics were not considered. In terms of the orbital elements, the Legendre associated functions for zonal up to J5 and sectorial terms C22 and C31, where the values for the spherical harmonic coefficients are given in the Appendix C, can be written in the following form [13, 14]:

where

In this paper, are taken into only the terms due to J2 and C22 accout. The zonal perturbation due to the oblateness is defined by [3] , where . However, the disturbing potential is

Substituting the relation , using the Cayley’s tables [17] to express the true anomaly in terms of the mean anomaly, and with some algebraic manipulations, we get

For the sectorial perturbation, we get [13, 14]

where (R0 is the equatorial radius of the Moon;  km).

With some manipulations, we get

where the disturbing potential is written in the form .

#### 3. The Hamiltonian System

We find in the literature several papers that use the method of the average to calculate perturbations of long-period on artificial satellites of the Moon, such as [1823]. However, our objective here is to compute analytically secular and periodic perturbations up to the second order and, using this, to analyze the coupling terms relating the harmonic coefficients. The Lie-Hori [1] perturbation method is applied to eliminate short-period terms of the Hamiltonian.

In [24], a different approach is proposed for the canonical version of Hori method. The reference [24] showed that the ordinary differential equation with an auxiliary parameter as independent variable, introduced through Hori auxiliary system, can be replaced by a partial differential equation in time t.

In what follows, first the Lie-Hori [1] method will be shortly presented and then applied to the problem of the orbital motion of the satellite around the Moon.

Consider the mth order equation of the algorithm of the perturbation method proposed by Hori [1]:

where braces stand for the Poisson brackets, is the undisturbed Hamiltonian, and is a function obtained from the preceding orders, involving , , , and , . All of these functions are written in terms of the new set of canonical variables and defined through the following equations:

where is the original set of canonical variables, is the original Hamiltonian, the new Hamiltonian and is the generating function of the canonical transformation, . The transformation is such that the new canonical system has some advantages for the solution.

In order to determine the functions and , Hori introduces an auxiliary parameter through the following system of canonical equations [1]:

Accordingly, (3.1) reduces to

with written in terms of the general solution of the system (3.3), involving constants of integration. Equation (3.4) has two unknown functions: and .

The Poisson brackets are defined as

with respect to the classical Delaunay variables set l, g, h, L, G, H. Since only l, g, h, L are explicitly present in the Hamiltonian, the partial derivatives with respect to L, G, H are computed as , , = where the bracket indicates the derivative with respect to L occurring explicitly, and c, s, are , , .

To solve (3.4), we separated in secular and periodic part as it is done by the principle of the mean [7]. Using (3.3) and (3.4), Hori [1] writes the equations that define the integration theory based on average principle to determine the new Hamiltonian and the generating function as follows.

Order zero:

Order one:

Order two:

where, at each function, the index s represents the secular part of the function and the index p the periodic part of the function. See Hori [1] for a detailed discussion.

The Hamiltonian of the dynamical system associated to the problem of the orbital motion of the satellite around the Moon can be written in the following form:

where

The term is added to reduce the degree of freedom, since the mean longitude of the Earth is time-dependent [14]. Here, the term is taken as order zero as suggested by Breiter [25].

Now, doing

we get,

With the purpose of applications of the perturbation method, the terms of the Hamiltonian are written in the following form:

The disturbing terms are represented in the first order of the applied method. To eliminate the short-period terms of equation (3.13), the method of Lie-Hori [1] perturbation theory is applied. In this work long-period terms are calculated, substituting the result in the planetary equations of Lagrange [26]. The equations of motion are integrated and finally the results analyzed. With a simplified model for the disturbing potential it is possible to do analyses for the orbital motion of the satellite.

Applying the method of Hori [1] to our problem to eliminate the terms of short period, we get the following:

Order zero:

Order one:

where .

Order two:

#### 4. Critical Inclination

We consider now the problem of a lunar artificial satellite with low altitude taking into account the oblateness (J2) and the equatorial ellipticity (sectorial term C22) of the Moon. The first order long-period disturbing potential (order of the method of perturbation theory) obtained by the Hori method algorithm can be written as

where and .

We observe that, at the order considered, the disturbing potential has the terms due to the oblateness (), that are secular, and terms due to the equatorial ellipticity of the Moon (), that appear multiplied for .

Taking into account (4.1) a formula for the critical inclination is found. In fact, substituting (4.1) in the planetary equations of Lagrange [26] and solving the equation , we get

this formula was already obtained by De Saedeleer and Henrard [5] and was here derived in an independent way, observing that in [5] . Thus, when we consider the terms due to the oblateness (J2) and the equatorial ellipticity of the Moon (C22), the critical inclination depends on the longitude of the ascending node. Figure 1 represents the variation between the inclination and the longitude of the ascending node. Table 2 represents the values of the critical inclination for some values of the ascending node.

Table 2: Critical inclination for the potential of first order, where: ; .
Figure 1: Variation of the critical inclination with respect to the longitude of the ascending node where i-degree and h-rad.

Now, let us consider the second-order disturbing potential where is given in the Appendix A (order of the method of perturbation theory); the potential k2 presents:

(a)coefficients of second order (b)coupling terms between J2 and C22.

Plugging the equations for the potential in the planetary equations of Lagrange and solving the equation , we present a new formula to compute the critical inclination taking into account the J2 and C22 terms of the second-order disturbing potential. It is a function of two variables: the argument of the periapsis (g) and the longitude of the ascending node (h). When the sectorial term C22 is considered, the first order disturbing potential is a function of the longitude of the ascending node and of both longitude of the ascending node and argument of the periapsis to the second order potential. The new formula is given in the Appendix B.

We observe that, at the second order, the disturbing potential presents terms due to the oblateness () and to the equatorial ellipticity of the Moon (), that also appear multiplied by periodic functions. Here, terms of couplings between the oblateness and the equatorial ellipticity of the Moon (J2, C22) and terms of second order appear. Several scenarios can be considered. For instance frizzing orbits with particular values of h and g, let us say and , we get a value of 53.46° for the critical inclination taking into account equation given in the Appendix B. Therefore, the critical inclination taking into account equation (4.2) is 58.56° (see Table 2).

#### 5. Sun-Synchronous Lunar Orbit

Now, an approach is presented for a Lunar Sun-synchronous orbit. The Moon rotates with angular rate about 360° for 27, 32 days, while the Earth rotates with angular rate about 360° by day. The perturbation caused by the orbital precession has been studied historically for orbits centered in the Earth because of near polar orbits the precession is about one degree per day and to provide attractive Sun-synchronous orbits for many missions around the Earth. Considering Sun-synchronous orbits for lunar satellites we show that it is not possible to produce near polar orbits. The precession of the ascending node due to the nonsphericity of the Moon, when only the effect of the J2 is considered, in (3.3) is well known in Brouwer theory [27] that the precession of the longitude of the ascending node is given by

The Moon’s orbital period is about 27, 32 days and the Earth’s orbital period is about 365, 26 days. Then, for a Sun-synchronous orbit, we have, in lunar day [2]:

An inclination for a Sun-synchronous orbit was presented by Park and Junkins [2] using (5.1) and (5.2), with the following initial conditions:  km; . The calculated inclination is . This inclination is not feasible for producing near-polar orbits, because this orbit does not pass sufficiently near the poles. Considering the disturbing potential given by (4.1) and substituting in the Lagrange planetary equations [26] to calculate the variation of the longitude of the ascending node, we get

where , . This new equation (5.3) gives the precession of the longitude of the ascending node due to nonsphericity of the Moon when the effect of the J2 and the C22 are considered. In first approximation the periodic terms due to the are negligible, however, for the C22 term in the first approximation appears the periodic term . When , we obtain the classic solution given by (5.1).

In what follows, the variation of the longitude of the ascending node will be analyzed to obtain polar or near-polar orbits for some special cases. Using (5.2) and (5.3) we get a new formula to compute the inclination of Sun-synchronous orbits for Moon’s satellites of low altitude:

where n is given in rad/s, ; . Equation (5.4) gives the inclination depending on the semi-major axis, the eccentricity and on the longitude of the ascending node. For Sun-synchronous orbits, considering  km; ; the calculated inclination is . This inclination is not also ideal for near-polar orbits. Thus, the obtained results are still distant from a polar Sun-synchronous orbit, but it is important to consider the term due to the Moon’s equatorial ellipticity to get more realistic results. Figure 2 represents the variation of Sun-synchronous inclinations with respect to the longitude of the ascending node where  km and . It can be observed in Figure 2 that, fixing , we get an inclination of about 132°.

Figure 2: Variation of the inclination Sun-synchronous (IS in degree) with longitude of the ascending node (h-rad),  km, .

#### 6. Applications for Low Altitude Satellites

The disturbing potential (first order (k1) and second order (k2)) is substituted in the Lagrange’s planetary equations [26] and numerically integrated. Considering the disturbing potential due to the nonsphericity of the Moon (J2 and C22), numerical applications (long-period potential) are performed to analyze the variation of the eccentricity for different initial conditions of the inclination.

Figures 3 and 4 represent the comparison for different orders of the disturbing potential. The variation of the eccentricity for lunar satellites in low altitude is constant at first order. This happens because the coefficients J2 and C22 do not affect the variation rate of the eccentricity (as we can verify in Lagrange’s planetary equations). Therefore, it is important to insert more terms in the potential to get more realistic results as, for example, to study the lifetimes of low altitude Moon artificial satellites [28], considering the zonal terms J2, J3 and J5 and the sectorial terms C22 and C31 in the disturbing potential. When taking into account the disturbing potential of second order considering the effect of the nonuniform distribution of mass of the Moon (J2 and C22), the coefficients J2 and C22 affect the variation rate of the eccentricity (as we can verify in the Lagrange’s planetary equations). For the second order the results shows a small variation of the eccentricity for larger inclinations and an accentuated increase for small inclinations. This is due to the C22 term that affects the eccentricity of the satellite directly in second order.

Figure 3: Initial conditions:  km, and t-days.
Figure 4: Initial conditions:  km, and t-days.

The expression of the eccentricity is presented in the following form where the terms due to the oblateness () and the equatorial ellipticity () appear multiplied by periodic functions, terms of couplings between J2 and C22 and terms of second order of the type and . Figures 3 and 4 shows the temporal variation of the eccentricity for cases where initial conditions are obtained from the frozen orbits condition. For instance frizzing orbits with particular values of h and g, let us say and . Another factor that also contributes for the variation of the eccentricity, using the potential up to the second order, is the presence of the coefficients of second order terms () and the coupling terms between J2 and C22. The choice of the initial inclination is very important to assure a frozen orbit when it is taken into account the second order disturbing potential.

Figures 5 and 6 show the inclination suffering a periodic variation that depends of the longitude of the ascending node, as we can verify by (4.2). Figure 5 shows a variation of the inclination for a value below of the critical inclination around 8 degrees in 100 days, while Figure 6 presents a more accentuated variation for values of the inclination above of the critical inclination around 50 degrees in 200 days.

Figure 5: Initial conditions:  km, and t-days. Fixing a value we find a value of 58.56° for the critical inclination.
Figure 6: Initial conditions:  km, and t-days. Fixing a value we find a value of 58.56° for the critical inclination.

Considering the second-order disturbing potential, Figures 7 and 8 represent the variation rate for the inclination. The variation represented is due to the initial values given for the argument of the periapsis and for the longitude of the ascending node. The given initial values for g and h are for the condition of frozen orbits. For instance frizzing orbits with particular values of h and g, let us say and . The same comments done for the eccentricity, including those for the coupling terms (zonal and sectorial), are valid for the variation of the orbital inclination.

Figure 7: Initial conditions:  km, and t-days.
Figure 8: Initial conditions:  km, and t-days.

Numerical applications with the first and second order disturbing potential are done taking into account the nonsphericity of the Moon and perturbations from the third-body in elliptical orbit (Earth is considered) considering the term P2 of the Legendre polynomial and the eccentricity of the disturbing body up to the second order. Figures 9, 10 and 11 represent the sum of the potential of the first order with the second order (). Figure 9 shows that, if and C22 = 0 the small inclinations cause small oscillations in the variation of the eccentricity and Figure 10 shows the effect caused by the C22 term, when , where the small inclinations cause a large increase in the variation of the eccentricity, as well as we can visualize in Figure 11, where it is taken into account the , and terms. Considering the disturbing potential up to second order, the terms J2 and C22 that does not cause perturbation in the eccentricity in first order, appear as disturbing term of the type and and terms of coupling of the type J2C22 that affects the eccentricity of the satellite directly. The terms appear due to the perturbation method used.

Figure 9: Initial conditions:  km, Considered Perturbations:
Figure 10: Initial conditions:  km, Considered Perturbations:
Figure 11: Initial conditions:  km, Considered Perturbations:

The temporal variation of the eccentricity is strongly affected by the initial inclination (i0). As it can be observed by Figure 11, for the variation of the eccentricity presents great amplitude but, for , the variation has small amplitude.

Figures 12 and 13 represent the variation of the eccentricity for a lunar satellite of low altitude considering different terms in the disturbing potential. Figure 12 considered third body perturbation and Figure 13 considered nonsphericity of the Moon and third body (P2) perturbation. A comparison is done between the perturbations for the case where the inclination is 90° (polar orbit). Therefore we can conclude that, besides the term due to the J2, the sectorial term C22 should also be considered in the disturbing potential to get more realistic results. For lunar satellites of low altitude it is impracticable to consider real applications taking into account only the perturbation of the third body (P2). In fact, taking into account only the perturbation of the third body (P2), the eccentricity of the satellite increases causing escape from the Moon, or crashing to the Moon in 600 days (see Figure 12). We observe that, for a Moon’s artificial satellite orbiting in low altitude, the combination of the two perturbations help to control the variation of the eccentricity, mainly for inclinations of 90° (see Figure 13).

Figure 12: Initial conditions:  km, , and t-days. Perturbations
Figure 13: Initial conditions:  km, , and t-days. Perturbations

These results agree with those obtained by d’Avanzo et al. [29] when it is considered just the effects of the zonal harmonics and a first order potential.

#### 7. Conclusions

Using Lie-Hori method, the disturbing potential due to the nonsphericity of the gravitational field of the Moon is obtained up to first and second order. The disturbing potential is substituted in Lagrange’s planetary equations and they are numerically integrated. Analyses for the variations of the orbital elements are done. Terms of couplings between the oblateness and the equatorial ellipticity of the Moon (J2, C22) and terms of second order of type and are obtained. A formula is developed to compute the critical inclination when the effect of the C22 (equatorial ellipticity) term is considered in the Hamiltonian in first and second order. The critical inclination can be strongly affected by the coefficient due to the equatorial ellipticity of the Moon and by the longitude of the ascending node. The formula for the critical inclination for the second order is a function of two variables: the argument of the periapsis and the longitude of the ascending node. At the first order this formula is a function of the longitude of the ascending node only. For Lunar low altitude satellites (LLAS), it is important to take into account both, the terms due to the oblateness and terms with the equatorial ellipticity of the Moon to get more realistic results.

The variation of the longitude of the ascending node is analyzed looking Lunar Sun-synchronous orbits and near-polar orbits. A new formula is obtained to compute inclinations of Lunar Sun-synchronous orbits when the terms due to the oblateness of the Moon (J2 and C22) are taken into account. The presented formula to the inclination depends on the semi-major axis, eccentricity and on the longitude of the ascending node of the satellite. The term due to the effect of the C22 must be considered for the case of a lunar satellite to analyze the precession of the longitude of the ascending node.

For a LLAS, when it is considered only the nonsphericity of the Moon in the disturbing potential, and at the first order, the orbital eccentricity of the satellite is constant along the time. This happens since the coefficients and , at this order, do not affect the variation rate of the eccentricity directly, therefore it is important to insert more terms in the potential as, for example, the zonal terms J3, J5 and the sectorial term C31, to get more realistic results. At the second order, small variations are present. In this case, the small variations of the eccentricity are due to a combination of the following factors: (a) initial conditions (given to get frozen orbits) (b) to couplings terms between the oblateness and the equatorial ellipticity of the Moon (J2, C22) and (c) terms of second order of type and . For small inclinations, second-order terms (including coupled terms) are greater than 1st order terms. This happens because the coefficients J2 and C22, at second order, affect the variation of the eccentricity directly.

To study lifetime of LLAS, due to the characteristics of the mass distribution of the Moon, it is necessary to take into account up the second order of the disturbing potential and develop up to the second order the Lie-Hori algorithm. In fact, at first order, the coefficients do not affect the eccentricity directly while at the second order the coefficients J2 and C22 affects the eccentricity directly and thus contributing efficiently to more complete studies.

#### A. Appendix A

Let us consider the second-order disturbing potential (order of the method of perturbation theory), where :

#### B. Appendix B

Plugging the equations for the potential in the planetary equations of Lagrange and solving the equation , we present a new formula to compute the critical inclination taking into account the and terms of the second-order disturbing potential. We get

#### C. Appendix C

Model given by Akim and Golikov [30] for the spherical harmonic coefficients is

#### Acknowledgment

This work was accomplished with support of the FAPESP under the contract no. 2007/04413-7 and 2006/04997-6, SP-Brazil, and CNPQ (300952/2008-2).

#### References

1. G. I. Hori, “Theory of general perturbations with unspecified canonical variables,” Publications of the Astronomical Society of Japan, vol. 18, no. 4, pp. 287–296, 1966.
2. S.-Y. Park and J. L. Junkins, “Orbital mission analysis for a lunar mapping satellite,” Journal of the Astronautical Sciences, vol. 43, no. 2, pp. 207–217, 1995.
3. B. De Saedeleer and J. Henrard, “Orbit of a lunar artificial satellite: analytical theory of perturbations,” in Transits of Venus: New Views of the Solar System and Galaxy, IAU Colloquium, no. 196, pp. 254–262, Cambridge University Press, Cambridge, Mass, USA, 2005.
4. B. De Saedeleer, “Analytical theory of a lunar artificial satellite with third body perturbations,” Celestial Mechanics & Dynamical Astronomy, vol. 95, no. 1–4, pp. 407–423, 2006.
5. B. De Saedeleer and J. Henrard, “The combined effect of ${J}_{2}$ and ${C}_{22}$ on the critical inclination of a lunar orbiter,” Advances in Space Research, vol. 37, no. 1, pp. 80–87, 2006.
6. A. Deprit, “Canonical transformations depending on a small parameter,” Celestial Mechanics, vol. 1, pp. 12–30, 1969.
7. S. Ferraz Mello, Canonical Perturbation Theories, Springer, New York, NY, USA, 2007.
8. J. P. dos Santos Carvalho, R. V. de Moraes, and A. F. B. de Almeida Prado, “Semi-analytic theory of a moon artificial satellite considering lunar oblateness and perturbations due to a third-body in elliptic orbit,” in Proceedings of the 7th Brazilian Conference on Dynamics, Control and Applications, pp. 51–57, Presidente Prudente, Brazil, 2008.
9. J. P. dos Santos Carvalho, R. V. de Moraes, and A. F. B. de Almeida Prado, “Moon artificial satellites: lunar oblateness and earth perturbations,” in Proceedings of the International Conference on Mathematical Problems in Engineering, Aerospace and Sciences (ICNPAA '08), pp. 1095–1106, Cambrige Scientific, Genoa, Italy, 2009.
10. J. P. dos Santos Carvalho, R. V. de Moraes, and A. F. B. de Almeida Prado, “A study on resonance for a lunar artificial satellite,” in Proceedings of the 8th Brazilian Conference on Dynamics, Control and Applications, pp. 1–6, Bauru, Brazil, 2009.
11. J. P. dos Santos Carvalho, R. V. de Moraes, and A. F. B. de Almeida Prado, “Non-sphericity of the moon and critical inclination,” in Proceedings of the 32nd Congresso Nacional de Matemática Aplicada e Computacional, Cuiabá, Brazil, 2009.
12. W. M. Kaula, Theory of Satellite Geodesy: Applications of Satellites to Geodesy, Dover, New York, NY, USA, 2000.
13. M. Radwan, “Analytical approach to the motion of a lunar artificial satellite,” Astrophysics and Space Science, vol. 283, no. 2, pp. 133–150, 2003.
14. G. E. O. Giacaglia, J. P. Murphy, and T. L. Felsentreger, “A semi-analytic theory for the motion of a lunar satellite,” Celestial Mechanics, vol. 3, no. 1, pp. 3–66, 1970.
15. A. Elipe and M. Lara, “Frozen orbits about the moon,” Journal of Guidance, Control, and Dynamics, vol. 26, no. 2, pp. 238–243, 2003.
16. Z. Knežević and A. Milani, “Orbit maintenance of a lunar polar orbiter,” Planetary and Space Science, vol. 46, no. 11-12, pp. 1605–1611, 1998.
17. A. Cayley, “Tables of the developments of functions in the theory of elliptic motion,” Memoirs of the Royal Astronomicai Society, vol. 29, p. 191, 1861.
18. 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.
19. I. V. Costa and A. F. B. de Almeida Prado, “Orbital Evolution of a Satellite Perturbed by a Third Body,” in Advances in Space Dynamics, A. F. Bertachini and A. F. B. de Almeida Prado, Eds., p. 176, 2000.
20. R. C. Domingos, R. V. de Moraes, and A. F. B. de Almeida Prado, “Third-body perturbation in the case of elliptic orbits for the disturbing body,” Advances in the Astronautical Sciences, vol. 129, part 3, pp. 2639–2652, 2008.
21. A. F. B. de Almeida Prado, “Third-body perturbation in orbits around natural satellites,” Journal of Guidance, Control, and Dynamics, vol. 26, no. 1, pp. 33–40, 2003.
22. C. R. H. Solórzano, Third-body perturbation using a single averaged model, M.S. dissertation, National Institute for Space Research (INPE), São José dos Campos, Brazil, 2002.
23. C. R. H. Solórzano and A. F. B. de Almeida Prado, “Third-body perturbation using a single averaged model,” in Advances in Space Dynamics 4: Celestial Mechanics and Astronautics, H. K. Kuga, Ed., pp. 47–56, 2004.
24. S. Da Silva Fernandes, “Notes on Hori method for canonical systems,” Celestial Mechanics & Dynamical Astronomy, vol. 85, no. 1, pp. 67–77, 2003.
25. S. Breiter, “Second-order solution for the zonal problem of satellite theory,” Celestial Mechanics and Dynamical Astronomy, vol. 67, pp. 237–249, 1997.
26. J. Kovalevsky, Introduction to Celestial Mechanics, Bureau dês Longitudes, Paris, France, 1967.
27. D. Brouwer, “Solution of the problem of artificial satellite theory without drag,” The Astronomical Journal, vol. 64, pp. 378–397, 1959.
28. W. K. Meyer, J. J. Buglia, and P. N. Desai, “Lifetimes of Lunar Satellite Orbits,” Tech. Rep. N-TP-3394 94, p. 27771, NASA STI/Recon, 1994.
29. P. D'Avanzo, P. Teofilatto, and C. Ulivieri, “Long-term effects on lunar orbiter,” Acta Astronautica, vol. 40, no. 1, pp. 13–20, 1997.
30. E. Akim and A. Golikov, “Combined model of the lunar gravity field,” in Proceedings of the 12th International Symposium on Space Flight Dynamics (ESOC '97), pp. 357–360, Darmstadt, Germany, June 1997.