Planetary Satellite Orbiters: Applications for the Moon

Abstract

Low-altitude, near-polar orbits are very desirable as science orbits for missions to planetary satellites, such as the Earth's Moon. In this paper, we present an analytical theory with numerical simulations to study the orbital motion of lunar low-altitude artificial satellite. We consider the problem of an artificial satellite perturbed by the nonuniform distribution of the mass of the Moon (𝐽₂–𝐽₅, 𝐽₇, and 𝐶₂₂). The conditions to get frozen orbits are presented. Using an approach that considers the single-averaged problem, we found families of periodic orbits for the problem of an orbiter travelling around the Moon, where frozen orbits valid for long periods of time are found. A comparison between the models for the zonal and tesseral harmonics coefficients is presented.

1. Introduction

Frozen orbits around the Moon, natural satellites, or asteroids are of current interest because several space missions have the goal of orbiting around such bodies (see, for instance, [1–5] and references therein). The dynamics of orbits around a planetary satellite, taking into account the gravitational attraction of a third-body and the nonuniform distribution of mass (𝐽₂, 𝐽₃, and 𝐶₂₂) of the planetary satellite (Europa), is studied in [6–11].

This paper contains a continuation of the work developed in [1, 12, 13], where the problem of the third-body gravitational perturbation and the non-uniform distribution of mass (𝐽₂, 𝐶₂₂) of the planetary satellite (Moon) on a spacecraft is studied using a simplified dynamical model. Also in [14, 15], the problem of critical inclination orbits for lunar artificial satellites including the terms (𝐽₂) zonal, (𝐶₂₂) sectorial, and the lunar rotation due to the lunar potential, besides the Keplerian term is considered. In [16], a study done to control the increasing eccentricity of a polar lunar satellite for different altitudes due to the perturbations of the Earth is presented. The approach presented here is the use of low-thrust propulsion in order to keep the orbital eccentricity of the satellite at low values. A fourth-order analytical theory is presented for the accurate computation of quasiperiodic frozen orbits in [17].

In our research, we use a simplified dynamical model that considers the effects caused by the nonsphericity of the Moon (𝐽₂–𝐽₅, 𝐽₇, and 𝐶₂₂) using single-averaged analytical model over the short satellite period. The motion of equations is obtained from the Lagrange planetary equations, and then they are integrated numerically, using the software Maple. Emphasis is given to the case of frozen orbits, defined as orbits which have constant values of eccentricity, inclination, and argument of the periapsis when the single-averaged system is considered. An approach is used for the case of frozen orbits similar to the one used in [1, 2, 5]. The basic idea is to eliminate the terms due to the short time periodic motion of the spacecraft to show smooth curves for the evolution of the mean orbital elements for a long time period. The interesting set of results allow us to find good conditions close to frozen orbits.

In other words (see [5, 11, 18, 19]), a better understanding of the physical phenomenon can be obtained and it allows the study of long-term stability of the orbits in the presence of disturbances that cause low changes in the orbital elements.

The paper has four sections. Section 2 is devoted to the terms due to the nonsphericity of the Moon and to the Hamiltonian of the system. Applications taking into account the long-period disturbing potential using the averaging method are presented in Section 3. Section 4 shows the conclusions.

2. Oblateness of the Moon

To analyze the motion of a lunar low-altitude artificial satellite, it is necessary to take into account the Moon's oblateness [1, 14, 20]. Besides that, the Moon is much less flattened than the Earth but it also causes perturbations in space vehicles. Table 1 shows the orders of magnitude for the zonal and sectorial harmonics for both celestial bodies to make a comparison. The 𝐶₂₀ term describes the equatorial bulge of the Moon, often referred to as oblateness. The coefficient 𝐶₂₂ measures the ellipticity of the equator.

Table 1 

Orders of magnitude for 𝐽₂ and 𝐶₂₂ [14].

The force function, the negative of the total energy as used in physics, is given by 𝐹=𝐾=𝑉−𝑇 here, 𝑉 is the opposite of the potential energy, and 𝑇 is the kinetic energy. The force function can be written as 𝐾=𝜇/𝑟+𝑅−𝑇or 𝐾=𝜇²/2𝑎+𝑅. The function 𝑅, comprising all terms of 𝑉 except the central term, is known as the disturbing potential. The space vehicle is a point of mass in a three-dimensional orbit with orbital elements 𝑎 (semimajor axis), 𝑒 (eccentricity), 𝑖 (inclination), 𝑔 (argument of the periapsis), ℎ (longitude of the ascending node), and 𝑛 (mean motion) given by the third Kepler law 𝑛²𝑎³=𝐺𝑚₀, where 𝐺 is the gravitational constant and 𝑟 is the radius vector of the body 𝑚₀ (mass of the planetary satellite). The term due to the unperturbed potential is given by𝐻₀=𝜇²2𝑎.(2.1)

Considering the lunar equatorial plane as the reference plane, the disturbing potential can be written in the form [21]𝑉_𝑀=−𝜇𝑟⁵_𝑛=2𝑅𝑟^𝑀^𝑛𝐽_𝑛𝑃_𝑛(sin𝜙)−𝑅^𝑀𝑟²𝐶₂₂𝑃₂₂(sin𝜙)cos2𝜆,(2.2) where 𝜇 is the lunar gravitational constant, 𝑅_𝑀 is the equatorial radius of the Moon (𝑅_𝑀=1738 km), 𝑃_𝑛 represent the Legendre polynomials, 𝑃_𝑛𝑚 represent the associated Legendre polynomials, the angle 𝜙 is the latitude of the orbit with respect to the equator of the Moon, and the angle 𝜆 is the longitude measured from the direction of the longest axis of the Moon, where 𝜆=𝜆′−𝜆₂₂, since 𝜆′ is the longitude measured from any fixed direction, and 𝜆₂₂ is the longitude of the Moon's longest meridian from the same fixed direction. However, 𝜆₂₂ contains explicitly the time [20]. Using spherical trigonometry, we have sin(𝜙)=sin(𝑖)sin(𝑓+𝑔), where 𝑓 is the true anomaly of the satellite.

The following assumptions 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 with respect to the ecliptic and the inclination of the lunar orbit with respect 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), and (d) the mean longitude of the Earth 𝜆_⊕ is equal to 𝜆₂₂.

Since the variables Ω and 𝜆_⊕ appear only as a combination of Ω−𝜆_⊕, where 𝜆_⊕=𝑛_𝑀𝑡+const, 𝑛_𝑀 being the lunar mean motion. The degrees 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 ℎ̇=Ω̇−𝑛_𝑀=−𝜕(𝐾+𝑛_𝑀𝐻)/𝜕𝐻. This term 𝑛_𝑀𝐻 is used here in the same way it is used in [20]. 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 𝐽₂. The rest are of higher order [2]. Opposite to the situation of the Earth, the first harmonics of the lunar potential are almost of the same order (see Table 1). This fact complicates the choice of the harmonic where the potential can be truncated, and this fact makes its choice a little arbitrary. The influences of the Earth and of the nonsphericity of the Moon on the stability of lunar satellites were also shown by [22], but the sectorial harmonics were not considered.

In terms of the orbital elements, the Legendre-associated functions for the zonal up to 𝐽₅ and sectorial (𝐶₂₂) terms can be written in the form [20, 21]𝑃₂(sin∅)=123𝑠²𝑃sin²(𝑓+𝑔)−1₃,(sin∅)=52s³sin³(𝑓+𝑔)−3𝑃2𝑠sin(𝑓+𝑔),₄(sin∅)=358𝑠⁴sin⁴(𝑓+𝑔)−354𝑠²sin²(𝑓+𝑔)+3𝑃₅8,(sin∅)=638𝑠⁵sin⁵(𝑓+𝑔)−354𝑠²sin³(𝑓+𝑔)+15𝑃8𝑠sin(𝑓+𝑔),₂₂(sin∅)cos2𝜆=6𝜉²cos²𝑓+𝑥²sin²𝑓+𝜉𝑥sin2𝑓−31−𝑠²sin²(𝑓+𝑔),(2.3) where we use the shortcuts 𝜉=cos𝑔cosℎ−𝑐sin𝑔sinℎ, 𝜒=−sin𝑔cosℎ−𝑐cos𝑔sinℎ, 𝑠=sin𝑖, and 𝑐=cos𝑖.

With this, the zonal perturbation due to the oblateness is [14] 𝐻₂₀=𝜖𝜇4𝑟³1−3𝑐²−3𝑠²cos(2𝑓+2𝑔),with𝜖=𝐽₂𝑅^2𝑀.(2.4)

Replacing the relation 𝜇=𝑛²𝑎³ and using Cayley's tables [23] to express the true anomaly in terms of the mean anomaly, after some algebraic manipulations, we get𝐻₂₀=3𝑛²8𝜖5𝑒²+17𝑒₂−2𝑐²+2−5𝑒²cos(2𝑔+3𝑙)+7𝑒𝑠²cos(2𝑔+3𝑙)+9𝑐₂𝑠²cos(2𝑔+4𝑙)−𝑒𝑠²cos(2𝑔+𝑙)−1392+23𝑒cos𝑙+31𝑒²+𝑒²cos2𝑙.(2.5)

For the sectorial perturbation [20, 21], we get 𝐻₂₂=𝛿3𝜇𝑟³2𝜉²cos²𝑓+2𝜒²sin²𝑓+2𝜉𝜒sin2𝑓−1+𝑠²sin²(𝑓+𝑔),(2.6)

where 𝛿=𝐶₂₂𝑅^2𝑀. With some manipulations, we get𝐻₂₂=−4516𝛿𝑛²𝑒²−25(𝑐−1)²−1cos(2𝑙+2𝑔−2ℎ)3(𝑐+1)²𝑒²−25cos(2𝑙−2𝑔−2ℎ)+(𝑐+1)²𝑒²−2−15cos(2𝑙+2𝑔+2ℎ)3𝑒²−25(𝑐−1)²cos(2𝑙−2𝑔+2ℎ)−175𝑒²−7(𝑐−1)²cos(4𝑙+2𝑔−2ℎ)5𝑒(𝑐−1)²cos(3𝑙+2𝑔−2ℎ)+1715𝑒²+7(𝑐+1)²cos(4𝑙−2𝑔−2ℎ)15𝑒(𝑐+1)²cos(3𝑙−2𝑔−2ℎ)−75𝑒(𝑐+1)²−17cos(3𝑙+2𝑔+2ℎ)5𝑒²(𝑐+1)²cos(4𝑙+2𝑔+2ℎ)+715𝑒(𝑐−1)²+17cos(3𝑙−2𝑔+2ℎ)15𝑒²(𝑐−1)²cos(4𝑙−2𝑔+2ℎ)−115𝑒(𝑐+1)²+1cos(𝑙−2𝑔−2ℎ)5𝑒(𝑐−1)²cos(𝑙+2𝑔−2ℎ)15𝑒(𝑐+1)²−1cos(𝑙+2𝑔+2ℎ)15𝑒(𝑐−1)²cos(𝑙−2𝑔+2ℎ)−23𝑠²𝑒²−2+95cos(2𝑙−2𝑔)5𝑒²cos(2𝑙−2ℎ)−7+95𝑒cos(3𝑙−2𝑔)+15𝑒cos(𝑙−2𝑔)5𝑒²cos(2𝑙+2ℎ)+253𝑒²+2cos2ℎ−175𝑒²+6cos(4𝑙−2𝑔)5𝑒cos(𝑙+2ℎ)+65𝑒cos(𝑙−2ℎ),(2.7) where the disturbing potential is written in the form 𝑅=𝐻₂₀+𝐻₂₂.

2.1. The Hamiltonian System

The Hamiltonian of the dynamical system associated with the problem of the orbital motion of the satellite around the Moon taking into account the non-uniform distribution of the mass of the planetary satellite (𝐽₂–𝐽₅, 𝐽₇, and 𝐶₂₂) can be written in the following form:𝐾=𝐻₀+⁵_𝑖=2𝐻_𝑖0+𝐻₇₀+𝐻₂₂,(2.8) where𝐻₀=𝜇²2𝑎+𝑛^𝑀𝐻.(2.9)

Now, using ∑^4𝑖=1𝐻_(𝑖+1)0=∑^4𝑖=1𝜖_𝑖𝐻_𝑖, 𝐻₇₀=𝜖₅𝐻₅ and 𝐻₂₂=𝛿𝐻₆, we get𝐾=𝐻₀+⁵_𝑖=1𝜖_𝑖𝐻_𝑖+𝛿𝐻₆,(2.10) where 𝜖₁=𝐽₂𝑅^2𝑀, 𝜖₂=𝐽₃𝑅^3𝑀, 𝜖₃=𝐽₄𝑅^4𝑀, 𝜖₄=𝐽₅𝑅^5𝑀, 𝜖₅=𝐽₇𝑅^7𝑀, 𝛿=𝐶₂₂𝑅^2𝑀.

We arrange the Hamiltonian 𝐾 as follows𝐾=𝐻⁰⁰+𝐻¹¹,(2.11) where𝐻⁰⁰=𝜇²2𝑎+𝑛^𝑀𝐻𝐻,(2.12)₁₁=⁵_𝑖=1𝜖_𝑖𝐻_𝑖+𝛿𝐻₆.(2.13)

The motion of the spacecraft is studied under the single-averaged analytical model. The average is taken over the mean motion of the satellite. The standard definition for average used in this research is [24] ⟨𝐹⟩=1/2𝜋∫^02𝜋(𝐹)𝑑𝑙 where 𝑙 is the mean anomaly of the satellite. The single-average method is applied in (2.13) to eliminate the terms of short period, to analyze the effect of the disturbing potential on the orbital elements. Periodic terms will be calculated, replacing the result in the Lagrange planetary equations [25] and numerically integrated using the Maple software. The long-period disturbing potential is given in the appendix.

3. Applications: The Long-Period Disturbing Potential Using the Averaging Method

With a simplified model for the disturbing potential, it is possible to analyse the orbital motion of the satellite.

To study the lifetimes of low-altitude artificial satellites of the Moon [26] we took into account the zonal (𝐽₂–𝐽₅) and sectorial (𝐶₂₂, 𝐶₃₁) terms in the disturbing potential. However, [5] shows that the effect of the 𝐽₇ term is an order of magnitude larger than the 𝐽₂–𝐽₅ terms (see Table 2, model given by [27]), and therefore it should be considered in the potential. In [5], an analytic model is presented to find frozen orbits for lunar satellites, where the potential only due to 𝐽₂ and 𝐽₇ zonal harmonics is taken into account. In this paper, we present an approach taking into account the 𝐽₂–𝐽₅, 𝐽₇, and 𝐶₂₂ terms.

Table 2 

Coefficients given by [27].

3.1. Frozen Orbits Solutions

The definition of frozen orbit is written in the form𝑑𝑔𝑑𝑡=0,𝑑𝑑𝑒𝑡=0,𝑑𝑑𝑡𝑖=0.(3.1)

From this, we can obtain frozen orbits, say, for 𝑔=𝜋/2 or 𝑔=3𝜋/2, and ℎ=𝜋/2 or ℎ=3𝜋/2 (see [1, 13]).

Replacing by the potential in (A.1) in the Lagrange's planetary equations [25] and solving the equation 𝑑𝑔/𝑑𝑡=0, we get an equation that depends on three variables, 𝑒, 𝑖, and 𝑎, that can be represented as a 3-D surface as shown in Figures 1, 2, 3, and 4. Points on this surface correspond to equilibria of the system in equations 𝑑𝑔/𝑑𝑡=0 and 𝑑𝑒/𝑑𝑡=0, that is, frozen orbits. Figures 1 and 3 take into account the coefficients 𝐽₂–𝐽₅ and 𝐽₇, while Figures 2 and 4 take into account the 𝐽₂–𝐽₅, 𝐽₇, and 𝐶₂₂ terms. Comparing the models given by [27] (see Table 2) and [28] (see Table 3), we can analyze that the results are very close and the model given by [28] presents better results to obtain a family of frozen orbits for inclinations above 69°. In [30], an analytical theory based on Lie-Deprit transformation is used to locate family of frozen orbits for asteroids and natural satellites.

Table 3 

Model given by [28] (LP165P).

Figure 1 

𝑒 versus 𝑖-rad versus 𝑎-km. Akim, [27], 𝐽2–𝐽5, 𝐽7. Green color is the surface when 𝑔=3𝜋/2, and the red is obtained for 𝑔=𝜋/2.

Figure 2 

𝑒 versus 𝑖-rad versus 𝑎-km. Akim, [27], 𝐽2–𝐽5, 𝐽7 and 𝐶22. Green color is the surface when 𝑔=3𝜋/2, and the red is obtained for 𝑔=𝜋/2.

Figure 3 

𝑒 versus 𝑖-rad versus 𝑎-km. LP165P, Konopliv et al. [28], 𝐽2–𝐽5, 𝐽7. Green color is the surface when 𝑔=3𝜋/2, and the red is obtained for 𝑔=𝜋/2.

Figure 4 

𝑒 versus 𝑖-rad versus 𝑎-km. LP165P, Konopliv et al. [28], 𝐽2–𝐽5, 𝐽7 and 𝐶22. Green color is the surface when 𝑔=3𝜋/2, and the red is obtained for 𝑔=𝜋/2.

For inclinations larger than 69°, we found a family of frozen orbit, as we can see in Figures 1, 2, 3, and 4, mainly considering the model LP165P for the coefficients of the lunar gravity field [28]. We also observed that the 𝐶₂₂ term contributed efficiently to obtain of such orbits.

The curves obtained for constant-a sections are very similar; thus, in order to go into more detail, we define a constant value for the semimajor axis 𝑎=1861 km, which indeed corresponds to low orbits. For the chosen value of 𝑎, we plot the curve 𝑖 versus 𝑒 (see Figures 5 and 6). Figure 5 does not take into account the 𝐶₂₂ term in the potential, while Figure 6 considers the sectorial term. We observed an increase of the eccentricity for inclinations between 28.65° and 45.84° (comparing Figures 5 and 6). For orbits with inclination larger than 69° degree, we find a family of frozen orbits (see Figure 6). The region where the eccentricity is frozen is for small values of the initial eccentricity, when compared with Figure 5.

Figure 5 

𝑒 versus 𝑖-rad, 𝑎=1861 km. LP165P, Konopliv et al. [28], 𝐽2–𝐽5 and 𝐽7. Green color is obtained for 𝑔=3𝜋/2, and the red is obtained for 𝑔=𝜋/2.

Figure 6 

𝑒 versus 𝑖-rad, 𝑎=1861 km. LP165P, Konopliv et al. [28], 𝐽2–𝐽5, 𝐽7, and 𝐶22. Green color is obtained for 𝑔=3𝜋/2, and the red is obtained for 𝑔=𝜋/2.

Taking into account that Figure 6 represents a family of frozen orbits for inclinations of 90°, an example to illustrate the behavior of the frozen orbits for a case where the semimajor axis is 1861 km is shown in Figures 7, 9, and 10. Figures 7–10 are obtained by numerical integration (using the software Maple) of the Lagrange planetary equations for different values of the eccentricity. Figure 8 shows frozen orbits for a case where the semimajor axis is 1838 km.

Figure 7 

𝑒 versus 𝑔-degree. Model LP165P. Initial conditions: 𝑖=90°, 𝑎=1861 km, 𝑒=0.001, black color; 𝑒=0.002, red color; 𝑒=0.003, blue color; 𝑒=0.004 green color. Perturbations: 𝐽2–𝐽5, 𝐽7, and 𝐶22.

Figure 8 

𝑒 versus 𝑔-degree. Model LP165P. Initial conditions: 𝑖=90°, 𝑎=1838 km. Perturbations: 𝐽2–𝐽5, 𝐽7, and 𝐶22.

Figure 9 

𝑒 versus 𝑔-degree. Model LP165P. Initial conditions: 𝑖=90°, 𝑎=1861 km. Perturbations: 𝐽2–𝐽5, 𝐽7, and 𝐶22.

Figure 10 

𝑒 versus 𝑔-degree. Model LP165P. Initial conditions: 𝑖=90°, 𝑎=1861 km. Perturbations: 𝐽2–𝐽5 and 𝐶22.

We found frozen orbits with initial eccentricity equal to 0.02 that can be applied for missions of long period without the needs for doing orbit corrections for a period of very long time (larger than 4000 days). The initial orbit with appropriate characteristics can extend the time life and reduce maintenance costs. The set of initial conditions for frozen orbits with long lifetime is given by 𝑎=1861 km, 𝑒=0.02, 𝑖=90°, 𝑔=90°, and ℎ=270°.

In the literature, in general, frozen polar orbits are found with periapsis 𝑔=270° [3, 31–33], however, for orbits with inclinations near 90°, the argument of periapsis moves to 𝑔=90° [2, 31, 34]. These authors take into account only the zonal terms. In [34], there is an approach that compares the disturbing potential when a model of the gravitational field of high order taking into account the terms 𝐽₇, 𝐽₂₀, and 𝐽₅₀ is considered. For small inclinations, the three models show similar results for low-altitude orbits (see [34]), while for high inclination orbits the 𝐽₂₀ truncation only approximates roughly the 𝐽₅₀ case, giving frozen orbits that, in general, show lower eccentricities and shifting the zero eccentricity solutions to lower inclinations by several degrees. In this case, the potential truncated in the 𝐽₇ term provides a qualitative description only correct in the case of low inclinations. In [35], a numerical study considering models for the lunar gravity field of high order is applied to near-circular, low-altitude orbits. The behavior of the periapsis altitude is analyzed for some harmonics, for example, 𝐽₇ and 𝐽₉.

Here, we truncate the potential of the lunar gravity field up to the 𝐽₇ term, to analyze the effect caused by this harmonic in particular [5]. When we consider the potential taking into account the terms 𝐽₂–𝐽₅ and 𝐶₂₂ (zonal and sectorial), the frozen orbit (𝑖=90°, see Figure 10) librates around the equilibrium point for 𝑔=270° which is in agreement with the literature (see [3, 31, 32]). But, when we take into account the potential up to the 𝐽₇ term, we found frozen orbits that librate around the equilibrium point for 𝑔=90°, which is in agreement with [5]. In [34], a frozen orbit considering the following initial conditions: 𝑎=1861 km, 𝑒=0.0036, 𝑖=84°, and 𝑔=90°, is found, and it corresponds to a near-polar orbit.

Then, we develop the disturbing potential up to the 𝐽₉ term (𝐽₂–𝐽₉, 𝐶₂₂) to analyze what happens with the value of 𝑔 for a low-altitude, near-circular, near-polar orbit. Figure 11 shows the behavior for the diagram 𝑒 versus 𝑔 for a polar orbit which librates around the equilibrium point for the value of 𝑔=270° (see [33] for more details with respect to the potential truncated up to the 𝐽₉ term.). In this case, the initial eccentricity varies from 0.001 to 0.004. The numerical values used for the harmonics up to 𝐽₉ is given in Table 3. Thus, we conclude that, because of the characteristics of the lunar gravity field, the potential truncated up to the 𝐽₇ produces the effect of the frozen orbit obtained that librates around 𝑔=90° and not 𝑔=270°.

Figure 11 

𝑒 versus 𝑔-degree. Model LP165P. Initial conditions: 𝑖=90°, 𝑎=1861 km. Perturbations: 𝐽2–𝐽9 and 𝐶22.

The idea here is to consider a simplified potential to analyze the effect of the harmonics on the orbital elements 𝑒 and 𝑔. Therefore, this study did not realize the full potential integrations, which can be viewed in [33, 34].

Another important factor that we must analyze is the value of the argument of the periapsis, that strongly affects the determination of the frozen orbit. As we previously showed, the values of the periapsis to assure the condition of frozen orbit are 𝑔=90° or 𝑔=270°, equally for the longitude of the ascending node that is ℎ=90° or ℎ=270°. Since that, in general, we find frozen orbits for the value of 𝑔=90°. Figures 12–15 show the behavior of the eccentricity as a function of time for different values of the inclinations. Note that Figure 12 presents an orbit with lesser variation of the eccentricity for inclination of 𝑖=90°, considering 𝑔=90°. When we compare Figure 12 with Figure 13, we get that, for the value of 𝑔=270°, the orbit presents great variations of the eccentricity. We also analyze the case for the initial eccentricity of the order of 10⁻³ and verify that the behavior of the orbit does not suffer the same effect as occurred for the case of the eccentricity of order 10⁻². Here, the value of 𝑔=90° or 𝑔=270° has smaller effect on the amplitude of the eccentricity (see Figures 14 and 15).

Figure 12 

LP165P Model, 𝑒=0.02, 𝑔=90°, and ℎ=270°.

Figure 13 

LP165P Model, 𝑒=0.02, 𝑔=270°, and ℎ=90°.

Figure 14 

LP165P Model, 𝑒=0.006, 𝑔=90°, and ℎ=270°.

Figure 15 

LP165P Model, 𝑒=0.006, 𝑔=270°, and ℎ=90°.

4. Conclusions

In this paper, the dynamics of orbits around a planetary satellite (Moon), taking into account the non-uniform distribution of mass (𝐽₂–𝐽₉ and 𝐶₂₂), was studied. Several analyses were performed, using the long-period disturbing potential obtained through the single-averaged model.

Considering the long-period disturbing potential, we found a family of frozen orbits for inclinations larger than 69° considering the model LP165P for the coefficients of the lunar gravity field. It was observed that the 𝐶₂₂ term contributed efficiently to obtain such orbits.

Low-altitude, near-polar orbits are very desirable for scientific missions to study natural satellites, such as the Moon. Due to the characteristics of the lunar gravity field, the potential truncated up to the 𝐽₇ produces the effect of the frozen orbit that librates around 𝑔=90° and not 𝑔=270°. Considering the long-period disturbing potential using the single-averaged method, we found near-circular, near-polar frozen orbits where this choice reduces the maintenance cost considerably.

In general, we conclude that, for the case of the Moon, it is necessary to consider a gravity model of high order, including the zonal and sectorial terms. Because the question of the coefficients are almost of the same order, this fact complicates the choice of the harmonic where the potential can be truncated and fact makes its choice a little arbitrary.

To study the lifetime of a lunar artificial satellite considering an analytical model, the zonal (𝐽₂–𝐽₉) and sectorial (𝐶₂₂) terms of the harmonic coefficients should be taken into account. We observe that the frozen orbits could be strongly affected by terms of the potential containing the coefficient due to the equatorial ellipticity of the Moon (𝐶₂₂) and by terms containing the argument of the periapsis and the longitude of the ascending node.

Appendix

The long-period disturbing potential up to the 𝐽₇ term is 𝑘1=7516𝑎⁵𝑛²−8807799×1−6318340033038129𝑐²[phone]+𝑐⁴𝜖₄𝑒⁴+−1113614665𝜖⁴𝑐⁴+18432419419𝑎²𝜖₃[phone]𝜖⁴𝑐²−942083774771𝑎²𝜖₃−104448×𝑒₂2097095𝜖⁴−2562933𝜖⁴𝑐⁴+2048419419𝑎²𝜖₃[phone]𝜖⁴𝑐²−20483774771𝑎²𝜖₃−768419419𝜖⁴𝑠³𝑒³sin(3𝑔)[phone]𝑠⁵×𝜖₄−1178377696865+𝑐²𝑒²[phone]𝜖⁴−662453605𝜖⁴𝑐²−368647665515𝑎²×𝑒₅𝜖₃sin(5𝑔)−4940𝑒²+27𝜖²𝑐²−17𝑎³+12𝑠²𝑒²cos(2𝑔)25𝛿𝑎⁵𝑠²𝑒²+23cos(2ℎ)[phone]𝜖⁴[phone]𝑠⁷𝑒⁷sin(7𝑔)1474560𝑐⁶−4009815930930133𝑐⁴−339439251020694389+1303793410913353𝑐²𝜖⁴𝑒⁶+19649635120𝜖⁴𝑐⁶+−5173160𝑎²𝜖₃−22567935120𝜖⁴𝑐⁴+1217651024𝜖⁴+2527120𝑎²𝜖₃𝑐²−637480𝑎²𝜖₃−42631024𝜖⁴𝑒⁴+27027320𝜖⁴𝑐⁶+−44140𝑎²𝜖₃−623764𝜖⁴𝑐⁴×147+20𝑎²𝜖₃+170164𝜖⁴+𝜖₁𝑎⁴𝑐²−2140𝑎²𝜖₃−15𝜖¹𝑎⁴−63×𝑒₂64𝜖⁴+3003320𝜖⁴𝑐⁶+−2110𝑎²𝜖₃−69364𝜖⁴𝑐⁴+18964𝜖⁴+25𝜖¹𝑎⁴+75𝑎²𝜖₃𝑐²−225𝜖¹𝑎⁴−764𝜖⁴−110𝑎²𝜖₃𝑠𝑒sin(𝑔)−147×𝑐₄32−67𝑐²+335𝜖²𝑒⁴+821𝜖²𝑐⁴+−641225𝑎²𝜖−1649𝜖²𝑐²+643675𝑎²𝜖+8×𝑒₂245𝜖²+8105𝜖²𝑐⁴+−16245𝜖²−1283675𝑎²𝜖𝑐²[phone]𝑎²𝜖+81225𝜖²𝑎³,(A.1)

where the values for the zonal (𝐽₂–𝐽₇) and tesseral (𝐶₂₂) harmonics coefficients taking into account the two models are given in Table 2 [27] and Table 3 [28].

Acknowledgment

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