#### 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 form *relegation*. 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.

#### 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 -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 a 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 Section 2 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. This *disturbing* potential is written in the usual form
where is the Earth's gravitational parameter, is the Earth's equatorial radius, are Legendre polynomials, , and are harmonic coefficients, and 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 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 where the small parameter , which establishes the order of each perturbation term, may be either physical or formal.

The averaged Hamiltonian, in mean elements, is to be constructed by the perturbation method. The new Hamiltonian terms are one of the sides of “Deprit's triangle”: a triangular array of terms derived from the recurrence where notes the Poisson bracket operator and terms define the generating function of the transformation from mean to osculating elements.

The transformed Hamiltonian (3) is stepwise constructed by solving at each step the *homological* equation
where terms are known from previous evaluations of Deprit's triangle, terms are chosen at our convenience, and terms must be solved from a partial differential equation. The operator is known as the Lie derivative generated by the Keplerian flow .

Finally, the transformation from mean to osculating elements is in turn computed from Deprit's triangle, by the simple expedient of replacing terms in (4) by , where notes any of the osculating elements and for .

##### 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 where are usual Keplerian elements, is commonly known as the eccentricity function, and notes the Earth’s rotation rate, which is assumed to be constant. The first order term is where is the true anomaly.

The term , with given in (1), carries effects of the second order of , which are contributed by other zonal and tesseral harmonics. Namely, it takes the form of a trigonometric series whose coefficients are functions of and , and whose arguments are combinations of the type , with , , and 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
where , are polynomial variables with integer coefficients and and are coefficients. Expressions of the form of (9) are commonly known as *Poisson series* of the type [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 where is the modulus of the angular momentum vector and its projection on the Earth's rotation axis, the zeroth-order Hamiltonian is written Therefore, the Lie derivative is where is the orbital mean motion. At each step of the perturbation algorithms we use (12) to compute the generating function term 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 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 or . Then, at each step of the perturbation algorithm, one can feel satisfied with finding an approximate solution which is computed as follows. We describe the subsynchronous orbits case , but the same philosophy can be applied to the supersynchronous case.

Replacing (13) into the homological equation, (6), we get with , which we reorganize as where 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 which is derived from the preservation of the angular momentum of Keplerian motion.

A rule for ceasing iterations may be choosing such that , 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 , for an orbit with km, we get , and hence .

Typically, terms grow in size with each iteration 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 which is rewritten as where

Therefore, the homological equation, (6), may be written as 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): which is taken as the first step in the relegation iterations.

Equation (13) is now written as Correspondingly, the analog of the homological equation (15) is where and, for ,

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 -dailies, namely, those terms of the Poisson series that only depend on the longitude of the ascending node in the rotating frame .

#### 3. Semianalytical Theory for a 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 . Then, from (1), the term in the initial Hamiltonian (2) is simply 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 the truncation, is just the main problem, which only considers the disturbing effect of the 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 , 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 between and the second order tesseral harmonics , , , and . Therefore, is chosen the same as , which makes at the first order of Deprit's triangle (4).

Next, we use the Keplerian relation , to express as a Poisson series in , , and , 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 and , the first term of the generating function is also made of 84 terms. Computing in (25) for the first iteration of the relegation involves multiplication of the whole series by , an operation that enlarges the resulting Poisson series by 16 new trigonometric terms, while the first approximation of the generating function is made of a total of 108 trigonometric terms. Similarly, the computation of from (26) involves multiplication of Poisson series by and , which again enlarges the size of the Poisson series, leading to a new approximation 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 function , that is, a Poisson series made of 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 of , , , , , and kB, for the iteration, respectively (throughout this paper, we refer to binary prefixes, that is 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 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 , into

Note that after carrying out the relegation transformation , 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 instead of Delaunay variables, where is radial velocity, , , , and .

##### 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 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, Since and have been removed by the perturbation theory, both and are integrals of the mean elements system defined by (30).

Then, the semianalytical integration of the 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.

Orbit 1 of Table 1 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, 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 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 , the mean distance to the node (not presented in Figure 1). Both the errors in and disclose a secular trend, which is of about 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 mas times orbital period, roughly indicating an along-track error of about 76 cm times orbital period.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Figure 2 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.

**(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 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 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 Table 1 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 , third order of effects will appear only at the third iteration of (24), where . 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 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 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 of . 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 in this example.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(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 of . Hence, expansions of both the tesseral Hamiltonian and the generating function of the tesseral averaging are also limited to , but the transformation equations for and , as well as the nonsingular element , are downgraded to the order of , 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 Table 1 (), 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, 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 (), 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 .

Finally, we provide a comparison for a Galileo-type orbit, with km, , and 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 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.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(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 in , the advantage goes clearly to the side of the relegation. Since errors in 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 , , , and and the transformation equations after two iterations for and , making a text file of 260 kB—provides the same errors as the eccentricity expansions approach, maintaining the clear advantage for . Therefore, we might say that for Galileo orbits, the efficiency is at least times better than the eccentricity expansions for all the elements, except for , for which the efficiency of the relegation is at least 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 term in (27), the initial iteration in the solution of (24) gives a term 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 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 the -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.

#### Acknowledgment

Part of this research has been supported by the Government of Spain (Projects AYA 2009-11896 and AYA 2010-18796 and grants from Gobierno de La Rioja fomenta 2010/16).