Third-Body Perturbation Using a Single Averaged Model : Application in Nonsingular Variables

The Lagrange’s planetary equations written in terms of the classical orbital elements have the disadvantage of singularities in eccentricity and inclination. These singularities are due to the mathematical model used and do not have physical reasons. In this paper, we studied the third-body perturbation using a single averaged model in nonsingular variables. The goal is to develop a semianalytical study of the perturbation caused in a spacecraft by a third body using a single averaged model to eliminate short-period terms caused by the motion of the spacecraft. This is valid if no resonance occurs with the moon or the sun. Several plots show the time histories of the Keplerian elements of equatorial and circular orbits, which are the situations with singularities. In this paper, the expansions are limited only to second order in eccentricity and for the ratio of the semimajor axis of the perturbing and perturbed bodies and to the fourth order for the inclination.


Introduction
Garofalo [1] and Pines [2] used several methods for removing the singularities caused by the use of the classical orbital elements.Cohen and Hubbard [3] described a new set of nonsingular elements, which are relatively simple functions of the classical elliptic orbital elements.They also give the differentials of the coordinates and velocity components in terms of the differentials of the nonsingular elements.Broucke and Cefola [4] studied the equinoctial orbit for the two-body problem, showing that the associated matrices are free from singularities for zero eccentriciticies and 0 • and 90 • inclinations.Giacaglia [5] removed the singularities from the geopotential expansion and its derivates, but allowed a mixed set of orbital elements to remain in the expansions.However, its development retains much of the original inclination and eccentricity functions, allowing their calculation by existing recursive relations.Nacozy and Dallas [6] removed singularities from the geopotential and its derivates for zero eccentricity and inclination and formulated the geopotential expansion entirely in terms of nonsingular orbital elements.Walker et al. [7] introduced a set of modified equinoctial orbit elements for perturbation analysis of all types of orbits.Wytrzyszczak [8] used nonsingular elements for the description of the motion of satellites with small eccentricity and inclination.
This work (Wytrzyszczak [8]) is an extension of the research performed by Nacozy and Dallas.Those researches present rich contributions and have a strong analytical approach, containing several derivations of equations.In the present paper, we developed the equations made by Giacaglia [5] and applied an averaged technique to obtain new equations for the evolution of the orbital elements in nonsingular variables.Using those equations, we made numerical simulations to understand the evolution of circular and equatorial orbits, plotting the results in terms of the standard orbital elements, in order to complement previous studies of this problem that used singular variables.It is aligned with papers that have emphasis in numerical results, like the ones by Broucke [9] that developed a semianalytical study up to the second-order theory, Prado [10] that extended this theory to the fourth order, and Sol órzano and Prado [11] that considered a single averaged model for the third-body perturbation.

Mathematical models
The model used in the present paper can be formulated in a very similar way to the planar restricted three-body problem.
There are three bodies involved in the dynamics: one body with mass m 0 fixed in the origin of the reference system, a second massless body m in a three-dimensional orbit around m 0 , and a third body (m ) in a circular orbit around m 0 (see Figure 2.1).
The motion of the spacecraft (the second massless body) is Keplerian and threedimensional, with its orbital elements perturbed by the third body.The main body m 0 is fixed in the center of the reference system X-Y -Z.The perturbing body m is in a circular orbit that stays in the referential plane (equatorial plane).The canonical system of units is used here, which means that the period of the disturbing body (moon) is 2π, which is equivalent to 27.322 days, and the semimajor axis of the moon is normalized to 1 (see Table 2.1).
The Lagrange planetary equations have the disadvantage of singularities in the eccentricity and inclination due to the presence of those quantities in the denominators of the equations.These singularities are mathematical and not physical.The presence of the eccentricity in the denominator causes problem in the evaluation of the variations of the argument of periapsis.The equations that depend on the inclination have a special behavior for small inclinations.Also, the longitude of the ascending node is not defined for inclination equal to zero.The circular orbits do not have argument of periapsis.It is possible to express these equations in terms of nonsingular elements, with the goal of eliminating the singularities.However, the set of nonsingular elements is not unique.In this paper, the equations used for the motion of the spacecraft are the ones determined by Giacaglia [5] using nonsingular elements and considering that i = π and e < 1.The following set of nonsingular elements are used: = Ω + ω. (2.1) The meaning of the symbols are as follows: (i) λ : mean longitude, (ii) M : mean anomaly, (iii) ω : argument of periapsis, (iv) Ω : longitude of the ascending node, (v) e : eccentricity, (vi) i : inclination, (vii) : longitude of pericentre, (viii) M 0 : initial mean anomaly, (ix) n : mean motion, (x) t : time, (xi) f : true anomaly.As a convention, primed variables are used for the perturbing body.
The Lagrange planetary equations in terms of the nonsingular elements are as follows (Giacaglia [5]): (2. 2) The disturbing function due to the third body is (Giacaglia [5]) Using equatorial coordinates for the satellite (α,δ) and the perturbing body (α ,δ ), we have (Giacaglia [5]) cos(S) = sin(δ)sin(δ ) + cos(δ)cos(δ )cos(α − α ), It is known that ζ 0 = 1 and ζ z = 2 for z = 0. From the Rodrigue formula, it is possible to obtain the associated Legendre function P ᐉz by Then, the harmonic coefficients are defined as where (Giacaglia [5]) In the disturbing function, F ᐉzp (i) and J(c) ᐉzp are the inclination functions.Those equations are in polynomial form, so its accuracy depends on the degree considered for polynomial.H ᐉpq are the Hansen coefficients that are functions of L ᐉpq (γ), which are functions introduced by Giacaglia [5], but that also depend on the power series for the eccentricity function.Equation (2.8) in orbital coordinates (Giacaglia [5]) can be written as ) (2.10) It is possible to relate (2.8) and (2.9) and also (2.7) and (2.8): (2.12) Using the Hansen coefficients H ᐉpq (Giacaglia [5]) obtained, r a ᐉcos sin where The functions L ᐉpq (Giacaglia [5]) are power series of γ (γ = √ 1 − e 2 ).The disturbing function can be written as (2.15) (2.18) Taking into account that the perturbing body is on the reference plane and that its orbit is circular (α = M 0 + n t, δ = 0), it is possible to write that   (Giacaglia [5]).(Giacaglia [5]). where Using Table 2.2 and the inclination functions F ᐉzp (i), J(c) ᐉzp (see Table 2.3), it is possible to have all the terms of the expansion of (2.20), where  (2.24) 10 Mathematical Problems in Engineering The definition of average is After writing (2.24) in terms of the nonsingular elements and using (2.25) to eliminate the short-periodic terms of the spacecraft motion (mean longitude), the final result for the disturbing function is (2.26)

Numerical results
This section shows several simulations for the case of equatorial and/or circular orbits.We studied the effect of the perturbations of the moon on high-altitude earth satellites, with semimajor axis of 0.070 and 0.110 canonical units.For the equatorial orbits, Figure 3.1 shows that the behavior of the eccentricity is constant for the time scale of the plot.They remain invariant for several values of the eccentricity.Orbits with semimajor axis of 0.110 (see Figure 3.2) show the same behavior.Our simulations, for several initial inclinations show a constant behavior (see Figure 3.3).The simulations for near equatorial orbits using several values for the initial eccentricity is shown in Figure 3.4 for a semimajor axis of 0.341 canonical units and an initial inclination near to 5 degrees.Looking at these results, we see that the behavior of the eccentricity shows oscillations with a very small amplitude that does not affect the stability of the orbit.In general, the inclination shows its oscillatory behavior.This fact is also reported in the literature (Broucke [9], Prado [10], and Sol órzano and Prado [11]).

Conclusions
The present research considered a semianalytical method to determine the effects of the disturbance of the third body (moon) on a spacecraft in equatorial and/or circular orbits using a single averaged nonsingular model.With the help of the inclination functions and the Hansen coefficients, the disturbing function was obtained.These inclination and eccentricity functions are polynomials in "s" and "e," so the accuracy is governed by the degree of these polynomials.Then, we used the Lagrange planetary equations in nonsingular variables and expanded the disturbing function until the second order in the eccentricity.The inclination was expanded until the fourth order, these terms appear on the expansion of the polynomials F ᐉzp(i) , J(c) ᐉzp .However, due to the fact that the series was truncated for ᐉ = 2, we leave out those terms.The results show small variations for the orbital elements.However, it is important to remember that the expansion of the disturbing function is not convergent for large values of the eccentricity.This limit is not important in our case because the moon moves in circular orbit and the spacecraft is in a low eccentric orbit.The results showed that the equatorial orbits remain with the eccentricity constant, for both values of the semimajor axis tested.It is also visible that, when the inclination leaves the zero value, the evolution of the eccentricity and inclination shows sinusoidal variations.In the ideal case of orbits that starts with zero eccentricity, its eccentricity remains always zero.This occurs because the right-hand side of the equation for the time derivate of the eccentricity is zero (Lagrange's equation) because it is a polynomial in the eccentricity with no independent term.Another property of those orbits is that the inclination is also constant for the same reason.However, these same behaviors are expected for orbits with small eccentricities.For the case of equatorial orbits, the inclination and eccentricity remain constant, and the orbits remain in the equatorial plane.The problem of the critical inclination was studied by several authors (Aoki [12,13] and Jupp [14,15]).However, Broucke [9] and Prado [10] studied the presence of a critical value for the inclination between the perturbed and the perturbing bodies.This critical inclination is related to the stability of near-circular orbits.Thus, if the inclination is higher than the critical value, the eccentricity increasses and the near-circular orbits become very elliptic.Alternatively, if the inclination is lower than this critical value, the orbit stays nearly circular.This critical value is 0@39,7 • and it represents the limit inclination that allows the existence of orbits with eccentricity, inclination, and argument of periapsis constant under the second-order model.
μ = m /(m 0 + m), m : mass of disturbing body; m 0 mass of the main body, (ii) a : semimajor axis of the disturbing body, (iii) S: angle between the position vectors of the perturbing and the perturbed body, (iv) r : geocentric distance of the perturbing body, (v) r: geocentric distance of the perturbed body, (vi) δ: declination, (vii) α: right ascension.C. R. H. Sol órzano and A. F. B. A.Prado 5