{'fraction_non_alphanumeric': 0.031286210892236384, 'fraction_numerical': 0.01086826227674112, 'mean_word_length': 4.8806390977443606, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Note on the analytical integration of circumterrestrial orbits

The acclaimed merits of analytical solutions based on a fictitious time developed in the 1970's were partially overvalued due to a common misuse of classical analytical solutions based on the physical time that were taken as reference. With the main problem of the artificial satellite theory as a model, we carry out a more objective comparison of both kinds of theories. We find that the proper initialization of classical solutions notably balances the performance of the two distinct approaches in what respects to accuracy. Besides, extension of both kinds of satellite theories to higher orders show additional pros and cons of each different perturbation approach, thus providing complementary information to prospective users on which kind of analytical solution may better support their needs.


Introduction
The advent of the perturbation theories of Kozai [21] and Brouwer [7] made a breakthrough in the analytical computation of circumterrestrial orbits. Both theories were computed in closed from of the eccentricity, thus extending their range of application from the lower eccentricities to highly elliptic orbits, and successfully included the effects of higher order zonal harmonics of the geopotential as well as the second order effects of the dominant, zonal harmonic of the second degree -customarily denoted J 2 . Brouwer and Kozai analytical solutions are given in the form of collections of truncated power series of J 2 with time as the argument. Several of the series represent the secular terms of the solution, which are expected to provide the average orbit evolution, and the remaining series comprise the periodic corrections that are needed for computing ephemeris. In the original proposals of Kozai and Brouwer, the former were accurate up to the second order of J 2 effects, whereas the later were truncated to O(J 2 ).
A crucial point in the use of these kinds of truncated power series is the initialization of the constants of the theory [8,43,42,16,18]. In particular, the accurate determination of the "mean" mean motion is mandatory to obtain the in-track error that would be expected from the truncation of the secular terms. This point was clearly highlighted by developers of perturbation solutions (see remarks in [22]). However, it created some confusion among users, who reported unexpected along-track errors in the propagation of initial conditions when Brouwer's theory was compared with alternative available solutions (see [4] and references therein). The trouble resulted from Brouwer's computations of the periodic correctionswhich are needed for the initialization of the constants of the theory from given initial conditions-to a lower order than the one of the secular terms. On the other hand, this inaccuracy did not happen when the value of the "mean" mean motion was known to the required accuracy by other means, like the usual fit to observations. Therefore, the problem was understood like an inconsistency in the use of Brouwer's solution.
Increasing the truncation order of the periodic corrections to solve the inconsistency involves heavy additional analytical computations [22,15], which commonly discouraged users from their implementation. However, since the increased accuracy was fundamentally needed for the osculating to mean conversion of the semimajor axis, simple ways of patching Brouwer's solution were soon devised [32]. In this regard, the alternative proposed by Breakwell and Vagners [6] is remarkable, because it does not require additional computations to those already carried out by Brouwer. Indeed, the smart use of the energy equation to calibrate the mean semimajor axis provided a computationally inexpensive way of reducing in-track errors to a comparable level of those in the radial and cross-track directions (see [24], for instance).
Conversely, because the accurate initialization of the mean semimajor axis is clearly associated with the correct treatment of the total energy, the pronounced loss of accuracy in the in-track direction does not happen to perturbation solutions formulated in the extended phase space, in which the total energy is one more integration variable [37,38,41]. In this scenario, the analytical solution is advantageously computed in a fictitious-time scale, whereas the physical time is incorporated into the perturbation scheme like the conjugate coordinate of the total energy. In particular, the annoying complications introduced by the equation of the center when dealing with geopotential perturbations in closed form [20,12,33,1,28] are now totally avoided.
The merits of the extended phase space formulation in the implementation of an analytical orbit generator are unquestionable [39]. However, assertions like the accuracy of first order analytical solutions based on this approach compare with traditional second order perturbation solutions based on the physical time [40] are true only when the latter is incorrectly initialized, and hence must be somewhat downgraded. To show that we take the main problem of the artificial satellite theory as a demonstration model, and make an independent implementation of both kinds of theories. Needless to say that our implementation of the traditional one, based on the physical time, incorporates the calibration of the mean semimajor axis in the initialization part of the orbit theory. While we still found higher accuracy of the secular terms in the case of the fictitious-time approach, in no way it reaches the accuracy of a second-order traditional (correctly initialized) perturbation solution. In addition, we show that, due to the unavoidable truncation of the time equation, additional errors arise in the fictitious-time approach that are of comparable magnitude to the obtained position errors. Moreover, the need of inverting the time equation, or resorting to iterative procedures for finding the fictitious time corresponding to a given physical time [5], which is needed for usual ephemeris generation, must also be taken as a shortcoming of the fictitious-time approach because it notably increases the computational burden of usual ephemeris evaluation.
To further assess the relative merits of each perturbation approach, we extended both kinds of satellite theories to the second order that is customarily used in accurate analytic ephemeris prediction programs [9,10]. In the case of the extended phase space formulation, spurious higher-order terms are generated at each order of the perturbation algorithm, a fact that makes the series defining the analytical solution notably longer than those resulting from the traditional, physical-time approach. Therefore, a post-processing is mandatory to keep the length of the series comprising the perturbation solution based on the fictitious time of comparable size to those of the traditional solution based on the physical time.
2 Solution to the J 2 -problem in the extended phase space The topic of perturbation solutions of the artificial satellite problem using the traditional approach is satisfactorily covered in the literature. In particular, for the main problem, explicit expressions of higher-order solutions have been provided in different places (see [19,23,29], for instance). Therefore, we focus on the less studied formulation of perturbation solutions in the extended phase space, which, as customary nowadays, we recompute using the method of Lie transforms [11] rather than the von Zeipel perturbation algorithm used in the original developments [40].
We approach the perturbation solution in the true anomaly Delaunay-similar elements [37,38,17]. These elements are the action-angle variables of the unperturbed problem -the Kepler problem in the extended phase space, with the true anomaly as the independent variable-in which the complete Hamiltonian reduction of perturbed problems is naturally achieved. The fact that the Delaunay-similar elements are harmed by singularities is not at all of concern. Indeed, once the perturbation solution is computed, we can always reformulate it at our will in our preferred choice of singular or non-singular, canonical or non-canonical variables [31,15].

The main problem in Delaunay-similar variables
The main problem admits Hamiltonian formulation. Using polar variables, the Hamiltonian is written in the form (see [14] for instance) where µ, R ⊕ , and J 2 , are physical parameters that characterize the gravity field, and stand for the Earth's gravitational parameter, equatorial radius, and oblateness coefficient, respectively; r is the radius, R the radial velocity, the polar angle θ is commonly dubbed as the argument of the latitude, Θ is the total angular momentum (per mass unit), P 2 is the Legendre polynomial of degree 2, and s ≡ sin I, with I denoting orbital inclination. The third component of the specific angular momentum N = Θ cos I is an integral that stems from the axial symmetry of the main problem.
Hamiltonian (1) is conveniently formulated in the extended phase space [35,41] by adding the pair of variables q = t and Q = −H, and defining the new Hamiltonian H * ≡ H+Q, which, therefore, is constrained to the manifold H * = 0. If, besides, we define a new time scale τ , such that dt = Ψ(r, θ, R, Θ, q, Q) dτ , it can be checked that the Hamiltonian H ≡ H * Ψ remains constant in the new time scale, and that the solution of the system governed by H coincides with the system governed by H * when the initial conditions for τ = 0 are the same as the ones for t = 0 [38]. In particular, the choice allows us to replace Eq. (1) by the solutions of which only make sense in the manifold H = 0. When choosing which only differs from Θ on the order of J 2 , then Hamiltonian (3) is suitably formulated in the Delaunay-similar variables (φ, g, h, λ, Φ, G, H, Λ), where φ is the true anomaly, g is the argument of the perigee, h is the right ascension of the ascending node, λ is the time element, Φ is related to the Keplerian energy, G = Θ, H = N , and Λ = Q. They are computed from the polar variables following the sequence (see [40,17] for details) In these variables, Eq. (3) takes the form in which, now,

Perturbation approach
Due to the smallness of the Earth's J 2 coefficient, the complete Hamiltonian reduction of the main problem Hamiltonian (11) can be achieved, up to some truncation order, by the Lie transforms method [11,3,30]. We assume that the reader is enough familiarized with this perturbation method and only provide the results.
In our approach, we adhere to the tradition and carry out the sequential elimination of short-and long-period terms [7,40]. This is done analytically by changing first from original to prime variables, and then from prime variables to mean (double-prime) variables. The complete Hamiltonian reduction is then achieved by neglecting higher-order terms from the mean elements Hamiltonian, to obtain The analytical solution is trivial in mean elements. Indeed, the Hamilton equations of Eq. (13) immediately show that the momenta are constant in mean elements where H and Λ are integrals of the original main problem, whereas the angles evolve linearly with constant frequencies given by To avoid offending divisors in the perturbation series, the analytical solution is customarily formulated in the canonical set of non-singular Poicaré-similar elements [39]. However, due to the axial symmetry of the main problem model, we find convenience in formulating the solution in a non-canonical set of variables that replaces the troublesome elements by the argument of the latitude θ = φ + g and the semi-equinoctial elements defining the eccentricity vector in the orbital plane C = e cos g, S = e sin g. The solution in mean elements is thus and the original dynamics is recovered by plugging this solution into the transformation from mean to osculating elements, which in the current, Hamiltonian case is derived from a generating function.
We extended the solution of [40] to the second order. That is, the secular terms in Eqs. . Accordingly, the analytical solution comprises 33 perturbation series in addition to those defining the secular frequencies. Namely, 5 × 3 for the long-period corrections (first-order, second-order direct, and second-order inverse), and 6 × 3 for the short-period corrections (first-order, second-order direct, and second-order inverse). Remarkably, there is no need of integrating Eq. (2) to recover the physical time, which, on the contrary, due to the extended phase space formulation is obtained by making explicit q from Eq. (10).
Note that, due to the truncation of the periodic corrections, the initialization of the perturbation theory from a given set of initial conditions provides the needed constants in Eqs. (14)-(18) only up to O(J 2 2 ) effects, whereas the secular frequencies are expected to be accurate to the order of J 3 3 . While this is an inconsistency in traditional perturbation theories, which need to be patched to avoid large in-track errors, we will see that this uneven truncation is not of concern in the extended phase space formulation.

First-order solution
The generating function of the short-period elimination is where we defined the auxiliary variable Remarkably, Eq. (24) is formally the same as the first order term of the generating function of the elimination of the parallax transformation when it is written in Delaunay variables (cf. [27]). Up to the first order, the transformation of the short-period elimination is computed like ξ = ξ +J 2 δ 1,ξ , with ξ denoting any of the involved variables. The shortperiod corrections δ 1,ξ are obtained from the computation of the Poisson bracket {ξ, W 1 }, and reformulating the result in prime variables. The inverse transformation is ξ = ξ − J 2 δ 1,ξ where now δ 1,ξ is written in the original, non-primed variables.
For the long-period elimination we obtain the first order generating function which must be written in the prime variables. In Eq. (26) we abbreviated and introduced the auxiliary, non-dimensional functions which are O(J 2 ), and hence produce a small displacement of similar magnitude in the critical inclination value, which now happens for ∆ = 0. That is, sin I = 4/5 + O(J 2 ). The first-order transformation to mean (double-prime) elements is computed like before. Thus, ξ = ξ + J 2 δ 1,ξ , where the long-period corrections δ 1,ξ are obtained computing the Poisson bracket {ξ , V 1 }, and reformulating the result in the double-prime variables. Analogously, ξ = ξ − J 2 δ 1,ξ with δ 1,ξ now written in prime variables.
The secular frequencies of the first order theory are obtained from the Hamilton equations of the secular Hamiltonian where and are written in mean (double-prime) variables.

Second-order solution
The second order generating function takes the form Now where δ 2,ξ is obtained after reformulating in prime variables the result of {ξ, W 1 }, W 1 + {ξ, W 2 }. For the inverse transformation, where δ 2,ξ is no longer the opposite of δ 2,ξ , but the result, in the original variables, The secular frequencies of the second-order theory are computed from the Hamiltonian (29), which is now complemented with the third-order term where the coefficients q k,i,j are given in Table 1.
The second-order term of the generating function of the transformation of the long-period elimination is where the inclination polynomials b l,k,i,j are listed in Table 2. The transformation of the long-period elimination is obtained analogously to Eq. (33) -resp. (34)and following formulas, replacing corresponding functions and variables by those of the long-period case.
Remarkably, at this order there is an explosion in the number of terms of the transformation of the long-period elimination. This is clearly shown by comparison with Brouwer's 2nd order generating function of the long-period elimination, Still, if third-order corrections are not of concern, after computing the second-order corrections stemming from Eq. (36), they can be notably simplified by making δ = υ = 0. Other terms of the different series that compose the solution in the extended phase space are certainly amenable of analogous simplifications.

Accuracy tests
A number of tests has been conducted to check the accuracy of the analytical solution based on the fictitious time and compare it with the precision provided by 4,0 −10368c 6 (13s 2 − 10) 6912c 6 (5s 2 − 2) −3888c 6  the traditional approach based on the physical time. For the later, the different analytical solutions in the literature for the same truncation order are essentially different arrangements of the same solution, which, therefore enjoy analogous accuracy [24]. However, one formulation may be preferred over others depending on the user needs [25]. For our comparisons we choose the approach based on the sequential elimination of the parallax [12,27], the elimination of the perigee [2,26], and the Delaunay normalization [13], which seems to be a popular option in the literature [9,36,10,23]. Remark that, as opposite to common implementations by different authors, we always improve the initialization of the constants of the traditional perturbation solution with the calibration of the mean semimajor axis from the energy equation [6].
For the different orbital configurations tested, we always obtained similar comparative results between the traditional and fictitious-time modalities. Therefore, we only illustrate them for the PRISMA orbit [34], which is a low-eccentricity, sunsynchronous orbit at around 700 km altitude. In particular, in the computation of the initial conditions we used the orbital parameters a = 6878.14 km, e = 0.001, I = 97.42 • , Ω = 168.2 • , ω = 20 • , and M = 30 • . The true, reference orbit was then computed from the numerical integration of the main problem in Cartesian coordinates. To guarantee that all the stored digits of the reference orbit are exact in the floating-point number representation, the integration was carried out in extended precision. In addition, the integration of the fictitious time variation, given by Eq. (2), was incorporated to the differential system in order to check the accuracy provided by each theory exactly at the same physical time.
Results showing the Root Sum Square (RSS) of the errors provided by the firstorder truncation of each perturbation solution are shown superimposed in Fig. 1 for a propagation interval of 10 days. The small secular trend in the errors of the traditional solution, of about 1 meter per day, is mostly nullified in the case of the solution based on the extended phase formulation, where the much smaller secular trend remains buried under the periodic oscillations of the errors. The nature of the position errors of each analytical solution is best illustrated when they are projected in the radial, along-track, and cross-track directions. Indeed, by simple inspection of Fig. 2, we check that the errors in the radial and cross-track directions are periodic in nature and of comparable magnitude in both theories -as it should be expected for perturbation solutions truncated to the same order. On the contrary, while the amplitude of the periodic components of the errors in the along-track direction remain in the case of traditional theory of the same level of accuracy as the other components of the errors, they reduce by about one order of magnitude, from meters to decimeters, in the case of the extended phase space formulation. The better performance of the extended phase space solution regarding along-track errors applies also to their evident secular component (bottom plots of Fig. 2), which are mostly responsible from the different behavior of both theories, and clearly illustrate the benefits of including the total energy among the integration variables regarding the stabilization of the secular errors in the along-track direction.
Undeniably, the greater secular trend of the errors of the traditional solution seems to make the extended phase space solution preferable in much longer propagation times. Still, for actual ephemeris computation, the fictitious time in which the extended phase space solution must be evaluated corresponding to the desired physical time will not be known in advance. On the contrary, its accurate computation requires a root-finding procedure that is commonly approached by Newton- Raphson iterations [5]. Moreover, due to the unavoidable truncation of the perturbation solution, the fictitious time corresponding to a given physical time cannot be determined exactly. As shown in Fig. 3, the nature of these errors is mostly periodic, with an amplitude of about one half of a millisecond. For the semimajor axis of the PRISMA orbit, this error in the timing may yield an additional in-track error in the meter level. These additional inaccuracies derived from the O(J 2 ) truncation of the series that provide the physical time as a function of the fictitious one, which were not presented in Fig. 1, are now illustrated in Fig. 4 where the along-track errors of the first-order theory in the extended phase space (in black) are computed taking the physical time as the argument, from which the needed fictitious time is then determined. This process requires the evaluation of the analytical theory in each iteration, thus notably increasing the computational burden, and, therefore, degrading the performance of the perturbation theory regarding computing time. When the second-order theory is used we obtain notable accuracy gains in both cases, as expected from the characteristics of perturbation solutions. Again, the amplitude of the periodic errors in the radial and cross-track directions is similar for both perturbation theories, now in the mm level, and corresponding errors are not presented. On the contrary, while the amplitude of the periodic errors in the along-track direction remains also in the mm level in the case of the solution in extended phase space, it is not at all the case of the errors of the traditional theory, whose amplitude now reaches the cm level. This is shown in Fig. 5, where we also notice the secular rate of the along-track errors of the traditional theory of about one third of cm per day, which is much smaller than the rate of about one tenth of mm per day of the extended phase space solution. Nonetheless, these figures yield an analogous ratio to the first-order case.
Again, additional inaccuracies to those presented in Fig. 5 arise in the case of the extended phase space solution due to the need of computing the physical time from the fictitious one, or vice versa. As shown in Fig. 6, the amplitude of these errors is now reduced to about 1 microsecond, which is not enough to hidden a clear secular trend of about 0.1 µs/day. These additional errors of the extended phase space formulation, now reaching the cm level for the PRISMA orbit, result in a notably increase with respect to the mm level shown in the bottom plot of Fig. 5. Again, timing errors notably balance the accuracy performance of both kinds of theories, as illustrated in Fig. 7. Therefore, like in the first-order case, this source of errors may be a dominant part that cannot be ignored in the evaluation performance of perturbation solutions based on the formulation in the extended phase space.

Conclusions
Inclusion of the total energy among the variables of an analytic orbit generator has beneficial, radical effects in the propagation of along-track errors of circumterrestrial orbits, which are commonly reduced by one order of magnitude with respect to usual along-track errors obtained with traditional (properly initialized) perturbation solutions. This fact could make the extended phase formulation preferable for mission analysis and planning. However, this kind of solution has to deal unavoidably with the transformation between physical and fictitious times in ephemeris computation. This is not only a clear inconvenience from the algorithmic point of view, but, most notably, this change cannot be carried out exactly due to the truncation that is inherent to any perturbation solution. Additional inaccuracies stemming from this conversion, which show both secular and periodic trends, may be as important as those derived from the intrinsic components of the position errors. Therefore, they cannot be ignored in the objective assessment of accuracy performance. From the point of view of the construction of the analytical solution, the fictitious-time approach has clear advantages regarding the computation of the integrals involved in the process, yet a post-processing is needed to reduce the length of the series comprising the solution to a comparable size to those resulting from the traditional approach. These pros and cons do not present a clear advantage for one of the options, and make the choice between traditional and extended phase space orbital perturbation solutions mostly a matter of the particular needs of prospective users.