# Averaging Tesseral Effects: Closed Form Relegation versus Expansions of Elliptic Motion

Authors: Martin Lara; Juan F. San-Juan; Luis M. López-Ochoa
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Longitude-dependent terms of the geopotential cause nonnegligible short-period effects in orbit propagation of artificial satellites. Hence, accurate analytical and semianalytical theories must cope with tesseral harmonics. Modern algorithms for dealing analytically with them allow for closed formrelegation. Nevertheless, current procedures for the relegation of tesseral effects from subsynchronous orbits are unavoidably related to orbit eccentricity, a key fact that is not enough emphasized and constrains application of this technique to small and moderate eccentricities. Comparisons with averaging procedures based on classical expansions of elliptic motion are carried out, and the pros and cons of each approach are discussed.


---

## Body

## 1. Introduction

Soon after Brouwer’s dream of a closed form solution to the artificial satellite theory (AST) to higher orders was achieved for the main problem of AST [1–3], the progress in mathematical methods made the closed form solution also possible for the general case of a tesseral potential. This notable success was feasible thanks to the invention of two algorithms, which are specifically suited for canonical simplifications of Hamiltonian problems, namely, the elimination of the parallax [4] and the relegation algorithm [5]. The former operates the tesseral averaging in closed form for low Earth orbits (LEO) [6–8] while the latter does the same in all other cases [9, 10]. However, many aerospace engineering applications of Earth artificial satellites are related to subsynchronous orbits, a case in which both algorithms limit their practical application to orbits with small or moderate eccentricities. Indeed, high-eccentricity orbits lead to impact with the surface of the Earth in the case of LEO, on the one hand, while the number of iterations required by the relegation algorithm makes it impracticable for high-eccentricity orbits, on the other [9].As reported in the literature, the relegation of tesseral effects from subsynchronous orbits is based on multiplications at each iteration by the ratio of the Earth’s rotation rate to the mean motion of the orbit and hence converges very slowly, except for low-altitude orbits. Another drawback of the standard algorithm, is that it assumes nonresonance conditions and finds a gradual deterioration in the solution when approaching synchronous condition [9]. On the contrary, usual series expansions of elliptic motion (see, [11–14], e.g.) allow for the effective isolation and handling of tesseral resonances [15–19].Both approaches, tesseral relegation and series expansions, have been implemented in analytical and semianalytical theories used in renowned laboratories. Thus, relegation is reported to be implemented in the Analytic orbit propagator program (AOPP) of the Naval Research Laboratory [20], whereas the Draper Semianalytical Satellite Theory (DSST) [21], these days maintained at MIT Lincoln Laboratory, or the numeric-analytical theory THEONA developed at Keldysh Institute of Applied Mathematics [22], makes recourse to series expansions, which are computed by means of efficient recurrences of different special functions (see [23], e.g.). However, to our knowledge, the literature lacks comparisons in the relative efficiency of these two radically different approaches when applied to the tesseral geopotential.Therefore, the question emerges whether averaging techniques based on expansions of elliptic motion are competitive or not, with algorithms devised for dealing in closed form with short-period effects, due to longitude-dependent terms of the geopotential. In the case of LEO, the advantage clearly favors the closed form solution, whose, based on the different perturbation orders of the Keplerian and the Coriolis term, tesseral effects can be averaged in closed form [24, 25]. Furthermore, in this particular case, the elimination of the parallax works wonders succeeding in the closed-form elimination of all tesseral terms, except for those related to m-dailies [6, 26]. This is a typical case in which the complexity of relegation makes it less efficient than normalization [5] or canonical simplification in the present context [27]. On the other hand, Keplerian and Coriolis terms can no longer be assumed to be of different orders when approaching resonance conditions, which is the common case of orbits in the medium Earth orbits (MEOs) region and above.Trying to give an answer to this question in the case of MEO, we carry out several comparisons in which we restrict the test model to a2×2 truncation of the geopotential. Waiting for more effective implementations of the relegation algorithm, which makes it useful also for the canonical simplification of tesseral resonance problems, we limit our comparisons to the case of nonresonant motion. In particular, we explore the performances of our algorithm in a set of selected cases borrowed from [9]. However, we also consider sample orbits with the same semimajor axis as Galileo orbits and different eccentricities, a case that better supports our claim that the new formulation of the tesseral relegation algorithm is valid for orbits with any semimajor axis. Besides, because both approaches, series expansion and tesseral relegation, after averaging may well end up being Brouwer’s long-term Hamiltonian (see [11, page 569]), the comparison criteria will be based on the size of the transformation equations from mean to osculating elements required by the different theories for obtaining analogous precision in the results, together with the relative time spent in their evaluation.The paper is organized as follows. For completeness, in Section2 we recall the basics of the relegation algorithm for the particular case of the geopotential, to which we contribute a slightly different, but definitely enlightening point of view, which justifies the application of the relegation technique to low-eccentricity orbits independently of the ratio of the orbital period to the Earth's rotation period [28]. In this section, a brief mention to the standard procedure for averaging longitude-dependent terms when using classical expansion of the elliptic motion is also included. Next, in Section 3 we briefly describe the whole procedure of computing the long-term Hamiltonian and the transformations equations from mean to osculating elements, which we implement for the simplest case of a second-order truncation of the geopotential. While application of the tesseral relegation to higher degree and order truncations does not involve any modifications to the basic algorithm, it produces a notable increase in the size of the series used—as is also the case when using the expansions of elliptic motion approach. Finally, comparisons between both procedures, relegation and series expansions, are presented in Section 4 for a variety of test cases, and the pros and cons of each method are discussed.


## 2. Averaging the Geopotential

The motion of artificial satellites close to the Earth departs slightly from Keplerian motion because of the noncentralities of the geopotential. Thisdisturbing potential is written in the usual form
(1)V=μr∑j≥2‍(αr)j[∑JjPj,0(sinφ)-∑k=1j‍(Cj,kcoskλ+Sj,ksinkλ)×Pj,k(sinφ)∑],
where μ is the Earth's gravitational parameter, α is the Earth's equatorial radius, Pj,k are Legendre polynomials, Jj=-Cj,0,Cj,k, and Sj,k are harmonic coefficients, and (φ,λ,r) are spherical coordinates—latitude, longitude, and distance, respectively—in an Earth-centered, Earth-fixed (ECEF) system of coordinates.Satellite motion under the influence of the disturbing potential in (1) is generally nonintegrable, even in the simplest case of the J2 problem [29, 30], with the remarkable exception of the equatorial case [31]. However, analytical or semianalytical integration of perturbed Keplerian motion is customarily approached by perturbation methods. Namely, short-period effects are analytically averaged to formulate the problem in mean elements; next, either the mean elements equations are numerically integrated with very long stepsizes, or a new averaging of long-period effects provides the secular evolution of the problem.


### 2.1. Lie Transforms: Deprit’s Triangle

In particular, the Lie transforms technique [32] allows for automatic computation of the perturbation theory by machine to higher orders, and these days is taken as a standard procedure. Briefly, in a Hamiltonian framework, the initial Hamiltonian in osculating variables is cast in the form
(2)ℋ=∑m≥0‍ϵmm!Hm,0,
where the small parameter ϵ, which establishes the order of each perturbation term, may be either physical or formal.The averaged Hamiltonian, in mean elements,(3)ℋ′=∑m≥0‍ϵmm!H0,m
is to be constructed by the perturbation method. The new Hamiltonian terms H0,m are one of the sides of “Deprit's triangle”: a triangular array of terms derived from the recurrence
(4)Hn,q=Hn+1,q-1+∑0≤m≤n(nm)‍{Hn-m,q-1;Wm+1},
where {-;-} notes the Poisson bracket operator and terms Wm+1 define the generating function
(5)𝒲=∑m≥0‍ϵmm!Wm+1
of the transformation from mean to osculating elements.The transformed Hamiltonian (3) is stepwise constructed by solving at each step the homological equation
(6)ℒ0(Wm)=H~m-H0,m,
where terms H~m are known from previous evaluations of Deprit's triangle, terms H0,m are chosen at our convenience, and terms Wm must be solved from a partial differential equation. The operator ℒ0(-)≡{-;H0,0} is known as the Lie derivative generated by the Keplerian flow H0,0.Finally, the transformation from mean to osculating elements is in turn computed from Deprit's triangle, by the simple expedient of replacing termsHi,j in (4) by ξi,j, where ξ0,0=ξ notes any of the osculating elements and ξm,0≡0 for m>0.


### 2.2. Removing Tesseral Effects from the Geopotential

Usual relations between spherical coordinates and orbital elements [12] permit to express the Hamiltonian of the geopotential in terms of orbital elements. Thus, the zeroth order of (2) corresponds to the Kepler problem in the ECEF frame
(7)H0,0=-μ2a-nEμaηcosi,
where (a,e,i,ω,Ω,M) are usual Keplerian elements, η=1-e2 is commonly known as the eccentricity function, and nE notes the Earth’s rotation rate, which is assumed to be constant. The first order term is
(8)H1,0=-μ2rα2r2J2[1-32sin2i+32sin2icos(2f+2g)],
where f is the true anomaly.The termH2=2(V-H1), with V given in (1), carries effects of the second order of J2, which are contributed by other zonal and tesseral harmonics. Namely, it takes the form of a trigonometric series whose coefficients are functions of e and i, and whose arguments are combinations of the type jf+kν+mω, with j, k, and m integers, where ν=Ω-ϑ is the longitude of the ascending node in the rotating frame and ϑ is Greenwich Sidereal Time. These series are particular instances of multivariate Fourier series of the form
(9)∑j1,j2,…jm≥0‍p1j1p2j2⋯pmjm×∑k1,k2,…kl≥0‍(Qj1,j2,…jmk1,k2,…klcosβ+Sj1,j2,…jmk1,k2,…klsinβ),
where β=k1β1+k2β2+⋯+klβl, p are polynomial variables with integer coefficients and Q and S are coefficients. Expressions of the form of (9) are commonly known as Poisson series of the type (m,l) [33].The Hamiltonian (2) must be formulated in canonical variables, whereas Keplerian elements are not canonical. The natural alternative is using Delaunay variables, the canonical counterpart of Keplerian elements, although other sets of canonical variables as nonsingular Poincaré elements are customarily used [20] (for a survey on orbital elements sets, see [34]). In what follows, it may be convenient using the orbital elements notation because of the natural insight that it provides to orbital motion problems. When this happens, orbital elements must be understood as functions of canonical variables instead of variables by themselves.Using Delaunay-type canonical variables(10)ℓ=M,g=ω,h=ν,L=μa,G=Lη,H=Gcosi,
where G is the modulus of the angular momentum vector and H its projection on the Earth's rotation axis, the zeroth-order Hamiltonian is written
(11)H0,0=-μ22L2-nEH.
Therefore, the Lie derivative is
(12)ℒ0=n∂∂ℓ-nE∂∂h,
where n=μ2/L3 is the orbital mean motion. At each step m of the perturbation algorithms we use (12) to compute the generating function term Wm from (6).


#### 2.2.1. Series Expansion

In spite of the simplicity of (12), the fact that the dependence of the geopotential on the mean anomaly does not happen explicitly but through the true anomaly, which is an implicit  function ℓ=ℓ(f,e) of the mean anomaly involving the solution of Kepler equation, complicates finding a closed form solution to the homological equation. Hence, expansions of the elliptic motion are customarily used to make explicit the dependence of the tesseral Hamiltonian on the mean anomaly (see, [11–14, 35], e.g.). Once the mean anomaly shows explicitly in the geopotential, the solution of the homological equation for the Lie derivative (12) is straightforward, not only allowing for the elimination of short-period terms related to tesseral harmonics in the case of nonresonant motion [36] but also for resonance problems [19].


#### 2.2.2. Relegation

Alternatively, one may realize that, except for synchronous orbits, eitherδ≡n/nE<1 or δ≡nE/n<1. Then, at each step m of the perturbation algorithm, one can feel satisfied with finding an approximate solution
(13)Wm=∑j≥0‍δjWm,j
which is computed as follows. We describe the subsynchronous orbits case δ≡nE/n<1, but the same philosophy can be applied to the supersynchronous case.Replacing (13) into the homological equation, (6), we get
(14)∑j≥0‍nδj(∂Wm,j∂ℓ-δ∂Wm,j∂h)=F(f,g,h,L,G,H)
with F≡H~m-H0,m, which we reorganize as
(15)∑j≥0‍δj∂Wm,j∂ℓ=∑j≥0‍δjFj,
where
(16)F0=1nF,Fj=∂Wm,j-1∂h(j>0).
Then, (15) is solved iteratively by equating equal powers of δ. At each iteration the corresponding differential equation is solved by quadrature. The corresponding solution may be reached in closed form with the help of the differential relation
(17)a2ηdℓ=r2df,
which is derived from the preservation of the angular momentum of Keplerian motion.A rule for ceasing iterations may be choosingj such that δj+1=𝒪(ϵ), which would imply neglecting from the generating function terms of higher order than the corresponding one of the perturbation theory. Thus, for example, in the usual case ϵ=𝒪(J2)=𝒪(10-3), for an orbit with a=10500 km, we get δ≈1/8, and hence j=2.Typically, termsWm,j grow in size with each iteration j making the relegation process unpractical for a relative low number of iterations [9]. This fact limits the practical application of the relegation algorithm to orbits relatively close to the Earth's surface, where nonimpact orbits must have small or moderate eccentricities. However, it has been recently shown that a slight modification of the relegation algorithm allows to extend its application to orbits with any semimajor axis, but constrained to the case of low-eccentricities [28]. For completeness, we summarize the procedure in what follows. Furthermore, we state the algorithm directly in Delaunay variables instead of the polar-nodal variables used in the original approach, thus relieving the relegation algorithm from the intricacies related to the specific algebra of the parallactic functions.Indeed, using (17), the Lie derivative is written as
(18)ℒ0=na2ηr2∂∂f-nE∂∂h,
which is rewritten as
(19)ℒ0=nη3(∂∂f-η3δ∂∂h+eℛ+e2𝒮),
where
(20)ℛ=2cosf∂∂f,𝒮=cos2f∂∂f.Therefore, the homological equation, (6), may be written as
(21)∂Wm∂f-η3δ∂Wm∂h=η3nF-2ecosf∂Wm∂f-e2cos2f∂Wm∂f,
showing that, independently of the value of the ratio δ (or, equivalently, for orbits with any semimajor axis), an analogous approach to the classical relegation algorithm can be taken also in the case of orbits with small eccentricities. Now, the integration of the homological equation is carried out in the true anomaly instead of the mean one, so there is no problem in solving the main part of the partial differential equation (21):
(22)∂Wm∂f-η3δ∂Wm∂h≈η3nF(f,g,h,L,G,H)
which is taken as the first step in the relegation iterations.Equation (13) is now written as
(23)Wm=∑j≥0‍ejWm,j.
Correspondingly, the analog of the homological equation (15) is
(24)∑j≥0‍ej(∂Wm,j∂f-η3δ∂Wm,j∂h)=∑j≥0‍ejFj,
where
(25)F0=η3nF,F1=-2cosf∂Wm,0∂f,
and, for j>1,
(26)Fj=-2cosf∂Wm,j-1∂f-cos2f∂Wm,j-2∂f.Note that this new formulation of the relegation algorithm shows explicitly that the effectiveness of the tesseral relegation is constrained to the case of low-eccentricity orbits. The limitation of the tesseral relegation to the case of small or moderate eccentricities seems to be an intrinsic limitation of the tesseral relegation for subsynchronous orbits and, although not enough emphasized, appears also in the classical formulation of the tesseral relegation because of the occurrence of Hamiltonian “type (A)” terms, in which the argument of the latitude appears explicitly [9]. An additional advantage of the present formulation of the tesseral relegation is that it does not need any special treatment of Hamiltonian terms related to m-dailies, namely, those terms of the Poisson series that only depend on the longitude of the ascending node in the rotating frame h.


## 2.1. Lie Transforms: Deprit’s Triangle

In particular, the Lie transforms technique [32] allows for automatic computation of the perturbation theory by machine to higher orders, and these days is taken as a standard procedure. Briefly, in a Hamiltonian framework, the initial Hamiltonian in osculating variables is cast in the form
(2)ℋ=∑m≥0‍ϵmm!Hm,0,
where the small parameter ϵ, which establishes the order of each perturbation term, may be either physical or formal.The averaged Hamiltonian, in mean elements,(3)ℋ′=∑m≥0‍ϵmm!H0,m
is to be constructed by the perturbation method. The new Hamiltonian terms H0,m are one of the sides of “Deprit's triangle”: a triangular array of terms derived from the recurrence
(4)Hn,q=Hn+1,q-1+∑0≤m≤n(nm)‍{Hn-m,q-1;Wm+1},
where {-;-} notes the Poisson bracket operator and terms Wm+1 define the generating function
(5)𝒲=∑m≥0‍ϵmm!Wm+1
of the transformation from mean to osculating elements.The transformed Hamiltonian (3) is stepwise constructed by solving at each step the homological equation
(6)ℒ0(Wm)=H~m-H0,m,
where terms H~m are known from previous evaluations of Deprit's triangle, terms H0,m are chosen at our convenience, and terms Wm must be solved from a partial differential equation. The operator ℒ0(-)≡{-;H0,0} is known as the Lie derivative generated by the Keplerian flow H0,0.Finally, the transformation from mean to osculating elements is in turn computed from Deprit's triangle, by the simple expedient of replacing termsHi,j in (4) by ξi,j, where ξ0,0=ξ notes any of the osculating elements and ξm,0≡0 for m>0.


## 2.2. Removing Tesseral Effects from the Geopotential

Usual relations between spherical coordinates and orbital elements [12] permit to express the Hamiltonian of the geopotential in terms of orbital elements. Thus, the zeroth order of (2) corresponds to the Kepler problem in the ECEF frame
(7)H0,0=-μ2a-nEμaηcosi,
where (a,e,i,ω,Ω,M) are usual Keplerian elements, η=1-e2 is commonly known as the eccentricity function, and nE notes the Earth’s rotation rate, which is assumed to be constant. The first order term is
(8)H1,0=-μ2rα2r2J2[1-32sin2i+32sin2icos(2f+2g)],
where f is the true anomaly.The termH2=2(V-H1), with V given in (1), carries effects of the second order of J2, which are contributed by other zonal and tesseral harmonics. Namely, it takes the form of a trigonometric series whose coefficients are functions of e and i, and whose arguments are combinations of the type jf+kν+mω, with j, k, and m integers, where ν=Ω-ϑ is the longitude of the ascending node in the rotating frame and ϑ is Greenwich Sidereal Time. These series are particular instances of multivariate Fourier series of the form
(9)∑j1,j2,…jm≥0‍p1j1p2j2⋯pmjm×∑k1,k2,…kl≥0‍(Qj1,j2,…jmk1,k2,…klcosβ+Sj1,j2,…jmk1,k2,…klsinβ),
where β=k1β1+k2β2+⋯+klβl, p are polynomial variables with integer coefficients and Q and S are coefficients. Expressions of the form of (9) are commonly known as Poisson series of the type (m,l) [33].The Hamiltonian (2) must be formulated in canonical variables, whereas Keplerian elements are not canonical. The natural alternative is using Delaunay variables, the canonical counterpart of Keplerian elements, although other sets of canonical variables as nonsingular Poincaré elements are customarily used [20] (for a survey on orbital elements sets, see [34]). In what follows, it may be convenient using the orbital elements notation because of the natural insight that it provides to orbital motion problems. When this happens, orbital elements must be understood as functions of canonical variables instead of variables by themselves.Using Delaunay-type canonical variables(10)ℓ=M,g=ω,h=ν,L=μa,G=Lη,H=Gcosi,
where G is the modulus of the angular momentum vector and H its projection on the Earth's rotation axis, the zeroth-order Hamiltonian is written
(11)H0,0=-μ22L2-nEH.
Therefore, the Lie derivative is
(12)ℒ0=n∂∂ℓ-nE∂∂h,
where n=μ2/L3 is the orbital mean motion. At each step m of the perturbation algorithms we use (12) to compute the generating function term Wm from (6).


### 2.2.1. Series Expansion

In spite of the simplicity of (12), the fact that the dependence of the geopotential on the mean anomaly does not happen explicitly but through the true anomaly, which is an implicit  function ℓ=ℓ(f,e) of the mean anomaly involving the solution of Kepler equation, complicates finding a closed form solution to the homological equation. Hence, expansions of the elliptic motion are customarily used to make explicit the dependence of the tesseral Hamiltonian on the mean anomaly (see, [11–14, 35], e.g.). Once the mean anomaly shows explicitly in the geopotential, the solution of the homological equation for the Lie derivative (12) is straightforward, not only allowing for the elimination of short-period terms related to tesseral harmonics in the case of nonresonant motion [36] but also for resonance problems [19].


### 2.2.2. Relegation

Alternatively, one may realize that, except for synchronous orbits, eitherδ≡n/nE<1 or δ≡nE/n<1. Then, at each step m of the perturbation algorithm, one can feel satisfied with finding an approximate solution
(13)Wm=∑j≥0‍δjWm,j
which is computed as follows. We describe the subsynchronous orbits case δ≡nE/n<1, but the same philosophy can be applied to the supersynchronous case.Replacing (13) into the homological equation, (6), we get
(14)∑j≥0‍nδj(∂Wm,j∂ℓ-δ∂Wm,j∂h)=F(f,g,h,L,G,H)
with F≡H~m-H0,m, which we reorganize as
(15)∑j≥0‍δj∂Wm,j∂ℓ=∑j≥0‍δjFj,
where
(16)F0=1nF,Fj=∂Wm,j-1∂h(j>0).
Then, (15) is solved iteratively by equating equal powers of δ. At each iteration the corresponding differential equation is solved by quadrature. The corresponding solution may be reached in closed form with the help of the differential relation
(17)a2ηdℓ=r2df,
which is derived from the preservation of the angular momentum of Keplerian motion.A rule for ceasing iterations may be choosingj such that δj+1=𝒪(ϵ), which would imply neglecting from the generating function terms of higher order than the corresponding one of the perturbation theory. Thus, for example, in the usual case ϵ=𝒪(J2)=𝒪(10-3), for an orbit with a=10500 km, we get δ≈1/8, and hence j=2.Typically, termsWm,j grow in size with each iteration j making the relegation process unpractical for a relative low number of iterations [9]. This fact limits the practical application of the relegation algorithm to orbits relatively close to the Earth's surface, where nonimpact orbits must have small or moderate eccentricities. However, it has been recently shown that a slight modification of the relegation algorithm allows to extend its application to orbits with any semimajor axis, but constrained to the case of low-eccentricities [28]. For completeness, we summarize the procedure in what follows. Furthermore, we state the algorithm directly in Delaunay variables instead of the polar-nodal variables used in the original approach, thus relieving the relegation algorithm from the intricacies related to the specific algebra of the parallactic functions.Indeed, using (17), the Lie derivative is written as
(18)ℒ0=na2ηr2∂∂f-nE∂∂h,
which is rewritten as
(19)ℒ0=nη3(∂∂f-η3δ∂∂h+eℛ+e2𝒮),
where
(20)ℛ=2cosf∂∂f,𝒮=cos2f∂∂f.Therefore, the homological equation, (6), may be written as
(21)∂Wm∂f-η3δ∂Wm∂h=η3nF-2ecosf∂Wm∂f-e2cos2f∂Wm∂f,
showing that, independently of the value of the ratio δ (or, equivalently, for orbits with any semimajor axis), an analogous approach to the classical relegation algorithm can be taken also in the case of orbits with small eccentricities. Now, the integration of the homological equation is carried out in the true anomaly instead of the mean one, so there is no problem in solving the main part of the partial differential equation (21):
(22)∂Wm∂f-η3δ∂Wm∂h≈η3nF(f,g,h,L,G,H)
which is taken as the first step in the relegation iterations.Equation (13) is now written as
(23)Wm=∑j≥0‍ejWm,j.
Correspondingly, the analog of the homological equation (15) is
(24)∑j≥0‍ej(∂Wm,j∂f-η3δ∂Wm,j∂h)=∑j≥0‍ejFj,
where
(25)F0=η3nF,F1=-2cosf∂Wm,0∂f,
and, for j>1,
(26)Fj=-2cosf∂Wm,j-1∂f-cos2f∂Wm,j-2∂f.Note that this new formulation of the relegation algorithm shows explicitly that the effectiveness of the tesseral relegation is constrained to the case of low-eccentricity orbits. The limitation of the tesseral relegation to the case of small or moderate eccentricities seems to be an intrinsic limitation of the tesseral relegation for subsynchronous orbits and, although not enough emphasized, appears also in the classical formulation of the tesseral relegation because of the occurrence of Hamiltonian “type (A)” terms, in which the argument of the latitude appears explicitly [9]. An additional advantage of the present formulation of the tesseral relegation is that it does not need any special treatment of Hamiltonian terms related to m-dailies, namely, those terms of the Poisson series that only depend on the longitude of the ascending node in the rotating frame h.


## 2.2.1. Series Expansion

In spite of the simplicity of (12), the fact that the dependence of the geopotential on the mean anomaly does not happen explicitly but through the true anomaly, which is an implicit  function ℓ=ℓ(f,e) of the mean anomaly involving the solution of Kepler equation, complicates finding a closed form solution to the homological equation. Hence, expansions of the elliptic motion are customarily used to make explicit the dependence of the tesseral Hamiltonian on the mean anomaly (see, [11–14, 35], e.g.). Once the mean anomaly shows explicitly in the geopotential, the solution of the homological equation for the Lie derivative (12) is straightforward, not only allowing for the elimination of short-period terms related to tesseral harmonics in the case of nonresonant motion [36] but also for resonance problems [19].


## 2.2.2. Relegation

Alternatively, one may realize that, except for synchronous orbits, eitherδ≡n/nE<1 or δ≡nE/n<1. Then, at each step m of the perturbation algorithm, one can feel satisfied with finding an approximate solution
(13)Wm=∑j≥0‍δjWm,j
which is computed as follows. We describe the subsynchronous orbits case δ≡nE/n<1, but the same philosophy can be applied to the supersynchronous case.Replacing (13) into the homological equation, (6), we get
(14)∑j≥0‍nδj(∂Wm,j∂ℓ-δ∂Wm,j∂h)=F(f,g,h,L,G,H)
with F≡H~m-H0,m, which we reorganize as
(15)∑j≥0‍δj∂Wm,j∂ℓ=∑j≥0‍δjFj,
where
(16)F0=1nF,Fj=∂Wm,j-1∂h(j>0).
Then, (15) is solved iteratively by equating equal powers of δ. At each iteration the corresponding differential equation is solved by quadrature. The corresponding solution may be reached in closed form with the help of the differential relation
(17)a2ηdℓ=r2df,
which is derived from the preservation of the angular momentum of Keplerian motion.A rule for ceasing iterations may be choosingj such that δj+1=𝒪(ϵ), which would imply neglecting from the generating function terms of higher order than the corresponding one of the perturbation theory. Thus, for example, in the usual case ϵ=𝒪(J2)=𝒪(10-3), for an orbit with a=10500 km, we get δ≈1/8, and hence j=2.Typically, termsWm,j grow in size with each iteration j making the relegation process unpractical for a relative low number of iterations [9]. This fact limits the practical application of the relegation algorithm to orbits relatively close to the Earth's surface, where nonimpact orbits must have small or moderate eccentricities. However, it has been recently shown that a slight modification of the relegation algorithm allows to extend its application to orbits with any semimajor axis, but constrained to the case of low-eccentricities [28]. For completeness, we summarize the procedure in what follows. Furthermore, we state the algorithm directly in Delaunay variables instead of the polar-nodal variables used in the original approach, thus relieving the relegation algorithm from the intricacies related to the specific algebra of the parallactic functions.Indeed, using (17), the Lie derivative is written as
(18)ℒ0=na2ηr2∂∂f-nE∂∂h,
which is rewritten as
(19)ℒ0=nη3(∂∂f-η3δ∂∂h+eℛ+e2𝒮),
where
(20)ℛ=2cosf∂∂f,𝒮=cos2f∂∂f.Therefore, the homological equation, (6), may be written as
(21)∂Wm∂f-η3δ∂Wm∂h=η3nF-2ecosf∂Wm∂f-e2cos2f∂Wm∂f,
showing that, independently of the value of the ratio δ (or, equivalently, for orbits with any semimajor axis), an analogous approach to the classical relegation algorithm can be taken also in the case of orbits with small eccentricities. Now, the integration of the homological equation is carried out in the true anomaly instead of the mean one, so there is no problem in solving the main part of the partial differential equation (21):
(22)∂Wm∂f-η3δ∂Wm∂h≈η3nF(f,g,h,L,G,H)
which is taken as the first step in the relegation iterations.Equation (13) is now written as
(23)Wm=∑j≥0‍ejWm,j.
Correspondingly, the analog of the homological equation (15) is
(24)∑j≥0‍ej(∂Wm,j∂f-η3δ∂Wm,j∂h)=∑j≥0‍ejFj,
where
(25)F0=η3nF,F1=-2cosf∂Wm,0∂f,
and, for j>1,
(26)Fj=-2cosf∂Wm,j-1∂f-cos2f∂Wm,j-2∂f.Note that this new formulation of the relegation algorithm shows explicitly that the effectiveness of the tesseral relegation is constrained to the case of low-eccentricity orbits. The limitation of the tesseral relegation to the case of small or moderate eccentricities seems to be an intrinsic limitation of the tesseral relegation for subsynchronous orbits and, although not enough emphasized, appears also in the classical formulation of the tesseral relegation because of the occurrence of Hamiltonian “type (A)” terms, in which the argument of the latitude appears explicitly [9]. An additional advantage of the present formulation of the tesseral relegation is that it does not need any special treatment of Hamiltonian terms related to m-dailies, namely, those terms of the Poisson series that only depend on the longitude of the ascending node in the rotating frame h.


## 3. Semianalytical Theory for a2×2 Geopotential

We illustrate the relegation procedure and investigate its performances by constructing a semianalytical theory for a truncation of the geopotential to degree and order 2. We limit to the second-order terms of the perturbation theory, including the transformation equations up to the second order of the small parameterϵ, which is taken of the same order as the second zonal harmonic J2. Then, from (1), the term H2,0=2(V-H1,0) in the initial Hamiltonian (2) is simply
(27)H2,0=-2μα2r3∑k=12‍(C2,kcoskλ+S2,ksinkλ)P2,k(sinφ),
where spherical coordinates are assumed to be expressed as functions of Delaunay variables.The construction of the theory consists of three different steps. First, tesseral effects are averaged using the relegation algorithm described above. This averaging reduces the Hamiltonian of the geopotential to the zonal part, which, because of the2×2 truncation, is just the main problem, which only considers the disturbing effect of the J2 harmonic. Next, short-period effects are partially removed via the elimination of the parallax [4]. Finally, remaining short periodics are averaged via a trivial Delaunay normalization [37]. It deserves to mention that, in view of we limit the perturbation theory to the second order of J2, the elimination of the parallax is not required for the closed form normalization [38], and, therefore, the parallax elimination and Delaunay transformation may be reduced to a single transformation. However, carrying out the elimination of the parallax still results convenient because the important reduction in the number of terms in the total transformation equations, compare [2].


### 3.1. Tesseral Relegation

Because the perturbation theory is limited to the second-order, there is no coupling betweenJ2 and the second order tesseral harmonics C2,1, S2,1, C2,2, and S2,2. Therefore, H0,1 is chosen the same as H1,0, which makes W1=0 at the first order of Deprit's triangle (4).Next, we use the Keplerian relationr=aη2/(1+ecosf), to express H2,0 as a Poisson series in f, g, and h, which is made of 84 trigonometric terms. Since the solution of each iteration of the homological equation (24) amounts to changing sines by cosines and vice versa, and dividing for adequate linear combinations of n and nE, the first term of the generating function W2,0 is also made of 84 terms. Computing F1 in (25) for the first iteration of the relegation involves multiplication of the whole series ∂W2,0/∂f by cosf, an operation that enlarges the resulting Poisson series by 16 new trigonometric terms, while the first approximation of the generating function W2=∑j=0,1W2,j is made of a total of 108 trigonometric terms. Similarly, the computation of F2 from (26) involves multiplication of Poisson series by cosf and cos2f, which again enlarges the size of the Poisson series, leading to a new approximation W2=∑j=0,2W2,j consisting of 132 trigonometric terms.Further iterations of the relegation are found to increase the generating function by the same amount, namely, by 24 new trigonometric terms at each new step. Thus, the fifth iteration of the relegation results in a total generating functionW2=∑j=0,5W2,j, that is, a Poisson series made of 84+5×24=204 trigonometric terms. We stop iterations at this order.It deserves mentioning that, in addition to the linear growing of the length of the Poisson series, the complexity of the coefficient of each trigonometric term notably increases with iterations of the relegation, as a result of the combination of the different series resulting from the procedure. Thus, if we take as a rough “complexity indicator” the size of the text file required for storing each generating function, we find out that the complexity of each new approximation to the generating function almost doubles with each iteration. Specifically, we obtain files of6.8, 13.9, 27.2, 48.1, 83.2, and 135.6 kB, for the 0,…,5 iteration, respectively (throughout this paper, we refer to binary prefixes, that is 1kB=210 bytes). Note that these sizes are the result of a standard simplification with a general purpose algebra system, without further investigating the structure of the series.


### 3.2. Elimination of the Parallax

After tesseral relegation, we get the well-known Hamiltonian of the main problem(28)ℋ=-μ2a-μ2rα2r2J2[1-32sin2i+32sin2icos(2f+2g)].
Note that since longitude-dependent terms were removed, the problem is now formulated in inertial space.The standard parallax elimination [4] transforms (28), up to the second order of J2, into
(29)ℋ=-μ2a-μ2aα2r2J2η2(1-32sin2i)-μ2aα4a2r2J22η6[54-218sin2i+2116sin4i+(38cos2i-1564sin4i)e2-(2116sin2i-4532sin4i)e2cos2g].Note that after carrying out the relegation transformation(ℓ,g,h,L,G,H)→ℛ(ℓ′,g′,h′,L′,G′,H′), the variables in which we express (28) are no longer osculating elements, and hence we should have used the prime notation. Analogously, the elimination of the parallax performs a new transformation from “prime” elements to “double prime” elements, in which should be expressed (29). In both cases, we relieved (28) and (29) from the prime notation because there is no risk of confusion. Recall also that the elimination of the parallax is carried out in polar-nodal variables (r,θ,ν,R,Θ,N) instead of Delaunay variables, where R is radial velocity, Θ=G, N=H=Θcosi, ν=h, and θ=f+g.


### 3.3. Long-Term Hamiltonian

A new transformation from double-prime to triple-prime variables normalizes (29) by removing the remaining short-period terms related to the radius. Then, disregarding the prime notation, we get the long-term Hamiltonian in mean (triple-prime) elements
(30)𝒦=-μ2a-μ2aα2a2J2η3(1-32s2)-μ2aα4a4J22η7[154-152s2+10532s4+38(2-3s2)2η-(34-34s2-1532s4)η2-(218s2-4516s4)e2cos2g]
which is easily verified to agree with Brouwer's original long-term Hamiltonian, compare [11, page 569].


### 3.4. Semianalytical Integration

The evolution equations are the Hamilton equations of the long-term Hamiltonian in (30), namely,
(31)d(ℓ,g,h)dt=∂𝒦∂(L,G,H),d(L,G,H)dt=-∂𝒦∂(ℓ,g,h).
Since ℓ and h have been removed by the perturbation theory, both L and H are integrals of the mean elements system defined by (30).Then, the semianalytical integration of the2×2 truncation of the geopotential performs the following steps:  (1)
choose osculating elements;(2)
compute prime elements related to relegation;(3)
compute second-prime elements related to the elimination of the parallax;(4)
compute mean elements related to normalization;(5)
numerically integrate (31), which allows for very long stepsizes.At each step of the numerical integration, osculating elements are computed from mean ones by recovering short-period effects. This requires undoing sequentially the transformations from triple to double primes, from double-primes to primes, and from primes to osculating.


## 3.1. Tesseral Relegation

Because the perturbation theory is limited to the second-order, there is no coupling betweenJ2 and the second order tesseral harmonics C2,1, S2,1, C2,2, and S2,2. Therefore, H0,1 is chosen the same as H1,0, which makes W1=0 at the first order of Deprit's triangle (4).Next, we use the Keplerian relationr=aη2/(1+ecosf), to express H2,0 as a Poisson series in f, g, and h, which is made of 84 trigonometric terms. Since the solution of each iteration of the homological equation (24) amounts to changing sines by cosines and vice versa, and dividing for adequate linear combinations of n and nE, the first term of the generating function W2,0 is also made of 84 terms. Computing F1 in (25) for the first iteration of the relegation involves multiplication of the whole series ∂W2,0/∂f by cosf, an operation that enlarges the resulting Poisson series by 16 new trigonometric terms, while the first approximation of the generating function W2=∑j=0,1W2,j is made of a total of 108 trigonometric terms. Similarly, the computation of F2 from (26) involves multiplication of Poisson series by cosf and cos2f, which again enlarges the size of the Poisson series, leading to a new approximation W2=∑j=0,2W2,j consisting of 132 trigonometric terms.Further iterations of the relegation are found to increase the generating function by the same amount, namely, by 24 new trigonometric terms at each new step. Thus, the fifth iteration of the relegation results in a total generating functionW2=∑j=0,5W2,j, that is, a Poisson series made of 84+5×24=204 trigonometric terms. We stop iterations at this order.It deserves mentioning that, in addition to the linear growing of the length of the Poisson series, the complexity of the coefficient of each trigonometric term notably increases with iterations of the relegation, as a result of the combination of the different series resulting from the procedure. Thus, if we take as a rough “complexity indicator” the size of the text file required for storing each generating function, we find out that the complexity of each new approximation to the generating function almost doubles with each iteration. Specifically, we obtain files of6.8, 13.9, 27.2, 48.1, 83.2, and 135.6 kB, for the 0,…,5 iteration, respectively (throughout this paper, we refer to binary prefixes, that is 1kB=210 bytes). Note that these sizes are the result of a standard simplification with a general purpose algebra system, without further investigating the structure of the series.


## 3.2. Elimination of the Parallax

After tesseral relegation, we get the well-known Hamiltonian of the main problem(28)ℋ=-μ2a-μ2rα2r2J2[1-32sin2i+32sin2icos(2f+2g)].
Note that since longitude-dependent terms were removed, the problem is now formulated in inertial space.The standard parallax elimination [4] transforms (28), up to the second order of J2, into
(29)ℋ=-μ2a-μ2aα2r2J2η2(1-32sin2i)-μ2aα4a2r2J22η6[54-218sin2i+2116sin4i+(38cos2i-1564sin4i)e2-(2116sin2i-4532sin4i)e2cos2g].Note that after carrying out the relegation transformation(ℓ,g,h,L,G,H)→ℛ(ℓ′,g′,h′,L′,G′,H′), the variables in which we express (28) are no longer osculating elements, and hence we should have used the prime notation. Analogously, the elimination of the parallax performs a new transformation from “prime” elements to “double prime” elements, in which should be expressed (29). In both cases, we relieved (28) and (29) from the prime notation because there is no risk of confusion. Recall also that the elimination of the parallax is carried out in polar-nodal variables (r,θ,ν,R,Θ,N) instead of Delaunay variables, where R is radial velocity, Θ=G, N=H=Θcosi, ν=h, and θ=f+g.


## 3.3. Long-Term Hamiltonian

A new transformation from double-prime to triple-prime variables normalizes (29) by removing the remaining short-period terms related to the radius. Then, disregarding the prime notation, we get the long-term Hamiltonian in mean (triple-prime) elements
(30)𝒦=-μ2a-μ2aα2a2J2η3(1-32s2)-μ2aα4a4J22η7[154-152s2+10532s4+38(2-3s2)2η-(34-34s2-1532s4)η2-(218s2-4516s4)e2cos2g]
which is easily verified to agree with Brouwer's original long-term Hamiltonian, compare [11, page 569].


## 3.4. Semianalytical Integration

The evolution equations are the Hamilton equations of the long-term Hamiltonian in (30), namely,
(31)d(ℓ,g,h)dt=∂𝒦∂(L,G,H),d(L,G,H)dt=-∂𝒦∂(ℓ,g,h).
Since ℓ and h have been removed by the perturbation theory, both L and H are integrals of the mean elements system defined by (30).Then, the semianalytical integration of the2×2 truncation of the geopotential performs the following steps:  (1)
choose osculating elements;(2)
compute prime elements related to relegation;(3)
compute second-prime elements related to the elimination of the parallax;(4)
compute mean elements related to normalization;(5)
numerically integrate (31), which allows for very long stepsizes.At each step of the numerical integration, osculating elements are computed from mean ones by recovering short-period effects. This requires undoing sequentially the transformations from triple to double primes, from double-primes to primes, and from primes to osculating.


## 4. Numerical Experiments

### 4.1. Preliminary Tests

In order to check the efficiency of our new variant of the relegation algorithm, we apply it to test some of the representative cases provided in [9]. In particular, in view of the center and bottom plots of Figure 2 of [9] deal with the case of low-altitude orbits, which are amenable to a more advantageous treatment of tesseral effects by using the standard parallax elimination [6, 8], we focus on the examples related to the top plot of Figure 2 of [9], namely, a set of orbits in the MEO region whose initial conditions are detailed in Table 1.Table 1
Initial osculating elements of selected test cases from [9].
Orbit
a (km)
e
i (deg)
Ω
ω (deg)
M
1
12159.596
0.01
5
0
270
0
2
13375.556
0.10
5
0
270
0
3
18520.000
0.35
5
0
270
0
4
18520.000
0.35
100
0
270
0Orbit 1 of Table1 is the orbit with the smallest eccentricity, and hence the better convergence of the relegation algorithm is expected. Indeed, while the propagation of the mean elements alone results in rough errors of about 200 meter for the semimajor axis, 𝒪(10-4) for the eccentricity, just a few arc seconds (as) for the inclination, tens of as for the Right Ascension of the ascending node (RAAN), and of just a few deg for both the argument of the perigee and the mean anomaly, Figure 1 shows that these errors are greatly improved when short-period zonal and tesseral effects are recovered, even without need of iterating the relegation algorithm. Thus, errors fall down to just one meter for the semi-major axis, to 𝒪(10-6) for the eccentricity, to one milliarcsecond (mas) in the case of inclination, below one as for the RAAN, and to tenths of deg both for the argument of the perigee and mean anomaly. The latter are shown to combine to just a few mas when using the nonsingular element F=ω+M, the mean distance to the node (not presented in Figure 1). Both the errors in F and Ω disclose a secular trend, which is of about 0.9 mas times orbital period for Ω (center plot of Figure 1), part of which should be attributed to the truncation of the perturbation theory to the second order. In the case of the mean distance to the node, the secular trend is of about 13 mas times orbital period, roughly indicating an along-track error of about 76 cm times orbital period.Errors in the orbital elements of orbit 1 of Table1 when tesseral short-period effects are recovered without iterations. Abscissas are orbital periods, corresponding to one month propagation.
(a)
(b)
(c)
(d)
(e)
(f)Figure2 shows that a single iteration of (24) clearly improves the errors in semimajor axis, eccentricity, and inclination, as well as the errors in the argument of the perigee and mean anomaly, which generally enhance by one order of magnitude. On the contrary, improvements in the short-period terms are outweighed by the dominant secular trend in the case of Ω, as displayed in the fourth row of Figure 2. Now, the secular trend of the errors in the mean distance to the node reduces to ~2.7 mas/cycle, roughly proportional to an along-track error of ~16 cm times orbital period.Errors in the orbital elements of orbit 1 of Table1 when tesseral short-period effects are recovered with a single iteration of (24). Abscissas are orbital periods, corresponding to one month propagation.
(a)
(b)
(c)
(d)
(e)
(f)Improvements in the propagation of orbit 1 obtained with further iterations of the relegation algorithm are very small, except for the eccentricity, whose errors reach𝒪(10-8) with a second iteration of the relegation. New iterations of (24) show no improvement for this example, in full agreement with the scaling by the eccentricity in (23). Indeed, performing a single iteration means neglecting terms factored by the square of the eccentricity, which is e2=0.0001<J2 for orbit 1. Therefore, a second iteration of the relegation would introduce third-order effects, which are not considered in the perturbation theory.Other orbits in Table1 have a notably larger eccentricity than orbit 1, a fact that makes more iterations of the relegation algorithm necessary. Thus, in view of orbit 2 has an eccentricity e=0.1, third order of J2 effects will appear only at the third iteration of (24), where e3=𝒪(J2). This resulted exactly to be the case where one month propagation of orbit 2 with the semianalytical theory produced errors of about 3 cm in the semi-major axis, of 𝒪(10-8) for the eccentricity, below the mas for the inclination, and below the as for both the argument of the perigee and the mean anomaly. The secular trend in Ω is now of half mas per orbit and about 2 mas per orbit in the case of the mean distance to the node, roughly equivalent to an along track error of ~14 cm per orbit.The moderate eccentricity of orbits 3 and 4 makes necessary the fifth iteration of the procedure, which takes account of terms of the order ofe5≈0.005. Results for one month propagation of orbit 4 with the semianalytical theory are presented in Figure 3, for the first iteration of the algorithm, and in Figure 4, for the fifth. Note the secular growing of the errors in the mean anomaly (bottom plot of Figure 4), roughly amounting to an along-track error of two meters per orbit. The rate of convergence of the relegation is poor for orbits 3 and 4 of Table 1, in agreement with the similar rate of reduction of successive powers of the eccentricity for the moderate value of e in this example.Errors in the orbital elements of orbit 4 of Table1 after a one iteration of (23). Abscissas are orbital periods, corresponding to one month propagation.
(a)
(b)
(c)
(d)
(e)
(f)Errors in the orbital elements of orbit 4 of Table1 after five iterations of the relegation algorithm. Abscissas are orbital periods, corresponding to one month propagation.
(a)
(b)
(c)
(d)
(e)
(f)


### 4.2. Comparisons with Expansions of Elliptic Motion

Once the performances of the new relegation algorithm have been checked, it rests to test its relative efficiency when compared with the usual approach based on expansions of elliptic motion. Having in mind that comparisons must limit to the case of small or moderate eccentricities, we decided to limit the expansion of the true anomaly as a function of the mean anomaly to the order ofe12. Hence, expansions of both the tesseral Hamiltonian and the generating function of the tesseral averaging are also limited to 𝒪(e12), but the transformation equations for ℓ and g, as well as the nonsingular element F=ℓ+g, are downgraded to the order of e11, compare [39]. With these limits for expansions in the eccentricity, we find out that the transformation equations of the tesseral averaging take up to 286 kB when stored in a text file. This size places the transformation equations of the expansions of the elliptic motion approach between those of the first (204 kB) and the second iteration of the relegation algorithm (434 kB). Note that zonal terms are subsequently averaged in closed form in both approaches, and hence the efficiency of each approach is just related to the way in which the tesseral averaging has been carried out.Thus, in view of the respective sizes of the files containing the transformation equations, we consider the tesseral averaging based on expansions of the elliptic motion to be more efficient than the tesseral relegation, in those cases in which it produces errors that are, at least, comparable to those of a second iteration of the relegation. On the contrary, in those cases in which the errors produced by the expansions of the elliptic motion approach are comparable or worse than those obtained with a single iteration of the relegation, we say that the tesseral relegation is more efficient than the expansions in the eccentricity procedure.It is not a surprise to check that tesseral averaging based on expansions of the elliptic motion is superior for orbits of moderate eccentricity. Thus, when applied to orbit 4 of Table1 (e=0.35), the eccentricity expansions approach reaches errors comparable to but slightly better than those of Figure 4, where the five iterations of the relegation used in the construction of Figure 4 require 2 135 kB of text storage, 7.5 times bigger than the eccentricity expansions file. Thus, in this case, we say that the relegation efficiency is about seven and a half times lesser than the efficiency of the eccentricity expansion approach. On the contrary, when applied to orbit 1 of Table 1 (e=0.01), the eccentricity expansions approach gets the same errors as those obtained with two iterations of the relegation, which are very similar to those of Figure 2. Therefore, we say in this case that the eccentricity expansion approach is of comparable efficiency to the relegation, with relative sizes of text storage 434kB/286kB≈1.5.Finally, we provide a comparison for a Galileo-type orbit, witha=29600 km, e=0.001, and i=56 deg. Note that the efficiency of the classical relegation approach of [9] will be very poor for orbits of these characteristics because of the small ratio nE/n≈7/12 that must be used as the relegation coefficient. Besides, Galileo orbits have the smallest eccentricity of all the cases tested. Results for this last comparison are presented in Figure 5, for a single iteration of the relegation approach, and Figure 6, for the eccentricity expansions case. We note that the errors are of a similar magnitude in both cases, with better results for the mean distance to the node when using relegation and better results for the eccentricity when using eccentricity expansions.Errors in the orbital elements of an almost-circular, Galileo-type orbit when tesseral short-period effects are recovered with one iteration of the relegation. Abscissas are orbital periods, corresponding to one month propagation.
(a)
(b)
(c)
(d)
(e)Same as Figure5 when using expansions truncated to 𝒪(e12).
(a)
(b)
(c)
(d)
(e)If we take a second iteration of the relegation, then errors in the eccentricity of both cases match, and in view of the respective magnitudes of the errors inF, the advantage goes clearly to the side of the relegation. Since errors in F are always better when using relegation, we say that the relegation is more efficient than the expansions in the eccentricity approach for Galileo orbits. Besides, in order to provide an efficiency estimator analogous to the previous cases, we checked that a “mixed relegation”—which includes the transformation equations of a single iteration for ℓ, g, h, and H and the transformation equations after two iterations for L and G, making a text file of 260 kB—provides the same errors as the eccentricity expansions approach, maintaining the clear advantage for F. Therefore, we might say that for Galileo orbits, the efficiency is at least 1.1 times better than the eccentricity expansions for all the elements, except for F, for which the efficiency of the relegation is at least 1.4 times better than the eccentricity expansions.However, in view of tesseral relegation for subsynchronous orbits seems to be useful only in the case of the lower eccentricity orbits, we further investigate this case. Thus, it is important to note that, because the term1/r3 in (27), the initial iteration in the solution of (24) gives a term W2,0 that involves eccentricity polynomials (not expansions) of degree 3. Since each iteration of the relegation involves multiplication by the eccentricity, we simplify the computation of the generating function by neglecting those powers of e that are higher than 3 in successive iterations of the procedure. Proceeding in this way, the transformation equations of the tesseral relegation can be shortened to accommodate in a file of 143.7 kB for a single iteration of the procedure and only 185.7 kB for two iterations. Therefore, we conclude that, for the tesseral model investigated, relegation may double the efficiency of tesseral averaging for the lower eccentricity orbits.


## 4.1. Preliminary Tests

In order to check the efficiency of our new variant of the relegation algorithm, we apply it to test some of the representative cases provided in [9]. In particular, in view of the center and bottom plots of Figure 2 of [9] deal with the case of low-altitude orbits, which are amenable to a more advantageous treatment of tesseral effects by using the standard parallax elimination [6, 8], we focus on the examples related to the top plot of Figure 2 of [9], namely, a set of orbits in the MEO region whose initial conditions are detailed in Table 1.Table 1
Initial osculating elements of selected test cases from [9].
Orbit
a (km)
e
i (deg)
Ω
ω (deg)
M
1
12159.596
0.01
5
0
270
0
2
13375.556
0.10
5
0
270
0
3
18520.000
0.35
5
0
270
0
4
18520.000
0.35
100
0
270
0Orbit 1 of Table1 is the orbit with the smallest eccentricity, and hence the better convergence of the relegation algorithm is expected. Indeed, while the propagation of the mean elements alone results in rough errors of about 200 meter for the semimajor axis, 𝒪(10-4) for the eccentricity, just a few arc seconds (as) for the inclination, tens of as for the Right Ascension of the ascending node (RAAN), and of just a few deg for both the argument of the perigee and the mean anomaly, Figure 1 shows that these errors are greatly improved when short-period zonal and tesseral effects are recovered, even without need of iterating the relegation algorithm. Thus, errors fall down to just one meter for the semi-major axis, to 𝒪(10-6) for the eccentricity, to one milliarcsecond (mas) in the case of inclination, below one as for the RAAN, and to tenths of deg both for the argument of the perigee and mean anomaly. The latter are shown to combine to just a few mas when using the nonsingular element F=ω+M, the mean distance to the node (not presented in Figure 1). Both the errors in F and Ω disclose a secular trend, which is of about 0.9 mas times orbital period for Ω (center plot of Figure 1), part of which should be attributed to the truncation of the perturbation theory to the second order. In the case of the mean distance to the node, the secular trend is of about 13 mas times orbital period, roughly indicating an along-track error of about 76 cm times orbital period.Errors in the orbital elements of orbit 1 of Table1 when tesseral short-period effects are recovered without iterations. Abscissas are orbital periods, corresponding to one month propagation.
(a)
(b)
(c)
(d)
(e)
(f)Figure2 shows that a single iteration of (24) clearly improves the errors in semimajor axis, eccentricity, and inclination, as well as the errors in the argument of the perigee and mean anomaly, which generally enhance by one order of magnitude. On the contrary, improvements in the short-period terms are outweighed by the dominant secular trend in the case of Ω, as displayed in the fourth row of Figure 2. Now, the secular trend of the errors in the mean distance to the node reduces to ~2.7 mas/cycle, roughly proportional to an along-track error of ~16 cm times orbital period.Errors in the orbital elements of orbit 1 of Table1 when tesseral short-period effects are recovered with a single iteration of (24). Abscissas are orbital periods, corresponding to one month propagation.
(a)
(b)
(c)
(d)
(e)
(f)Improvements in the propagation of orbit 1 obtained with further iterations of the relegation algorithm are very small, except for the eccentricity, whose errors reach𝒪(10-8) with a second iteration of the relegation. New iterations of (24) show no improvement for this example, in full agreement with the scaling by the eccentricity in (23). Indeed, performing a single iteration means neglecting terms factored by the square of the eccentricity, which is e2=0.0001<J2 for orbit 1. Therefore, a second iteration of the relegation would introduce third-order effects, which are not considered in the perturbation theory.Other orbits in Table1 have a notably larger eccentricity than orbit 1, a fact that makes more iterations of the relegation algorithm necessary. Thus, in view of orbit 2 has an eccentricity e=0.1, third order of J2 effects will appear only at the third iteration of (24), where e3=𝒪(J2). This resulted exactly to be the case where one month propagation of orbit 2 with the semianalytical theory produced errors of about 3 cm in the semi-major axis, of 𝒪(10-8) for the eccentricity, below the mas for the inclination, and below the as for both the argument of the perigee and the mean anomaly. The secular trend in Ω is now of half mas per orbit and about 2 mas per orbit in the case of the mean distance to the node, roughly equivalent to an along track error of ~14 cm per orbit.The moderate eccentricity of orbits 3 and 4 makes necessary the fifth iteration of the procedure, which takes account of terms of the order ofe5≈0.005. Results for one month propagation of orbit 4 with the semianalytical theory are presented in Figure 3, for the first iteration of the algorithm, and in Figure 4, for the fifth. Note the secular growing of the errors in the mean anomaly (bottom plot of Figure 4), roughly amounting to an along-track error of two meters per orbit. The rate of convergence of the relegation is poor for orbits 3 and 4 of Table 1, in agreement with the similar rate of reduction of successive powers of the eccentricity for the moderate value of e in this example.Errors in the orbital elements of orbit 4 of Table1 after a one iteration of (23). Abscissas are orbital periods, corresponding to one month propagation.
(a)
(b)
(c)
(d)
(e)
(f)Errors in the orbital elements of orbit 4 of Table1 after five iterations of the relegation algorithm. Abscissas are orbital periods, corresponding to one month propagation.
(a)
(b)
(c)
(d)
(e)
(f)


## 4.2. Comparisons with Expansions of Elliptic Motion

Once the performances of the new relegation algorithm have been checked, it rests to test its relative efficiency when compared with the usual approach based on expansions of elliptic motion. Having in mind that comparisons must limit to the case of small or moderate eccentricities, we decided to limit the expansion of the true anomaly as a function of the mean anomaly to the order ofe12. Hence, expansions of both the tesseral Hamiltonian and the generating function of the tesseral averaging are also limited to 𝒪(e12), but the transformation equations for ℓ and g, as well as the nonsingular element F=ℓ+g, are downgraded to the order of e11, compare [39]. With these limits for expansions in the eccentricity, we find out that the transformation equations of the tesseral averaging take up to 286 kB when stored in a text file. This size places the transformation equations of the expansions of the elliptic motion approach between those of the first (204 kB) and the second iteration of the relegation algorithm (434 kB). Note that zonal terms are subsequently averaged in closed form in both approaches, and hence the efficiency of each approach is just related to the way in which the tesseral averaging has been carried out.Thus, in view of the respective sizes of the files containing the transformation equations, we consider the tesseral averaging based on expansions of the elliptic motion to be more efficient than the tesseral relegation, in those cases in which it produces errors that are, at least, comparable to those of a second iteration of the relegation. On the contrary, in those cases in which the errors produced by the expansions of the elliptic motion approach are comparable or worse than those obtained with a single iteration of the relegation, we say that the tesseral relegation is more efficient than the expansions in the eccentricity procedure.It is not a surprise to check that tesseral averaging based on expansions of the elliptic motion is superior for orbits of moderate eccentricity. Thus, when applied to orbit 4 of Table1 (e=0.35), the eccentricity expansions approach reaches errors comparable to but slightly better than those of Figure 4, where the five iterations of the relegation used in the construction of Figure 4 require 2 135 kB of text storage, 7.5 times bigger than the eccentricity expansions file. Thus, in this case, we say that the relegation efficiency is about seven and a half times lesser than the efficiency of the eccentricity expansion approach. On the contrary, when applied to orbit 1 of Table 1 (e=0.01), the eccentricity expansions approach gets the same errors as those obtained with two iterations of the relegation, which are very similar to those of Figure 2. Therefore, we say in this case that the eccentricity expansion approach is of comparable efficiency to the relegation, with relative sizes of text storage 434kB/286kB≈1.5.Finally, we provide a comparison for a Galileo-type orbit, witha=29600 km, e=0.001, and i=56 deg. Note that the efficiency of the classical relegation approach of [9] will be very poor for orbits of these characteristics because of the small ratio nE/n≈7/12 that must be used as the relegation coefficient. Besides, Galileo orbits have the smallest eccentricity of all the cases tested. Results for this last comparison are presented in Figure 5, for a single iteration of the relegation approach, and Figure 6, for the eccentricity expansions case. We note that the errors are of a similar magnitude in both cases, with better results for the mean distance to the node when using relegation and better results for the eccentricity when using eccentricity expansions.Errors in the orbital elements of an almost-circular, Galileo-type orbit when tesseral short-period effects are recovered with one iteration of the relegation. Abscissas are orbital periods, corresponding to one month propagation.
(a)
(b)
(c)
(d)
(e)Same as Figure5 when using expansions truncated to 𝒪(e12).
(a)
(b)
(c)
(d)
(e)If we take a second iteration of the relegation, then errors in the eccentricity of both cases match, and in view of the respective magnitudes of the errors inF, the advantage goes clearly to the side of the relegation. Since errors in F are always better when using relegation, we say that the relegation is more efficient than the expansions in the eccentricity approach for Galileo orbits. Besides, in order to provide an efficiency estimator analogous to the previous cases, we checked that a “mixed relegation”—which includes the transformation equations of a single iteration for ℓ, g, h, and H and the transformation equations after two iterations for L and G, making a text file of 260 kB—provides the same errors as the eccentricity expansions approach, maintaining the clear advantage for F. Therefore, we might say that for Galileo orbits, the efficiency is at least 1.1 times better than the eccentricity expansions for all the elements, except for F, for which the efficiency of the relegation is at least 1.4 times better than the eccentricity expansions.However, in view of tesseral relegation for subsynchronous orbits seems to be useful only in the case of the lower eccentricity orbits, we further investigate this case. Thus, it is important to note that, because the term1/r3 in (27), the initial iteration in the solution of (24) gives a term W2,0 that involves eccentricity polynomials (not expansions) of degree 3. Since each iteration of the relegation involves multiplication by the eccentricity, we simplify the computation of the generating function by neglecting those powers of e that are higher than 3 in successive iterations of the procedure. Proceeding in this way, the transformation equations of the tesseral relegation can be shortened to accommodate in a file of 143.7 kB for a single iteration of the procedure and only 185.7 kB for two iterations. Therefore, we conclude that, for the tesseral model investigated, relegation may double the efficiency of tesseral averaging for the lower eccentricity orbits.


## 5. Conclusions

The relegation algorithm has been suitably formulated in Delaunay variables, contrary to the standard approach in polar-nodal variables. The use of Delaunay variables relieves the relegation algorithm from unnecessary subtleties related to the particular algebra of the parallactic functions and discloses the fundamental differences between the cases of super- and subsynchronous orbit relegation. Specifically, while the relegation factor for super-synchronous orbits is the ratio of one day to the satellite's orbital period, in the subsynchronous case the relegation factor unavoidably depends on the eccentricity, a fact that limits application of the subsynchronous relegation to orbits with small and moderate eccentricities. Besides, it deserves mentioning that the present formulation of the tesseral relegation does not need any special treatment of them-daily terms, albeit, of course, allowing for separate treatment if preferred.Computations in this paper show that the tesseral relegation algorithm for subsynchronous orbits may be competitive with classical expansions of elliptic motion, but only for orbits with small eccentricities, a case in which relegation doubles the efficiency of the series expansions approach. This fact constrains the practical application of the relegation but in no way invalidates it, because low-eccentricity orbits are used in a variety of aerospace engineering applications.



---

*Source: 570127-2013-04-17.xml*