Efficient Formulation of the Periodic Corrections in Brouwer’s Gravity Solution

Abstract

The periodic terms of Brouwer’s gravity solution are reconstructed in a nonsingular set of variables which are derived from the well-known polar-nodal variables. This change does not affect the essence of the solution, which still keeps all the benefits of the action-angle variables approach and yields two major improvements. Namely, the periodic corrections of Brouwer’s solution are now valid for any eccentricity below one and any inclination except the critical inclination and, besides, are significantly simpler than the nonsingular corrections in Lyddane’s reformulation of Brouwer’s theory.

1. Introduction

Concern in space situational awareness by an increasing number of satellite operators and in particular the necessity of timely scheduling collision avoidance maneuvers motivates current interest in improving the capabilities of orbit prediction programs.

Satellite short-term prediction is customarily carried out with SGP4 [1], an analytical solution that has its roots in Brouwer’s celebrated gravity solution to the artificial satellite problem [2] and which is optimized for the propagation of satellite ephemeris using the element sets in the two-line format specified by the North American Aerospace Defense Command [3, 4]. However, it has been claimed that SGP4 may lack sufficient capabilities for conjunction analysis tasks [5–7]. Besides, terms that may be missing in SGP4 could be responsible for detected noteworthy along-track errors in the SGP4 predictions for GPS satellites [8, 9]. The known limitations of predicted ephemeris from two-line elements motivate the current development of new algorithms as, for instance, those in the software STELA of the Centre National d’Etudes Spatiales [10].

Brouwer found his solution using a perturbation approach, the so-called von Zeipel method [11, 12], which splits the satellite motion into secular terms, long-period corrections, related to the evolution of the argument of the perigee, and short-period corrections, related to the satellite’s mean motion. For the secular terms, Brouwer’s theory includes gravitational effects up to the second order of , the second degree zonal harmonic coefficient of the spherical harmonics expansion of the geopotential, which for the Earth is of the order of one thousandth. But in the case of periodic corrections the theory is limited to first order effects of . Therefore, the short-period corrections are only related to the contribution of , whereas the long-period corrections of Brouwer’s gravitational solution include first order corrections due to the few first zonal harmonics [2].

Brouwer developed his original theory in Delaunay variables, the canonical counterpart of classical Keplerian orbital elements, which, like them, are singular for circular orbits and equatorial orbits. This fact may cause troubles in the computation of the periodic corrections for both low-eccentricity and low-inclination orbits, but the problem is easily solved by reformulating Brouwer’s gravitational solution in nonsingular variables as, for instance, Poincaré’s canonical variables [13]. However, the periodic corrections either when formulated in Delaunay variables or in Poincaré variables are made of long trigonometric series, a fact that very soon motivated efforts in improving their evaluation [14]. Among the different efforts in improving the evaluation of Brouwer’s gravitational theory, the use of polar-nodal variables for computing the periodic corrections was advocated by different authors [15–17]. These variables, also called Hill or Whittaker variables, are valid for orbits with any eccentricity below 1, but still they are singular in the case of equatorial orbits, a fact that may cause trouble in the evaluation of the periodic corrections of almost equatorial orbits.

It is worthy to remember that Brouwer’s gravitational theory breaks at the critical inclination of 63.4 degrees. Indeed, since the secular terms in Brouwer’s perturbation approach are computed by averaging periodic effects, resonant inclination orbits and, in particular, the critical inclination are excluded from the field of applicability of Brouwer’s solution (see [18] and references therein).

The efficient evaluation of the periodic corrections is even more critical when taking into account second order corrections, which notably improve the performance of the perturbation theory [19] but for which the trigonometric series are significantly longer [15], and hence the advantages of using polar-nodal variables are more evident. Besides, the benefits of formulating the periodic corrections in polar-nodal variables are not limited to the case of geopotential perturbations, and this set of canonical variables has been revealed very efficient in the evaluation of periodic corrections due to third-body perturbations [20].

To avoid the troubles related to the evaluation of the long-period corrections of low-inclination orbits, Aksnes’ suggestion of computing the corrections for the satellite’s latitude and (true) longitude for these orbits [17] is followed in the present research. Indeed, without limiting to the case of low-inclination orbits, the periodic corrections are rewritten in a set of noncanonical variables which are directly constructed from the polar-nodal ones. These variables are nonsingular and provide a more efficient evaluation of the periodic corrections than the corresponding corrections in Poincaré variables which are used in Lyddane’s modifications to Brouwer’s gravitational solution.

The paper is organized as follows. For completeness, the construction of short-period corrections in polar-nodal variables of Brouwer’s gravitational solution is recalled in Section 2, while the construction of long-period corrections in polar-nodal variables is illustrated in Section 3. Then, the new set of nonsingular, noncanonical variables is introduced in Section 4, and the long-period and short-period corrections are reformulated in the nonsingular variables. The transformations of the nonsingular elements from and to Cartesian variables are free from singularities and are also documented in Section 4.

2. Short-Period Elimination

Since this research deals only with perturbations of gravitational origin, the problem of disturbed Keplerian motion can take benefit from Hamiltonian formulation. Thus, the motion of a massless particle in the gravitational field of the Earth is derived from the Hamiltonian where represents the integrable Keplerian Hamiltonian and is the disturbing function, which comprises the noncentralities of the geopotential. From the usual solution of Laplace’s equation in spherical coordinates, the forces model is further limited to the zonal harmonics case, in which the disturbing function is written as follows: where is the Earth’s gravitational parameter, is the Earths’ equatorial radius, is the radial distance from the Earth’s center of mass, is latitude, are Legendre polynomials of degree , and are corresponding zonal harmonic coefficients.

The problem of small inclinations in Brouwer’s solution is related to the effects of odd zonal harmonics, so to illustrate this case it is enough to consider the impact of , in addition to the main problem, and hence the zonal gravitational potential in (2) is further truncated to the degree . Besides, because of the different orders of the harmonic coefficients, where , it is found convenient to make the Hamiltonian perturbative arrangement in which where the relation has been used, with being the orbital inclination, the argument of the perigee, and the true anomaly; and are abbreviations for the sine and cosine of the inclination, respectively, is the semimajor axis, and, in consequence with the Hamiltonian formulation, all the symbols, that is to say, , , , , and , are assumed to be functions of some set of canonical variables.

In particular, Brouwer finds a transformation from “old” to “new” (or primes) variables, such that the Hamiltonian in the new variables only depends on momenta, whereas the angles have been averaged. Therefore, it relies on the action-angle variables of the Kepler problem, the so-called Delaunay variables, namely, the mean anomaly and its conjugate momentum (the Delaunay action), the argument of the perigee and its conjugate momentum (the total angular momentum), where is the orbital eccentricity, and the argument of the node and its conjugate momentum (the polar component of the angular momentum). By using this canonical set it is simple to see that is cyclic in (3) and, therefore, is an integral of the zonal problem.

The Hamiltonian reduction of (3) by perturbation methods is thoroughly documented in the literature, and hence results are provided without giving details on the method. In particular, the computations carried out were based on the implementation of the Lie transforms method known as Deprit’s triangle algorithm, which is nowadays considered standard for Hamiltonian perturbations. Readers interested in the Lie transforms method can find all the required details in the original papers of Hori [21] and Deprit [22], as well as in modern celestial mechanics textbooks like [23, 24] or other specialized books as [12]. Note that, following tradition, the notation in prime variables is avoided in what follows when there is no risk of confusion.

At the first order of the perturbation approach, Deprit’s triangle gives where notes the Poisson bracket of two functions and of the canonical variables, which in this case are the Delaunay variables. In order to obtain Brouwer’s solution, the new Hamiltonian term is chosen as the averaging of over the mean anomaly where is the eccentricity function and, for the sake of abbreviating notation, the function has been introduced, which is given by where is the semilatus rectum.

The corresponding term of the generating function is computed from (5) by quadrature where is the mean motion. In view of the differential relation , (9) can be integrated in closed form of the eccentricity to give where is the equation of the center (cf. Equation of [2], keeping in mind the different sign convention in the Hamilton equations).

Up to the first order, the transformation equations of the averaging are computed from where, here, and .

Corresponding transformation equations in Delaunay variables can be expressed as Fourier series which involve sine and cosine functions of 10 different arguments of the form with and (cf. the first order terms in Equations and of [15], for instance). However, important simplifications can be achieved by using the function instead of wholly expanding the transformation equations as Fourier series. In this way, the number of trigonometric arguments is reduced to just four: , , , and (cf. Equations and of [2]).

Alternatively, as pointed out by Izsak [16], the generating function can be expressed in the canonical set of polar-nodal variables , which stand for the radial distance, the argument of latitude, the argument of the node, the radial velocity, the total angular momentum, and the polar component of the angular momentum. Rewriting (10) in polar-nodal variables as is straightforward, leading to where the functions are the projections of the eccentricity vector in the orbital frame when written in polar-nodal variables, which are trivially derived from (12) and from its time derivative

Note that (13) differs from Equation of [16] in the sign; however, both equations are equivalent because of the different sign convention used in the derivation of Hamilton equations.

The first order transformation equations of the short-period averaging in polar-nodal variables are obtained, again, from (11), where, now, and . In this case, the equation of the center is , and, in particular, the partial derivatives are needed in the computation of the Poisson brackets.

Hence, it is easily obtained that which must be evaluated in prime variables for direct corrections and in original variables for inverse corrections . Remarkably, now the evaluation of the corrections only requires dealing with sine and cosine functions of the single argument . Note that the evaluation of the equation of the center is required in (18) and (19). It is done using Kepler equation , where , and is unambiguously computed from and .

The second order of Deprit’s triangle gives and the new Hamiltonian term is chosen as the average of plus the computable Poisson brackets to give Second order corrections to the orbital elements are of the order of the square of and are normally omitted. Therefore, there is no need for solving from (24) and the short-period transformation limits to the first order corrections in (17)–(22), whose simple inspection shows that they are free from singularities either for equatorial or circular orbits.

3. Long-Period Elimination

After the mean anomaly has been averaged, the long-period Hamiltonian is where , , , which are expressed in prime elements, although primes have been dropped for alleviating notation.

In the new notation, the first order of Deprit’s triangle in (5) is rewritten as Because in (26) does not depend on , the new first order Hamiltonian term is chosen, , and hence from (27). However, this does not mean to make null the first order term of the long-period generating function. Quite on the contrary, since the generating function of the long-period averaging does not depend on , then necessarily vanishes in (27). Therefore, the term can only be determined at the next order of the perturbation algorithm.

The Poisson bracket vanishes likewise, and the second order of Deprit’s triangle in (24) is simplified in this case to where the term is chosen as the average of over the argument of the perigee. That is,

It follows that the computation of from (28) by a quadrature is which is trivially solved to give where, for conciseness, the notation has been introduced.

Next, the long-period generating function is obtained by rewriting (31) in polar-nodal variables, and the first order transformation equations in polar-nodal variables are obtained again from (11), with and where now is replaced by . In this way the long-period corrections have been obtained, where the inclination polynomials , , are given in . Equations (34)–(39) must be evaluated in second prime variables for direct corrections and in prime variables from the inverse corrections .

Note that the term in denominators of (34)–(38) prevents application of the long-period corrections to orbits with the critical inclination of 63.4 degrees. This singularity is not related to the variables used and simply reflects the fact that inclination resonances are out of the range of applicability of Brouwer’s gravitational solution (see [18] and references therein).

Finally, it is worthy to mention that after computing the double-prime Delaunay variables from the secular terms, the Kepler equation must be solved to find first and then , in order to compute corresponding double-prime polar-nodal variables.

Inclination Polynomials. Consider

4. The Case of Low Inclinations

Due to the contribution of the odd zonal harmonic , it happens that and are singular for equatorial orbits. However, as the simple inspection of (35) and (36) may suggest, the trouble in the case of low inclinations is easily remedied by computing the long-period corrections to the nonsingular, noncanonical elements , where

Indeed, by simple differentiation, where neither nor are affected by singularities, as easily checked in (38) and (35), respectively. Alternatively, because (11) applies to any function of the canonical variables [22, 25], the corrections in (41) can be computed from corresponding Poisson brackets, namely, Straightforward manipulations lead to the explicit equations in the nonsingular variables where the coefficients are given in , with taken from and    from (40),   , and and are given in (14). Note the almost symmetric form of the corrections and in (44) and (45), respectively.

Coefficients in (44)-(45). Consider

In the case of the Earth, , and hence terms of the order of and higher can be neglected for the lower inclination orbits, because they only produce higher order effects. Therefore, the corrections are straightforwardly derived from (34)–(39). Note that (49)–(51) have been previously provided by Aksnes [17].

4.1. Transformation from Cartesian Variables

The direct transformation from nonsingular to Cartesian variables is obtained by means of the usual rotations applied to the projections of the position and velocity vectors in the orbital frame. Thus, where , are the usual rotation matrices:

After replacing and , , in (55), the transformation from nonsingular to Cartesian variables can be obtained from the sequence where and . Remark that is an integral of the zonal problem and, therefore, its value is always known from given initial conditions.

The inverse transformation, from Cartesian to nonsingular variables, is obtained from the sequence where the computation of , which is unambiguously determined from (70), requires the previous computation of and from (63).

Note that (63) are singular for equatorial retrograde orbits, a case in which . However, this drawback is easily remedied, and the case of almost equatorial, retrograde orbits is effectively addressed by using the variable instead of . Then, the corrections in (49)–(54) still apply, yet the conversion from nonsingular to Cartesian coordinates is slightly modified. Indeed, and in (58) and (61) must be replaced by and , respectively, whereas changing in (63) by allows for computing and from this equation in both cases of direct and retrograde inclinations.

4.2. Short-Period Corrections in Nonsingular Variables

In spite of the fact that there is no trouble in the evaluation of the short-period corrections in the case of low-inclination orbits, it may be convenient to compute (17)–(22) also in nonsingular variables. In this case, and, except for higher order effects, the short-period corrections to the lower inclination orbits may be written in nonsingular elements as

Working in real arithmetic, a state , which would correspond to an exactly (instantaneous) equatorial orbit, would be rarely obtained. If this case occurs either in the original or double prime space, then . In consequence, it makes no sense to speak of the node or the argument of latitude. However, periodic corrections still exist for and . Indeed, while short-period corrections and vanish for equatorial orbits, as derived from (73)-(74), corresponding long-period corrections do not and (50) and (51) result in

5. Conclusions

Soon after Brouwer’s solution was announced, the reformulation in polar-nodal variables of both the short-period and long-period corrections was suggested as a way of simplifying their evaluation. Indeed, as odd as it may seem to introduce short-period terms in the computation of long-period corrections, this artifact prevents the usual deterioration of the corrections in the case of low-eccentricity orbits, yet the case of low-inclination orbits must be treated separately. However, the elementary inspection of the long-period corrections in polar-nodal variables reveals a simple set of (noncanonical) elements that may be used for dealing properly with that case. The new formalism which is nonsingular yields significantly less computational effort than Lyddane’s nonsingular variables approach and can be extended to reformulate third-body periodic corrections in a compact form. The latter is under development and will be published elsewhere.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

Partial support from Projects ESP2013-41634-P and ESP2014-57071-R of the Ministry of Economic Affairs and Competitiveness of Spain is recognized.