An Adaptive Analytic Continuation Method for Computing the Perturbed Two-Body Problem State Transition Matrix

In this work, the Taylor series based technique, Analytic Continuation is implemented to develop a method for the computation of the gravity and drag perturbed State Transition Matrix (STM) incorporating adaptive time steps and expansion order. Analytic Continuation has been developed for the two-body problem based on two scalar variables f and gp and their higher order time derivatives using Leibniz rule. The method has been proven to be very precise and efficient in trajectory propagation. The method is expanded to include the computation of the STM for the perturbed two-body problem. Leibniz product rule is used to compute the partials for the recursive formulas and an arbitrary order Taylor series is used to compute the STM. Four types of orbits, LEO, MEO, GTO and HEO, are presented and the simulations are run for 10 orbit periods. The accuracy of the STM is evaluated via RMS error for the unperturbed cases, symplectic check for the gravity perturbed cases and error propagation for the gravity and drag perturbed orbits. The results are compared against analytical and high order numerical solvers (ODE45, ODE113 and ODE87) in terms of accuracy. The results show that the method maintains double-precision accuracy for all test cases and 1-2 orders of magnitude improvement in linear prediction results compared to ODE87. The present approach is simple, adaptive and can readily be expanded to compute the full spherical harmonics gravity perturbations as well as the higher order state transition tensors.


Introduction
In Astrodynamics, the State Transition Matrix (STM) of the two-body problem works as a sensitivity of the current states to the initial conditions. Hence, it computes the propagation of error of the initial states over time. Computing the STM is ubiquitous to spaceflight dynamics, navigation and control, [20,21,26,32]. Goodyear presented an exact analytical solution of the two-body Cartesian STM and the method is valid for all types of orbits for the attractive force and for the hyperbolic and rectilinear orbits for the repulsive force, [10]. However, it is complicated to introduce third body effect in this method and the method requires transcendental function evaluation, [21]. This method has been simplified to increase the numerical efficiency for the Keplerian elliptical orbit, [6,19]. However, gravity perturbation was not incorporated in this simplified method, [6]. Recently a decoupled direct method is developed by Hatten and Russel using fixed step size Runge-Kutta method, [13]. In this method, first and second order STM is decoupled from the state propagation to make the procedure computationally efficient. To formulate the STM with Spherical Harmonic Gravity model, Read et al. applied Modified Chebyshev Picard Iteration method, [25]. The method is shown to be well suited for parallel implementation for additional speed of computation. Markley presented a J 2 perturbed approximate Cartesian STM, [21]. This method is also capable of handling the n-body problem. A major drawback of this method is that it is restricted by short time intervals. It can compute the J 2 perturbation term for only up to 6 second time step.
A geometric method is developed by Gim and Alfriend to derive the STM for the relative motion including gravitational perturbation, [9]. This method uses the relationship between the relative states and the differential orbital elements. To implement the STM in practical engineering applications, Yamanaka and Ankersen linearized the differential equations of relative motion on elliptical orbits, [35]. This linearized method is valid for unperturbed elliptical orbits. Koenig et al. derived the STM using classical Keplerian orbit elements, [18]. They introduced J 2 and differential drag perturbation in their work.
In low earth orbit (LEO), there is another inevitable perturbation, which is atmospheric drag. It depends on the shape, size and orientation of the orbiting body. Drag perturbation has been extensively studied, [5,7,17]. Mavraganis and Michalakis reviewed drag along with radiation pressure exerted by the major body on the minor one, [22]. Vallado and Finkleman discussed the sources of uncertainty in the drag models and evaluated the interrelation among different parameters of the problem, [34]. The literature is rich with numerical techniques that study the orbit problem in the presence of perturbations, [33]. More recently, parallel implementation of Modified Chebyshev Picard iteration on various gravity fields using multiple cores and Graphic Processing Units (GPUs) are discussed, [1]. Atallah et al. also focused on the comparison among different numerical integrators for perturbed orbit propagation in terms of accuracy and efficiency, [2]. The Analytic Continuation method is a numerical integration method that is based on Taylor series expansion and Leibniz product rule. This method has already been proven to be very precise in solving unperturbed and J 2 − J 6 gravity perturbed two-body problem, [16,31]. Later on, exponential atmospheric drag model is introduced, [14]. Recently, the method has been further developed to handle the full spherical harmonics gravity potential up to 250 degree and order, [15].
In this paper, the Analytic Continuation Method is used to derive the State Transition Matrix for the perturbed two-body problem. In "STM Derivation Using Adaptive Analytic Continuation", the procedure to derive the drag and gravity perturbed adaptive STM is discussed. In "Numerical Simulations", first results for the unperturbed case are compared with Lagrange's F and G solution. Then the method is expanded for the J 2 − J 6 gravity perturbation terms and compared against ODE45 (based on Runge-Kutta (4,5) integration method, [4,8]), ODE113 (based on variable order Adams-Bashforth-Moulton predictor corrector method, [27,28]) and ODE87 (based on 8-7th order Dormand and Prince formulas, [11,12,24]). Atmospheric drag model is included and error propagation of the states using Analytic Continuation STM are checked against the results using ODE87. Finally, discussion of the results and concluding remarks are presented in "Discussion" and "Conclusion", respectively.

STM Derivation Using Adaptive Analytic Continuation
In the Analytic Continuation technique, two scalar variables are defined as, f = r.r and g p = f −p/2 , where r is the position vector and p is an integer. In the two-body problem, the acceleration vector is defined as, where, r(t) is the position vector at current time and μ is the standard gravitational parameter. Next, Leibniz product rule is implemented to derive the recursive formulas for computing the higher order time derivatives of the variables: Finally, the higher order time derivatives of the position vector are substituted into the Taylor series expansion to obtain position and velocity of the next time step as shown in Eqs. 5 and 6, To increase efficiency of computation, adaptive time step and expansion order are applied to the method, [15]. As in embedded Runge-Kutta methods, the adaptation is based on defining absolute and relative tolerances, δ a and δ r , respectively. The expansion order is then determined from, where, δ is defined as, min(δ a , |r|δ r ) in the first step and max(δ a , δ a × f ac) in the subsequent steps. Here, f ac for canonical units is defined as, The time step is then determined using, where, dT a and dT b are calculated from δ×(N−1)!
, respectively which represent the last two terms in the Taylor series. The adaptation parameters δ f ac , k inc and h f ac are tuned towards accuracy or speed and to increase the method robustness to handling different types of orbits.
According to the definition of the STM for the two-body problem, [3,26] ∂r(t+dT ) ∂r (1) (t) ∂r (1) ∂r (1) (t+dT ) ∂r (1) where, φ 11 is the sensitivity of the next position to the current position, φ 12 is the sensitivity of the next position to the current velocity, φ 21 is the sensitivity of the next velocity to the current position and φ 22 is the sensitivity of the next velocity to the current velocity. The same approach of the Taylor series expansion of the state variables are followed to expand the elements of the STM, [29] (1) (t + dT ) ∂r (1) where, dT is the time step between the current position and the next position and ψ (m) r and ψ (m) v are defined as ∂r (m) (t) ∂r(t) and ∂r (m) (t) ∂r (1) (t) , respectively.
The partial derivatives of r(t), r (1) (t), f and g p with respect to r(t) and r (1) (t) can be written as, ∂r(t) = ∂r (1) (t) ∂r (1) ∂r (1) (t) = ∂r (1) ∂r (1) ∂r (1) By substitution, the recursive formulas to calculate partial derivatives of the higher order time derivatives of the variables, with respect to position and velocity are given by, where, χ = r(t) or r (1)

Zonal Perturbations
For the J 2 − J 6 zonal perturbations, the higher order partial derivatives are obtained. In this section we show the detailed derivation for the J 2 perturbation and refer the reader to the appendix for the detailed derivation of J 3 − J 6 . The J 2 acceleration is given by, [16,26] where, r eq is equatorial radius of earth and J 2 = 1082.63 × 10 −6 , [26].
In order to obtain the recursions for the higher order partials, we define a constant, C J 2 = − 3 2 J 2 μr 2 eq , and 2 new variables, B p and C α , as, The higher order time derivatives of the newly defined variables are given by, Hence, the higher order time derivatives of the J 2 perturbation acceleration is computed as, The partial derivatives of Eqs. 21 and 22 with respect to position and velocity are then expressed as, where, χ is r(t) or r (1)

(t).
Finally, Eq. 25 is added to Eq. 11 through Eq. 14 to obtain the J 2 perturbed STM as,

Atmospheric Drag Perturbation
A canon ball drag model is re-derived using analytic continuation method, [14]. The acceleration for the drag model is given by, [34] a dr = r (2) where, ρ(r) is the atmospheric density, C D is the drag coefficient, A is the reference area, and v rel is the relative velocity of the satellite with respect to Earth defined as, From Eq. 31, the higher order time derivatives of the relative velocity is, To calculate the higher order time derivatives of the relative velocity magnitude, ||v rel || (n) , two new variables similar to f and g p , Eqs. 3 and 4, are introduced as, Applying Leibniz product rule, the higher order time derivatives of the variables in Eq. 33 are expressed as, ρ(r) is defined according to the Exponential Model of atmospheric density, [33], as where, ρ 0 is the nominal density, R eq is the equatorial radius of Earth, h ellp is the altitude, h 0 is the base altitude, and H is the scale height. From Eqs. 30 and 36, the equation ofρ is given by, where, β is a constant. The first order time derivative ofρ is, Applying Leibniz product rule on Eq. 38 the higher order time derivative of ρ is derived as,ρ The higher order time derivatives of the acceleration due to drag is calculated using, [14] a (n) The recursive formula for the partial derivatives of the drag perturbation acceleration with respect to position and velocity is then expressed as,  (1) are then added to Eqs. 26 to 29 to include drag perturbation in computing the STM.
Algorithms 1 -3 in the Appendix show the steps for the full implementation of the analytic continuation method to compute the J 2 − J 6 and drag perturbed STM for the two-body problem.

Numerical Simulations
In this section, numerical results are presented for the derived STM via the Analytic Continuation technique for 4 orbits as shown in Table 1. Each orbit is propagated for 10 orbit periods.
First, Analytic Continuation is utilized to compute the STM for the unperturbed orbits and compared versus the closed-form solution by Battin, [3], via computing the RMS error. Next, J 2 − J 6 perturbation is added and the STM computed with Analytic Continuation is compared against MATLAB ODE45, ODE113 and ODE87 in terms of accuracy via the symplectic check. Finally, an exponential drag model is added to the J 2 − J 6 perturbed STM and the results are compared against ODE87 via error propagation. All codes are written and compiled using MATLAB R2019a and canonical units. To compare the numerical results, three different methods are used; calculation of RMS error, the symplectic check and error propagation of the states.

RMS Error for the Unperturbed STM
The RMS Error for each element of the STM in the time domain of 10 orbit periods is computed from the difference between the elements of the unperturbed STMs computed using Battin's method, [3], and the Analytic Continuation technique as shown in Eq. 42.
where, M ijk and L ijk are (i,j)th terms of the STM at k-th time step from Battin's method, [3] and Analytic Continuation method, respectively, E ij is the RMS Error of the (i,j)th term of the unperturbed STM and N is the total number of steps. The RMS errors for the 36 elements of the STM are calculated using Eq. 42 for the four orbit test cases and the results are shown in Tables 2, 3, 4 and 5, [30].

Symplectic Error for the Gravity Perturbed STM
The symplectic nature of the J 2 − J 6 perturbed STM is examined and elements of the error matrix are plotted versus orbital periods. The matrix [φ] is called symplectic, if it satisfies Eq. 43, [26] [ where, [φ] is the STM and [J ] is a skew-symmetric matrix defined by, The absolute error in the symplectic nature of the J 2 −J 6 perturbed STM is calculated by Eq. 45, [25] [ where, [E sym. ] is the symplectic error matrix. For the J 2 − J 6 perturbed cases, the results of the 10 orbit periods of the elements of the [E sym. ] matrix of every step using Analytic Continuation method are compared with the results using ODE87, ODE113 and ODE45. In this case, the initial and final time is provided to the solvers and the flexibility is given to the methods to select the time steps for the calculation. The results of the four different orbits are plotted as scattered points versus orbital period as shown in Figs. 1, 2, 3, 4, 5, 6, 7 and 8. Next, the gravity perturbed STM is generated up to 1,000 orbit periods using Analytic Continuation and compared against the results of ODE87, ODE113 and ODE45. The comparison results are shown in Figs. 9, 10, 11 and 12 via the average of the 36  elements of the symplectic error matrix, E sym. , versus orbital period. The Analytic Continuation is then used to propagate the STMs of the four orbit types for 10,000 orbit periods. The average symplectic check results at per orbit is used and the results are shown in Fig. 13. When performing these symplectic error calculations, both relative and absolute tolerance are set to 10 −13 to get the highest possible accurate results for the STM generation from the MATLAB ODE suite.

Error Propagation for Gravity and Drag Perturbed STM
To check the error propagation of the states using the STM, first the trajectory, x i , is computed with a nominal set of initial conditions, Table 1. A perturbation is introduced to the initial states as shown in Table 6. Linear prediction is used to compute the trajectory,x i , at the desired time steps using the STM as shown in Eq. 46. The "true" trajectory is propagated separately using the perturbed initial conditions. The error is then computed as the difference between the true states, x ⊗i , and the predicted statesx i , [9]. The error propagation compares J 2 − J 6 gravity and drag perturbed Analytic Continuation results for the four orbit types in Table 1 against ODE87 for 10 orbit periods. The results of the error propagation of the position and velocity are shown in Figs. 14, 15, 16, 17, 18, 19, 20, and 21.

Discussion
As shown in Tables 2 to 5 the unperturbed STMs computed via the Analytic Continuation technique maintained 11 -16 digits of accuracy in the RMS error across all 36 elements for all types of orbits when compared against the analytic solution by Battin, [3]. The check is performed for 10 orbit periods and the RMS error is computed for each element at each time step.       accuracy, ODE87 performed the best relative to ODE113 and ODE45, and ODE113 marginally outperformed ODE45. However, for ODE87, ODE113 and ODE45, the accumulation of truncation and round-off errors are evident in the loss of accuracy as the simulation time increases, whereas for Analytic Continuation, the ability to adapt to arbitrary order in the Taylor series enabled it to maintain double precision throughout the simulation. As the eccentricity, and consequently the nonlinearity, [23], is increased, the accuracy of ODE87, ODE113 and ODE45 degrades over the 10 orbit periods, whereas Analytic Continuation maintains double precision for MEO, GTO and HEO as shown in Figs. 3 to 8. The Analytic Continuation method is then used to generate STMs up to 1,000 orbit periods and compared against ODE87, ODE113 and ODE45 using the average symplectic error of the STMs, as shown in Figs. 9 to 12. Results for 1000 orbits and beyond (Fig. 13) show that the Analytic Continuation method is able to maintain double precision in the symplectic check whereas the other methods lost between 10 -12 digits of accuracy at 1000 orbits. It has to be noted that all numerical integration methods, including Analytic Continuation, were set to the lowest absolute and relative tolerance which explains the similarity in results of ODE45, ODE113 and ODE87 at 1000 orbits. If tolerances were to be set to their default values the lower order methods (particularly ODE45) will lose an additional 4-5 orders of magnitude when compared with ODE87. Finally, in the case of drag perturbation, linear prediction is used to verify the results where Analytic Continuation is compared against ODE87 for the four types of orbits. From the results presented in Figs. 14 to 21, Analytic Continuation outperforms ODE87 in linear prediction error in all test cases. For the LEO case, Analytic Continuation is two orders of magnitude more accurate compared to ODE87 in position, Fig. 14, and three order of magnitude in velocity, Fig. 15. For MEO and GTO, Analytic continuation achieved one order of magnitude more accuracy in position and velocity, Figs. 16 to 19. Finally, for HEO, Analytic Continuation is shown to be 3 -4 times more accurate as shown in Figs. 20 and 21 for position and velocity, respectively. We attribute the significant improvement in LEO results to the higher effect drag perturbation has on LEO versus the other orbits which take advantage of the more accurate STM produced by Analytic Continuation to produce better linear prediction results.

Conclusion
In this work, the State Transition Matrix of the perturbed two body problem is derived using the higher order Analytic Continuation technique. All the higher order partial derivatives are computed recursively by utilizing Leibniz rule and scalar variables transformations. The recursive relations shown for J 2 − J 6 and drag perturbations can be readily generalized for arbitrary order perturbations, [15]. Simulation results are presented for four different types of orbits: LEO, MEO, GTO and HEO, for up to 10,000 orbit periods and compared with ODE45, ODE113 and ODE87. For all cases, Analytic Continuation STM maintains machine precision that is enabled by the adaptation scheme combined with the ability of computing arbitrary higher order partials. The cases presented include unperturbed orbits where the RMS error for each element of the STM is computed and compared against Lagrange's F and G solution. For the gravity perturbed cases, the accumulation of truncation and round-off errors in ODE45, ODE113 and ODE87 is observed whereas Analytic Continuation maintained double-precision accuracy. Finally, with the introduction of drag, the STM computed via Analytic Continuation is shown to outperform ODE87 (in terms of accuracy) for all the test cases in linear prediction results.
The present method for computing the STM is simple, highly accurate and can be readily expanded to include higher order perturbations. The same recursions can also be expanded to compute the higher order state transition tensors. Areas that will be explored in future research. The partial derivatives of the higher order time derivatives in Eqs. 53 to 56 are recursively computed as,