A Second‑Order Closed‑Form J 2 Model for the Draper Semi‑Analytical Satellite Theory

A second-order closed-form semi-analytical solution of the main problem of the artificial satellite theory ( J 2 contribution) consistent with the Draper Semi-analytic Satellite Theory (DSST) is presented. This paper aims to improve the computational speed of the numerical-based approach, which is only available in the GTDS-DSST version. The short-period terms are removed by means of an extension of the Lie-Deprit method using Delaunay variables. The averaged equations of motion are given explicitly and transformed to the non-singular equinoctial elements. Finally, the second-order terms in the equations of motion are included in the C/C++ version of the DSST orbit propagator.


Introduction
The Draper Semi-analytic Satellite Theory (DSST) orbit propagator can be found in two forms, as an option within the Massachusetts Institute of Technology version of the Goddard Trajectory Determination System (GTDS) computer program [1,2], and as the DSST Standalone orbit propagator package [3][4][5][6][7].
The original implementations of the DSST, both in GTDS and in the Standalone versions, were done in Fortran 77 (F77). Between 2012 and 2015, the DSST was re-implemented in Java and included in the Orekit flight dynamics library [5,6]. During the same time frame, the University of La Rioja provided web access to the F77 DSST Standalone via a friendly and intuitive interface [8]. In 2016, Setty [7] extended the F77-DSST-Standalone force models, state-transition matrix, and dynamic-parameter partial derivatives to provide the F77-DSST Standalone with orbit determination capability. More recently, the migration of the F77-DSST Standalone code to C/C++ has been done at the University of La Rioja [9,10].
The theory underlying DSST makes use of the non-singular equinoctial elements and is based on the generalized method of averaging [11][12][13]. The Lagrange Planetary equations are used for handling the conservative terms, 50 × 50 gravity fields (M-daily, tesseral resonance and approximate solution of the second-order effect in J 2 , J 2 2 , based on eccentricity expansions), lunar-solar point masses, and the solid Earth tides, whereas the Gauss equations are used for the atmospheric drag and the solar radiation pressure (SRP).
In particular, the J 2 2 effect was formulated by Zeis [14] in DSST. This author developed an approximate second-order solution as a power series of the eccentricity in J 2 for mean-element and short-periodic motions. This analytical model was implemented with the help of the computer algebra system, Macsyma. Latter, Fisher [15] described the computational procedure required to obtain the J 2 2 effect in closed form. This author proposed two alternatives methods for completing this effect: first, using numerical quadratures and, second, symbolically expanding and manipulating all the mathematical expressions. Unfortunately, time and computer constraints prevented obtaining the full second-order effect. Recently, Folcik and Cefola [16] ported this code to Maxima, the open-source descendant of Macsyma, to calculate the integrant of the mean element rate expressions. The averaged process was carried out using Gauss-Kronrod numerical quadrature. These authors also performed a detailed numerical comparison between Zeis and their numerical model. Finally, they concluded with the necessity of developing an analytical closed-form model to replace the numerical quadrature process and, so, improve the speed of evaluating the theory.
To obtain the closed-form averaged equations and the mean-to-osculating transformation, we propose an alternative method to build fully second-order closedform analytical expressions for the J 2 2 contribution based on the Lie-transform method [17][18][19][20][21] and canonical variables consistent with the DSST orbit propagator.

3
The main difference between the Lie transform method and the generalized method of averaging is regarded more as a matter of the algorithm than the underlying idea. In fact, both methods are connected through an integration constant [22,23]. It is worth noting that the Lie transform method allows us to obtain the transformations from osculating-to-mean and mean-to-osculating elements simultaneously from the generating function.
The main characteristic observed in the application of the generalized method of averaging in the zonal problem in DSST can be found in the calculation of the zonal harmonic short-periodic generator in DSST [24] to which must be added an integration constant to guarantee that the short-period terms do not contain any long-period terms. The equivalent approach is followed by Kozai [25] using Von Zeipel method [26].
In this work, we present a second-order semi-analytical theory for the J 2 problem to remove the short-period terms from the equations of motion using Hamiltonian formalism, an extension of the Lie-Deprit method [27], which is implemented in MathATESAT [28], and Delaunay variables. The resultant theory is equivalent to the elimination of mean anomaly in Kozai [25]. Then, the equations of motion are expressed in equinoctial elements. After that, the transformed equations of motion are analytically and numerically validated. Finally, the Mathematica expressions are migrated to C and included in DSST C.

On the J Problem and its Normalization
Extensive investigations dealing with the J 2 problem, or classically called the main problem of artificial satellite theory, have been done since Brouwer's work [29]. In this section, we present a semi-analytical theory equivalent to the first canonical transformation described in Kozai [25] using an extension of the Lie-Deprit method.
The J 2 problem is defined as a Kepler problem perturbed by Earths oblateness. The Hamiltonian of this dynamical system can be written in a cartesian coordinate system ( , ) as is the satellite velocity, is the gravitational constant, the equatorial radius of the Earth, J 2 the oblateness coefficient and P 2 the second-degree Legendre polynomial.
The first step to carry out the analytical transformation consists of expressing the Hamiltonian (1) in terms of Delaunay variables ( , g, h, L, G, H) [30,31]. This set of canonical action-angle variables can be defined in terms of the orbital elements such as

3
The Journal of the Astronautical Sciences (2022) 69:1292-1318 where M, , Ω , a, e, i are the mean anomaly, argument of the perigee, longitude of the ascending node, semi-major axis, eccentricity and inclination, respectively. Then, the transformed Hamiltonian is given by where where = J 2 is a small parameter, s = sin i and f is the true anomaly. To avoid any ambiguity in the notation between Delaunay variable and the equinoctial elements (a, h, k, p, q, ) , we use Ω instead of h to refer to the argument of the node in Delaunay variables. Next, we normalize the Hamiltonian (3) by applying the Lie transform , the so-called Delaunay Normalization [32], using the Extended Lie-Deprit method [27]. This perturbation method is based on transforming the analytic Hamiltonian function into the new one which satisfies some specific prerequisites. The double index notation in the Hamiltonian is introduced to simplify the automatization of this algorithm. As well as the Lie-Deprit method [18], the extended method looks for a generating function W n = ∑ n≥0 n n! W n+1 of so that the terms H n , K n and W n verify the partial differential equation, called the Homological Equation, where L H 0 is the Lie operator associated with H 0 , a linear operator given in terms of a Poisson bracket. (2) The right-hand side of the homological equation, H 0,n , is computed from H n , (W i ) 1≤i≤n−1 and (H p,q ) p+q≤n−1 , where the latter are obtained by means of the recursive formula with i ≥ 0 and j ≥ 0 and { , } represents the Poisson bracket (for more details, see [18]). On the other hand, C n is an arbitrary integration function which can depend on In this case, the Lie operator is defined from the zero order Hamiltonian H 0 as L H 0 = n ∕ � , where n = 2 ∕L �3 is the mean motion. On the other hand, the solution of the homological equation is obtained as where C n , in accordance with Kozai's work, is chosen to remove the long-period terms contained in the generating function and is given by This semi-analytical theory inherits the intrinsic singularities of the Delaunay variables for the eccentricity, e = 0 , and for the inclination, i = 0, . However, the critical inclination is not a problem in the present theory as DSST because perigee's motion remains in the equation of motion. In the remainder of this paper, we shall drop the primes on the transformed variables to alleviate the notation.
After this summary, we briefly outline the operations that accomplished the elimination of the short period terms using the extended Lie-Deprit method, which up to first-order reads Using Eq. (10), we obtain the first-order transformed Hamiltonian as the average over the fastest angle :

3
The The direct and inverse transformations of the Lie transform are obtained from the generating functions W 1 and W 2 using the classical recurrent algorithms [18][19][20]. (16)

3
It is worth noting that both transformations expressed in Delaunay variables suffer from the inclusion of the eccentricity in the denominator. Their explicit expressions are not provided here due to a large number of terms. Finally, Hamilton's equations, the partial derivatives of the Hamiltonian (19) with respect to the new variables, provide the equations of motion as The first-order coefficients Δ In order to include this theory in DSST, the transformations and equations of motion are formed in terms of the non-singular equinoctial element set.

Mean Equinoctial Variational Equations
In this section, we derive the mean equinoctial variational equations from Hamilton's equations (21). The equinoctial elements [33] are defined in terms of the orbital elements as follows: where I is called the retrograde factor. The factor I takes the value 1 for the direct equinoctial elements and -1 for the retrograde equinoctial elements, although direct orbits can include the range of inclination [0, ).
Finally, we will obtain Eqs (28) in equinoctial elements taking into account Eqs. (25) and the transformation from the equinoctial elements to the orbital elements given by where is an auxiliary angle, which is defined by: The zero-order terms in equinoctial elements are is the mean motion. The first-order terms only depend on a =a,

3
The Journal of the Astronautical Sciences (2022) 69:1292-1318 the momenta L, G and , s, c allowing us to obtain the Δ 1 terms in a straightforward manner as where C 1 = 3 2 n∕(4a 2 4 ) and = 5s 2 + 2Ic − 4 . The expressions of , c and s in equinoctial elements are As it can be seen, Eq. (31) are valid for both direct and retrograde equinoctial elements and agree with known results [14,34]. The expressions of the second-order terms to equinoctial elements are not directly calculated, The second-order Hamilton equations introduce the contribution of the argument of the perigee g through terms in sin 2g and cos 2g . Using Eq. (25), these trigonometric functions are converted to equinoctial elements and yield After a laborious simplification process, we obtained that the Δ 2 terms were not entirely non-singular. These terms depended on the factor (p 2 + q 2 ) −2I . Therefore, it was not possible to derive a unique set of expressions valid for both direct and retrograde orbits.
For the case of direct orbits (I = 1) and the range of inclination [0, ) , we obtain the following non-singular expressions  Table 1.
Finally, the second-order terms have been tested reproducing the solution provided by Zeis [14]. This approach has been derived from the beginning considering the first-order approximation in the eccentricity of the transformed Hamiltonian K 2 (19) which yields Following the same process used to build the Δ 2 terms (33), the corresponding second-order terms in the Zeis case are (34)  where C Z 2 = 3 4 n∕(4a 4 ) . In this simplified case, the expressions do not include any singularity and are valid for both direct and retrograde orbits and agree with known results [14,34].

Numerical Experiments
Validation of any analytical or semi-analytical theories is not a trivial task since the accuracy of the theory mainly depends on its correct initialization [35,36]. In the semi-analytical theory presented in this paper, the osculating-to-mean and mean-toosculating transformations are expressed in the singular Delaunay variables, which introduce an additional problem to its validation. Both issues do not allow us to reproduce or reach the accuracy of the results presented in Fig. 6 in the work of Folcik and Cefola [16]. These authors use orbital determination techniques to adjust the initial conditions properly, whereas the mean-to-oscillating transform is formulated in the non-singular equinoctial variables.
Finally, we compared the first-order, the truncate second-order (Zeis's solution), and the closed-form second-order semi-analytical theories (SA) with the non-transformed equations of motion for three families of orbits. The first family, low-altitude near circular orbits, has a 200 km perigee height by 210 km apogee, the second 200 km perigee height by 500 km apogee, and the third 200 km perigee height by 1000 km apogee. The inclination varies from zero to 180 degrees. Both transformed and non-transformed equations of motion were integrated using a highly accurate eighth-order Runge-Kutta method [37] (NUM) for two days. Then, their respective outputs, SA i and NUM i , were then compared in terms of their RMS of position differences, as given in Eq. (38).
The equations of motion are formulated in equinoctial elements, whereas the changes from mean-to-osculating and osculating-to-mean remain in the Delaunay variables. Figure 1 describes the evaluation process of the semi-analytical theories. The initial osculating elements t 0 are transformed to initial mean elements ̄ t 0 using the osculating-to-mean transformation. Then, the mean equations of motion are numerically integrated. Finally, the osculating elements t are obtained through the mean-to-osculating transformation. It is worth noting that the non-transformed equations of motion are integrated starting with the same initial conditions t 0 . Figure 2 shows the position error RMS between the first-order semi-analytical theory and the numerical integration of the non-transformed problem for the three families of orbits. As can be seen, the errors are symmetric for the 90 degrees of inclination. For the 200×210 km family, the maximum error is 10 km, obtained at an inclination of 90 degrees, whereas the minimum is 4.37 km, at an inclination around 52 and 128 degrees. The maximum errors are reduced as the apogee height is increased to 500 km and 1000 km to 2.2 and 2 km for inclinations 0 and 180 degrees, respectively. The minimum errors are approximately 0.47 km at inclinations around 74 and 106 degrees. Figure 3 shows the position error RMS for the truncate (blue) and the closedform (red) second-order semi-analytical theories for the family near-circular orbits ( e = 0.000759 ). The position errors remain the symmetry to the 90 degrees of inclination. As can be seen, the truncate and the closed-form solutions are very close to each other, which is consistent with the small eccentricity case. The difference between both lines is very small. The maximum error of the two theories   Figure 4 depicts the difference between the closed-form and the truncate position errors plotted in Fig. 3. The accuracy provided by the truncate solution is only better between 65 and 180, whereas the closed-form solution is in the rest. We assume that this behaviour mainly lies in the small eccentricity of this family and in conserving the transformations from osculating-to-mean and mean-to-osculating in Delaunay variables. Figure 5 shows the position error RMS for the truncate (dashed line) and the closed-form (continuous line) second-order semi-analytical theories for two families of orbits with 500 km (red) and 1000 km (blue) apogee height, respectively. As can be seen, the position errors are symmetrical to the 90 degrees of inclination. The loss of accuracy produced by conserving the singular Delaunay variables in the transformations is visible in the family with 500 km apogee height as inclination varies from 0 to 50 and 130 to 180 degrees. The maximum errors of both theories are reached for the inclination of 0 and 180 degrees, with around 180 m for the closed-form and 160 m for the truncate. The minimum errors or the closed and truncated theories are around 63 m and 45 m for an inclination of 55 and 45 degrees, respectively. However, this behaviour changes for the family with 1000 km  Finally, Fig. 6 simulates the current status of the different implementations of the J 2 2 effect in GTDS-DSST and DSST Standalone using the closed-form and truncate mean equations of motion and the families of 500 km (red) and 1000 km (blue) apogee height. The second-order osculating-to-mean transformation was applied to the initial conditions. On the other hand, the change from mean to osculating is done using the first-order expressions. For both families of orbits, the accuracy of the closed-form equations of motion slightly predominates over the accuracy given by the truncate equations. It is not clear for the family of 500 km because both lines are close to each other, as can be seen in Fig. 7. However, this behaviour can be appreciated in the family of 1000 km apogee height, especially for inclinations from 40 to 130 degrees. It is worth noting that the asymmetry errors at 0 and 180 degrees inclination detected in [16] do not appear in our analytical expressions.

Conclusion and Future Work
In this paper, a second-order closed-form semi-analytical solution of the J 2 problem has been developed. An extension of the Lie-Deprit method and Delaunay variables is used to remove the short-period terms from the equations of motion using the Hamiltonian formalism. This semi-analytical theory has been implemented using MathATESAT a Mathematica-based tool. The averaged equations of motion are given explicitly and transformed into the non-singular equinoctial elements. This semi-analytical theory is equivalent to others in which Lagrange Planetary equations, as well as General Averaging Method and equinoctial elements, are used, and hence, consistent with the Draper Semianalytic Satellite Theory. These analytical expressions improve the computational speed of the numerical-based approach, which is only available in the GTDS-DSST version.
Additionally, the second-order terms have been tested reproducing the solution provided by Zeis, validated numerically, and included in DSST C.
Finally, we are currently working on expressing the second-order mean-to-osculating transformation in equinoctial elements, calculating the partial derivatives necessary for an orbit determination system and integrating GTDS-DSST and DSST Standalone. It is worth noting that these expressions are essential to reach the complete accuracy of the second-order theory. Fig. 7 Difference between the closed-form and truncate position errors plotted in Fig. 6 Appendix A W 2 The second-order generating function is where S 1 = s 2 5s 2 − 4 , S 2 = s 2 3s 2 − 2 and P i represent polynomials in given by