A Study of Single-and Double-Averaged Second-Order Models to Evaluate Third-Body Perturbation Considering Elliptic Orbits for the Perturbing Body

The equations for the variations of the Keplerian elements of the orbit of a spacecraft perturbed by a third body are developed using a single average over the motion of the spacecraft, considering an elliptic orbit for the disturbing body. A comparison is made between this approach and the more used double averaged technique, as well as with the full elliptic restricted three-body problem. The disturbing function is expanded in Legendre polynomials up to the second order in both cases. The equations of motion are obtained from the planetary equations, and several numerical simulations are made to show the evolution of the orbit of the spacecraft. Some characteristics known from the circular perturbing body are studied: circular, elliptic equatorial, and frozen orbits. Different initial eccentricities for the perturbed body are considered, since the effect of this variable is one of the goals of the present study. The results show the impact of this parameter as well as the differences between both models compared to the full elliptic restricted three-body problem. Regions below, near, and above the critical angle of the third-body perturbation are considered, as well as different altitudes for the orbit of the spacecraft.


Introduction
Most of the papers on this topic consider the third-body perturbation due to the Sun and due to the Moon in a satellite around the Earth.This is the most immediate application of the third-body perturbation.One of the first studies was made by Kozai [1] that developed the most important longperiod and secular terms of the perturbing potential of the lunisolar perturbations, written as a function of the orbital elements of the Sun, the Moon, and the satellite.Moe [2], Musen [3], and Cook [4] studied long-period effects of the Sun and the Moon on artificial satellites of the Earth.They applied Lagrange's planetary equations to study the variation of the orbital elements of the satellite and its rate of variation.This idea was expanded by Musen et al. [5] that included the parallactic term in the perturbing potential.Kozai [6] studied the secular perturbations in asteroids that are in orbits with high inclination and eccentricity, assuming that the bodies are perturbed by Jupiter.Blitzer [7] made estimates for the lunisolar perturbation for the secular terms.Around the same time, Kaula [8] obtained the disturbing function considering the perturbations of the Sun and the Moon.
Later, Giacaglia [9] calculated the disturbing function due to the Moon using equatorial elements for the satellite and ecliptic elements for the Moon.All terms were expressed in closed forms.Kozai [10] worked again on that problem and expressed the perturbing function as a function of the polar geocentric coordinates of the Moon, the Sun, and the orbital elements of the satellite.The short-period terms are shown in an analytical form, and the secular and long-period terms are obtained by numerical integration.
Hough [11] considered the effects of the perturbation of the Sun and the Moon in orbits near the inclinations of 63.4 ∘ and 116.6 ∘ (critical inclinations when considering the geopotential of the Earth) and showed that those effects are important only in high altitudes.Ash [12] also studied this problem using the double-averaged technique for a high-altitude satellite around the Earth.Collins and Cefola [13] used the same technique for the prediction of orbits for a long time span.
In the last years, several researches, based on both analytical and numerical approaches, studied the third-body perturbation.The majority of them concentrated on studying the effects of a perturbation caused by a third body using a double-averaged technique, like those of Šidlichovskyý [14], Kwok [15,16], Broucke [17], and Prado [18].The disturbing function is expanded in Legendre's polynomials, and the average is made over the periods of the perturbing and the perturbed body, using an approximated model expanded to some order.Some others studied this problem based on a single-averaged analytical model, where the disturbing function is expanded in Legendre's polynomials but the average is made over the short period of the perturbed body only; see, for example, [19] and references therein.In all of the researches cited before, the perturbing body was in a circular orbit.There are also several papers considering the thirdbody perturbation in constellations of satellites (see [20,21]) and in other planetary systems not involving the Earth, like those of Carvalho (see [22][23][24][25]), Folta and Quinn [26], Kinochita and Nakai [27], Lara [28], Paskowitz and Scheers (see [29,30]), Russel and Brinckerhoff [31], and Scheers et al. [32].In particular, Roscoe et al. [33] used the analytical equations derived in Prado [20] to study the dynamics of a constellation of satellites perturbed by a third body.
Using the double-averaged analytical model in Domingos et al. [34], the analytical expansion to study the third-body perturbation was extended to the case where the perturbing body is in an elliptical orbit.Liu et al. [35] included the inclination of the perturbing body's equatorial plane with respect to its orbital plane.
The main idea behind the single-averaged technique is that it eliminates only the terms due to the short-time periodic motion of the perturbed body.The results are expected to show smooth curves for the evolution of the mean orbital elements for a long-time period when compared with the full restricted three-body problem.In other words, 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 slow changes in the orbital elements.
So, an interesting point would be a study to show the differences between those two averaged models when compared with the full elliptic restricted three-body problem.The idea of the present paper is to study this problem, but assuming that the perturbing body is in an elliptical orbit.Studies under this assumption are available in Domingos et al. (see [34,36,37]).In particular, in the present paper, a study of the effects of the eccentricity of the perturbing body in the orbit of the perturbed body is made for several different conditions regarding the orbit of the primaries as well as the orbit of the spacecraft.
So, our goal is to make a more complete study of this problem and to perform some tests to verify the differences in the results obtained by those two approximated techniques.A comparative investigation is made to verify the differences between the single-averaged analytical model with the model based on the double-averaged technique, as well as a comparison with the full restricted elliptic three-body problem that can provide some insights about their applications in celestial mechanics.
The assumptions used to develop the single-averaged analytical model are the same ones of the restricted elliptic three-body problem (planet-satellite-spacecraft).Our analysis for the evolution of the mean orbital elements will be based only on gravitational forces.The equations of motion are obtained from Lagrange's planetary equations, and then we numerically integrated those equations.Different initial eccentricities for the perturbing body are considered.
The set of results obtained in this research performs an analysis of several well-known characteristics of the thirdbody perturbation, like (i) the stability of near-circular orbits, to investigate under which conditions this orbit remains near circular; (ii) the behavior of elliptic equatorial orbits; and (iii) the existence of frozen orbits, which are orbits that keep the eccentricity, inclination, and argument of periapsis constants.
A detailed study considering the "critical angle" of the third-body perturbation, which is an inclination that makes a near circular orbit that has an inclination below this value to remain near-circular, is made.This work is structured as follows.In Section 2, we present the equations of motion used for the numerical simulations.Section 3 is devoted to the analysis of the numerical results for near circular orbits.The theory developed here is used to study the behavior of a spacecraft around the Earth, where the Moon is the disturbing body.The choice of the Moon as the disturbing body is made based on the fact that its effect on the spacecraft is larger than the effect of the Sun.Several plots show the time histories of the Keplerian elements of the orbits involved.Our conclusions are presented in Section 4.

Dynamical Model
For the determination of the equations of motion, we started by assuming the existence of a central body, with mass  0 , that is fixed in the center of the reference system --.The perturbing body, with mass   , is in an elliptic orbit with semimajor axis   , eccentricity   , inclination   , argument of periapsis   , longitude of the ascending node Ω  , and mean motion   .The spacecraft with negligible mass  is in a generic orbit with orbital elements:  (semimajor axis),  (eccentricity),  (inclination),  (argument of periapsis), Ω (longitude of the ascending node), and it has mean motion .This system is shown in Figure 1.
Using the traditional expansion in Legendre's polynomials (assuming that   ≫ ), the second-order term of the disturbing potential averaged over the eccentric anomaly of the spacecraft is given by the following (see [17,18,34]): where (Murray and Dermott [38]) For the case of elliptic orbits of the perturbing body, the parameters  and  are written as follows (see [34]): where  = Ω −   −   .The mean anomaly   of the perturbing body is given as Thus, the variations in the orbital elements of the perturbed body are obtained.To do this, we derived Lagrange's planetary equations that describe the variations of the mean orbital elements of the spacecraft.The semimajor axis is constant, since the mean anomaly  was eliminated from the perturbing function.Those equations can be written as shown below: where Those equations show some characteristics of this system compared with similar researches (see [18,34]) which are listed as follows.
(i) Circular orbits exist, but they are not planar.It is immediate from the inspections of the time derivatives of the Keplerian elements.If the initial eccentricity is zero, then the time derivative of the eccentricity is also zero and the orbit remains circular, but the time derivative of the inclination is not zero, so the orbit does not stay planar.
(ii) Elliptic equatorial orbits are not stable.The equatorial case ( = 0) has a singularity in the equation for the time derivative of the inclination, so this model is not able to make prediction for the behavior of the inclination.Regarding the eccentricity, it is visible that this case has a nonzero value for the time derivative of the eccentricity, so those orbits do not keep the eccentricity constant.The only exception is when the initial orbit is circular, as shown above, where the eccentricity remains zero.In fact, the equation for the time derivative of the eccentricity for equatorial orbits can be simplified to the following: (iii) Those facts explained above also show that there are no frozen orbits under this analytical model, which would be orbits where the time derivatives of the inclination, eccentricity, and argument of periapsis are all zero.[34]) made a similar study, but using the double-averaged technique and showed that circular and equatorial orbits exist in the double-averaged model, since orbits with zero initial inclination or zero initial eccentricity remain with zero eccentricity and zero inclination.The difference in the equations of motion when considering the circular and the elliptic motion of the primaries under the doubleaveraged models is the existence of the term (1 + (3/2)  2 + (15/8)  4 ), which is an extra term that depends on the eccentricity of the primaries   and this term does not destroy those properties.The same is true for the frozen orbits, since the appearance of a multiplicative nonzero constant does not change this property.

Numerical Results
In this section, we show the effects of the nonzero eccentricities of the perturbing body in the orbit of the perturbed body.We numerically investigate the evolutions of the eccentricity and the inclination for a spacecraft within the elliptic restricted three-body problem Earth-Moon-spacecraft, as well as using the single-and the double-averaged models up to the second order.The spacecraft is in an elliptic three-dimensional orbit around the Earth, and its motion is perturbed by the Moon.To justify our study and to make it a representative sample of the range of possibilities, we used three values for the semimajor axis of the spacecraft: 8000 km, 26000 km, and 42000 km.They represent a Low-Altitude Earth Orbit, a Medium-Altitude Earth Orbit (the one used by the GPS satellites), and a High-Altitude Earth Orbit (the one used by geosynchronous satellites), respectively.The spacecraft is assumed to be in an orbit with the following initial Keplerian elements:  0 = 0.01,  0 = Ω 0 = 0.
We focus our attention On the stability of near-circular orbits.Results of numerical integrations showed that this stability depends on  0 (the initial inclination between the perturbed and perturbing bodies).When  0 is above the critical value (  ), the orbit becomes very elliptic.If  0 <   , then the orbit remains near circular (see [17,18,34]).The value of   is 0.684719203 radians or 39.23152048 degrees.Thus, for our numerical integrations, the initial inclinations for the orbit of the spacecraft received values near the critical inclination (in the interval  0 = 38 ∘ and 41 ∘ ), below this critical value ( 0 = 30 ∘ , 20 ∘ , 10 ∘ , and zero), and then above this critical value ( 0 = 40 ∘ , 50  , 60 ∘ , 70 ∘ , and 80 ∘ ).The eccentricity of the perturbing body is assumed to be in the range 0 ≤   ≤ 0.5, in general, but some results are omitted when the approximations are found to be very poor.It means that we generalized the Earth-Moon system in order to measure the effects of the eccentricity of the orbit of the primaries.
All of the cases considering the Low-Altitude Earth Orbit were simulated for 80000 canonical units of time, which correspond to 12800 orbits of the disturbing body.This time shows the characteristics of the system, since the third-body perturbations are weaker at  0 = 8000 km.The simulations for the other values of the semimajor axis were made only for 20000 canonical units of time, which correspond to 3200 orbits of the disturbing body, because it appeared to be a good number to evidence the characteristics of the system.
Our numerical results are summarized in Figures 2-6.Such figures refer to the orbital evolution of the eccentricity and inclination as a function of time for the initial inclinations cited above and a direct comparison between the models.

Low-Inclination Orbits.
Figure 2 shows the results of the evolution of the eccentricity as a function of time for the case where the initial inclination assumes the values  0 = 30 ∘ (red line), 20 ∘ (blue line), 10 ∘ (pink line), and zero (orange line).The first, second, and third columns show the results for the cases of Low-, Medium-, and High-altitude orbits, respectively.Figure 2 uses a solid line to represent the full elliptic restricted three-body problem, a dashed line to represent the singleaveraged model, and a dashed-dotted line to represent the double-averaged model.Regarding the overall behavior, the results are according to the expected from the literature (see [17,18]), and the increase of the initial inclination causes oscillations with larger amplitudes in the eccentricity.It is also clear that the averaged models are very accurate for the case where the perturbing body is in a circular orbit.This accuracy decreases with the increase of the eccentricity of the primaries for all values of the initial inclinations.
The use of second-order averaged models is not recommended for eccentricities of the primaries of 0.3 or larger.A larger expansion is required in this situation.In general, both averaged models have results that are much closer to each other than closer to those of the full model.It means that both averages tend to give the same errors, at least for eccentricities of the primaries up to 0.2.For eccentricity of 0.3 and above, the averaged models begin to show different behaviors, and those differences increase with this eccentricity.Regarding the comparison between both averaged models, the results depend on the value of the initial inclination.
The increase of the semimajor axis makes the spacecraft to stay more time at shorter distances of the Moon, so the shown the corresponding eccentricity of the perturbing body (  ).In the first column, the results are for the case of Low-altitude orbits.The second column is for the case of Medium-altitude orbits.The third column is for the case of High-altitude orbits.The lines are as follows: full elliptic restricted three-body problem (solid line) and single-averaged (dashed line) and double-averaged models (dashed-dotted line).For the colors,  0 = 30 ∘ (red line), 20 ∘ (blue line), 10 ∘ (pink line), and zero (orange line).The increase of the semimajor axis and the eccentricity of the primaries reduce the average distance from the spacecraft to the Moon, making the perturbations stronger.
third-body perturbations are stronger.This fact accelerates the dynamics of the system, and the period of the oscillation of the eccentricity is shorter when compared with the lowaltitude orbits.It is visible that high-altitude orbits have oscillations of eccentricity increased in amplitude with the increase of the eccentricity of the primaries and the averaged models are not able to predict this property.They have results that decrease in quality with the increase of the eccentricity of the primaries.The amplitude of the oscillations increased about ten times with respect to the previous lower orbits  shown the corresponding eccentricity of the perturbing body (  ).In the first column, the results are for the case of Low-altitude orbits.The second column is for the case of Medium-altitude orbits.The third column is for the case of High-altitude orbits.The lines are as follows: full elliptic restricted three-body problem (solid line) and single-averaged (dashed line) and double-averaged models (dashed-dotted line).For the colors,  0 = 41 ∘ (red line), 40 ∘ (orange line), 39 ∘ (blue line), and 38 ∘ (pink line).The increase of the semimajor axis and the eccentricity of the primaries reduce the average distance from the spacecraft to the Moon, making the perturbations stronger.
due to the increase of the perturbation effects.The period of the oscillations is also reduced as an effect of the increase of the perturbations.The increase of the eccentricity of the primaries has the same effect of increasing the amplitudes and reducing the period of oscillations.The evolutions of the inclinations are not shown here because they remain constant for all of the situations considered in Figure 2.

Near Critical Inclination
Orbits.Figures 3 and 4 show the results of the evolutions of the inclination and the eccentricity, respectively, as a function of time for initial inclinations  0 = 41 ∘ (red line), 40 ∘ (orange line), 39 ∘ (blue line), and 38 ∘ (pink line).The same representations for the models are used here.The general behaviors are also according to the expected from the same literature cited above, and the increase of the initial inclination causes oscillations with larger amplitudes.Both averaged models have results that have good accuracy for values of   up to 0.2.For eccentricities of 0.3 and above, the averaged models start to have results that are not so good, and those deviations from the real cases increase with the Figure 4: These figures show the evolution of the eccentricity as a function of time for near critical inclination spacecraft orbits.The first, second, and third rows correspond to the results for eccentricity of the perturbing body (  ) 0.0, 0.1, and 0.2.In the first column, the results are for the case of Low-altitude orbits.The second column is for the case of Medium-altitude orbits.The third column is for the case of High-altitude orbits.The lines are as follows: full elliptic restricted three-body problem (solid line) and single-averaged (dashed line) and double-averaged models (dashed-dotted line).For the colors,  0 = 41 ∘ (red line), 40 ∘ (orange line), 39 ∘ (blue line), and 38 ∘ (pink line).The increase of the semimajor axis and the eccentricity of the primaries reduce the average distance from the spacecraft to the Moon, making the perturbations stronger.
eccentricity of the primaries, so they are omitted here.It is clear that, for values of the eccentricity of the primaries up to 0.1, both averaged models are very accurate.When   = 0.2, the single-averaged model gives more accurate results for the evolution of the eccentricity of the perturbed body, while the double-averaged model gives more accurate results for the evolution of the inclination.The increase of   also does not change the dynamics of the system in a noticeable form for low-and medium-altitude orbits.For high-altitude orbits, there is a noticeable acceleration of the dynamics.
It can be noticed that near the time 20000 canonical units, medium-altitude orbits have presented some deviations between the results of the full and the averaged models, in particular for the higher inclinations.This is a result of the increase of the semimajor axis from low to medium orbits, which places the spacecraft in an orbit that is much more perturbed by the third body.Again, those models begin to give results that are also not so accurate when the eccentricity of the primaries increases.It is visible that the two-second order averaged models have results that are very closer to each other with similar deviations from the full model.So, in this range of inclinations and altitudes, both averaged models have about the same quality of results.The increases of the semimajor axis accelerate the dynamics also in this situation.Note that the time scale of the axis is different for low-altitude earth orbits and the deviations of the averaged models occur very early.
The third columns of Figures 3 and 4 show that at higher altitude both averaged models do not fit perfectly even for the circular case for the perturbing body (  = 0).It can be noticed that, near the time 8000 canonical units, some deviations between the full and the averaged models occur, in particular for the higher inclinations.This is a result of the increase of the semimajor axis, which places the spacecraft in an orbit that is much more perturbed by the third body because it is closer to the Moon.The period of the oscillations is also reduced by this stronger perturbation.Again, those models begin to give results that are also not so accurate when ).In the first column, the results are for the case of Low-altitude orbits.The second column is for the case of Medium-altitude orbits.The third column is for the case of High-altitude orbits.The lines are as follows: full elliptic restricted three-body problem (solid line), single-averaged (dashed line), and double-averaged models (dashed-dotted line).For the colors,  0 = 80 ∘ (red line), 70 ∘ (blue line), 60 ∘ (green line), 50 ∘ (pink line), and 40 ∘ (orange line).The increase of the semimajor axis and the eccentricity of the primaries reduce the average distance from the spacecraft to the Moon, making the perturbations stronger.
the eccentricity of the primaries increases, not predicting correctly the oscillations performed by the full model.At high altitudes, both averaged approximations begin to lose quality for   = 0.2.It is visible that the two second-order averaged models have results that are very closer to each other with similar deviations from the full model.So, in this range of inclinations and altitudes, both averaged models have about the same quality of results for values of the eccentricity of the primaries below 0.2.

Above Critical Inclination
Orbits.Figures 5 and 6 show the evolutions of the inclinations and the eccentricities as a function of time for initial inclinations above the critical value.The general behaviors are the expected ones (see [17,18]), and the inclinations remain constant for a long time, return very fast to the critical value, and then increase again to its original values.The eccentricity shows an opposite behavior and increases when the inclination decreases and vice versa.The increase of the initial inclinations also causes oscillations with larger amplitudes.Both averaged models have results that are excellent for   = 0.0, and they still have good accuracy for eccentricities of the primaries up to 0.2.
The single-averaged model for low-altitude orbits, in an average over the time, gives results that are closer to those of the full model.This fact is very clear for eccentricities of the primaries 0.2.It is also noted that, for   = 0.4, the values given by the full model for the inclination and eccentricity lie in the middle of the curves described by the averaged models, showing a decrease to the critical value later than the decrease predicted by the single-averaged model but before the predictions made by the double-averaged model.This fact is due to the truncation of the expansion.The results show that the effect of this truncation depends on the initial conditions of the simulations, so the effects of the neglected terms are different for different orbits.For this eccentricity of primaries (  = 0.4), the second-order averaged models are no longer accurate enough to study this problem.The accuracy of the averaged models decreases with the increase of the initial inclination.Note that the averaged models are still good for   = 0.4 in the case  0 = 40 ∘ , but this is not so Figure 6: These figures show the evolution of the eccentricity as a function of time for high-inclination orbits.The first, second, and third rows correspond to the results for eccentricity of the perturbing body (  ) 0.0, 0.2, and 0.4.In the first column, the results are for the case of Lowaltitude orbits.The second column is for the case of Medium-altitude orbits.The third column is for the case of High-altitude orbits.The lines are as follows: full elliptic restricted three-body problem (solid line), single-averaged (dashed line), and double-averaged models (dasheddotted line).For the colors,  0 = 80 ∘ (red line), 70 ∘ (blue line), 60 ∘ (green line), 50 ∘ (pink line), and 40 ∘ (orange line).The increase of the semimajor axis and the eccentricity of the primaries reduce the average distance from the spacecraft to the Moon, making the perturbations stronger.
good for higher values of the initial inclinations, with the worst results happening for  0 = 80 ∘ .Another result that comes from those figures is that the increase of   causes an acceleration of the dynamics of the system.Note that the time required by the inclination and eccentricity to suffer the largeamplitude oscillation becomes shorter when the primaries are in a more elliptic orbit.For medium altitudes, when   = 0.4, the accuracy of the averaged models increases with the increase of the initial inclinations.Note that the averaged models are still good for the case of initial inclination 80 ∘ , and the situation with the worst results occurs for 60 ∘ .Inclined orbits are less perturbed by a third body because the mean distance between the spacecraft and the Moon is larger than in planar orbits.
Note that the strong changes in inclination and eccentricity occur earlier as a result of the stronger perturbations resulting from the increase of the semimajor axis of the orbit.
The increase of the eccentricity of the primaries has the same effect, showing that an elliptic orbit for the perturbing body causes more perturbation on the spacecraft, when compared with circular orbits with the same semimajor axis, because the distance between the bodies decreases at the periapsis.The correspondent increase of that distance at the apoapsis does not compensate the previous aspect.For values of   up to 0.2, the region where the second-order averaged models work better and the differences between the full model and the averaged models are larger in higher orbits.This shows that the increase of the perturbation effects also increases the differences between the models.
For high-altitude orbits, both averaged models have results that are excellent for   = 0.0 and that still have good accuracy for eccentricities of the primaries up to 0.2, at least for times up to 17000 canonical units.After that, the approximations are no longer good enough to study this problem under the second-order averaged hypotheses.In general, in this situation, the single-averaged model gives better results.For   = 0.1 and higher, the averaged models do not predict very well the time for the second jump of the inclination and eccentricity for all of the inclinations studied.The eccentricity of the primaries still accelerates the system, reducing the period of the oscillations, in the same way made by the higher altitudes.Note that, at this high altitude, the system shows two jumps for the inclination and eccentricity.

Conclusions
The equations of motion of a spacecraft perturbed by a third body using a single-averaged technique are developed, considering an elliptic orbit for the disturbing body.
Looking at the overall behavior, the results are according to the literature: short oscillations in the inclination and eccentricity for initial inclinations below the critical value, which increases fast in amplitude around the critical value, then the typical behavior of having the inclinations remaining constant for a long time, and then returning very fast to the critical value, to increase fast again to its original values.The eccentricity shows an opposite behavior and increases when the inclination decreases and vice versa.
The results also showed that circular orbits exist, but frozen orbits do not exist under this model.The circular orbits, in general, do not keep the inclination constant, as happened in the double-averaged model.
The increase of the semimajor axis causes an increase of the third-body perturbation, and this fact accelerates the dynamics of the system.It was noticed that the period of the oscillations of the inclination and the eccentricity are shorter when compared with lower-altitude orbits.The increase of the eccentricity of the primaries, considering the semimajor axis constant, has the same effect of accelerating the dynamics.So, elliptic orbits for the disturbing body have the effect of increasing the perturbation, if all other elements remain the same.
It was also showed that inclined orbits are less perturbed by a third body, since the mean distance between the spacecraft and the Moon is larger than in planar orbits.
A comparison with the double-averaged technique up to the second order and with the full restricted elliptic problem was made, and it showed some other facts.The second-order averaged models are very accurate when the perturbing body is in a circular orbit.This accuracy decreases with the increase of the eccentricity of the primaries.The use of second-order averaged models is not recommended for eccentricities of the primaries of 0.3 or larger, and a better expansion, including more terms, is required in this situation.In general, both second-order averaged models have results that are much closer with each other than closer to those of the full model.It means that both averages tend to give similar errors, when compared with the full model.

Figure 1 :
Figure 1: Illustration of the dynamical system.

Figure 2 :
Figure2: These figures show the evolution of the eccentricity as a function of time for low-inclination orbits.On the top of each figure is shown the corresponding eccentricity of the perturbing body (  ).In the first column, the results are for the case of Low-altitude orbits.The second column is for the case of Medium-altitude orbits.The third column is for the case of High-altitude orbits.The lines are as follows: full elliptic restricted three-body problem (solid line) and single-averaged (dashed line) and double-averaged models (dashed-dotted line).For the colors,  0 = 30 ∘ (red line), 20 ∘ (blue line), 10 ∘ (pink line), and zero (orange line).The increase of the semimajor axis and the eccentricity of the primaries reduce the average distance from the spacecraft to the Moon, making the perturbations stronger.

Figure 3 :
Figure3: These figures show the evolution of the inclination as a function of time for near critical inclinations.On the top of each figure is shown the corresponding eccentricity of the perturbing body (  ).In the first column, the results are for the case of Low-altitude orbits.The second column is for the case of Medium-altitude orbits.The third column is for the case of High-altitude orbits.The lines are as follows: full elliptic restricted three-body problem (solid line) and single-averaged (dashed line) and double-averaged models (dashed-dotted line).For the colors,  0 = 41 ∘ (red line), 40 ∘ (orange line), 39 ∘ (blue line), and 38 ∘ (pink line).The increase of the semimajor axis and the eccentricity of the primaries reduce the average distance from the spacecraft to the Moon, making the perturbations stronger.

Figure 5 :
Figure5: These figures show the evolution of the inclination as a function of time for high-inclination spacecraft orbit.On the top of each figure is shown the corresponding eccentricity of the perturbing body (  ).In the first column, the results are for the case of Low-altitude orbits.The second column is for the case of Medium-altitude orbits.The third column is for the case of High-altitude orbits.The lines are as follows: full elliptic restricted three-body problem (solid line), single-averaged (dashed line), and double-averaged models (dashed-dotted line).For the colors,  0 = 80 ∘ (red line), 70 ∘ (blue line), 60 ∘ (green line), 50 ∘ (pink line), and 40 ∘ (orange line).The increase of the semimajor axis and the eccentricity of the primaries reduce the average distance from the spacecraft to the Moon, making the perturbations stronger.