Semianalytic Integration of High-Altitude Orbits under Lunisolar Effects

The long-term effect of lunisolar perturbations on high-altitude orbits is studied after a double averaging procedure that removes both the mean anomaly of the satellite and that of the moon. Lunisolar effects acting on high-altitude orbits are comparable in magnitude to the Earth’s oblateness perturbation. Hence, their accurate modeling does not allow for the usual truncation of the expansion of the third-body disturbing function up to the second degree. Using canonical perturbation theory, the averaging is carried out up to the order where second-order terms in the Earth oblateness coefficient are apparent. This truncation order forces to take into account up to the fifth degree in the expansion of the lunar disturbing function. The small values of the moon’s orbital eccentricity and inclinationwith respect to the ecliptic allow for some simplification. Nevertheless, as far as the averaging is carried out in closed form of the satellite’s orbit eccentricity, it is not restricted to low-eccentricity orbits.


Introduction
In an increasingly saturated space about the Earth, aerospace engineers confront the mathematical problem of accurately predicting the position of Earth's artificial satellites.This is required not only for the correct operation of satellites but also for preserving the integrity of space assets.Thus, operational satellites are threatened by the not so remote possibility of a collision with a defunct satellite 1 but most probably by the impact with other uncontrolled man-made space objects as spent rocket stages or collision fragments-all of them commonly called space debris.
Precise predictions require the integration of complete force models including both gravitational and nongravitational effects, like a high-degree and order geopotential, ephemerides-based lunisolar perturbations, drag, solar radiation pressure including eclipses, and so forth see 2, 3 , for instance .The most accurate integration is expected from numerical methods, although precision ephemeris can also be obtained by means of semianalytical integration 4 .In fact, both approaches, numerical and semi-analytical, do not need to enter a competition.Thus, for instance, while semi-analytical methods may be efficient in keeping a running catalog of hundreds of thousands of space objects within a reasonable accuracy, if the probability of impact with an operational satellite is detected to surpass a certain threshold, then a more accurate numerical integration will help control engineers in deciding whether a collision avoidance maneuver is required or, on the contrary, the integrity of the satellite is not in risk 5 .
In a semi-analytical approach the highest frequencies of the motion, which normally have small amplitudes, are filtered analytically via averaging procedures.This filtering allows the numerical integration of the averaged system to proceed with very long step sizes.Then, the short-period terms are recovered analytically, if desired, at any step of the numerical integration 6-8 .
Averaging techniques are also used for exploring questions affecting stability, such as those derived from tesseral resonances or third-body perturbations, in a reduced-phase space 9 .In this respect, much attention has been recently paid to the long-term evolution of classical GNSS constellations, either for operational or disposal orbits 10 .
While the noncentralities of the Earth gravitational potential play a key role in the motion of low altitude satellites, third-body perturbations have also a decisive influence in the long-term evolution of medium-and high-altitude Earth orbits.The third-body disturbing function is commonly given by a series expansion in Legendre polynomials.Often, the series is truncated to the first term in the expansion 11 , but this early truncation is not always accurate enough to represent the real dynamics 4, 12, 13 .Nevertheless, recursions in the literature allow to extend the Legendre polynomial expansion to any desired order either in classical or nonsingular elements 14-16 .
Because of the physical characteristics of the Earth gravitational potential, where the second-order zonal coefficient J 2 clearly dominates over all other harmonics, second-order effects of J 2 may be important and must be taken into account when the effect of higherorder harmonics is studied.Correspondingly, the truncation in the expansion of third-body perturbations must include terms of magnitude comparable to J 2 -squared.Because the thirdbody disturbing function is expanded in the ratio semi-major axis of the satellite's orbit to semi-major axis of the third-body's orbit, the degree at which the expansion must be truncated depends on the altitude of the satellites to be studied.
In this paper we study the effect of lunisolar perturbations on high-altitude orbits about a noncentral Earth, which is assumed to be oblate although without equatorial symmetry.More specifically, we are interested in the semi-analytical integration of satellites on altitudes of classical global navigation satellite system GNSS constellations such as GPS, Glonass, or Galileo.Note, however, that the assumption of an axisymmetric geopotential prevents to tackle the tesseral resonance problem that commonly suffer this kind of orbits.With respect to previous research 17 , where we approached the general case of third-body perturbations via double averaging, we release here the common simplification of assuming the third-body in circular orbit.Also, we focus on the case of Earth's artificial satellites dealing with more real lunisolar perturbations.We do that because recent results 18 seem to contradict the claim that taking the third-body in circular orbit does not produce any noticeable degradation in the long-term propagation of real Earth orbits 19 .Besides, for the actual values of the orbits of the sun and moon, neglecting terms in the eccentricity is not consistent with a higher-order expansion of the lunisolar disturbing function.Hence, both the eccentricity of the orbits of the sun and moon are taken into account in the present work.For the latter, the moon, the orbit is assumed to remain with constant inclination with respect to the ecliptic plane, over which the longitude of the ascending node moves with linear motion.The argument of the perigee of the moon is also assumed to evolve linearly, while we take the apparent orbit of the sun to be purely Keplerian.
We use canonical perturbation theory by Lie transforms 20-22 .The order of each term of the disturbing function is determined by a virtual small parameter that is taken proportional to the ratio of the satellite's orbit semimajor axis to the moon's orbit semi-major axis.Then, we check that the magnitude of second order terms due to the Earth oblateness is comparable to that of the fifth degree in the expansion of the moon's third-body potential, and hence we truncate the moon's third-body potential up to the fifth degree.On the contrary, the usual truncation up to the second degree is enough for modeling the sun disturbing function assumed that we limit the theory to the order of J 2 -squared terms.The Hamiltonian model also takes into account the asymmetries in the Earth potential caused by the J 3 term, whose influence in the dynamics of the lower eccentricity orbits is clearly noted.
The initial Hamiltonian is of three degrees of freedom and time dependent.Note, however, that the time appears in different scales related to the mean motion of the moon, that of the sun, and the rate of variation of the moon's argument of the perigee and longitude of the ascending node.In our averaging procedure, we eliminate the mean anomaly of the satellite and that of the moon to obtain a two-degrees-of-freedom Hamiltonian, which still depends on time albeit only through very slowly varying quantities.As far as we do not find resonances between the mean motions of the satellite and the moon, the averaging can be done in two steps: the mean anomaly of the satellite is removed first by means of a classical Delaunay normalization 23 ; then, a similar transformation is used to eliminate the mean anomaly of the moon.The splitting of the averaging has the advantage of simplifying computations.Specifically, the generating function of each transformation is obtained from the solution of simple quadratures.
The averaging is carried out not only in closed form of the satellite's orbit eccentricity, thus making the simplified Hamiltonian useful for studying the long-term evolution of any orbit, but also in closed form of the eccentricity of the moon's orbit.Nevertheless, because of the small values of the eccentricity of the orbits of both perturbing bodies as well as the small inclination of the orbit of the moon with respect to the ecliptic, further simplifications can be done by neglecting terms of order higher than the truncation order of the theory.
Numerical experiments using the doubly averaged Hamiltonian demonstrate that the inclusion in the model of the orbital eccentricities of the sun and moon does not cause qualitative differences with respect to the circular orbit approximation.Besides, the recovery of the short-period effects by means of the analytical transformation equations of the averaging provides a quite reasonable precision in the long-term integrations.

Model
We study the motion of an Earth artificial satellite of negligible mass whose Keplerian motion is slightly distorted by the noncentralities of the geopotential due to the Earth's oblateness and latitudinal asymmetry, as represented by the harmonic coefficient J 2 and J 3 , respectively, and under the point-mass attraction of the sun and moon.Thus, the motion of the satellite is derived from the potential where μ is the Earth's gravitational parameter, α its equatorial radius, J i is the zonal harmonic coefficient of degree i, r is distance from the origin, and z is the satellite's coordinate in the direction of the symmetry axis of the Earth.In the mass-point approximation, the third-body disturbing potential V , ∈ , , is where μ is the third-body's gravitational parameter, r and r are the radius vector of the satellite and of the third-body, respectively, of corresponding modulus r and r .If the disturbing body is far away from the perturbed body, 2.2 can be expanded in power series of the ratio r/r : where β m / m m , n is the mean motion of the third-body in its orbit of semi-major axis a , P j are Legendre polynomials, and ψ is the angle encompassed by r and r .The absence of the Keplerian term −μ /r as a summand in 2.3 has no effect in the restricted problem approximation, in which the mass of the satellite has negligible effects on the third-body's orbit.
In our model we assume an Earth-centered frame with the origin defined by the intersection of the Earth's equator and the ecliptic; we both of which consider fixed planes.Whereas we take the sun's apparent orbit about the Earth to be purely Keplerian, we assume that the moon moves with Keplerian motion on a precessing ellipse that evolves with constant rate of the argument of perigee.Besides, the moon's orbital plane is assumed to have constant inclination with respect to the ecliptic; yet it is also assumed to be precessing over the ecliptic with linear variation of the longitude of the ascending node.
Calling a, e, i, ω, Ω, and M to semimajor axis, eccentricity, inclination, argument of the perigee, longitude of the ascending node, and mean anomaly, respectively, we use the approximate values We recall that J and N are referred to the ecliptic.Besides, β 1 for the sun and β 1/28.8245 for the moon.
Remark that this simple model is not, of course, valid for computing precise ephemeris of the moon, who may be clearly out of the predicted position because of the amplitude of periodic oscillations in N and J, and also in e and i , comper 24 .Nevertheless, it will be enough for our purposes of investigating the qualitative features that lunisolar perturbations introduce in the long-term behavior of Earth's artificial satellites.

Perturbation Approach
By using canonical variables we can study the problem in Hamiltonian form.Besides, in order to apply perturbation theory, we arrange the Hamiltonian as a power series in a small parameter: We are interested in high-altitude orbits in the range of 20 000 to 35 000 km, the socalled upper mean earth orbit MEO region.In order to take account of all relevant lunisolar perturbations, we take a virtual small parameter proportional to the ratio semimajor axis of the satellite semimajor axis of the moon.For the altitudes of interest, this ratio is δ ∼ O 10 −1 .Then, we find that the Earth's oblateness coefficient is a fourth-order quantity, and the J 3 effect appears at the eighth order.With respect to the moon the consecutive terms of the Legendre polynomials expansion of the moon disturbing function appear at consecutive orders of the Hamiltonian, starting by the fifth order.We include up to the fifth degree, whose effect may be comparable to second-order effects of J 2 .For the sun, it is enough to take the first term in the Legendre polynomials expansion.Then, the zero-order term is the Keplerian Besides, H i,0 0 for i < 4, and where f is the true anomaly and η √ 1 − e 2 is the usual eccentricity function.The equations of motion are obtained from the Hamilton equations: Recall that the Cartesian coordinates of the satellite are expressed in orbital elements by means of simple rotations.Thus, where θ f ω is the argument of latitude and R 1 and R 3 are the usual rotation matrices about the xand z-axes, respectively.Since the Hamiltonian must be expressed in canonical variables, we assume hereafter that the orbital elements of the satellite are always expressed as function of the Delaunay variables , g, h, L, G, H , given by the known relations M, g ω, h Ω, L √ μ a, G L η, and H G cos i.Nevertheless, we use the orbital elements notation because of its immediate physical insight.Delaunay's variables are singular for zero eccentricity and/or zero inclination.In spite of that, the Lie transforms theory is naturally computed in Delaunay variables.Once the generating function of the averaging is computed, it can be applied to any function of the Delaunay's variables.Specifically nonsingular variables are assumed to be functions of Delaunay variables to obtain the transformation equations of the averaging in nonsingular variables 25 .
We note that the sun and moon coordinates are assumed to be known functions of time.Therefore, 3.1 is a time-dependent Hamiltonian of three degrees of freedom.In our perturbation approach, we avoid dealing with time by moving to an "extended" phase space see, for instance, 26 .Since the time-dependent variables evolve in quite different time scales, we find convenient to introduce four new pairs of canonical conjugate variables: l, L for the mean anomaly of the sun and its conjugate momentum λ, Λ the same for the moon, w, W for the moon argument of perigee and its conjugate momentum, and N, Γ for the moon longitude of the ascending node and its conjugate momentum.The specific values of the introduced canonical momenta are irrelevant for our purposes, and we find convenient to make the ordering Once we have set the Hamiltonian order, we will apply perturbation theory by the Lie transforms in order to filter the short-period effects in the potential equation 2.1 .More precisely, we base our computations on Deprit's algorithm, which is specifically designed for automatic computation by machine 21, 27 .

The Delaunay Normalization
We compute the canonical transformation T : , g, h, L, G, H → , g , h , L , G , H that removes the mean anomaly of the satellite from the Hamiltonian up to the eighth order.Note that the mean anomaly does not appear explicitly but through its dependence on the true anomaly.This fact introduces some subtleties in the averaging procedure, but the usual differential relations between true, mean, and eccentric anomalies allow to carry out the normalization on closed form of the eccentricity.
After normalization, we get the new Hamiltonian where the H 0,m terms are expressed in the prime variables.We find H 0,m H m,0 for m < 4, and where m/2 is an integer division.The eccentricity coefficients E, the inclination ones p, and the third-body direction coefficients t and d are given in Tables 1, 2, and 3. Note that except for J 2 -squared terms the new Hamiltonian is obtained by the simple average over the mean anomaly of all the terms in 3.1 .Besides, we introduced Kozai's arbitrary constant 28 in the solution of the fourth-order generating function to keep the prime variables as close as possible to the average value of corresponding osculating ones.
The semimajor axis of the satellite remains constant after averaging, as well as its canonical partner, the Delaunay action L .Then, the time evolution of the mean anomaly decouples from the two-degrees-of-freedom, time-dependent system The numerical integration of this system can be done with longer step sizes than the original one because of the filtering of short periodic effects via averaging.At each step of the numerical integration, the osculating elements can be recovered analytically using the transformation equations computed also by the Lie transforms method.

Elimination of the Mean Anomaly of the Moon
A new canonical transformation T : , g , h , L , G , H → , g , h , L , G , H removes the mean anomaly from the Hamiltonian.The algorithm starts now by setting the new Hamiltonian where K m,0 H 0,m .Then, we express the direction of the moon x , y , z in terms of the moon's orbital elements by the composition of rotations: where θ f ω .The components of the moon direction vector are then replaced in coefficients t , d of Table 1.Finally, application of the Lie triangle provides the new Hamiltonian where the K 0,m terms are expressed in the double prime variables, as well as the generating function of the transformation.
The new averaging is similar to the preceding Delaunay normalization in the sense that we base on the differential relations between the true and mean anomalies to perform the averaging.As before, up to the computed order in the perturbation theory, there is no coupling between the different Hamiltonian terms, which, therefore, reduce to their averaging over the mean anomaly of the moon λ.Thus, K 0,m K m,0 for m < 5, and

3.16
Table 4: Averaged moon direction coefficients t j,k .Here, c ≡ cos J, s ≡ sin J, and e ≡ e .

3.18
To get some further simplification, we note that e ∼ i ∼ O δ .Thus, in our eighthorder theory, we neglect higher-order terms factored by e m sin n J δ k such that m n k > 8.
Under these simplifying assumptions, the coefficients t m,i and d m,i are presented in Tables 4  and 5, respectively.
As in the preceding averaging, adequate arbitrary constants have been introduced in the computation of the generating function to guarantee its average to zero.
After the double averaging, the Hamiltonian only depends on long-period terms related to the sun's apparent motion and very long-period terms related to the precession of the nodes and recession of the line of apsides of the moon's orbit.The numerical integration of corresponding Hamilton equations is, now, very much faster and efficient.

Numerical Experiments
For the numerical tests we used a higher-order Runge-Kutta method.Specifically, the numerical integration was performed with the Dormand and Prince implementation of the Runge-Kutta method coded in FORTRAN 77 by Hairer et al. 29 .To illustrate the significance of recovering the short-period effects up to higher orders, we first show in Figure 1 a sequence of the errors obtained in the semimajor axis after one month semi-analytical propagation.For this example, we take the initial osculating elements a 28 560 km, e 0.2, i 56 deg., ω 0, Ω 72 deg., and 0; besides we assumed that both the mean anomalies of the moon and the sun are zero at the origin of time.Corresponding elements in the single-and double-averaged phase space depend on the truncation order of the perturbation theory and are presented in the top and bottom parts of Table 6, respectively.In reference to Figure 1, the top plot shows a direct comparison between the numerical integration of the Hamilton equations of the original, nonaveraged problemwhose disturbing potential is given in 2.1 -and that of the long-term Hamiltonian after the double averaging up to the eight order of the small parameter.Then, from top to bottom, we show the errors obtained at each step of the integration when recovering short-period terms by computing the transformation equations of the averaging up to the fourth, fifth, sixth, seventh, and eighth orders.For the latter, the amplitude of the periodic errors is reduced to several centimeter; note that, in addition to short-period errors related to the orbital period of the satellite, we can appreciate a two-week modulation related to the moon's mean motion.Remark that the amplitude of the periodic errors roughly divides by ten with each order of the transformation equations, which is consistent with the assumed magnitude of the virtual small parameter δ ∼ 0.1.

Mathematical Problems in Engineering
The short-period effects can be ignored in the study of the long-term orbital behavior, where the simple propagation of the double-averaged equations appears very fast and efficient.Sample propagations are shown in figures below, which show the important effects that have the initial right ascension of the ascending node and argument of perigee of the satellite's orbit in the long term and specifically in the eccentricity and inclination evolution of the satellite's orbit.Thus, Figure 2 shows the notably different evolution of the satellite's eccentricity and inclination for different initial nodes; for the other initial conditions we have taken the same as in the preceding short-term propagations but now assuming directly that they are mean elements.In fact, the more relevant parameter is the difference between the node of the satellite's orbit and that of the moon's orbit, as illustrated in Figure 3 where the initial longitude of the ascending node of the moon's orbit over the ecliptic has been taken as N 180 instead of N 0 deg.
Figure 4 shows that the effect of the initial argument of the perigee of the satellite's orbit is also important in the long term, eccentricity evolution.The effect is almost negligible in the orbital inclination in the long-term and it is not presented.
Finally, we must mention that further tests demonstrated that there are no important qualitative differences in the long term when using lower-order truncations of the theory, resulting in faster numerical integration of the mean elements.So the fifth-order truncation should be the preferred long-term Hamiltonian.Besides, we also checked for a variety of orbits that making e e 0 does not either introduce qualitative differences in the longterm, in agreement with 19 .

Conclusions
Modeling lunisolar perturbations on high-altitude Earth orbits requires to retain high degrees in the Legendre polynomials expansion of the third-body disturbing function of the moon.In consequence the eccentricity of the third-body's orbit cannot be neglected in the case of either the moon or the sun.
The long-term behavior of high-altitude Earth orbits is approached in a semianalytical way via averaging procedures, in which we take advantage of the different scales in which appears the time to do the averaging in the extended phase space.In addition, up to secondorder terms in the Earth's oblateness coefficient, the averaging has been computed in closed form of the eccentricity, and, therefore, the semianalytical integration can be applied to any orbit.
Sample numerical propagations of test cases show that the more relevant parameter on the long-term behavior is the difference between the right ascension of the ascending node of the satellite's orbit and the longitude of the ascending node of the moon's orbit over the ecliptic, having an apparent effect manifested by the almost-secular growing of the orbit eccentricity and also by very-long-period oscillations of the inclination with an amplitude of several degrees.Also the initial argument of the perigee of the satellite's orbit has notable effects in the satellite's orbit, but only in what to the long-term evolution of the eccentricity.

3 . 10 where
s and c stand for sine and cosine of the inclination, respectively; 2j k,i t m,i cos φ d m,i sin φ , φ 2j k ω iΩ, k m mod 2, ∈ , , 3.11

Figure 1 :
Figure 1:Semianalytical theory versus numerical integration.From top to bottom, errors in semimajor axis of the secular terms and after computing the fourth, fifth, sixth, seventh, and eighth order transformation equations, respectively.Abscissas are days.

Figure 3 :
Figure 3: The same as Figure 2 but now N 180 instead of N 0 deg.

Table 1 :
Direction coefficients t j,k and d j,k and eccentricity coefficients E j,k .The unit vector x, y, z defines the direction of the center of mass of the third body.

Table 5 :
Averaged moon direction coefficients d j,k .Here, c cos J, s sin J, and e ≡ e .

Table 6 :
Initial conditions in the single-top and double-averaged phase space bottom after different truncation orders n of the semi-analytical theory.Distances are km, and angles are deg.