¹¹institutetext: M. Lara ²²institutetext: Scientific Computing Group–GRUCACI, University of La Rioja, Madre de Dios 53, 26006 Logroño, La Rioja, Spain
Tel.: +34-941-299440
Fax: +34-941-299460
²²email: [email protected] ³³institutetext: M. Lara ⁴⁴institutetext: Space Dynamics Group, Polytechnic University of Madrid–UPM, Plaza Cardenal Cisneros 3, 28040 Madrid, Spain

Solution to the main problem of the artificial satellite by reverse normalization

Martin Lara

( draft of April 16, 2020)

Abstract

The non-linearities of the dynamics of Earth artificial satellites are encapsulated by two formal integrals that are customarily computed by perturbation methods. Standard procedures begin with a Hamiltonian simplification that removes non-essential short-period terms from the Geopotential, and follow with the removal of both short- and long-period terms by means of two different canonical transformations that can be carried out in either order. We depart from the tradition and proceed by standard normalization to show that the Hamiltonian simplification part is dispensable. Decoupling first the motion of the orbital plane from the in-plane motion reveals as a feasible strategy to reach higher orders of the perturbation solution, which, besides, permits an efficient evaluation of the long series that comprise the analytical solution.

Keywords:

main problem Hamiltonian mechanics normalization canonical perturbation theory floating-point arithmetic

^†^†journal: Nonlinear Dynamics

1 Introduction

The dynamics of close Earth satellites under gravitational effects are mostly driven by perturbations of the Keplerian motion induced by the Earth oblateness. For this reason, the approximation obtained when truncating the Legendre polynomials expansion of the Geopotential to the only contribution of the zonal harmonic of the second degree, whose coefficient is usually denoted $J_{2}$ , is traditionally known as the main problem of artificial satellite theory. The $J_{2}$ -truncation of the gravitational potential is known to give rise to nonintegrable dynamics Danby1968 ; IrigoyenSimo1993 ; Broucke1994 ; CellettiNegrini1995 that comprise short- and long-period effects, as well as secular terms Kozai1959 ; Kaula1961 . However, due to the smallness of the $J_{2}$ coefficient of the Earth, the full system can be replaced by a separable approximation, which is customarily obtained by removing the periodic effects by means of perturbation methods Nayfeh2004 .

When written in the action-angle variables of the Kepler problem, also called Delaunay variables, the main problem Hamiltonian immediately shows that the right ascension of the ascending node is a cyclic variable. In consequence, its conjugate momentum, the projection of the angular momentum vector along the Earth’s rotation axis, is an integral of the main problem, which, therefore, is just of two degrees of freedom. Then, following Brouwer Brouwer1959 , the main problem Hamiltonian is normalized in two steps. First, the short period effects are removed by means of a canonical transformation that, after truncation to some order of the perturbation approach, turns the conjugate momentum to the mean anomaly (the Delaunay action) into a formal integral. Next, a new canonical transformation converts the total angular momentum into another formal integral. The main problem Hamiltonian is in this way completely reduced to a function of only the momenta in the new variables, whose Hamilton equations are trivially integrable.

In spite of the normalization procedure is standard from the point of view of perturbation theory, it happens that not all the action-angle variables of the Kepler problem appear explicitly in the Geopotential, as is usually advisable in the construction of perturbation solutions Kinoshita1972 ; LaraFerrer2013 ; LaraFerrer2015 ; Lara2018 . In particular, the mean anomaly remains implicit in the gravitational potential through its dependence on the polar angle. Unavoidably, Kepler’s equation must be solved to show the mean anomaly explicit, a fact that entails expanding the elliptic motion in powers of the eccentricity Hansen1855 ; BrouwerClemence1961 . This standard procedure is quite successful when dealing with orbits of low eccentricities, like is typical in a variety of astronomy problems Delaunay1860 ; Tisserand1889 . However, it involves handling very long Fourier series in the case of orbits with moderate eccentricities DepritRom1970 ; Kinoshita1977SAOSR , and hence the application of this method to different problems of interest in astrodynamics is de facto prevented.

On the contrary, the integrals appearing in the solution of the artificial satellite problem can be solved in closed form of the eccentricity Brouwer1959 ; Kozai1959 . It only requires the help of the standard relation between the differentials of the true and mean anomalies that is derived from the preservation of the total angular momentum of the Keplerian motion. Regrettably, the closed form approach soon finds difficulties in achieving higher orders of the short-period elimination, which stem from the impossibility of obtaining the antiderivative of the equation of the center (the difference between the true and mean anomalies) in closed form of the eccentricity in the realm of trigonometric functions Jefferys1971 ; OsacarPalacian1994 . Nonetheless, the difficulties are overcome by the artifact of grouping the equation of the center with other functions appearing in the procedure previously to approaching their integration Kozai1962 ; Aksnes1971 ; Deprit1982 . Alternatively, the application of a preliminary Hamiltonian simplification, the elimination of the parallax Deprit1981 ; LaraSanJuanLopezOchoa2013b , eases the consequent removal of short-period effects to some extent CoffeyDeprit1982 ; Healy2000 .

From the point of view of the perturbation approach, removing the short-period effects in the first place seems the more natural in view of the degeneracy of the Kepler problem. Indeed, the Kepler Hamiltonian in action-angle variables, on which the perturbation approach hinges on, only depends on the Delaunay action GoldsteinPooleSafko2001 ; LaraTossa2016 . However, the order in which the formal integrals are sequentially introduced in the perturbation solution is not relevant in a total normalization procedure, whose result is unique Arnold1989 . In fact, it happens that relegating the transformation of the Delaunay action into a formal integral to the last step of the perturbation approach provides clear simplifications in dealing with the equation of the center AlfriendCoffey1984 . In this way the task of extending the solution of the main problem to higher orders is notably simplified.

Converting the total angular momentum into a formal integral of the main problem requires making cyclic the argument of the perigee, up to some truncation order of the perturbation approach, in the transformed Hamiltonian. However, as it appeared in the literature, the transformation called by their authors the elimination of the perigee AlfriendCoffey1984 is not the typical normalization procedure, although it operates analogous results. Indeed, on the one hand, the elimination of the perigee is only applied to a Hamiltonian obtained after the elimination of the parallax, to which simplification it could seem to be unavoidably attached. On the other hand, when extended to higher orders, it removes more terms than those needed in a partial normalization. Last, the fact that the technique was originally devised in the canonical set of polar variables, to which the argument of the perigee does not pertain, neither helps in grasping the essence of the transformation. Reimplementation of the procedure in the usual set of action-angle variables makes the process of converting the argument of the perigee into a cyclic variable much more evident LaraSanJuanLopezOchoa2013c , but it still bears the same differences with respect to a classical normalization procedure.

We disregard the claimed benefits of Hamiltonian simplification procedures and compute the solution to the main problem of the artificial satellite theory by standard normalization. It is called reverse normalization because we exchange the order in which the formal integrals are traditionally introduced when solving the artificial satellite problem. More precisely, the total angular momentum is transformed into a formal integral in the first place, in this way decoupling the motion of the orbital plane from the satellite’s motion on that plane.¹¹1The advantages of decoupling the motion of the instantaneous orbital plane from the in-plane motion are well known, and are commonly pursued in the search for efficient numerical integration methods, q.v. Lara2017if and references therein. Then, a second canonical transformation converts the mean anomaly into a cyclic variable, in this way achieving the total reduction of the main problem Hamiltonian.

The procedure for making the argument of the perigee cyclic in the first place follows an analogous strategy to the one devised in the classical elimination of the perigee transformation AlfriendCoffey1984 . However, in our approach it is applied directly to the original main problem Hamiltonian, and differs from the original technique, as well as from an analogous procedure carried out in SanJuanetal2013 , in which the parallactic terms (inverse powers of the radius with exponents higher than 2) are not removed from the new, partially normalized Hamiltonian. In spite of that, we did not find trouble in dealing with the equation of the center in closed form in the subsequent Delaunay normalization Deprit1982 , a convenience that might had been anticipated from the discussions in Lara2019CMDA .

The Hamiltonian reduction has been approached in Delaunay variables. Unfortunately, these variables share the deficiencies of their partner Keplerian elements, which are singular for circular orbits and for equatorial orbits. Because of that, the secular terms are reformulated in a set of non-singular variables that replaces the mean anomaly, the argument of the perigee, and the total angular momentum, by the mean argument of the latitude and the projections of the eccentricity vector in the nodal frame, which are sometimes denoted semi-equinoctial variables Konopliv1990 . For the periodic corrections, we find convenience in using polar variables, which in the particular case of the main problem are free from small divisors of any kind except for those related to the critical inclination resonance, and are known to produce compact expressions that yield faster evaluation Izsak1963 ; LaraVilhenaSanchezPrado2015 .

We extended the complete normalization to the order six of the perturbation approach, which, to our knowledge, is the maximum order that has been reported in the literature (yet limited to partial normalization cases) Healy2000 ; SanJuanetal2013 . The aim of computing such a high order is not to enter a competition. On the contrary, we did it simply because the particular value of the Earth’s $J_{2}$ coefficient is $\sim 10^{-3}$ , and hence the computed solution should be exact to the numerical precision of standard floating-point arithmetic Kahan1997 already at the fifth order. We checked that it is exactly the case, and extending the computations to the sixth order only served us to verify that we don’t find observable improvements in our tests. In order to compute this high-order solution we resorted to the practicalities of standard commercial software, in which Deprit’s perturbation algorithm based on the Lie transforms method Deprit1969 is easily implemented. On account of the current computational power, it should be quite feasible to extend the perturbation solution, if desired, to even higher orders, although we didn’t try that.

On the other hand, making cyclic the argument of the perigee with the standard normalization in action-angle variables seems to be a particularly efficient procedure from the computational point of view, resulting in a generating function whose size is astonishingly smaller than alternative proposals in the literature. This fact, combined with the reduction in the number of transformations required by the perturbation solution to just two, as opposite to the three transformations needed when the elimination of the parallax is carried out in the first place, might make this latter simplification dispensable, as well as the discussions in DepritMiller1989 questionable, at least in what respects to the computation of explicit analytical solutions as opposite to partial normalizations to be integrated semianalytically.

Tests carried out on typical Earth orbits of different types confirm that the computed solution bears exactly the expected characteristics of a perturbation solution of perturbed Keplerian motion. In particular, we verified the degree of convergence of successive approximations. We confirmed too that, as expected, the quality of the analytical solution degrades in the vicinity of the critical inclination resonance —a case that, of course, would require a specific treatment CoffeyDepritMiller1986 ; Jupp1988 ; Lara2015IR . In addition, we checked that the computation of the initialization constants of the analytical solution from corresponding initial conditions gets a clear benefit of using a higher order of the periodic corrections than the order needed for ephemeris evaluation, yet this additional accuracy is not needed in all the variables and can be limited to the computation of the formal integral given by the Delaunay action BreakwellVagners1970 ; DepritRom1970 ; HautesserresLaraCeMDA2017 .

2 Construction of the analytical solution

The solution of Laplace’s equation in spherical coordinates provides the gravitational potential in the form of an expansion in harmonic functions. In the case of the Earth, the mass distribution is almost symmetric with respect to the rotation axis. For this reason, the Geopotential is simplified in different applications to just the contribution of the zonal harmonics (see Danby1992 , for instance)

\mathcal{V}=\frac{\mu}{r}\sum_{m\geq{0}}\frac{R_{\oplus}^{m}}{r^{m}}J_{m}P_{m}(\sin\varphi),

(1)

where $r$ is distance from the origin, $\varphi$ is latitude, $\mu$ is the Earth’s gravitational parameter, $R_{\oplus}$ is the Earth’s mean equatorial radius, $J_{m}$ are zonal harmonic coefficients, and $P_{m}$ denote Legendre polynomials of degrees $m$ .

The flow stemming from the potential (1) admits Hamiltonian formulation. Besides, because $J_{2}\sim 10^{-3}$ whereas $J_{m}$ ( $m>2$ ) are $\mathcal{O}(J_{2}^{2})$ or smaller, the Hamiltonian

\mathcal{H}=-\frac{\mu}{2a}-\frac{\mu}{r}\frac{R_{\oplus}^{2}}{r^{2}}J_{2}\frac{1}{2}\left(1-3\sin^{2}\varphi\right),

(2)

is representative of the main characteristics of the dynamics of close Earth satellites, and is commonly known as the main problem of artificial satellite theory. The symbol $a$ in Eq. (2) stands for the orbit semi-major axis, and, from standard trigonometric relations, $\sin\varphi=\sin{I}\sin\theta$ , where $I$ and $\theta$ are used to denote the inclination and the argument of the latitude, respectively.

When using the classic set of Keplerian elements, the argument of the latitude is $\theta=f+\omega$ , where $f$ and $\omega$ denote the true anomaly and the argument of the periapsis, respectively. The radius is computed from the conic equation

r=\frac{p}{1+e\cos{f}},

(3)

where $p=a(1-e^{2})$ is the parameter of the conic and $e$ is the eccentricity of the orbit. Note, however, that, because we are using Hamiltonian formulation, the symbols appearing in Eq. (2) need to be expressed as functions of some set of canonical variables. While the Keplerian variables are not canonical, they are conveniently expressed in the set of Delaunay canonical variables $(\ell,g,h,L,G,H)$ , in which $\ell$ is the mean anomaly, $g=\omega$ , $h$ is the argument of the ascending node, $L=\sqrt{\mu{a}}$ is known as the Delaunay action, $G=L\sqrt{1-e^{2}}$ is the total angular momentum, and $H=G\cos{I}$ is the projection of the angular momentum along the Earth’s rotation axis. That the latter is an integral of the main problem becomes evident after checking that $h$ is an ignorable variable in Hamiltonian (2), which, in consequence, is of just two degrees of freedom.

Besides, the main problem Hamiltonian is conservative, and, therefore, Eq. (2) remains constant (the total energy) for a given set of initial conditions. On the other hand, the existence of a third integral has not been proved. Then, exact solutions of the main problem are not known. Alternatively, it is approached with the usual tools of non-linear dynamics, such as Poincaré surfaces of section or the computation of families of periodic orbits Danby1968 ; Broucke1994 ; SADSaM1999 ; Lara1999 . Still, in some regions of phase space, and for particular energy values, the third integral can be constructed formally with the help of perturbation methods, obtaining in this way useful analytical approximations to the main problem dynamics in these particular regions.

2.1 Perturbation approach

We rewrite the main problem in the form of a perturbation Hamiltonian

\mathcal{H}=\mathcal{H}(\ell,g,-,L,G,H)\equiv\sum_{m\geq 0}\frac{\varepsilon^{m}}{m!}K_{m,0},

(4)

in which $\varepsilon$ is a formal small parameter ( $\varepsilon=1$ ) used to manifest the strength of each summand of Eq. (4), and

$\displaystyle K_{0,0}$	$\displaystyle=$	$\displaystyle-\frac{\mu}{2a},$	(5)
$\displaystyle K_{1,0}$	$\displaystyle=$	$\displaystyle-\frac{\mu}{r}\frac{R_{\oplus}^{2}}{r^{2}}\frac{1}{4}J_{2}\left[2-3s^{2}+3s^{2}\cos(2f+2\omega)\right],$	(6)
$\displaystyle K_{m,0}$	$\displaystyle=$	$\displaystyle 0,\quad m\geq 2.$	(7)

The symbol $s$ in Eq. (6) stands for the usual abbreviation of the sine of the inclination.

The aim is to find a transformation of variables

(\ell,g,h,L,G,H;\varepsilon)\mapsto(\ell^{\prime},g^{\prime},h^{\prime},L^{\prime},G^{\prime},H^{\prime}),

given also by an expansion in powers of the small parameter, such that, up to some truncation order $m=k$ , the Hamiltonian in the new variables is transformed into a separable Hamiltonian. Namely,

\mathcal{H}^{\prime}=\sum_{m=0}^{k}\frac{\varepsilon^{m}}{m!}K_{0,m}(-,-,-,L^{\prime},G^{\prime},H)+\mathcal{O}(\varepsilon^{k+1}).

(8)

The desired transformation is derived from a scalar generating function

\mathcal{W}=\sum_{m\geq 0}\frac{\varepsilon^{m}}{m!}W_{m+1}.

(9)

which, in our case, is computed using Deprit’s perturbation algorithm based on Lie transforms Deprit1969 . The procedure is summarized in finding a particular solution of the homological equation

\mathcal{L}_{0}(W_{m})\equiv\{W_{m};K_{0,0}\}=\widetilde{K}_{0,m}-K_{0,m},

(10)

in which $\{W_{m};K_{0,0}\}$ stands for the computation of the Poisson bracket of $W_{m}$ and the zeroth order term of the Hamiltonian. Terms $\widetilde{K}_{0,m}$ are known from the original Hamiltonian as well as from previous computations to the step $m$ , which are carried out using Deprit’s fundamental recursion

{F}_{n,q}={F}_{n+1,q-1}+\sum_{m=0}^{n}{n\choose{m}}\{{F}_{n-m,q-1};W_{m+1}\}.

(11)

Equation (11) can be used to formulate any function $\mathcal{F}$ of the canonical set of original variables like a function of the new variables, the Hamiltonian (4) being just one among the different possibilities. Finally, the term $K_{0,m}$ is chosen at our will, with the only condition of making the homological equation solvable in $W_{m}$ .

When using Delaunay variables the homological equation of the main problem Hamiltonian is solved by indefinite integration. Indeed, Eq. (5) turns into $K_{0,0}=-\frac{1}{2}\mu^{2}/L^{2}$ , and hence the Lie derivative operator $\mathcal{L}_{0}$ takes the extremely simple form

\mathcal{L}_{0}(W_{m})=n\frac{\partial{W}_{m}}{\partial\ell},

in which $n=\sqrt{\mu/a^{3}}=\mu^{2}/L^{3}$ is the mean motion of the unperturbed problem. Then, from Eq. (10),

W_{m}=\frac{1}{n}\int(\widetilde{K}_{0,m}-K_{0,m})\,\mathrm{d}\ell+C_{m},

(12)

where the functions $C_{m}\equiv{C}_{m}(-,g,L,G;H)$ play the role of integration “constants” that verify $\mathrm{d}C_{m}/\mathrm{d}\ell=0$ . Equation (12) is solved in closed form of the eccentricity with the help of the differential relation between the true and mean anomalies provided by the preservation of the angular momentum of the Kepler problem. That is, $G=r^{2}\mathrm{d}f/\mathrm{d}t$ , from which

a^{2}\eta\,\mathrm{d}\ell=r^{2}\mathrm{d}f,

(13)

where $\eta=\sqrt{1-e^{2}}=G/L$ is usually known as the eccentricity function. Hence,

W_{m}=\frac{1}{n}\int(\widetilde{K}_{0,m}-K_{0,m})\frac{r^{2}}{a^{2}\eta}\,\mathrm{d}f+C_{m}.

(14)

The transformation that makes the main problem Hamiltonian separable, up to the truncation order, is split into two different canonical transformations. With the first one

(\ell,g,h,L,G,H;\varepsilon)\mapsto(\ell^{\prime},g^{\prime},h^{\prime},L^{\prime},G^{\prime},H^{\prime}),

we make the argument of the perigee cyclic, thus converting the total angular momentum into a formal integral. That is, the motion of the satellite in the orbital plane, whose inclination remains constant in the new variables, is decoupled from the motion of the orbital plane. Therefore, the reduced system representing the motion on the orbital plane could be integrated separately. Rather, we carry out a second canonical transformation

(\ell^{\prime},g^{\prime},h^{\prime},L^{\prime},G^{\prime},H^{\prime};\varepsilon)\mapsto(\ell^{\prime\prime},g^{\prime\prime},h^{\prime\prime},L^{\prime\prime},G^{\prime\prime},H^{\prime\prime}),

that makes ignorable the mean anomaly in the transformed Hamiltonian in double-prime variables. In this way, the Delaunay action is converted into a formal integral too, and the complete Hamiltonian reduction is achieved up to the truncation order of the perturbation solution, whose Hamilton equations are thus trivially integrable.

The computation of the secular terms from the normalized Hamiltonian is only part of the perturbation solution. To be complete, it requires the correct initialization of the constants of the perturbation solution from corresponding initial conditions using the inverse transformation (from original to double prime variables), on the one hand, and the recovery of periodic effects using the direct transformation (from double prime to original variables) to obtain the ephemeris corresponding to the secular terms propagation. Both the direct and inverse transformations are obtained from standard application of the Lie transforms method. Because there are no specific artifices related to the computation of the current solution in what respects to that part, we do not provide details on their computation and an interested reader is referred to the literature Deprit1969 .

2.2 Normalization of the total angular momentum

At the first order, we check from Eq. (11) that the known terms involved in the homological equation consist only of $\widetilde{K}_{0,1}=K_{1,0}$ . Therefore, we chose the new Hamiltonian term

K_{0,1}=-\frac{\mu}{r}\frac{R_{\oplus}^{2}}{r^{2}}\frac{1}{4}J_{2}\left(2-3s^{2}\right),

(15)

which is the part of Eq. (6) that is free from the argument of the perigee, as desired. Then, Eq. (14) turns into

W_{1}=-\frac{3}{4}GJ_{2}\frac{R_{\oplus}^{2}}{p^{2}}s^{2}\int\frac{p}{r}\cos(2f+2\omega)\,\mathrm{d}f+C_{1},

which is solved by standard integration after replacing the radius with the conic equation (3). We obtain

W_{1}=-\epsilon G\frac{1}{2}s^{2}\sum_{i=1}^{3}3^{\left\lfloor 2-\frac{1}{2}i\right\rfloor}e^{\left|i-2\right|}\sin(if+2\omega)+C_{1},

(16)

where $C_{1}$ is hold arbitrary by now, the symbol $\lfloor\;\rfloor$ stands for integer part, and we abbreviated

\epsilon=\frac{1}{4}J_{2}\frac{R_{\oplus}^{2}}{p^{2}},

(17)

in which $p=G^{2}/\mu$ , and, therefore $\epsilon$ is a function of the total angular momentum. Recall that the symbols $p$ , $s$ , $e$ , $\omega$ , and $f$ in Eqs. (15), (16), and (17) are functions of the Delaunay variables.

At the second order we compute $K_{0,2}$ from Eq. (11), to obtain

K_{0,2}=\{K_{0,1},W_{1}\}+K_{1,1},

(18)

in which $K_{1,1}$ is computed using again Eq. (11). Namely,

K_{1,1}=\{K_{0,0},W_{2}\}+\{K_{1,0},W_{1}\}+K_{2,0}.

(19)

On account of $K_{2,0}=0$ from Eq. (7), the known terms hitherto of the homological equation (10) are

\widetilde{K}_{0,2}=\{K_{0,1},W_{1}\}+\{K_{1,0},W_{1}\},

whose computation only involves partial differentiation. Before approaching the solution of Eq. (14), we first check that the term $\widetilde{K}_{0,2}$ is made of trigonometric coefficients $T_{i}$ whose arguments always involve the true anomaly as argument, except for the terms

$\displaystyle T_{0}$	$\displaystyle=$	$\displaystyle\epsilon^{2}\frac{3}{4}\frac{\mu}{r}\frac{p}{r}s^{2}\left[e^{2}(23s^{2}-16\right)+8\left(4s^{2}-3)\right],$
$\displaystyle T_{1}$	$\displaystyle=$	$\displaystyle-\frac{3}{2}\epsilon^{2}G^{2}\frac{1}{r^{2}}(15s^{2}-14)s^{2}e^{2}\cos 2\omega,$
$\displaystyle T_{2}$	$\displaystyle=$	$\displaystyle 6\epsilon{G}\frac{1}{r^{2}}(5s^{2}-4)\frac{\partial{C}_{1}}{\partial{g}}.$

Integration of these 3 terms in Eq. (14) would yield corresponding terms that grow unbounded with $f$ . The term $T_{0}$ is free from the argument of the periapsis, and hence it can be cancelled by choosing the new Hamiltonian $K_{0,2}$ having a corresponding term $T_{0}$ . On the contrary, the terms $T_{1}$ and $T_{2}$ depend on the argument of the perigee, a fact that prevents their appearance in the new Hamiltonian under our requirement of making $g$ a cyclic variable. Nevertheless, since we had left $C_{1}$ arbitrary, both terms will cancel each other if $C_{1}$ is determined from the partial differential equation $T_{1}+T_{2}=0$ . We readily check that the particular solution

C_{1}=\epsilon{G}\frac{15s^{2}-14}{8(5s^{2}-4)}s^{2}e^{2}\sin 2\omega,

(20)

matches this purpose. Remark that the appearance of the divisor $5s^{2}-4$ prevents application of the solution to the resonant cases that happen at the so-called critical inclinations in which $\sin^{2}I=4/5$ . That is, the supplementary inclinations $I=63.435^{\circ}$ and $I=116.565^{\circ}$ .

After computing the partial derivatives of Eq. (20) with respect to $L$ and $G$ , which also appear among the coefficients of the terms remaining in $\widetilde{K}_{0,2}$ , the known terms of the homological equation become fully determined. Then, we can safely choose $K_{0,2}$ so that it comprises those terms of $\widetilde{K}_{0,2}$ that are free from the argument of the periapsis. We obtain,

K_{0,2}=\epsilon^{2}\frac{\mu}{r}\frac{p^{2}}{r^{2}}\frac{3s^{2}}{8(5s^{2}-4)^{2}}\sum_{j=0}^{2}\frac{p^{j}}{r^{j}}\sum_{k=0}^{\lfloor 1-\frac{1}{2}j\rfloor}e^{2k}\gamma_{2,j,k},

(21)

in which the inclination polynomials $\gamma_{2,j,k}\equiv\gamma_{2,j,k}(s)$ are provided in Table 1.

Table 1: Inclination polynomials

\gamma_{2,j,k}

in Eq. (21).

_0,0	$-8\left(200s^{6}-455s^{4}+345s^{2}-88\right)$
_0,1	$375s^{6}-930s^{4}+780s^{2}-224$
_1,0	$5\left(805s^{6}-1878s^{4}+1464s^{2}-384\right)$
_2,0	$-825s^{6}+1990s^{4}-1616s^{2}+448$

The second order term of the generating function is then readily computed from Eq. (14), yielding

	$\displaystyle W_{2}$	$\displaystyle=$	$\displaystyle\frac{\epsilon^{2}G}{32(5s^{2}-4)^{2}}\sum_{l=1}^{2}\sum_{\begin{subarray}{c}k=l-2\\ k\neq 0\end{subarray}}^{2l+2}\sum_{j=0}^{1}\Gamma_{2,j,k,l}e^{2j+k^{*}}$		(22)
			$\displaystyle\times s^{2l}\sin(kf+2l\omega)+C_{2},$		(22)

with $k^{*}\equiv{k}\bmod 2$ , and the inclination polynomials $\Gamma_{2,j,k,l}$ are given in Table 2.

Table 2: Non-zero inclination polynomials

\Gamma_{2,j,k,l}

in Eqs. (22)–(23).

_1,-1,1	$-12(5s^{2}-4)(7s^{2}-6)(15s^{2}-14)$
_0,1,1	$-48(5s^{2}-4)(195s^{4}-340s^{2}+148)$
_1,1,1	$24(5s^{2}-4)^{2}(15s^{2}-14)$
_1,1,2	$-3(225s^{4}-430s^{2}+208)$
_0,2,1	$-96(5s^{2}-4)^{2}(9s^{2}-8)$
_1,2,1	$-24(5s^{2}-4)(65s^{4}-116s^{2}+52)$
_1,2,2	$-60(50s^{4}-87s^{2}+38)$
_0,3,1	$-64(5s^{2}-4)^{2}(8s^{2}-7)$
_1,3,1	$4(3s^{2}-2)(5s^{2}-4)(15s^{2}-14)$
_0,3,2	$-4(5s^{2}-4)(135s^{2}-122)$
_1,3,2	$-8(75s^{4}-135s^{2}+61)$
_1,4,1	$-12(5s^{2}-4)^{2}(7s^{2}-6)$
_0,4,2	$24(5s^{2}-4)^{2}$
_1,4,2	$-12(5s^{2}-4)(25s^{2}-23)$
_0,5,2	$24(5s^{2}-4)^{2}$
_1,5,2	$-3(5s^{2}-4)(15s^{2}-14)$
_1,6,2	$6(5s^{2}-4)^{2}$
_0,0,2	$(15s^{2}-14)^{2}(15s^{2}-13)$
_0,0,1	$8(5s^{2}-4)^{2}(1215s^{4}-1997s^{2}+824)$
_1,0,1	$-2(5s^{2}-4)(15s^{2}-14)(45s^{4}+36s^{2}-56)$

Analogously as we did with $C_{1}$ , the arbitrary function $C_{2}$ in which Eq. (22) depends upon is determined at the next step of the perturbation approach in such a way that the appearance of unbounded terms in the solution of the homological equation of the third order are avoided. Thus, we compute $\widetilde{K}_{0,3}$ by successive evaluations of the fundamental recursion (11). After identifying the problematic terms of $\widetilde{K}_{0,3}$ , they are cancelled by computing

C_{2}=\frac{\epsilon^{2}G}{64(5s^{2}-4)^{3}}\sum_{l=1}^{2}\sum_{j=0}^{2-l}\Gamma_{2,j,0,l}e^{2(j+l)}s^{2l}\sin 2l\omega,

(23)

which contributes three additional inclination functions that are also listed in Table 2. In this way, the second order term of the generating function given by Eq. (22) becomes fully determined.

The needed partial derivatives of $C_{2}$ appearing in $\widetilde{K}_{0,3}$ are then computed, and the third order term of the new Hamiltonian is chosen, as before, to comprise those terms of $\widetilde{K}_{0,3}$ that are free from the argument of the periapsis. We obtain

K_{0,3}=\epsilon^{3}\frac{\mu}{r}\frac{p^{2}}{r^{2}}\frac{3s^{2}}{32(5s^{2}-4)^{3}}\sum_{j=0}^{4}\frac{p^{j}}{r^{j}}\sum_{k=0}^{\lfloor 2-\frac{1}{2}j\rfloor}e^{2k}\gamma_{3,j,k},

(24)

where the inclination polynomials $\gamma_{3,j,k}$ are displayed in Table 3.

Table 3: Inclination polynomials

\gamma_{3,j,k}

in Eq. (24).

_0,0	$-8(5s^{2}-4)(313525s^{8}-899030s^{6}+933656s^{4}$
	$-409296s^{2}+61824)$
_0,1	$4(1551625s^{10}-5675700s^{8}+8148960s^{6}-5706408s^{4}$
	$+1930272s^{2}-248064)$
_0,2	$-2(40500s^{10}-99525s^{8}+64840s^{6}+18788s^{4}$
	$-33936s^{2}+9408)$
_1,0	$2(5s^{2}-4)(2631475s^{8}-7558270s^{6}+7872692s^{4}$
	$-3470616s^{2}+530304)$
_1,1	$-3457125s^{10}+12282750s^{8}-17085020s^{6}$
	$+11554040s^{4}-3756000s^{2}+459648$
_2,0	$-2(5s^{2}-4)(1584375s^{8}-4536150s^{6}+4716436s^{4}$
	$-2082712s^{2}+321408)$
_2,1	$138375s^{10}-128250s^{8}-351900s^{6}+612440s^{4}$
	$-326368s^{2}+56448$
_3,0	$8(5s^{2}-4)(93300s^{8}-259915s^{6}+264982s^{4}-116928s^{2}$
	$+18816)$
_4,0	$-20s^{2}(5s^{2}-4)(15s^{2}-14)(45s^{4}+36s^{2}-56)$

Then, after solving Eq. (14), we find that the third order term of the generating function can be arranged in the form

	$\displaystyle W_{3}$	$\displaystyle=$	$\displaystyle\frac{\epsilon^{3}G}{8960(5s^{2}-4)^{4}}\sum_{l=1}^{3}\sum_{\begin{subarray}{c}k=l-4\\ k\neq 0\end{subarray}}^{2l+4}\sum_{j=0}^{2}\Gamma_{3,j,k,l}e^{2j+k^{*}}$		(25)
			$\displaystyle\times s^{2l}\sin(kf+2lg)+C_{3}$		(25)

in which the inclination polynomials $\Gamma_{3,j,k,l}$ that do not vanish are listed in Table 4. The constant $C_{3}$ is determined at the next order of the perturbation approach in the same way as we did previously. We obtain

C_{3}=\frac{\epsilon^{3}G}{1536(5s^{2}-4)^{5}}\sum_{l=1}^{3}\sum_{j=0}^{3-l}\Gamma_{3,j,0,l}e^{2(j+l)}s^{2l}\sin 2l\omega,

which contributes six additional inclination polynomials $\Gamma_{3,j,0,l}$ , that are also listed in Table 4.

Table 4: Non-zero inclination polynomials

\Gamma_{3,j,k,l}

in Eq. (25).

_2,-3,1	$35\left(15s^{2}-14\right)\left(87375s^{10}-335550s^{8}+505080s^{6}-371184s^{4}+132096s^{2}-17920\right)$
_2,-2,1	$105\left(15s^{2}-14\right)\left(399375s^{10}-1863400s^{8}+3389440s^{6}-3023632s^{4}+1328128s^{2}-230400\right)$
_1,-1,1	$-840\left(5s^{2}-4\right)\left(228125s^{10}-549325s^{8}+255940s^{6}+324664s^{4}-352992s^{2}+93824\right)$
_2,-1,1	$210\left(15s^{2}-14\right)\left(100875s^{10}-275600s^{8}+228220s^{6}-6408s^{4}-70272s^{2}+23296\right)$
_2,-1,2	$-105\left(5s^{2}-4\right)\left(15s^{2}-14\right)\left(13725s^{6}-37680s^{4}+34228s^{2}-10304\right)$
_0,1,1	$-1680\left(5s^{2}-4\right)^{2}\left(486525s^{8}-1594290s^{6}+1955772s^{4}-1064576s^{2}+216960\right)$
_1,1,1	$840\left(5s^{2}-4\right)\left(1531125s^{10}-6503075s^{8}+10982780s^{6}-9224760s^{4}+3855648s^{2}-641920\right)$
_2,1,1	$-420\left(15s^{2}-14\right)\left(61875s^{10}-138825s^{8}+51640s^{6}+92200s^{4}-89088s^{2}+22400\right)$
_1,1,2	$1680\left(5s^{2}-4\right)\left(240750s^{8}-775475s^{6}+932445s^{4}-495822s^{2}+98320\right)$
_2,1,2	$105\left(5s^{2}-4\right)\left(226125s^{8}-787950s^{6}+1015020s^{4}-572056s^{2}+118944\right)$
_2,1,3	$-840\left(15s^{2}-14\right)\left(1125s^{6}-3300s^{4}+3235s^{2}-1058\right)$
_0,2,1	$-1680\left(5s^{2}-4\right)^{3}\left(41615s^{6}-97838s^{4}+76016s^{2}-19488\right)$
_1,2,1	$-1680\left(5s^{2}-4\right)\left(666875s^{10}-2586600s^{8}+4014940s^{6}-3117320s^{4}+1210368s^{2}-187904\right)$
_2,2,1	$210\left(5398125s^{12}-27480750s^{10}+57999400s^{8}-64973520s^{6}+40757888s^{4}-13579264s^{2}+1878016\right)$
_1,2,2	$-3360\left(5s^{2}-4\right)^{2}\left(42850s^{6}-108830s^{4}+92099s^{2}-25958\right)$
_2,2,2	$840\left(5s^{2}-4\right)\left(123750s^{8}-359475s^{6}+378010s^{4}-167724s^{2}+25624\right)$
_2,2,3	$-105\left(639375s^{8}-2259750s^{6}+2991200s^{4}-1757840s^{2}+387104\right)$
_0,3,1	$-1120\left(5s^{2}-4\right)^{2}\left(270650s^{8}-828285s^{6}+945816s^{4}-477232s^{2}+89664\right)$
_1,3,1	$280\left(5s^{2}-4\right)\left(634500s^{10}-2623725s^{8}+4300340s^{6}-3496152s^{4}+1412352s^{2}-227584\right)$
_2,3,1	$-70\left(15s^{2}-14\right)\left(76875s^{10}-202950s^{8}+167700s^{6}-16544s^{4}-37376s^{2}+12544\right)$
_0,3,2	$-560\left(5s^{2}-4\right)^{2}\left(605775s^{6}-1524950s^{4}+1277728s^{2}-356256\right)$
_1,3,2	$560\left(5s^{2}-4\right)^{2}\left(9150s^{6}-6435s^{4}-9741s^{2}+7154\right)$
_2,3,2	$35\left(5s^{2}-4\right)\left(104625s^{8}-68850s^{6}-286620s^{4}+378936s^{2}-127232\right)$
_1,3,3	$-280\left(5s^{2}-4\right)\left(52875s^{6}-129225s^{4}+102900s^{2}-26456\right)$
_2,3,3	$-140\left(15s^{2}-14\right)\left(7875s^{6}-21225s^{4}+19120s^{2}-5756\right)$
_1,4,1	$-840\left(5s^{2}-4\right)\left(516875s^{10}-1956050s^{8}+2950940s^{6}-2217472s^{4}+829440s^{2}-123392\right)$
_2,4,1	$420\left(5s^{2}-4\right)\left(173625s^{10}-668250s^{8}+1023720s^{6}-780120s^{4}+295872s^{2}-44800\right)$
_0,4,2	$-6720\left(5s^{2}-4\right)^{3}\left(730s^{4}-1153s^{2}+444\right)$
_1,4,2	$-3360\left(5s^{2}-4\right)^{2}\left(66050s^{6}-166215s^{4}+139230s^{2}-38808\right)$
_2,4,2	$420\left(5s^{2}-4\right)^{2}\left(20925s^{6}-38700s^{4}+19984s^{2}-1976\right)$
_1,4,3	$840\left(5s^{2}-4\right)\left(10125s^{6}-32150s^{4}+33500s^{2}-11456\right)$
_2,4,3	$-2100\left(5s^{2}-4\right)\left(3525s^{6}-9240s^{4}+7996s^{2}-2280\right)$
_1,5,1	$-168\left(5s^{2}-4\right)\left(115000s^{10}-434875s^{8}+663700s^{6}-512080s^{4}+199936s^{2}-31488\right)$
_2,5,1	$-105s^{2}\left(5s^{2}-4\right)\left(15s^{2}-14\right)\left(825s^{6}-1990s^{4}+1616s^{2}-448\right)$
_0,5,2	$-336\left(5s^{2}-4\right)^{3}\left(15425s^{4}-24050s^{2}+9112\right)$
_1,5,2	$-84\left(5s^{2}-4\right)^{2}\left(551625s^{6}-1390850s^{4}+1167040s^{2}-325728\right)$
_2,5,2	$105\left(5s^{2}-4\right)^{2}\left(15s^{2}-14\right)\left(225s^{4}+288s^{2}-364\right)$
_0,5,3	$3360\left(5s^{2}-4\right)^{2}\left(1575s^{4}-2795s^{2}+1256\right)$
_1,5,3	$840\left(5s^{2}-4\right)\left(8250s^{6}-24975s^{4}+25080s^{2}-8336\right)$
_2,5,3	$-420\left(5s^{2}-4\right)\left(15s^{2}-14\right)^{2}\left(15s^{2}-13\right)$
_2,6,1	$35\left(5s^{2}-4\right)\left(171375s^{10}-616950s^{8}+871680s^{6}-600128s^{4}+199168s^{2}-25088\right)$
_1,6,2	$-280\left(5s^{2}-4\right)^{3}\left(8265s^{4}-12874s^{2}+4872\right)$
_2,6,2	$-280\left(5s^{2}-4\right)^{2}\left(11475s^{6}-29280s^{4}+24828s^{2}-6992\right)$
_0,6,3	$560\left(5s^{2}-4\right)^{3}\left(1335s^{2}-1166\right)$
_1,6,3	$560\left(5s^{2}-4\right)^{2}\left(8325s^{4}-14910s^{2}+6764\right)$
_2,6,3	$70\left(5s^{2}-4\right)\left(21375s^{6}-61950s^{4}+59960s^{2}-19328\right)$
_1,7,2	$-60\left(5s^{2}-4\right)^{3}\left(8385s^{4}-13226s^{2}+5080\right)$
_0,7,3	$10080\left(5s^{2}-4\right)^{3}\left(90s^{2}-79\right)$
_1,7,3	$12600\left(5s^{2}-4\right)^{2}\left(105s^{4}-191s^{2}+88\right)$
_2,8,2	$-840\left(3s^{2}-2\right)\left(5s^{2}-4\right)^{3}\left(15s^{2}-14\right)$
_1,8,3	$4200\left(5s^{2}-4\right)^{3}\left(96s^{2}-85\right)$
_2,8,3	$210\left(5s^{2}-4\right)^{2}\left(525s^{4}-990s^{2}+472\right)$
_1,9,3	$7560\left(5s^{2}-4\right)^{3}\left(10s^{2}-9\right)$
_2,10,3	$315\left(5s^{2}-4\right)^{3}\left(15s^{2}-14\right)$
_0,0,3	$2\left(15s^{2}-14\right)^{3}\left(825s^{4}-1445s^{2}+634\right)$
_0,0,2	$-6\left(5s^{2}-4\right)^{2}\left(2171250s^{8}-7719525s^{6}+10225470s^{4}-5983260s^{2}+1305248\right)$
_1,0,2	$-3\left(5s^{2}-4\right)\left(15s^{2}-14\right)^{2}\left(1800s^{6}+2655s^{4}-8208s^{2}+3928\right)$
_0,0,1	$48\left(5s^{2}-4\right)^{2}\left(9060750s^{10}-34431275s^{8}+51858720s^{6}-38675200s^{4}+14258176s^{2}-2072064\right)$
_1,0,1	$-12\left(5s^{2}-4\right)\left(93223125s^{12}-421210500s^{10}+784654200s^{8}-771469840s^{6}+422629664s^{4}-122600960s^{2}+14780416\right)$
_2,0,1	$6\left(15s^{2}-14\right)\left(2328750s^{12}-8703375s^{10}+13317150s^{8}-10848180s^{6}+5157560s^{4}-1450624s^{2}+200704\right)$

The computation of additional orders finds similar structures. Thus, for instance, the fourth order term of the new Hamiltonian takes the form

K_{0,4}=\epsilon^{4}\frac{\mu}{r}\frac{p^{2}}{r^{2}}\frac{9s^{2}}{1280(5s^{2}-4)^{6}}\sum_{j=0}^{5}\frac{p^{j}}{r^{j}}\sum_{k=0}^{\lfloor 3-\frac{1}{2}j\rfloor}e^{2k}\gamma_{4,j,k},

(26)

and the fourth order term of the generating function takes the form

	$\displaystyle W_{4}$	$\displaystyle=$	$\displaystyle\frac{\epsilon^{4}G}{30105600(5s^{2}-4)^{6}}\sum_{l=1}^{4}\sum_{\begin{subarray}{c}k=l-6\\ k\neq 0\end{subarray}}^{2l+6}\sum_{j=0}^{6}\Gamma_{4,j,k,l}e^{2j+k^{*}}$		(27)
			$\displaystyle\times s^{2l}\sin(kf+2l\omega)+C_{4}.$		(27)

While the former involves 15 inclination polynomials $\gamma_{4,j,k}$ , the later comprises up to 124 non-vanishing inclination polynomials $\Gamma_{4,j,k,l}$ , of degree 9 in $s^{2}$ . Therefore, these and other polynomials resulting from following orders are not listed due to their length.

The normalization of $G$ has been extended up to the order six without major trouble, except for the increase of the computational burden of successive orders due to the notable growth of the length of the series to be handled, on the one hand, and the increasing size of the rational coefficients resulting from the integer arithmetic used, on the other. Thus, the new Hamiltonian terms $K_{0,5}$ and $K_{0,6}$ take analogous forms to the previous orders, with 24 inclination polynomials $\gamma_{5,j,k}$ , and 35 $\gamma_{6,j,k}$ . Similarly, the generating function terms $W_{5}$ and $W_{6}$ have been arranged in the same form as previous orders of the perturbation approach, with 254 non-vanishing coefficients $\Gamma_{5,j,k,l}$ and 429 $\Gamma_{6,j,k,l}$ , 15 of which correspond to the integration constant $C_{5}$ and 21 to $C_{6}$ .

At the end, since the generating function is known, the transformation equations from the original Delaunay variables to the new ones, and vice versa, are readily computed by standard application of the fundamental recursion (11) (see Deprit1969 for details). The procedure ends by replacing the original variables by the new ones in the computed terms $K_{0,m}$ .

We remark that the procedure described here is not equivalent to the combination of the elimination of the parallax and the elimination of the perigee into a single transformation, in spite of the fact that the total angular momentum is converted into a formal integral in both cases. On the contrary, the perigee has been removed here keeping as much short-period terms as possible in the new Hamiltonian. With this strategy, the size of the generating function of the first partial normalization is astonishingly smaller than the one that would be obtained with other alternatives in the literature. To check that, we fully expanded both the Hamiltonian and the generating function and reckoned the number of terms. This procedure is the one that has been traditionally used in the literature to assess the complexity of a perturbation solution DepritRom1970 ; CoffeyDeprit1982 . Thus, for instance, we find 5, 56, 367, 1152, 2627, and 4897 expanded terms for the 1st, 2nd …, 6th order terms of the generating function, respectively, whereas much longer expressions have been reported in the literature. For instance, the sixth order of the generating function reported in SanJuanetal2013 entails 39630 terms, which is almost one order of magnitude larger than the one obtained with the current approach. It is worth mentioning that further simplifications could be obtained in particular cases, like when constraining the application of the analytical solution to the case of low eccentricity orbits, in which case many of the terms involved can be neglected Lara2008 ; GaiasColomboLara2020 .

The comparisons are, nevertheless, inconclusive due to the diversity of approaches used in the computation of the variety of solutions reported in the literature, which involve representation in different variables, on the one hand, and yield distinct one-degree-of-freedom Hamiltonians, on the other. The partially normalized Hamiltonian obtained here is definitely longer than the one obtained after the classical elimination of the perigee. However, this is not of concern in the complete, as opposite to partial, normalization of the main problem. Indeed, after the following elimination of short-period terms either procedure should arrive to the same Hamiltonian. We didn’t find reported data for the second normalization after the elimination of the perigee, but one would expect that the figures should be balanced to some extent —the generating function of the short-period elimination probably being heavier in our case than in other prospective approaches. While the needed data for the thorough comparison is not available, our claims must restrict to the checked fact that our approach makes both transformations of manageable size —as we will show in the next section, where the short-period effects are removed in a Delaunay normalization. This feature makes the evaluation of higher orders of our perturbation solution certainly practicable.

2.3 Delaunay normalization

The partially normalized Hamiltonian is just of one degree of freedom, yet it is not separable. To get an explicit analytical solution we remove the short-period terms that are associated to the radius by means of a Delaunay normalization Deprit1982 . The partially reduced Hamiltonian from which we start takes again the form of Eq. (4), but it is now written in prime Delaunay variables $(\ell^{\prime},g^{\prime},h^{\prime},L^{\prime},G^{\prime},H^{\prime})$ . In this Hamiltonian the variables $g^{\prime}$ and $h^{\prime}$ are cyclic, and hence $H^{\prime}=H$ and $G^{\prime}$ remain constant for given initial conditions. The zeroth order term is the same as in Eq. (5), the term $K_{1,0}$ is given by Eq. (15), whereas terms $K_{m,0}$ with $m\geq 2$ are no longer void, and, on the contrary are given by Eqs. (21), (24), (26), for $m=2,3,4$ , respectively, and analogous equations for higher orders that, due to its length, are not printed in the paper. Moreover, the definition in Eq. (17) turns $\epsilon$ into a physical constant, rather than a function, when witten in the prime variables. In consequence $\epsilon$ can replace now the formal small parameter $\varepsilon$ used before in the Lie transforms procedure, assumed, of course, that the corresponding scaling of the Hamiltonian terms is properly made.

The homological equation to be solved at each step of the new Lie transformation is the same as before, either in the form given by Eq. (12) or Eq. (14), except for choosing a non-zero integration constant does not provide any advantage in the current case. In fact, since the new Hamiltonian must be free from terms depending on the mean anomaly, which are obtained by averaging, the homological equation can be further particularized. That is,

K_{0,m}=\langle\widetilde{K}_{0,m}\rangle_{\ell}=\frac{1}{2\pi}\int_{0}^{2\pi}\widetilde{K}_{0,m}\frac{r^{2}}{a^{2}\eta}\,\mathrm{d}f,

(28)

and hence

	$\displaystyle W_{m}$	$\displaystyle=$	$\displaystyle-\frac{\ell}{n}K_{0,m}+\frac{1}{n}\int\widetilde{K}_{0,m}\frac{r^{2}}{a^{2}\eta}\,\mathrm{d}f$		(29)
		$\displaystyle=$	$\displaystyle K_{0,m}\frac{\phi}{n}+\frac{1}{n}\int\left(\widetilde{K}_{0,m}\frac{r^{2}}{a^{2}\eta}-K_{0,m}\right)\mathrm{d}f,$		(29)

in which $\phi=f-\ell$ is the equation of the center, and the terms under the integral sign in the final form of the equation are purely periodic in $f$ .

Thus, at the first order of the Lie transforms procedure we find

K_{0,1}=\langle{K}_{1,0}\rangle_{\ell}\equiv\frac{\mu}{p}\eta^{3}\left(3s^{2}-2\right),

(30)

from which

W_{1}=G^{\prime}\left(3s^{2}-2\right)(e\sin{f}+\phi).

(31)

Recall that the symbols $p$ , $\eta$ , $e$ , etc. are now functions of the Delaunay prime variables.

At the second order, the known terms are given again by Eqs. (18) and (19), from which, using Eq. (28),

K_{0,2}=-\frac{\mu}{p}\frac{3}{4}\eta^{3}\sum_{j=0}^{2}\lambda_{2,j}\eta^{j},

(32)

with $\lambda_{2,0}=5(7s^{4}-16s^{2}+8)$ , $\lambda_{2,1}=4(3s^{2}-2)^{2}$ , and $\lambda_{2,2}=5s^{4}+8s^{2}-8$ . The second order term of the generating function is then trivially integrated from Eq. (29), to yield

	$\displaystyle W_{2}$	$\displaystyle=$	$\displaystyle-\frac{G^{\prime}\beta}{32(5s^{2}-4)^{2}}\sum_{j=1}^{3}\sum_{k=0}^{3-\lfloor\frac{1}{2}j\rfloor}\Lambda_{2,j,k}\eta^{k}e^{j}\sin{jf}$		(33)
			$\displaystyle-G^{\prime}\frac{3}{4}\phi\sum_{j=0}^{1}\Phi_{2,j}e^{2j},$		(33)

where $\beta=1/(1+\eta)$ , $\Phi_{2,0}=8(s^{2}-1)(5s^{2}-4)$ , $\Phi_{2,1}=8-8s^{2}-5s^{4}$ , and the inclination polynomials $\Lambda_{2,j,k}$ are provided in Table 5.

Table 5: Inclination polynomials

\Lambda_{2,j,k}

in Eq. (33).

_1,0	$15(3s^{2}-2)(805s^{6}-2448s^{4}+2400s^{2}-768)$
_1,1	$3(3s^{2}-2)(2225s^{6}-8160s^{4}+8928s^{2}-3072)$
_1,2	$3(825s^{8}-3030s^{6}+4064s^{4}-2368s^{2}+512)$
_1,3	$-3s^{2}(975s^{6}-2250s^{4}+1728s^{2}-448)$
_2,0	$6(1925s^{8}-6210s^{6}+7452s^{4}-3936s^{2}+768)$
_2,1	$6(125s^{8}-930s^{6}+1660s^{4}-1120s^{2}+256)$
_3,0	$2625s^{8}-7270s^{6}+7408s^{4}-3264s^{2}+512$
_3,1	$s^{2}(825s^{6}-1990s^{4}+1616s^{2}-448)$

Integrals cease to be trivial at the third order. Indeed, the formal computation of $K_{0,3}$ using the fundamental recursion (11) yields two different types of terms.

The first type consists of terms depending on the equation of the center, which can be reduced to the form

\frac{\mu}{r}\frac{p}{r}P(e)Q(s)\phi\sin{mf},

with $P$ and $Q$ denoting arbitrary eccentricity and inclination polynomials. Definite integration of these kinds of terms is carried out from expressions in Metris1991 , whereas the indefinite integration is achieved by parts, to get

ma^{2}\eta\int\frac{\sin{mf}}{r^{2}}\phi\,\mathrm{d}\ell=\frac{\sin{mf}}{m}-\phi\cos{mf}-\int\cos{mf}\mathrm{d}\ell,

in which the antiderivatives of cosines of multiples of the true anomaly are carried out after expressing them in terms of $r$ and $R=\mathrm{d}r/\mathrm{d}t$ , rather than in trigonometric functions, as discussed in Jefferys1971 .

The second type consists of terms free from the equation of the center. In these terms, the trigonometric functions of the true anomaly can be replaced by inverse powers of the radius, without involving the radial velocity. We found that the exponents of the inverse of $r$ range from 0 to 8 missing the exponent 1. Therefore, both definite and indefinite integrals of terms of the second type are trivially solved.

Proceeding in this way, we compute the third order Hamiltonian term

K_{0,3}=\frac{\mu}{p}\frac{9\eta^{3}}{16(5s^{2}-4)^{2}}\sum_{j=0}^{4}\lambda_{3,j}\eta^{j},

(34)

with the inclination polynomials $\lambda_{3,j}$ given in Table 6.

Table 6: Inclination polynomials

\lambda_{3,j}

in Eq. (34).

₀	$5(28700s^{10}-107205s^{8}+158960s^{6}-118492s^{4}$
	$+45152s^{2}-7168)$
₁	$60(3s^{2}-2)(5s^{2}-4)^{2}(7s^{4}-16s^{2}+8)$
₂	$-2(28675s^{10}-98005s^{8}+130852s^{6}-87164s^{4}$
	$+30176s^{2}-4608)$
₃	$20(3s^{2}-2)(5s^{2}-4)^{2}(5s^{4}+8s^{2}-8)$
₄	$-s^{2}(15s^{2}-14)(450s^{6}-925s^{4}+590s^{2}-112)$

The corresponding term of the generating function is

	$\displaystyle W_{3}$	$\displaystyle=$	$\displaystyle\frac{G^{\prime}\beta^{2}}{128(5s^{2}-4)^{3}}\sum_{j=1}^{6}\sum_{k=0}^{7-2\lfloor\frac{1}{2}j\rfloor}\Lambda_{3,j,k}\eta^{k-1}e^{j}\sin{jf}$
			$\displaystyle+\frac{3G^{\prime}}{16(5s^{2}-4)^{2}}\phi\sum_{j=0}^{3}\sum_{k=0}^{4}\Phi_{3,j,k}\eta^{k}e^{j}\cos{jf},$		(35)

in which the inclination polynomials $\Phi_{3,j,k}$ are given in Table 7, and those $\Lambda_{3,j,k}$ in Table 8.

Table 7: Non-zero coefficients

\Phi_{3,j,k}

in Eq. (35).

\Phi_{3,1,0}=\frac{15}{2}\Phi_{3,0,3}

\Phi_{3,1,2}=-\frac{3}{2}\Phi_{3,0,3}

\Phi_{3,2,0}=3\Phi_{3,0,3}

, and

\Phi_{3,3,0}=\frac{1}{2}\Phi_{3,0,3}

_0,0	$5(89100s^{10}-323615s^{8}+466320s^{6}-337684s^{4}$
	$+125216s^{2}-19456)$
_0,2	$-2(112125s^{10}-374775s^{8}+488460s^{6}-314932s^{4}$
	$+103840s^{2}-14848)$
_0,3	$8(3s^{2}-2)(5s^{2}-4)^{2}(5s^{4}+8s^{2}-8)$
_0,4	$-3s^{2}(15s^{2}-14)(450s^{6}-925s^{4}+590s^{2}-112)$

Table 8: Non-zero inclination polynomials

\Lambda_{3,j,k}

in Eq. (35).

_1,0	$-864(3s^{2}-2)^{3}(5s^{2}-4)^{3}$
_1,1	$3(22218875s^{12}-104346550s^{10}+202703740s^{8}-209869352s^{6}+123038240s^{4}-39033472s^{2}+5275648)$
_1,2	$12(10925500s^{12}-50711075s^{10}+97386820s^{8}-99715748s^{6}+57863024s^{4}-18199872s^{2}+2445312)$
_1,3	$3(27560125s^{12}-119080550s^{10}+212650740s^{8}-202245448s^{6}+109190304s^{4}-32208768s^{2}+4128768)$
_1,4	$12(3155125s^{12}-10820800s^{10}+13899620s^{8}-7620256s^{6}+944256s^{4}+613248s^{2}-167936)$
_1,5	$3(5410625s^{12}-17331450s^{10}+19448180s^{8}-6842968s^{6}-2742560s^{4}+2556544s^{2}-491520)$
_1,6	$-12(59625s^{12}-415275s^{10}+942920s^{8}-994980s^{6}+529776s^{4}-134592s^{2}+12288)$
_1,7	$-3s^{2}(77625s^{10}-568950s^{8}+1256420s^{6}-1222216s^{4}+550816s^{2}-94080)$
_2,0	$-1044(3s^{2}-2)^{3}(5s^{2}-4)^{3}$
_2,1	$24(131000s^{12}-1121875s^{10}+3061340s^{8}-3989664s^{6}+2758768s^{4}-983648s^{2}+143360)$
_2,2	$96(51625s^{12}-437800s^{10}+1183290s^{8}-1528682s^{6}+1049344s^{4}-372240s^{2}+54144)$
_2,3	$-12(5s^{2}-4)(16375s^{10}+64070s^{8}-257508s^{6}+297320s^{4}-145792s^{2}+26624)$
_2,4	$-12(263625s^{12}-1000750s^{10}+1526820s^{8}-1206712s^{6}+539616s^{4}-142592s^{2}+19968)$
_2,5	$-12s^{2}(162375s^{10}-576100s^{8}+787020s^{6}-506248s^{4}+145984s^{2}-12992)$
_3,0	$-656(3s^{2}-2)^{3}(5s^{2}-4)^{3}$
_3,1	$-934875s^{12}+605000s^{10}+4973120s^{8}-10412952s^{6}+8554272s^{4}-3281280s^{2}+491520$
_3,2	$-2080125s^{12}+3562000s^{10}+2881300s^{8}-11103360s^{6}+10155456s^{4}-4038912s^{2}+614400$
_3,3	$-(5s^{2}-4)(254925s^{10}-526480s^{8}+318084s^{6}-10672s^{4}-40064s^{2}+8192)$
_3,4	$3(28875s^{12}-147800s^{10}+254260s^{8}-160992s^{6}-11968s^{4}+54016s^{2}-16384)$
_3,5	$-12s^{2}(4500s^{10}-4125s^{8}-16075s^{6}+31670s^{4}-20664s^{2}+4704)$
_4,0	$-240(3s^{2}-2)^{3}(5s^{2}-4)^{3}$
_4,1	$6(5s^{2}-4)(50325s^{10}-157660s^{8}+180520s^{6}-90312s^{4}+18112s^{2}-1024)$
_4,2	$12(5s^{2}-4)(47475s^{10}-147000s^{8}+164852s^{6}-79056s^{4}+14208s^{2}-512)$
_4,3	$6s^{2}(5s^{2}-4)(44625s^{8}-136340s^{6}+149184s^{4}-67800s^{2}+10304)$
_5,0	$-48(3s^{2}-2)^{3}(5s^{2}-4)^{3}$
_5,1	$6s^{2}(5s^{2}-4)^{2}(180s^{6}-609s^{4}+530s^{2}-112)$
_5,2	$3s^{2}(5s^{2}-4)(1125s^{8}-7440s^{6}+11516s^{4}-6144s^{2}+896)$
_5,3	$-3s^{4}(5s^{2}-4)(15s^{2}-14)(45s^{4}+36s^{2}-56)$
_6,0	$-4(3s^{2}-2)^{3}(5s^{2}-4)^{3}$

In the process of carrying out the normalization to higher orders we need to deal with trigonometric series of notably increasing length. However, we only found terms of the same two types as before in the solution of the homological equation, which, in consequence, is analogously integrated. In this way, the second normalization has been extended up to the order 6 in the small parameter, in agreement with the order to which the perigee was previously eliminated. Corresponding Hamiltonian and generating function terms are analogously arranged in the form of Eqs. (34) and (35), respectively, yet the inclination polynomials comprise much more longer listings, as expected. Thus, for instance, at the fourth order the normalized Hamiltonian term

K_{0,4}=\frac{\mu}{p}\frac{9\eta^{3}}{64(5s^{2}-4)^{3}}\sum_{j=0}^{6}\lambda_{4,j}\eta^{j},

(36)

contributes 7 new inclination polynomials, which are polynomials of degree 7 in the square of the sine of the inclination, whereas the generating function term

	$\displaystyle W_{4}$	$\displaystyle=$	$\displaystyle\frac{G^{\prime}\beta^{3}}{20480(5s^{2}-4)^{6}}\sum_{j=1}^{6}\sum_{k=0}^{7-2\lfloor\frac{1}{2}j\rfloor}\Lambda_{4,j,k}\eta^{k-3}e^{j}\sin{jf}$		(37)
			$\displaystyle+\frac{3G^{\prime}}{256(5s^{2}-4)^{3}}\phi\sum_{j=0}^{5}\sum_{k=0}^{6}\Phi_{4,j,k}\eta^{k}e^{j}\cos{jf},$		(37)

is made of 20 non-zero coefficients $\Phi_{4,j,k}$ , of degree 7 in $s^{2}$ , and 72 non-zero coefficients $\Lambda_{4,j,k}$ , which are of degree 10 in $s^{2}$ . Note that, because of the denominators factoring the summations, the maximum degree amounts to 4 in any case, in agreement with the order of the perturbation. At the 5th order we find 9 coefficients of the type $\lambda$ and 41 nonvanishing coefficients of the type $\Phi$ , which are of degree 11 in $s^{2}$ although they are divided by $(5s^{2}-4)^{6}$ , as well as 134 nonvanishing trigonometric polynomials of the type $\Lambda$ , which are of degree 12 in $s^{2}$ yet they must be divided by $(5s^{2}-4)^{7}$ . Finally, the figures of the 6th order are 11 and 70 of degree 13 in $s^{2}$ for $\lambda$ and $\Phi$ , respectively, and 218 for $\Lambda$ of degree 16.

Once reached the desired order of the second normalization, the procedure ends by changing prime variables by double-prime variables in the new Hamiltonian terms $K_{0,m}$ .

In order to provide comparative figures of the computational burden of this normalization with other approaches that might be carried out, we expanded the series that comprise the solution and reckon the number of separate terms, as we already did in the normalization of the total angular momentum. We found that the 1st, 2nd, …6th-order terms of the generating function of the Delaunay normalization comprise 4, 48, 257, 931, 2266, and 4826 terms respectively. Note that different arrangements from the one chosen by us —Eqs. (35), (37), etc.— may provide different figures than those reported, yet should be of analogous magnitude. Following the same procedure with the completely reduced Hamiltonian, we reckon up to 2, 9, 29, 55, 106, and 152 coefficients, for the 1st, 2nd, …6th-order term, respectively.

2.4 Secular terms

After neglecting higher order effects of $J_{2}$ , the normalized Hamiltonian is

\mathcal{K}^{\prime\prime}=\mathcal{K}^{\prime\prime}(-,-,-,L^{\prime\prime},G^{\prime\prime},H^{\prime\prime})\equiv\sum_{m\geq 0}^{\tilde{m}}\frac{\epsilon^{m}}{m!}K_{m,0},

(38)

in which $\tilde{m}\leq 6$ for the different approximations provided by the computed perturbation solution, and terms $K_{m,0}$ are obtained by replacing prime by double-prime variables in corresponding terms given by Eqs. (30), (32), (34), as well as higher order Hamiltonian terms $K_{0,m}$ that have not been displayed. That is, the symbols $p$ , $\eta$ , and $s$ in these equations are assumed to be functions of the double-prime Delaunay momenta. Recall that $\epsilon$ is obtained from Eq. (17) by making $p=G^{\prime\prime 2}/\mu$ .

The corresponding solution to the flow stemming from Eq. (38) is

$\displaystyle\ell^{\prime\prime}$	$\displaystyle=$	$\displaystyle\ell^{\prime\prime}_{0}+n_{\ell}{t},$
$\displaystyle g^{\prime\prime}$	$\displaystyle=$	$\displaystyle g^{\prime\prime}_{0}+{n}_{\omega}{t}$
$\displaystyle h^{\prime\prime}$	$\displaystyle=$	$\displaystyle h^{\prime\prime}_{0}+n_{h}t,$

in which, from Hamilton equations,

n_{\ell}=\frac{\partial\mathcal{K}^{\prime\prime}}{\partial{L}^{\prime\prime}},\qquad n_{g}=\frac{\partial\mathcal{K}^{\prime\prime}}{\partial{G}^{\prime\prime}},\qquad n_{h}=\frac{\partial\mathcal{K}^{\prime\prime}}{\partial{H}^{\prime\prime}},

$L^{\prime\prime}=L^{\prime\prime}_{0}$ , $G^{\prime\prime}=G^{\prime}_{0}$ , and $H^{\prime\prime}=H_{0}$ , are the initialization constants of the analytical solution, and $\ell^{\prime\prime}_{0}$ , $g^{\prime\prime}_{0}$ , $h^{\prime\prime}_{0}$ , $L^{\prime\prime}_{0}$ , and $G^{\prime}_{0}$ , are computed applying consecutively the inverse transformation of the normalization of the angular momentum and the Delaunay normalization to corresponding initial conditions of the original problem.

On the other hand, Delaunay variables are singular for equatorial orbits, in which the argument of the node is not defined, as well as for circular orbits, in which the argument of the periapsis is not defined. In particular, the singularity for circular orbits reveals immediately by the appearance of the eccentricity in denominators of the transformation equations of both the mean anomaly and the argument of the perigee (see Eqs. (20) and (21) of Brouwer1959 , for instance). This fact not only makes singular the transformation of these elements for exactly circular orbits, but prevents convergence of the respective perturbation series for small values of the eccentricity.

Different sets of non-singular variables can be used to avoid these issues Lyddane1963 . In particular, troubles related with low eccentricities are commonly avoided by replacing the mean anomaly, the argument of the periapsis, and the total angular momentum with the non-canonical variables given by the mean distance to the node $F=\ell+g$ , also called mean argument of the latitude, and the components of the eccentricity vector in the nodal frame $C=e\cos\omega$ , $S=e\sin\omega$ , also called semi-equinoctial elements Konopliv1990 ; Cook1992 . In these variables, the secular terms of the main problem are given by DepritRom1970

$\displaystyle F$	$\displaystyle=$	$\displaystyle F_{0}+n_{F}t,$	(39)
$\displaystyle C$	$\displaystyle=$	$\displaystyle e\cos(g_{0}+n_{g}{t})=C_{0}\cos{n}_{g}{t}-S_{0}\sin{n}_{g}{t},$	(40)
$\displaystyle S$	$\displaystyle=$	$\displaystyle e\sin(g_{0}+n_{g}{t})=S_{0}\cos{n}_{g}{t}+C_{0}\sin{n}_{g}{t},$	(41)
$\displaystyle L$	$\displaystyle=$	$\displaystyle L_{0},$	(42)
$\displaystyle h$	$\displaystyle=$	$\displaystyle h_{0}+n_{h}t,$	(43)
$\displaystyle H$	$\displaystyle=$	$\displaystyle H_{0},$	(44)

in which $n_{F}=n_{\ell}+n_{g}$ , $C_{0}=e\cos{g}_{0}$ , $S_{0}=e\sin{g}_{0}$ , and $e=(1-G_{0}^{2}/L_{0}^{2})^{1/2}$ . The double-prime notation has been omitted for simplicity.

After standard partial differentiation, we obtain

$\displaystyle n_{F}$	$\displaystyle=$	$\displaystyle n+n\sum_{m=1}^{\tilde{m}}\frac{\epsilon^{m}}{(5s^{2}-4)^{m}}\,\sum_{i=0}^{2m-1}\Psi_{m,i}(s)\eta^{i},$	(45)
$\displaystyle n_{g}$	$\displaystyle=$	$\displaystyle n\sum_{m=1}^{\tilde{m}}\frac{\epsilon^{m}}{(5s^{2}-4)^{m}}\sum_{i=0}^{2m-2}\omega_{m,i}(s)\eta^{i},$	(46)
$\displaystyle n_{h}$	$\displaystyle=$	$\displaystyle n{c}\sum_{m=1}^{\tilde{m}}\frac{\epsilon^{m}}{(5s^{2}-4)^{m}}\sum_{i=0}^{2m-2}\Omega_{m,i}(s)\eta^{i},$	(47)

where the inclination polynomials $\Psi_{m,i}$ , $\omega_{m,i}$ , and $\Omega_{m,i}$ are provided in Tables 9, 10, and 11, respectively, up to the third order of the perturbation approach. In these tables we can check that, in spite of the general arrangement of the secular frequencies in Eqs. (45)–(47), the critical inclination divisors $5s^{2}-4$ start to appear only at the third order truncation.

Table 9: Inclination polynomials

\Psi_{i,j}

in Eq. (45).

_1,0	$-3(5s^{2}-4)^{2}$
_1,1	$-3(3s^{2}-2)(5s^{2}-4)$
_2,0	$\frac{15}{8}(5s^{2}-4)^{2}(77s^{4}-172s^{2}+88)$
_2,1	$\frac{9}{8}(5s^{2}-4)^{2}(155s^{4}-256s^{2}+104)$
_2,2	$\frac{3}{8}(5s^{2}-4)^{2}(189s^{4}-156s^{2}+8)$
_2,3	$\frac{15}{8}(5s^{2}-4)^{2}(5s^{4}+8s^{2}-8)$
_3,0	$-\frac{15}{32}(2439500s^{12}-11312175s^{10}+21772080s^{8}$
	$-22346500s^{6}+12956400s^{4}-4043136s^{2}+533248)$
_3,1	$-\frac{45}{32}(5s^{2}-4)(62300s^{10}-260365s^{8}+431504s^{6}$
	$-356508s^{4}+147552s^{2}-24576)$
_3,2	$\frac{3}{16}(1835625s^{12}-7723875s^{10}+13291500s^{8}$
	$-12015300s^{6}+6064176s^{4}-1644928s^{2}+192256)$
_3,3	$\frac{15}{16}(5s^{2}-4)(18175s^{10}-85105s^{8}+153172s^{6}-136540s^{4}$
	$+61408s^{2}-11264)$
_3,4	$\frac{3}{32}(213750s^{12}-1441125s^{10}+3537000s^{8}-4313100s^{6}$
	$+2835280s^{4}-967808s^{2}+135424)$
_3,5	$\frac{21}{32}s^{2}(5s^{2}-4)(15s^{2}-14)(450s^{6}-925s^{4}+590s^{2}-112)$

Table 10: Inclination polynomials

\omega_{i,j}

in Eq. (46).

_1,0	$-3(5s^{2}-4)^{2}$
_2,0	$\frac{15}{8}(5s^{2}-4)^{2}(77s^{4}-172s^{2}+88)$
_2,1	$9(3s^{2}-2)(5s^{2}-4)^{3}$
_2,2	$\frac{3}{8}(5s^{2}-4)^{2}(45s^{4}+36s^{2}-56)$
_3,0	$-\frac{15}{32}(2439500s^{12}-11312175s^{10}+21772080s^{8}$
	$-22346500s^{6}+12956400s^{4}-4043136s^{2}+533248)$
_3,1	$-\frac{45}{4}(5s^{2}-4)^{3}(168s^{6}-497s^{4}+460s^{2}-136)$
_3,2	$\frac{3}{16}(2150625s^{12}-9409875s^{10}+16968300s^{8}$
	$-16218180s^{6}+8729136s^{4}-2535808s^{2}+315136)$
_3,3	$-\frac{15}{4}(5s^{2}-4)^{3}(105s^{6}+39s^{4}-228s^{2}+104)$
_3,4	$\frac{3}{32}(438750s^{12}-1771125s^{10}+2865000s^{8}-2345100s^{6}$
	$+999760s^{4}-199808s^{2}+12544)$

Table 11: Inclination polynomials

\Omega_{i,j}

in Eq. (47).

_1,0	$-6(5s^{2}-4)$
_2,0	$\frac{15}{2}(5s^{2}-4)^{2}(7s^{2}-8)$
_2,1	$18(3s^{2}-2)(5s^{2}-4)^{2}$
_2,2	$\frac{3}{2}(5s^{2}-4)^{2}(5s^{2}+4)$
_3,0	$-\frac{15}{8}(215250s^{10}-823025s^{8}+1255040s^{6}-953760s^{4}$
	$+361088s^{2}-54464)$
_3,1	$-\frac{45}{4}(5s^{2}-4)^{3}(63s^{4}-124s^{2}+56)$
_3,2	$\frac{3}{8}(430125s^{10}-1553550s^{8}+2222340s^{6}-1570224s^{4}$
	$+546432s^{2}-74624)$
_3,3	$-\frac{15}{4}(5s^{2}-4)^{3}(45s^{4}+28s^{2}-40)$
_3,4	$\frac{3}{8}(50625s^{10}-168375s^{8}+215900s^{6}-130800s^{4}$
	$+35840s^{2}-3136)$

Computing the periodic corrections of the semi-equinoctial variables would require to carry out expansions of the true anomaly. Therefore, we rather formulate them in polar variables. In this way we avoid the trouble of small denominators in the $J_{2}$ -problem, yet, of course, the nodes of exactly equatorial orbits remain undefined. For simplicity we do not deal with this latter case, which, if desired, could be approached using different sets of non-singular variables.

3 Performance of the solution

The performance of the analytical solution is illustrated in three different cases. The first one is a low-altitude, almost-circular orbit with the orbital parameters of the PRISMA mission PerssonJacobssonGill2005 . The second is a high-eccentricity orbit with the typical characteristics of a geostationary transfer orbit (GTO), which we borrowed from Konopliv1991 . The third test has been specifically carried out to check the behavior of the analytical solution when approaching the critical inclination resonance. For this last case we selected an orbit with orbital parameters similar to the TOPEX orbit, which departs only $\sim 3^{\circ}$ from the inclination resonance condition Frauenholzetal1998 . Orbital elements corresponding to these three cases are presented in Table 12. The reference, true orbits have been propagated numerically in extended precision to assure that all the computed points are accurate with 15 significant digits in the decimal representation. We requested that precision because this is the maximum number of digits with which exact decimal operations are guaranteed in standard double-precision Kahan1997 . Because the analytical solution is evaluated in double-precision the reference numeric orbit is considered exact.

Table 12: Initial conditions of the test cases. Angles are in degrees.

type	$a$ (km)	$e$	$I$	$\Omega$	$\omega$	$\ell$
PRISMA	$6878.137$	$0.001$	$97.42$	$168.162$	$20$	$30$
TOPEX	$7707.270$	$0.0001$	$66.04$	$180.001$	$270$	$180$
GTO	$24460.00$	$0.73$	$30.$	$170.1$	$280$	$0$

Two different kind of errors are associated to the truncation order of the analytical solution. On the one hand, the truncation of the secular terms introduces an error in the computation of the frequencies of the analytical solution, Eqs. (45)–(47), which will make the perturbation solution to degrade with time. On the other hand, the truncation affects also the generating function from which the periodic corrections are derived. The latter fundamentally affects the amplitude of the periodic errors, and, therefore, the quality of the ephemeris provided by the analytical solution is not expected to deteriorate significantly with time by effect of this error. However, this is only true for the direct transformation, a case in which we will see that it is acceptable to use a lower order truncation than the truncation used for the secular terms. On the contrary, the truncation order of the inverse corrections has a direct effect in the precision with which the initialization constants are computed from a given set of initial conditions. Therefore, not computing the inverse transformation with the same accuracy as that of the secular terms will also contribute to the deterioration of the latter. In practice, the highest accuracy of the inverse transformation can be limited to the computation of the periodic corrections of the initial semimajor axis (or its partner canonical variable, the Delaunay action) because it affects directly the secular mean motion $n_{\ell}$ , while the other elements only affect the frequencies of the secular node and perigee, which are $\mathcal{O}(J_{2})$ when compared to $n_{\ell}$ .

3.1 Accuracy of the periodic corrections

First of all, we check the accuracy of the periodic corrections for increasing orders of the perturbation solution. If the periodic corrections were exact, transforming different states of the same orbit provided by the reference solution will result in the same, constant secular values of the momenta, and in an exactly linear growing of the secular angles. On the contrary, the secular values obtained from the different states of the true orbit oscillate periodically due to the truncation order of the solution. This is illustrated in Fig. 1 for the semimajor axis of the PRISMA test orbit, where the relative errors with respect to a constant reference value, which has been obtained as the arithmetic mean of the computed secular values, are shown for different truncations of the perturbation solution. We recall that the expected errors of a perturbation solution are of the order of the neglected terms, as follows from Eq. (8). Therefore, for a truncation to the order $m$ , we would expect errors of the order $\varepsilon^{m+1}$ .

Refer to caption — Figure 1: Relative errors of the semimajor axis of the PRISMA-type orbit for, from top to bottom, the 1st, 2nd, 3rd, and 4th order truncations of the perturbation solution. The latter shows the relative errors of the 5th order truncation superimposed. Abscissas are hours.

We found that the relative errors of the secular semimajor axis of the first order truncation (top plot of Fig. 1) are of the order of $J_{2}^{2}$ , as it should be the case, which amounts to 3 meters in absolute value. The relative errors of the second order truncation (second to top plot) are of the order of $J_{2}^{3}$ , or less than 3 millimeter in absolute value. The third and fourth order truncations of the solution yield relative errors of the order of $J_{2}^{4}$ and $J_{2}^{5}$ , respectively (second from bottom and bottom plots, respectively), or absolute errors of micrometers and hundredths of micrometers, respectively. The latter are very close to the 15 exact digits reachable working in double precision. Still, additional improvements are obtained when the truncation of the periodic corrections is extended to the fifth order of $J_{2}$ , now effectively reaching the numerical precision, as shown in the bottom plot of Fig. 1, in which the relative errors of the fifth order truncation (black dots) are superimposed to the previous case.

The behavior is analogous in the case of the GTO orbit, yet now the errors notable peak at perigee passages. This is illustrated in Fig. 2, where, from top to bottom, we present now the relative errors of the inclination $I=\arccos(H/G^{\prime})$ , rather than the total angular momentum $G^{\prime}$ , for the 1st, 2nd, …, 4th order truncations of the periodic corrections, respectively. The relative errors of the 4th order truncation in the bottom plot of Fig. 2 (black dots) are superimposed to the third order ones (gray line) for reference. Corresponding absolute errors are of the order of tens of milliarc seconds for the 1st order, hundredths of mas for the second, hundredths of microarc seconds for the third, and below the level of 1 thousandth of $\mu\mathrm{as}$ for the fourth order. Due to the larger semimajor axis, the GTO orbit is less perturbed in general, and we found that the numerical precision is achieved at the fourth order of the perturbation approach.

In the case of the TOPEX orbit the inclination of $\sim 66^{\circ}$ makes the coefficient $5s^{2}-4$ to become less than 2 tenths. However, this small divisor is not harmful at all in a first order truncation because, as follows from Eq. (20), it is multiplied by the square of the eccentricity, which is really small for the TOPEX orbit (about one thousandth, on average). Hence, the relative errors in the semimajor axis introduced by the periodic corrections are slightly better than those previously found for PRISMA, as shown in the top plot of Fig. 3. The improvements should be a consequence of the larger semimajor axis of the TOPEX orbit when compared with that of PRISMA.

The effects of the offending divisors are more apparent when using higher order truncations of the periodic corrections. Indeed, the coefficients of $\sin(f+2\omega)$ and $\sin(3f+4\omega)$ in Eq. (22), as well as the coefficient of $\sin{f}$ in Eq. (33), are multiplied by the factor $e/(5s^{2}-4)^{2}$ , which, while being bigger than before, it is neither of worry for TOPEX yet, as demonstrated in the second to the top plot in Fig. 3, where it can be checked that the relative errors in the semimajor axis remain of $\mathcal{O}(J_{2}^{3})$ . On the contrary, at the third order we find the small divisor alone in coefficients of $\sin 2m(f+\omega)$ , $m=1,2,3$ , in Eq. (25), in this way lessening to some extent the corrector effect of the terms factored by $J_{2}^{3}$ . This slight deterioration of convergence is noted in the second-from-bottom plot of Fig. 3, where we see that the relative errors in the semimajor axis are now about 10 times bigger than corresponding errors of PRISMA. Bigger offending coefficients are found in higher orders of the periodic corrections, which make that the relative errors of the fourth order truncation are clearly larger than corresponding ones of PRISMA. These relative errors are depicted in the bottom plot of Fig. 3 (gray line) superimposed with the errors of the fifth order truncation (black dots), showing that, at the end, the numerical precision is reachable also in this case with the fifth order truncation of the solution. The absolute error of the different truncations are now of the order of one meter, a few tenths of cm, sevral hundredths of mm, tenths of $\mu\mathrm{m}$ , and thousandths of $\mu\mathrm{m}$ , respectively.

3.2 Accuracy of the analytical solution

The reference orbits of the three test cases are now compared with the ones provided by the analytical solution. That is, starting from initial conditions in Table 12, we first apply the inverse periodic corrections to initialize the constants of the secular solution. Then, the secular terms are evaluated in Eqs. (39)–(44) at the same times $t_{i}$ as those in which the true solution has been obtained from the numerical integration. It follows the application of the direct periodic corrections to each secular term to get the ephemeris predicted by the analytical solution at the times $t_{i}$ . The accuracy of the different truncations of the analytical solution is assessed by computing the root-sum-square (RSS) of the difference between the position and velocity predicted by the analytical solution and those of the reference orbit and the times $t_{i}$ .

The results of the PRISMA test case are shown in Fig. 4, in which ordinates are displayed in a logarithmic scale to ease comparisons. The letter $S$ in the labels ( $S$ : $P$ ) stand for the truncation order of both the inverse corrections and the secular terms, whereas $P$ indicates the truncation order of the direct periodic corrections. As we already pointed out, the latter do not need to be computed to the same accuracy as the inverse transformation. We observe in Fig. 4 that a simple first order truncation of the analytical solutions —curve labeled (1:1)— starts with a RSS error of approximately 1 meter in position, but, due to the errors introduced by the first order truncation of the secular terms, the RSS errors reach more than 10 km after one month. The solution is notably improved when taking the second order terms of the inverse corrections and the secular terms into account. This is shown with the curve labeled (2:1), that only reaches about 30 m at the end of the one-month propagation. The improvements are obtained with only a slight increase of the computational burden, because the inverse corrections and the secular frequencies are evaluated only once. The (3:2) propagation starts from a much smaller RSS error, of less than 1 cm, that only grows to about 10 cm by the end of day 30. Errors fall clearly below the mm level in the case of the (4:3) truncation, and to just a few $\mu\mathrm{m}$ in the (5:4) case. No further improvement of the RSS errors is observed if the direct periodic corrections are taken also up to the fifth order, yet a slight improvement is achieved in that case in the preservation of the energy integral, that reaches in this last case the numerical precision.

The behavior of the analytical solution is analogous in the case of the high-eccentricity GTO orbit in what respects to the secular terms, yet the errors of the periodic corrections are now of larger amplitude, as clearly observed in Fig. 5. Now, the periodic oscillations of the (1:1) truncation may allow the RSS errors to grow to many tens of km. A secular trend in the RSS errors of about half meter per day of the (2:1) solution remains almost hidden along the 30 days propagation under periodic oscillations of about 30 m. A similar behavior is observed in the (3:2) propagation, where the amplitude of the oscillations is now reduced to the cm order. This amplitude is further reduced to hundredths of mm with the (4:3) truncation. Finally the (5:4) solution improves slightly the propagation, and the RSS errors remain of just a few $\mu\mathrm{m}$ along the whole propagation interval, thus showing the same quality as in the previous test case.

Analogous tests carried out for TOPEX are presented in Fig. 6. The results are similar to those previously presented in Fig. 4 for PRISMA, and the apparent discrepancies when using the (2:1) solution are only due to the logarithmic scale, which encompasses a different range in each figure. In fact, the periodic oscillations of the RSS errors are of the same amplitude in both cases ( $\sim 2$ m), yet the secular trend grows at a low rate of less than 2 cm per day in the case of TOPEX while it does almost two orders of magnitude faster for PRISMA ( $\sim 1$ m/day). This fact should be attributed to the less perturbed orbit of TOPEX due to its higher altitude. The (3:2) solution provides quite similar RSS errors in both cases, TOPEX and PRISMA, whereas the (4:3) solution performs worse for TOPEX, whose RSS errors grow now at a secular rate about 5 times faster than in the case of PRISMA. These behavior is in agreement with the slight deterioration of the solution previously observed when testing the periodic corrections due to the proximity to the critical inclination. Things are balanced with the (5:4) order solution —also in agreement with the previously observed behavior due to the increased precision of the analytical solution— for which the RSS errors reach only the micrometer level at the end of the 30 day propagation.

As it is confirmed by the examples carried out, at the end of a long-term propagation the secular trend of the errors prevails over their periodic oscillations. Therefore, in practice, perturbation solutions in which the periodic terms are recovered to a lower order than the order at which the secular frequencies are truncated make full sense. The computational burden of these kinds of solutions is notably alleviated, and hence are definitely much more practicable for ephemeris computation. This is further illustrated for TOPEX in Fig. 7, where it is shown that the (5:3) truncation of the analytical solution might replace the more accurate (5:4) truncation in a long-term propagation with substantial reductions in computing time and minimum degradation of accuracy, which, besides, remains much more uniform along the whole propagation interval. Note that the scale of the ordinates axis has been changed from meters in Fig. 6 to mm in Fig. 7.

4 Conclusions

The main problem of artificial satellite theory is a two degrees of freedom model that captures the bulk of the gravitational effects undergone by low Earth orbits. Still, it lacks of enough integrals to obtain a closed form solution. On the other hand, perturbation methods disclose as specially successful in providing analytical solutions that can replace the true dynamics within a high degree of approximation. Using them, it is shown that there are wide regions of phase space in which the main problem behaves as if it were separable.

In particular, the main problem Hamiltonian has been completely reduced by reverse normalization. That is, the motion of the orbital plane is first decoupled from the motion on that plane, and then the remaining short-period terms are removed by averaging. The normalization is achieved without need of resorting to Hamiltonian simplification procedures, in this way radically departing from a well stablished tradition in artificial satellite theory. The computed analytical solution is exact to all practical purposes in the sense that, working in double-precision floating point aritmethic, there are no substantial differences with respect to exact numerical integrations (computed with extended precision) for reasonably long time intervals. In agreement with the particular value of the Earth’s oblateness coefficient, this precision is achieved only when the perturbation approach is extended up to, at least, the fifth order of the small parameter.

Acknowledgements.

Support by the Spanish State Research Agency and the European Regional Development Fund under Projects ESP2016 -76585-R and ESP2017-87271-P (MINECO/ AEI/ERDF, EU) is recognized.

Conflict of Interest

: The author declares that he has no conflict of interest.

References

(1) K. Aksnes. A note on ‘The main problem of satellite theory for small eccentricities, by A. Deprit and A. Rom, 1970’. Celestial Mechanics, 4(1):119–121, September 1971.
(2) K. T. Alfriend and S. L. Coffey. Elimination of the perigee in the satellite problem. Celestial Mechanics, 32(2):163–172, February 1984.
(3) V. I. Arnol’d. Mathematical Methods of Classical Mechanics, volume 60 of Graduate Texts in Mathematics. Springer-Verlag, New York, 2nd edition, 1989.
(4) J. V. Breakwell and J. Vagners. On Error Bounds and Initialization in Satellite Orbit Theories. Celestial Mechanics, 2:253–264, June 1970.
(5) R. A. Broucke. Numerical integration of periodic orbits in the main problem of artificial satellite theory. Celestial Mechanics and Dynamical Astronomy, 58(2):99–123, February 1994.
(6) D. Brouwer. Solution of the problem of artificial satellite theory without drag. The Astronomical Journal, 64:378–397, November 1959.
(7) D. Brouwer and G. M. Clemence. Methods of Celestial Mechanics. Academic Press, New York and London, 1961.
(8) A. Celletti and P. Negrini. Non-integrability of the problem of motion around an oblate planet. Celestial Mechanics and Dynamical Astronomy, 61(3):253–260, March 1995.
(9) S. Coffey and A. Deprit. Third-Order Solution to the Main Problem in Satellite Theory. Journal of Guidance, Control and Dynamics, 5(4):366–371, 1982.
(10) S. L. Coffey, A. Deprit, and B. R. Miller. The critical inclination in artificial satellite theory. Celestial Mechanics, 39(4):365–406, December 1986.
(11) R. A. Cook. The long-term behavior of near-circular orbits in a zonal gravity field (AAS 91-463). In B. Kaufman, K. T. Alfriend, R. L. Roehrich, and R. R. Dasenbrock, editors, Astrodynamics 1991, volume 76 of Advances in the Astronautical Sciences, pages 2205–2221, P.O. Box 28130, San Diego, California 92198, USA, 1992. American Astronautical Society, Univelt, Inc.
(12) J. M. A. Danby. Motion of a Satellite of a Very Oblate Planet. The Astronomical Journal, 73(10):1031–1038, December 1968.
(13) J. M. A. Danby. Fundamentals of Celestial Mechanics. Willmann-Bell, Richmond VA, 2nd edition, 1992.
(14) C. E. Delaunay. La Théorie du Mouvement de la Lune, Premier volume, volume 28 of Mémoires de l’Academie des Sciences de l’Institut Impérial de France. Mallet-Bachellier, Paris, 1860.
(15) A. Deprit. Canonical transformations depending on a small parameter. Celestial Mechanics, 1(1):12–30, 1969.
(16) A. Deprit. The elimination of the parallax in satellite theory. Celestial Mechanics, 24(2):111–153, 1981.
(17) A. Deprit. Delaunay normalisations. Celestial Mechanics, 26:9–21, January 1982.
(18) A. Deprit and B. Miller. Simplify or Perish. Celestial Mechanics, 45:189–200, 1989.
(19) A. Deprit and A. Rom. The Main Problem of Artificial Satellite Theory for Small and Moderate Eccentricities. Celestial Mechanics, 2(2):166–206, June 1970.
(20) R. B. Frauenholz, R. S. Bhat, B. E. Shapiro, and R. K. Leavitt. Analysis of the TOPEX/Poseidon Operational Orbit: Observed Variations and Why. Journal of Spacecraft and Rockets, 35:212–224, March 1998.
(21) G. Gaias, C. Colombo, and M. Lara. Analytical Framework for Precise Relative Motion in Low Earth Orbits. Journal of Guidance Control Dynamics, on line:1–13, March 2020.
(22) H. Goldstein, C. P. Poole, and J. L. Safko. Classical Mechanics. Addison-Wesley, 3rd edition, 2001.
(23) P. A. Hansen. Expansions of the product of a power of the radius vector with the sinus or cosinus of a multiple of the true anomaly in terms of series containing the sinuses or cosinuses of the multiples of the true, eccentric or mean anomaly. Abhandlungen der Koniglich Sachsischen Gesellschaft der Wissenschaften, 2(3):183–281, 1855. English translation by J. C. Van der Ha, ESA/ESOC, Darmstadt, Germany, 1977.
(24) D. Hautesserres and M. Lara. Intermediary LEO propagation including higher order zonal harmonics. Celestial Mechanics and Dynamical Astronomy, 127:505–526, April 2017.
(25) L. M. Healy. The Main Problem in Satellite Theory Revisited. Celestial Mechanics and Dynamical Astronomy, 76(2):79–120, 2000.
(26) M. Irigoyen and C. Simó. Non integrability of theJ ₂ problem. Celestial Mechanics and Dynamical Astronomy, 55(3):281–287, March 1993.
(27) I. G. Izsak. A Second-Order Solution of Vinti’s Dynamical Problem. Smithsonian Contributions to Astrophysics, 6:81–107, 1963.
(28) W. H. Jefferys. Automated, Closed Form Integration of Formulas in Elliptic Motion. Celestial Mechanics, 3:390–394, September 1971.
(29) A. H. Jupp. The critical inclination problem - 30 years of progress. Celestial Mechanics, 43(1-4):127–138, 1988.
(30) W. Kahan. ”Lecture Notes on the Status of IEEE Standard 754 for Binary Floating-Point Arithmetic. Technical report, Electrical Engineering and Computer Science, University of California, Berkeley, CA, October 1997.
(31) W. M. Kaula. Analysis of Gravitational and Geometric Aspects of Geodetic Utilization of Satellites. Geophysical Journal, 5:104–133, July 1961.
(32) H. Kinoshita. First-Order Perturbations of the Two Finite Body Problem. Publications of the Astronomical Society of Japan, 24:423–457, 1972.
(33) H. Kinoshita. Third-Order Solution of an Artificial-Satellite Theory. SAO Special Report, 379, July 1977.
(34) A. Konopliv. A Perturbation Method and Some Applications. Celestial Mechanics and Dynamical Astronomy, 47:305, 1990.
(35) A. S. Konopliv. A Third Order of $J_{2}$ Solution with a Transformed Time”. Interoffice Memorandum IOM 314.3 - 970, Jet Propulsion Laboratory, 4800 Oak Grove Dr, Pasadena, CA 91109, USA, March 1991.
(36) Y. Kozai. The motion of a close earth satellite. The Astronomical Journal, 64:367–377, November 1959.
(37) Y. Kozai. Second-order solution of artificial satellite theory without air drag. The Astronomical Journal, 67:446–461, September 1962.
(38) M. Lara. SADSaM: A Software Assistant for Designing SAtellite Missions. Technical Report DTS/MPI/MS/MN/99-053, Centre National d’Études Spatiales, 18, avenue Edouard Belin - 31401 Toulouse Cedex 9, France, May 1999.
(39) M. Lara. Searching for Repeating Ground Track Orbits: A Systematic Approach. Journal of the Astronautical Sciences, 47(3-4):177–188, 1999.
(40) M. Lara. Simplified Equations for Computing Science Orbits Around Planetary Satellites. Journal of Guidance Control Dynamics, 31(1):172–181, January 2008.
(41) M. Lara. Analytical and Semianalytical Propagation of Space Orbits: The Role of Polar-Nodal Variables. In G. Gómez and J. J. Masdemont, editors, Astrodynamics Network AstroNet-II: The Final Conference, volume 44 of Astrophysics and Space Science Proceedings, pages 151–166. Springer, Cham, 2016.
(42) M. Lara. Note on the ideal frame formulation. Celestial Mechanics and Dynamical Astronomy, 129:137–151, September 2017.
(43) M. Lara and S. Ferrer. Closed form perturbation solution of a fast rotating triaxial satellite under gravity-gradient torque. Cosmic Research, 51(4):289–303, July 2013.
(44) M. Lara and S. Ferrer. Expanding the simple pendulum’s rotation solution in action-angle variables. European Journal of Physics, 36(5):055040, September 2015. Expanded version: arXiv eprint nlin.SI 1503.03358.
(45) M. Lara, J. F. San-Juan, and L. M. López-Ochoa. Delaunay variables approach to the elimination of the perigee in Artificial Satellite Theory. Celestial Mechanics and Dynamical Astronomy, 120(1):39–56, September 2014.
(46) M. Lara, J. F. San-Juan, and L. M. López-Ochoa. Proper Averaging Via Parallax Elimination (AAS 13-722). In S. B. Broschart, J. D. Turner, K. C. Howell, and F. R. Hoots, editors, Astrodynamics 2013, volume 150 of Advances in the Astronautical Sciences, pages 315–331, P.O. Box 28130, San Diego, California 92198, USA, January 2014. American Astronautical Society, Univelt, Inc.
(47) M. Lara, R. Vilhena de Moraes, D. M. Sanchez, and A. F. B. A. Prado. Efficient computation of short-period analytical corrections due to third-body effects (AAS 15-295). In R. Furfaro, S. Casotto, A. Trask, and S. Zimmer, editors, AAS/AIAA Spaceflight Mechanics Meeting 2015, volume 155 of Advances in the Astronautical Sciences, pages 437–455, P.O. Box 28130, San Diego, California 92198, USA, 2015. American Astronautical Society, Univelt, Inc.
(48) M. Lara. On inclination resonances in Artificial Satellite Theory. Acta Astronautica, 110:239–246, May 2015. Dynamics and Control of Space Systems.
(49) M. Lara. Nonlinear librations of distant retrograde orbits: a perturbative approach—the Hill problem case. Nonlinear Dynamics, 93(4):2019–2038, apr 2018.
(50) M. Lara. A new radial, natural, higher-order intermediary of the main problem four decades after the elimination of the parallax. Celestial Mechanics and Dynamical Astronomy, 131(9):20, September 2019.
(51) R. H. Lyddane. Small eccentricities or inclinations in the Brouwer theory of the artificial satellite. Astronomical Journal, 68(8):555–558, October 1963.
(52) G. Metris. Mean values of particular functions in the elliptic motion. Celestial Mechanics and Dynamical Astronomy, 52:79–84, March 1991.
(53) A. H. Nayfeh. Perturbation Methods. Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2004.
(54) C. Osácar and J. F. Palacián. Decomposition of functions for elliptical orbits. Celestial Mechanics and Dynamical Astronomy, 60(2):207–223, October 1994.
(55) S. Persson, B. Jacobsson, and E. Gill. PRISMA – Demonstration Mission for Advanced Rendezvous and Formation Flying Technologies and Sensors (paper IAC-05-B56B07). In Proceedings of the 56th International Astronautical Congress (IAC), October 17 - 21 2005, Fukuoka, Japan, pages 1–10. International Astronautical Federation (IAF), International Astronautical Federation (IAF), October 2005.
(56) J. F. San-Juan, D. Ortigosa, L. M. López-Ochoa, and R. López. Deprit’s Elimination of the Parallax Revisited. Journal of the Astronautical Sciences, 60:137–148, June 2013.
(57) F. Tisserand. Traité de mécanique céleste. Tome I: Perturbations des planètes d’aprés la méthode de la variation des constantes arbitraries. Gauthier-Villars et fils, Quai des Grands-Augustins, 55, Paris, 1889.