# Planetary Satellite Orbiters: Applications for the Moon

Authors: Jean Paulo dos Santos Carvalho; Rodolpho Vilhena de Moraes; Antônio Fernando Bertachini de Almeida Prado
Journal: Mathematical Problems in Engineering
 (2011)

---

## 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 (J2–J5, J7, and C22). 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.


---

## Body

## 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 (J2, J3, and C22) 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 (J2, C22) 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 (J2) zonal, (C22) 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 (J2–J5, J7, and C22) 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. Section2 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 C20 term describes the equatorial bulge of the Moon, often referred to as oblateness. The coefficient C22 measures the ellipticity of the equator.Table 1
Orders of magnitude forJ2 and C22 [14].
EarthMoonC20≡(-J2)-10-3-2×10-4C222×10-62×10-5The force function, the negative of the total energy as used in physics, is given byF=K=V-T here, V is the opposite of the potential energy, and T is the kinetic energy. The force function can be written as K=μ/r+R-Tor K=μ2/2a+R. The function R, comprising all terms of V 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 a (semimajor axis), e (eccentricity), i (inclination), g (argument of the periapsis), h (longitude of the ascending node), and n (mean motion) given by the third Kepler law n2a3=Gm0, where G is the gravitational constant and r is the radius vector of the body m0 (mass of the planetary satellite). The term due to the unperturbed potential is given by(2.1)H0=μ22a.Considering the lunar equatorial plane as the reference plane, the disturbing potential can be written in the form [21](2.2)VM=-μr[∑n=25(RMr)nJnPn(sinϕ)-(RMr)2C22P22(sinϕ)cos⁡2λ],
where μ is the lunar gravitational constant, RM is the equatorial radius of the Moon (RM=1738 km), Pn represent the Legendre polynomials, Pnm 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 λ=λ′-λ22, since λ′ is the longitude measured from any fixed direction, and λ22 is the longitude of the Moon's longest meridian from the same fixed direction. However, λ22 contains explicitly the time [20]. Using spherical trigonometry, we have sin(ϕ)=sin(i)sin(f+g), where f 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 λ22.Since the variablesΩ and λ⊕ appear only as a combination of Ω-λ⊕, where λ⊕=nMt+const, nM being the lunar mean motion. The degrees of freedom can be reduced by choosing h=Ω-λ⊕ as a new variable. A new term must then be added to the Hamiltonian in order to get ḣ=Ω̇-nM=-∂(K+nMH)/∂H. This term nMH 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 x-axis passes through this meridian.Regarding the Earth's potential, the dominant coefficient isJ2. 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 toJ5 and sectorial (C22) terms can be written in the form [20, 21](2.3)P2(sin∅)=12(3s2sin2(f+g)-1),P3(sin∅)=52s3sin3(f+g)-32ssin⁡(f+g),P4(sin∅)=358s4sin4(f+g)-354s2sin2(f+g)+38,P5(sin∅)=638s5sin5(f+g)-354s2sin3(f+g)+158ssin(f+g),P22(sin∅)cos⁡2λ=6(ξ2cos⁡2f+x2sin2f+ξxsin2f)-3(1-s2sin2(f+g)),
where we use the shortcuts ξ=cosgcosh-csingsinh, χ=-singcosh-ccosgsinh, s=sini, and c=cosi.With this, the zonal perturbation due to the oblateness is [14] (2.4)H20=ϵμ4r3(1-3c2-3s2cos⁡(2f+2g)),withϵ=J2RM2.Replacing the relationμ=n2a3 and using Cayley's tables [23] to express the true anomaly in terms of the mean anomaly, after some algebraic manipulations, we get(2.5)H20=3n28ϵ[(c2-13)(29+23ecos⁡l+13e2+e2cos⁡2l)((5e2-2)c2+2-5e2)cos⁡⁡(2g+3l)+7es2cos⁡⁡(2g+3l)+17e2s2cos⁡⁡(2g+4l)-es2cos⁡⁡(2g+l)+9(c2-13)(29+23ecos⁡l+13e2+e2cos⁡2l)].For the sectorial perturbation [20, 21], we get (2.6)H22=δ3μr3[2ξ2cos⁡2f+2χ2sin2f+2ξχsin2f-1+s2sin2(f+g)],whereδ=C22RM2. With some manipulations, we get(2.7)H22=-4516δn2[(e2-25)(c-1)2cos⁡(2l+2g-2h)-13(c+1)2(e2-25)cos⁡(2l-2g-2h)+(c+1)2(e2-25)cos⁡(2l+2g+2h)-13(e2-25)(c-1)2cos⁡(2l-2g+2h)-175e2(c-1)2cos⁡(4l+2g-2h)-75e(c-1)2cos⁡(3l+2g-2h)+1715e2(c+1)2cos⁡(4l-2g-2h)+715e(c+1)2cos⁡(3l-2g-2h)-75e(c+1)2cos⁡(3l+2g+2h)-175e2(c+1)2cos⁡(4l+2g+2h)+715e(c-1)2cos⁡(3l-2g+2h)+1715e2(c-1)2cos⁡(4l-2g+2h)-115e(c+1)2cos⁡(l-2g-2h)+15e(c-1)2cos⁡(l+2g-2h)15e(c+1)2cos⁡(l+2g+2h)-115e(c-1)2cos⁡(l-2g+2h)-23s2(e2-25)cos⁡(2l-2g)+95e2cos⁡(2l-2h)-75ecos⁡(3l-2g)+15ecos⁡(l-2g)+95e2cos⁡(2l+2h)+25(3e2+2)cos⁡2h-175e2cos⁡(4l-2g)+65ecos⁡(l+2h)+65ecos⁡(l-2h)],
where the disturbing potential is written in the form R=H20+H22.


### 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 (J2–J5, J7, and C22) can be written in the following form:(2.8)K=H0+∑i=25Hi0+H70+H22,
where(2.9)H0=μ22a+nMH.Now, using∑i=14H(i+1)0=∑i=14ϵiHi, H70=ϵ5H5 and H22=δH6, we get(2.10)K=H0+∑i=15ϵiHi+δH6,
where ϵ1=J2RM2, ϵ2=J3RM3, ϵ3=J4RM4, ϵ4=J5RM5, ϵ5=J7RM7, δ=C22RM2.We arrange the HamiltonianK as follows(2.11)K=H00+H11,
where(2.12)H00=μ22a+nMH,(2.13)H11=∑i=15ϵiHi+δH6.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] 〈F〉=1/2π∫02π(F)dl where l 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.


## 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 (J2–J5, J7, and C22) can be written in the following form:(2.8)K=H0+∑i=25Hi0+H70+H22,
where(2.9)H0=μ22a+nMH.Now, using∑i=14H(i+1)0=∑i=14ϵiHi, H70=ϵ5H5 and H22=δH6, we get(2.10)K=H0+∑i=15ϵiHi+δH6,
where ϵ1=J2RM2, ϵ2=J3RM3, ϵ3=J4RM4, ϵ4=J5RM5, ϵ5=J7RM7, δ=C22RM2.We arrange the HamiltonianK as follows(2.11)K=H00+H11,
where(2.12)H00=μ22a+nMH,(2.13)H11=∑i=15ϵiHi+δH6.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] 〈F〉=1/2π∫02π(F)dl where l 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 (J2–J5) and sectorial (C22, C31) terms in the disturbing potential. However, [5] shows that the effect of the J7 term is an order of magnitude larger than the J2–J5 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 J2 and J7 zonal harmonics is taken into account. In this paper, we present an approach taking into account the J2–J5, J7, and C22 terms.Table 2
Coefficients given by [27].
CoefficientsJ2+2.070×10-4J3+4.900×10-6J4+8.000×10-7J5-3.600×10-6J6-1.100×10-6J7-2.870×10-5C22+2.447305×10-5


### 3.1. Frozen Orbits Solutions

The definition of frozen orbit is written in the form(3.1)dgdt=0,dedt=0,didt=0.From this, we can obtain frozen orbits, say, forg=π/2 or g=3π/2, and h=π/2 or h=3π/2 (see [1, 13]).Replacing by the potential in (A.1) in the Lagrange's planetary equations [25] and solving the equation dg/dt=0, we get an equation that depends on three variables, e, i, and a, 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 dg/dt=0 and de/dt=0, that is, frozen orbits. Figures 1 and 3 take into account the coefficients J2–J5 and J7, while Figures 2 and 4 take into account the J2–J5, J7, and C22 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).
Coefficients [29]J2+2.032337×10-4J3+8.475900×10-6J4-9.591930×10-6J5+7.154090×10-7J6-1.357770×10-5J7-2.177470×10-5J8-9.674870×10-6J9+1.549600×10-5C22+2.235700×10-5Figure 1
e versus i-rad versus a-km. Akim, [27], J2–J5, J7. Green color is the surface when g=3π/2, and the red is obtained for g=π/2.Figure 2
e versus i-rad versus a-km. Akim, [27], J2–J5, J7 and C22. Green color is the surface when g=3π/2, and the red is obtained for g=π/2.Figure 3
e versus i-rad versus a-km. LP165P, Konopliv et al. [28], J2–J5, J7. Green color is the surface when g=3π/2, and the red is obtained for g=π/2.Figure 4
e versus i-rad versus a-km. LP165P, Konopliv et al. [28], J2–J5, J7 and C22. Green color is the surface when g=3π/2, and the red is obtained for g=π/2.For inclinations larger than 69°, we found a family of frozen orbit, as we can see in Figures1, 2, 3, and 4, mainly considering the model LP165P for the coefficients of the lunar gravity field [28]. We also observed that the C22 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 axisa=1861 km, which indeed corresponds to low orbits. For the chosen value of a, we plot the curve i versus e (see Figures 5 and 6). Figure 5 does not take into account the C22 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
e versus i-rad, a=1861 km. LP165P, Konopliv et al. [28], J2–J5 and J7. Green color is obtained for g=3π/2, and the red is obtained for g=π/2.Figure 6
e versus i-rad, a=1861 km. LP165P, Konopliv et al. [28], J2–J5, J7, and C22. Green color is obtained for g=3π/2, and the red is obtained for g=π/2.Taking into account that Figure6 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
e versus g-degree. Model LP165P. Initial conditions: i=90°, a=1861 km, e=0.001, black color; e=0.002, red color; e=0.003, blue color; e=0.004 green color. Perturbations: J2–J5, J7, and C22.Figure 8
e versus g-degree. Model LP165P. Initial conditions: i=90°, a=1838 km. Perturbations: J2–J5, J7, and C22.Figure 9
e versus g-degree. Model LP165P. Initial conditions: i=90°, a=1861 km. Perturbations: J2–J5, J7, and C22.Figure 10
e versus g-degree. Model LP165P. Initial conditions: i=90°, a=1861 km. Perturbations: J2–J5 and C22.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 bya=1861 km, e=0.02, i=90°, g=90°, and h=270°.In the literature, in general, frozen polar orbits are found with periapsisg=270° [3, 31–33], however, for orbits with inclinations near 90°, the argument of periapsis moves to g=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 J7, J20, and J50 is considered. For small inclinations, the three models show similar results for low-altitude orbits (see [34]), while for high inclination orbits the J20 truncation only approximates roughly the J50 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 J7 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, J7 and J9.Here, we truncate the potential of the lunar gravity field up to theJ7 term, to analyze the effect caused by this harmonic in particular [5]. When we consider the potential taking into account the terms J2–J5 and C22 (zonal and sectorial), the frozen orbit (i=90°, see Figure 10) librates around the equilibrium point for g=270° which is in agreement with the literature (see [3, 31, 32]). But, when we take into account the potential up to the J7 term, we found frozen orbits that librate around the equilibrium point for g=90°, which is in agreement with [5]. In [34], a frozen orbit considering the following initial conditions: a=1861 km, e=0.0036, i=84°, and g=90°, is found, and it corresponds to a near-polar orbit.Then, we develop the disturbing potential up to theJ9 term (J2–J9, C22) to analyze what happens with the value of g for a low-altitude, near-circular, near-polar orbit. Figure 11 shows the behavior for the diagram e versus g for a polar orbit which librates around the equilibrium point for the value of g=270° (see [33] for more details with respect to the potential truncated up to the J9 term.). In this case, the initial eccentricity varies from 0.001 to 0.004. The numerical values used for the harmonics up to J9 is given in Table 3. Thus, we conclude that, because of the characteristics of the lunar gravity field, the potential truncated up to the J7 produces the effect of the frozen orbit obtained that librates around g=90° and not g=270°.Figure 11
e versus g-degree. Model LP165P. Initial conditions: i=90°, a=1861 km. Perturbations: J2–J9 and C22.The idea here is to consider a simplified potential to analyze the effect of the harmonics on the orbital elementse and g. 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 areg=90° or g=270°, equally for the longitude of the ascending node that is h=90° or h=270°. Since that, in general, we find frozen orbits for the value of g=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 i=90°, considering g=90°. When we compare Figure 12 with Figure 13, we get that, for the value of g=270°, the orbit presents great variations of the eccentricity. We also analyze the case for the initial eccentricity of the order of 10−3 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−2. Here, the value of g=90° or g=270° has smaller effect on the amplitude of the eccentricity (see Figures 14 and 15).Figure 12
LP165P Model,e=0.02, g=90°, and h=270°.Figure 13
LP165P Model,e=0.02, g=270°, and h=90°.Figure 14
LP165P Model,e=0.006, g=90°, and h=270°.Figure 15
LP165P Model,e=0.006, g=270°, and h=90°.


## 3.1. Frozen Orbits Solutions

The definition of frozen orbit is written in the form(3.1)dgdt=0,dedt=0,didt=0.From this, we can obtain frozen orbits, say, forg=π/2 or g=3π/2, and h=π/2 or h=3π/2 (see [1, 13]).Replacing by the potential in (A.1) in the Lagrange's planetary equations [25] and solving the equation dg/dt=0, we get an equation that depends on three variables, e, i, and a, 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 dg/dt=0 and de/dt=0, that is, frozen orbits. Figures 1 and 3 take into account the coefficients J2–J5 and J7, while Figures 2 and 4 take into account the J2–J5, J7, and C22 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).
Coefficients [29]J2+2.032337×10-4J3+8.475900×10-6J4-9.591930×10-6J5+7.154090×10-7J6-1.357770×10-5J7-2.177470×10-5J8-9.674870×10-6J9+1.549600×10-5C22+2.235700×10-5Figure 1
e versus i-rad versus a-km. Akim, [27], J2–J5, J7. Green color is the surface when g=3π/2, and the red is obtained for g=π/2.Figure 2
e versus i-rad versus a-km. Akim, [27], J2–J5, J7 and C22. Green color is the surface when g=3π/2, and the red is obtained for g=π/2.Figure 3
e versus i-rad versus a-km. LP165P, Konopliv et al. [28], J2–J5, J7. Green color is the surface when g=3π/2, and the red is obtained for g=π/2.Figure 4
e versus i-rad versus a-km. LP165P, Konopliv et al. [28], J2–J5, J7 and C22. Green color is the surface when g=3π/2, and the red is obtained for g=π/2.For inclinations larger than 69°, we found a family of frozen orbit, as we can see in Figures1, 2, 3, and 4, mainly considering the model LP165P for the coefficients of the lunar gravity field [28]. We also observed that the C22 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 axisa=1861 km, which indeed corresponds to low orbits. For the chosen value of a, we plot the curve i versus e (see Figures 5 and 6). Figure 5 does not take into account the C22 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
e versus i-rad, a=1861 km. LP165P, Konopliv et al. [28], J2–J5 and J7. Green color is obtained for g=3π/2, and the red is obtained for g=π/2.Figure 6
e versus i-rad, a=1861 km. LP165P, Konopliv et al. [28], J2–J5, J7, and C22. Green color is obtained for g=3π/2, and the red is obtained for g=π/2.Taking into account that Figure6 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
e versus g-degree. Model LP165P. Initial conditions: i=90°, a=1861 km, e=0.001, black color; e=0.002, red color; e=0.003, blue color; e=0.004 green color. Perturbations: J2–J5, J7, and C22.Figure 8
e versus g-degree. Model LP165P. Initial conditions: i=90°, a=1838 km. Perturbations: J2–J5, J7, and C22.Figure 9
e versus g-degree. Model LP165P. Initial conditions: i=90°, a=1861 km. Perturbations: J2–J5, J7, and C22.Figure 10
e versus g-degree. Model LP165P. Initial conditions: i=90°, a=1861 km. Perturbations: J2–J5 and C22.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 bya=1861 km, e=0.02, i=90°, g=90°, and h=270°.In the literature, in general, frozen polar orbits are found with periapsisg=270° [3, 31–33], however, for orbits with inclinations near 90°, the argument of periapsis moves to g=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 J7, J20, and J50 is considered. For small inclinations, the three models show similar results for low-altitude orbits (see [34]), while for high inclination orbits the J20 truncation only approximates roughly the J50 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 J7 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, J7 and J9.Here, we truncate the potential of the lunar gravity field up to theJ7 term, to analyze the effect caused by this harmonic in particular [5]. When we consider the potential taking into account the terms J2–J5 and C22 (zonal and sectorial), the frozen orbit (i=90°, see Figure 10) librates around the equilibrium point for g=270° which is in agreement with the literature (see [3, 31, 32]). But, when we take into account the potential up to the J7 term, we found frozen orbits that librate around the equilibrium point for g=90°, which is in agreement with [5]. In [34], a frozen orbit considering the following initial conditions: a=1861 km, e=0.0036, i=84°, and g=90°, is found, and it corresponds to a near-polar orbit.Then, we develop the disturbing potential up to theJ9 term (J2–J9, C22) to analyze what happens with the value of g for a low-altitude, near-circular, near-polar orbit. Figure 11 shows the behavior for the diagram e versus g for a polar orbit which librates around the equilibrium point for the value of g=270° (see [33] for more details with respect to the potential truncated up to the J9 term.). In this case, the initial eccentricity varies from 0.001 to 0.004. The numerical values used for the harmonics up to J9 is given in Table 3. Thus, we conclude that, because of the characteristics of the lunar gravity field, the potential truncated up to the J7 produces the effect of the frozen orbit obtained that librates around g=90° and not g=270°.Figure 11
e versus g-degree. Model LP165P. Initial conditions: i=90°, a=1861 km. Perturbations: J2–J9 and C22.The idea here is to consider a simplified potential to analyze the effect of the harmonics on the orbital elementse and g. 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 areg=90° or g=270°, equally for the longitude of the ascending node that is h=90° or h=270°. Since that, in general, we find frozen orbits for the value of g=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 i=90°, considering g=90°. When we compare Figure 12 with Figure 13, we get that, for the value of g=270°, the orbit presents great variations of the eccentricity. We also analyze the case for the initial eccentricity of the order of 10−3 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−2. Here, the value of g=90° or g=270° has smaller effect on the amplitude of the eccentricity (see Figures 14 and 15).Figure 12
LP165P Model,e=0.02, g=90°, and h=270°.Figure 13
LP165P Model,e=0.02, g=270°, and h=90°.Figure 14
LP165P Model,e=0.006, g=90°, and h=270°.Figure 15
LP165P Model,e=0.006, g=270°, and h=90°.


## 4. Conclusions

In this paper, the dynamics of orbits around a planetary satellite (Moon), taking into account the non-uniform distribution of mass (J2–J9 and C22), 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 theC22 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 theJ7 produces the effect of the frozen orbit that librates around g=90° and not g=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 (J2–J9) and sectorial (C22) 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 (C22) and by terms containing the argument of the periapsis and the longitude of the ascending node.



---

*Source: 187478-2011-05-25.xml*