Extended analytical formulae for the perturbed Keplerian motion under low-thrust acceleration and orbital perturbations

This paper presents a collection of analytical formulae that can be used in the long-term propagation of the motion of a spacecraft subject to low-thrust acceleration and orbital perturbations. The paper considers accelerations due to: a low-thrust profile following an inverse square law, gravity perturbations due to the central body gravity field and the third-body gravitational perturbation. The analytical formulae are expressed in terms of non-singular equinoctial elements. The formulae for the third-body gravitational perturbation have been obtained starting from equations for the third-body potential already available in the literature. However, the final analytical formulae for the variation of the equinoctial orbital elements are a novel derivation. The results are validated, for different orbital regimes, using high-precision numerical orbit propagators.


Introduction
This paper is an extension of the work presented in  and Zuiani and Vasile (2015) and in Di Carlo et al. (2017b). In these works, analytical formulae were derived for the motion of a spacecraft subject to constant tangential acceleration, constant acceleration in the radial-transverse-normal reference frame, constant acceleration in an inertial reference frame, and orbital perturbations due to the second-order zonal harmonic of the central body gravitational perturbation, J 2 . The analytical formulae were obtained using a first-order expansion in the perturbing acceleration.
The paper is structured as follows: the theoretical background for the development of the analytical solutions is presented in Sect. 2. Sections 3, 4 and 5 present the mathematical derivation of the equations relative to the perturbations and low-thrust acceleration included in this work. Finally, Sect. 6 discusses the validation of the analytical formulae against highprecision numerical orbit propagators. Section 7 concludes the paper.

First-order analytical solution in non-singular equinoctial parameters
In this work, the motion of the spacecraft is described in non-singular equinoctial elements, in order to avoid singularities when the eccentricity or inclination of the orbit is zero. The set of equinoctial elements, as defined by Broucke and Cefola (1972), is a, P 1 = e sin ( + ω) , P 2 = e cos ( + ω) , where L = + ω + θ is the true longitude. The perturbing acceleration is expressed in a radial-transverse-normal reference frame (RTN). Any perturbing acceleration to the Keplerian orbital motion (including low-thrust propulsion actions) is, therefore, expressed in the RTN frame as: The Gauss' planetary equations, expressed in terms of equinoctial elements, equinoctial elements, are, as presented in Battin (1987): (L) = 1 + e cos θ = 1 + P 1 sin L + P 2 cos L .
The sixth Gauss' equation for dL/dt is, as presented in Kechichian (1997): Under the assumption that the perturbing acceleration is small compared to the local central body's Keplerian gravitational acceleration,  obtained the following approximation of Eq. (6): The validity of the approximation in Eq. (7) is shown in Fig. 1. Figure 1 shows the first and second terms of Eq. (6), representing, respectively, the effect of the central body's Keplerian gravitational acceleration and the effect of the perturbations f N . The results are relative to the SSO orbit defined in Table 2 and are obtained considering the perturbations due to J 2 and to a low-thrust acceleration aligned along the N direction. These accelerations are expected to be the biggest among the ones considered in this work; therefore, any other acceleration will have a smaller contribution to the second term of Eq. (6). Figure 1 shows that the second term of Eq. (6) is more than 3 order of magnitudes smaller than the first term. Combining Eqs.
(3) and (7), the variations of the equinoctial elements with the true longitude can be expressed as: A first-order analytical solution to Eq. (8) can be generated with the method of perturbations, as presented in Vallado (2007). The idea at the basis of the method of perturbations is that small disturbing forces cause small deviations from the known solution to the unperturbed problem. The small perturbing forces can be associated with small parameters which characterise the magnitude of the disturbing forces. If one calls X the state of the spacecraft expressed in terms of equinoctial elements, X = [a, P 1 , P 2 , Q 1 , Q 2 , L] T , and the reduced vector with no true longitude X = [a, P 1 , P 2 , Q 1 , Q 2 ] T , Eqs. (8) can be rewritten in compact form as where = f . One can then look for the first-order approximation of the solution to Eq. where and X 0 represents the vector of initial conditions X 0 = [a 0 , P 10 , P 20 , Q 10 , Q 20 ] T . An analytical solutions of Eq. (11) was previously presented in Zuiani and Vasile (2015) and in Di Carlo et al. (2017b) for some orbital perturbations and low-thrust profiles. In particular, in Zuiani and Vasile (2015) it was shown that the analytical formulas in osculating elements can be used to efficiently propagate moderately long spirals composed of several tens of orbital revolutions. In this case, an analytical formula can be derived also for the variation of time with the true longitude L (see Eq. (7)). However, in this work the time equation is integrated with a quadrature method, as proposed in Zuiani and Vasile (2015) because it does not increase the computational time and provides accurate results. For longer spirals, an averaged propagation of the orbital elements is implemented; the variation of the equinoctial elements is, in this case, computed as andF We note that if one takes X 0 =X 0 , then integral (11) can be used to compute integral (14) analytically. According to Verhulst (1990) (see page 161 and following), the approximation is O( 2 ) and remains small on a time scale that is proportional to 1/ . While the integrals in Eq. (11) are computed analytically, the time integral in Eq. (12) is computed numerically; the resulting averaged propagator is, therefore, a semi-analytical method. In the remainder of this paper, we will focus on the averaged formulation only. By using this approach, the authors have derived analytical formulae for the following accelerations and orbital perturbations: 1. second zonal harmonic of the Earth's gravitational perturbation, J 2 (as in Zuiani and Vasile 2015); 2. third, fourth and fifth zonal harmonics of the Earth's gravitational perturbation, J 3 , J 4 , J 5 (Sect. 3); 3. atmospheric drag (as in Di Carlo et al. 2017b); 4. solar radiation pressure, including eclipses (as in Zuiani and Vasile 2015); 5. third-body gravitational perturbation (Sect. 4); 6. constant tangential acceleration (as in Zuiani and Vasile 2015); 7. constant acceleration in a radial-transverse-normal reference frame (as in Zuiani and Vasile 2015); 8. acceleration with constant direction in a radial-transverse-normal reference frame, and with magnitude of the acceleration proportional to 1/r 2 , where r is the distance from the central body (Sun) (Sect. 5); 9. constant acceleration in an inertial reference frame (as in Zuiani and Vasile 2015).
In the following sections, we will present only the development of the analytical formulae for the zonal harmonic J 3 , J 4 and J 5 , the third-body gravitational attraction and the inversesquare low-thrust acceleration. In the examples that follow, we will demonstrate the validity of our formulae at computing an averaged solution as in Eq. (12), including all the effects presented in this paper and in previous works.
3 Analytical formulae for the effects of J 3 , J 4 and J 5 The potential due to the zonal terms of the Earth's gravity field is given in Vallado (2007): In Eq. (15), r is the distance of the considered point from the centre of mass of the Earth, δ is its declination, G is the gravitational constant, m ⊕ is the mass of the Earth, R ⊕ its radius, and P l are the Legendre polynomials of order l in sin δ as reported in Abramowitz and Stegun (1972): The expression for the coordinate z, z = r sin δ, is used to substitute sin δ = z/r in Eq. (15). The perturbing acceleration due to J l can be computed from the gradient of the associated potential U J l as whereî R is the versor of the RTN reference frame andk is the z-component versor of the Earth Centred Inertial (ECI) reference frame, as presented in Vallado (2007). The components of the perturbing acceleration due to J l can be expressed, in the RTN reference frame, as: The scalar products in the previous equations arê where u is the argument of latitude. The analytical formulae for the variation of the equinoctial orbital elements under the effect of J 3 , J 4 and J 5 perturbations are derived in the following subsections.
Note that for Earth's orbits, the value of the zonal coefficients are such that the terms J 3 and J 4 are of order J 2 2 . However, this is not true for other celestial bodies. Table 1 reports the value of the zonal coefficients from J 2 to J 5 , for Earth, Moon, Mars and Venus. Table 1 shows that for Moon and Venus, J 5 is larger than J 2 2 , while for Mars J 4 is larger than J 2 2 .

Analytical formulae for the effect of J 3
The potential due to J 3 is where μ ⊕ = Gm ⊕ is the Earth's gravitational parameter and P 3 (x) is: The derivatives of the potential with respect to r and z are: Using Eqs. (18), the components of the perturbing acceleration due to J 3 are: The terms in i and u in Eqs. (24) can be expressed in terms of the equinoctial elements. Moreover, using it is possible to obtain the three components of the perturbing acceleration expressed in terms of equinoctial elements: where S = 1 + Q 2 1 + Q 2 2 . The accelerations in Eqs. (26) are substituted in Eq. (8). After integration, the following expression for X J 3 1 (Eq. 10) is obtained: The contribution of J 3 to Eq. (10) can be obtained using X J 3 1 and J 3 = J 3 R 3 ⊕ a −3 0 . The integrals I J 3 ,a,1 , I J 3 ,a,2 , I J 3 ,P 1 , I J 3 ,P 2 , I J 3 ,P 1 ,P 2 , I J 3 ,Q 1 and I J 3 ,Q 2 in Eqs. (27) are reported in "Appendix A". They are computed analytically using the software Wolfram Mathematica 1

Analytical formulae for the effect of J 4
The potential due to J 4 is expressed as: Following the same approach used for J 3 , the components of the perturbing accelerations are: The resulting expression for X J 4 1 is: The contribution of J 4 to Eq. (10) can be obtained using X J 4 1 and J 4 = J 4 R 4 ⊕ a −4 0 . The integrals in Eqs. (30) are reported in Appendix A.

Analytical formulae for the effect of J 5
The potential due to J 5 is: The components of the corresponding perturbing acceleration are: The resulting expression for X J 5 1 is: The contribution of J 5 to Eq. (10) can be obtained using X J 5 1 and J 5 = J 5 R 5 ⊕ a −5 0 . The integrals in Eqs. (33) are reported in Appendix A.

Third-body perturbations
Similarly to the aspherical gravitational potential of the previous section, the third-body potential is commonly written as a sum of Legendre polynomials written using the Legendre polynomials P n (Eq. 35), as in Vallado (2007). Our work here is similar to that of Cefola et al. (1974); however, instead of using the Lagrange VOP equations, we obtain the corresponding accelerations and introduce them in the Gauss VOP equations. Furthermore, instead of just obtaining averaged solutions, we also obtain osculating solutions.
In this section, n refers to the order of the polynomial. The quantities r andr refer to the norm and normalised vector of the position of the satellite relative to the Earth. The third-body disturbing potential can be written as where R P n refers to the disturbing third-body potential of Legendre polynomial order n and is given by: In Eq. (35), R 3 andR 3 are the norm and normalised vector of the position of the third-body relative to Earth. We can now obtain f P n , the acceleration vector caused by this perturbation, by taking its gradient: where P n is the derivative of the Legendre polynomial. The exact third-body acceleration is f 3 rd = ∞ n=2 f P n . We approximate to order n = 5. To obtain these accelerations in the RTN frame, they are first calculated in the equinoctial frame and then rotated. Therefore, in the following equations, the vectorr will be expressed in the equinoctial reference frame. The equinoctial reference frame is defined as having the x axis pointing towards the satellite when L = 0, the z direction along the direction of the angular momentum, and y, in order to complete the right-handed system, pointing towards the satellite when L = π/2. In this frame, the quantityr in Eq. (36) iŝ and we use the same direction cosines as in Cefola et al. (1974): Introducing the equations above into the formula for the acceleration will give us the acceleration in the equinoctial reference frame, which can be converted to RTN using where R z [−L] is the rotation matrix that rotates about the z axis of an angle L. Defining the function F as the acceleration can be written as: The values of the coefficients c and s are presented in Appendix B.1. These formulas were obtained using Wolfram Mathematica. Using the following properties of the function F, cos L F l (L; c mn , s mn ) = 1 2 F l (L; c m−1;n + c m+1;n , s m−1;n + s m+1;n ) , we can insert the formula for the acceleration in the Gauss' equations, multiply by dt/d L (Eq. (7)) and obtain the following equations: where 3 rd is given by and the coefficients are defined as, using c m as short hand for c mn : Because the third-body moves along its orbit, the coefficients c and s, which depend solely on the direction cosines, vary as well as the value of R 3 . In order to capture this effect, one would have to express the motion of the third-body as a function of L. Alternatively, one can split the calculation of integral (11) into segments and change the position of the third-body on each segment. In this paper, however, we use the simplifying assumption that the thirdbody moves slowly with respect to the orbital motion of the spacecraft. Thus, the averaging integral is computed keeping the position of the third-body constant. By applying integral (11) to the right hand side of Eqs. (43), one obtains the system where The integrals I c m,n+l and I s m,n+l are defined as and their formulas are presented in Appendix B.2. The third-body perturbation is obtained by using X 3 rd 1 and 3 rd in Eq. (10).

Low-thrust propulsion following an inverse square law
For solar electric propulsion (SEP) applications or solar sails, it is assumed that the acceleration of the propulsion system decreases as the distance from the Sun increases. This is modelled by considering the following inverse square law expression for the magnitude of the acceleration where˜ is a reference acceleration, delivered by the propulsion system at a distance from the Sun equal tor . The distance of the spacecraft from the central body, r , can be expressed as in Eq. (25). The magnitude of the acceleration is, therefore, and the acceleration vector can be expressed as: In Eq. (51), α and β are the in-plane azimuth and out-of-plane elevation angle of the acceleration vector in the RTN reference frame. In this paper, the notation is consistent with Zuiani, Zuiani and Vasile (2015), who measured the in-plane angle of the perturbing acceleration from the radial direction, rather than from the classic transverse direction. In order to avoid confusion with the notation traditionally used in the literature, in this paper the in-plane thrust angle will be therefore identified as α , where α = π/2 − α and α is the angle measured from the transverse direction (see Fig. 2). Due to the introduction of the angles α and β, in the case of low-thrust acceleration, Eqs. (9) and (11) are replaced, respectively, by and Substituting Eq. (51) into Eqs. (8) and integrating from L 0 to L with constant azimuth angle α and elevation angle β results in the following components for X 1 due to the low-thrust acceleration: − cos α cos β (sin L − sin L 0 ) + sin α cos β (− cos L + cos L 0 ) + sin α cos β + P 20 Q 20 sin β I s1 − Q 10 P 20 sin β I c1 + P 10 sin α cos β I 11 cos α cos β (cos L 0 − cos L) + sin α cos β (sin L − sin L 0 ) + sin α cos β + P 10 Q 10 sin β I c1 − Q 20 P 10 sin β I s1 + P 20 sin α cos β I 11 The contribution of the low-thrust acceleration to Eq. (10) can be obtained using X LT 1 and LT =˜ r 2 /μ. The expression for I 11 is reported in Zuiani and Vasile (2015). The integrals I s1 and I c1 are and If e 0 ≈ 0, the following non-singular expressions are used:

Validation against high-precision numerical propagators
The equations presented in the previous sections and in  and Di Carlo et al. (2017b) constitute the core of a semi-analytical propagation tool called CALYPSO (Computational AnaLYtical Propagator of Space Orbits). CALYPSO includes the orbital perturbations and low-thrust acceleration summarised in Sect. 2.
In this section, we validate the averaged semi-analytical propagation in CALYPSO by comparing its results against the direct numerical integration of the Gauss' equations (Eq. (3)) and against the NASA open-source software General Mission Analysis Tool (GMAT) 2 . As a further validation, we compared our results on the effects of J 2 against the J 2 propagator of the AGI software Systems Tool Kit (STK) 3 . While the averaged propagator CALYPSO provides the mean orbital elements of the orbit, GMAT gives as output a set of osculating orbital elements. In order to be able to compare the results, the initial osculating orbital elements, used to define the initial state of the orbit in GMAT, are converted into mean orbital elements. These are then used as initial conditions for the averaged propagator. The method used to convert from osculating to mean elements is described in Sect. 6.1.
The validation against GMAT accounted for the perturbations due to J 2 , J 3 , J 4 and J 5 , the third-body effect of Sun and Moon, and solar radiation pressure with eclipses. GMAT was used with the Runge-Kutta 89 integrator with default values for accuracy and step size. 4 The direct numerical integration of the Gauss' planetary equations was used to validate only the equations resulting from the low-thrust acceleration presented in Sect. 5. The numerical integration was performed using MATLAB ode113, a variable-step, variable-order Adams-Bashforth-Moulton predict-evaluate-correct-evaluate solver of order 13. In the following, therefore, the validation of the low-thrust acceleration is presented separately from the validation of the other analytical formulae. For the low-thrust acceleration, in fact, the aim is only to verify that the proposed analytical solutions are accurate. Therefore, a comparison with the numerical integration of the Gauss' equations is deemed sufficient. For the orbital perturbations, instead, other factors, such as orbital transformation or ephemeris computation, might influence the final results. The perturbation due to the atmospheric drag is not included in the comparison against GMAT, because the analytical solution introduced in Di Carlo et al. (2017b) would require extracting specific data on the computation of the air density in GMAT which could not be directly accessed at the time of the comparison. For more information on the validation of the propagation with atmospheric drag perturbation against numerical integration of the Gauss' equations, the interested reader is referred to Di Carlo et al. (2017b).

Conversion from osculating to mean elements
This section describes the method used to convert from osculating to mean orbital elements. Corresponding osculating and mean orbital elements are used as initial conditions for GMAT and CALYPSO, respectively.
The formula that relates the mean orbital elementsX(t) to the osculating orbital elements X(t) is, as in Vallado (2007): where η contains the short-periodic variations and is defined such that its average over one orbital period is zero. We write it as where X 0 andX 0 are the initial osculating and mean orbital elements. The quantities X(X 0 , t) andX(X 0 , t) are the osculating and mean elements, obtained from the integration of Eq. (9) and Eq. (13), respectively, using X 0 andX 0 as initial conditions. An iterative process is used to estimateX 0 . LetX (k) 0 be the k-th estimate of the initial mean elements. Since the estimatē X (k) 0 is not exact, the difference X(X 0 , t) −X(X (k) 0 , t) does not have a zero average over an orbital period. Thus, the short-periodic variation η is estimated such that its average over a period is zero by subtracting the average of this quantity, that is, The (k + 1)-th estimate of the initial mean elements is then obtained by solving Eq. (58) with t = t 0 forX(t 0 ) =X 0 , using Eq. (60): The process starts by using the initial osculating elements X 0 as first estimate forX (0) 0 . The iterative process stops when the difference between consecutive estimates is less than 10 −6 . This condition is always met for k ≤ 2. For the purpose of obtaining the initial conditions of CALYPSO, the integral in Eq. (61) and the osculating elements X(X 0 , t) are computed numerically.

Numerical test set-up
Seven different types of orbits are considered in this test set: two low earth orbits (LEO), sunsynchronous orbit (SSO), medium earth orbit (MEO), geostationary transfer orbit (GTO), geostationary equatorial orbit (GEO) and high elliptic orbit (HEO). The initial osculating orbital elements of each orbit are reported in Table 2. These orbits were selected because they are the orbital regions where Earth orbiting spacecraft operate, and because different orbit perturbations are preponderant in different orbital regimes.
A 1000 kg mass spacecraft with initial propagation date 21 March 2030, 00:00, is considered. The spacecraft reflectivity coefficient for the solar radiation pressure coefficient is set to C r = 1.3. The area-to-mass ratio for the solar radiation pressure is 0.1 m 2 /kg. A cylindrical shadow model is considered for the solar radiation pressure. For the low-thrust acceleration, the thrust is assumed to be 0.1 N atr = R ⊕ . The azimuth and elevation angles are set to α = 90 deg and β = 30 deg. The variations of the equinoctial elements, for propagations of one year starting from the initial orbits described in Table 2, are shown in the next subsections. The analytical equations presented in this paper were developed with the aim of being used in the context of the optimisation of low-thrust trajectories (as presented in Di Carlo et al. (2017a, b) and Di Carlo et al. (2017c)), rather than for long-term propagations on times of the order of hundreds of year. Therefore, a propagation time of one year is deemed sufficient for validation purposes. Table 3 reports the computational time, in seconds, required by CALYPSO and by the numerical integration of the osculating Gauss' equations (Eq. (3)). Both integrations were performed using ode45 on MATLAB R2018b, with absolute and relative tolerance equal to 10 −10 , on a system with Intel(R) Core(TM) i7-8700 CPU 3.20 GHz with 8GB RAM. Different sets of perturbations were considered, as reported in Table 3; the last column of the table gives the reference to the figure where the results of CALYPSO are presented. It is important to stress that the implementation of CALYPSO was not optimised, and future work will focus on refactoring it for greater efficiency. At the same time, the advantage of CALYPSO in terms of computational time is already evident, especially at low altitudes (LEO1, LEO2, SSO). At higher altitudes on circular orbits, the advantage is less evident. This is due to the fact that a lower number of integration steps are required by the osculating numerical integration for high-altitude circular orbits, thus bringing its computational cost closer to the one of CALYPSO.
The following sections present the comparison of the results of the semi-analytical propagator CALYPSO with those of GMAT and of the numerical integration of the equations of Gauss, for the orbits defined in Table 2.

Results for the LEO 1 Orbit Type
In order to show the effect of the different zonal harmonics on orbit LEO 1, Figures from 3 to 6 represent the results of the propagation considering J 2 only (Fig. 3), J 2 and J 3 (Fig. 4), J 2 , J 3 and J 4 (Fig. 5) and J 2 , J 3 , J 4 and J 5 (Fig. 6). It is evident that the biggest contributions come from J 2 and J 3 ; J 4 and J 5 introduce minor changes. Figure 7 shows the results of the propagation considering the Earth's zonal harmonics J 2 , J 3 , J 4 , J 5 , Sun and Moon gravitational perturbation and solar radiation pressure with eclipses; the results of the semi-analytical propagator (in blue) are in agreement with those of GMAT (in red), thus demonstrating the validity of the method presented in Sects. 3 and 4 and in Zuiani and Vasile (2015). Figure 8 presents the results of the propagation with low-thrust Table 3 Computational times required by CALYPSO and by the numerical integration of the osculating Gauss' equations (Eq. (3)), for the propagation of the orbits reported in   Fig. 7 Validation against GMAT for orbit LEO 1: J 2 to J 5 , Sun, Moon, SRP and eclipses acceleration following an inverse square law 1/r 2 (Sect. 5). In this case the comparison is against a numerical integration of the Gauss' equations. Figure 9 shows the results of the propagation considering the Earth's zonal harmonics J 2 , J 3 , J 4 , J 5 , Sun and Moon gravitational perturbation and solar radiation pressure with eclipses for orbit LEO 2. Note the long-term effects in the semi-major axis that were captured by splitting the integral (14) to properly model eclipses. The error in Q 1 is due to the fact that nutation and precession are not modelled in the semi-analytical averaged propagator CALYPSO; the absence of a model for the nutation and precession affects the acceleration of the Earth's gravitational perturbations. The Earth's gravitational accelerations presented in Sect. 3 are, in fact, expressed in an Earth-fixed reference frame and should be transformed into an inertial reference frame before the propagation, as presented in Montenbruck and Gill (2000). The transformation to the inertial reference frame would take into account Earth's nutation and precession. These phenomena are not currently modelled in CALYPSO, but will be the subject of future works. The effect caused by the lack of model for nutation and precession will be further discussed for the MEO orbit in Sect. 6.6. Figure 10 presents the results of the propagation with low-thrust acceleration following an inverse square law (Sect. 5).  Figure 11 shows the results of the propagation with Earth's zonal harmonic perturbations due to J 2 , J 3 , J 4 , J 5 , Sun and Moon gravitational perturbation and solar radiation pressure with eclipses for the SSO. Figure 12 presents the results of the propagation with low-thrust acceleration following an inverse square law 1/r 2 (Sect. 5).

Results for the MEO orbit type
This section presents the comparison of the results of the averaged propagator CALYPSO against those of GMAT and of the direct numerical integration of the equations of Gauss for a MEO (Table 2). In addition, in this section comparisons against STK are also presented. Figure 13 shows the results of a propagation with perturbations due to the Earth's zonal harmonic J 2 , J 3 , J 4 , J 5 , Sun and Moon gravitational perturbations and solar radiation pressure with eclipses. Results in Fig. 13 show a 0.076 % relative error in the equinoctial parameter Q 2 at the end of the propagation. This value was obtained by comparing the values of Q 2 provided by CALYPSO and by GMAT at the end of the propagation period. The error in Q 2 is due to the fact that nutation and precession are not modelled in the semi-analytical averaged propagator CALYPSO, as explained in Sect. 6.3.
To further analyse the effect of the absence of a model for nutation and precession, Fig. 14 shows that the same kind of error can be also seen when considering only the perturbation due to the second Earth's zonal harmonic, J 2 . However, this error is not present in Fig. 15, which shows a comparison of the results produced by the semi-analytical propagator CALYPSO Validation against GMAT for orbit LEO 2: J 2 to J 5 , Sun, Moon, SRP and eclipses against those of STK using the option "J 2 propagator" 5 . In this case, the results of the two propagators for Q 2 are in agreement.
In order to further demonstrate that the difference in Q 2 in Fig. 13 is due only to the way in which the Earth's gravitational acceleration is modelled, Figs. 16 and 17 show the comparison of the results obtained considering the third-body gravitational perturbations only (Sun and Moon) and solar radiation pressure. Neither the third-body perturbations nor the solar radiation pressure shows an error in Q 2 .
The comparison considering solar radiation pressure (SRP) only (Fig. 17) shows some differences, mostly in Q 1 . Further investigation has revealed that these differences are mostly due to the variation of the SRP with the distance from the sun (r ). In Fig. 18, an osculating numerical propagation which intentionally neglects this variation is introduced. It is clear that our averaged semi-analytical propagator follows this curve much more closely than GMAT. We have verified that if this numerical propagator is made to model the variation of the SRP with r , the curve does approximate GMAT's results.
Finally, Fig. 19 presents the results of the propagation with low-thrust acceleration changing with 1/r 2 (Sect. 5).    Figure 20 shows the results of a propagation with perturbations due to the Earth's zonal harmonics J 2 , J 3 , J 4 , J 5 , Sun and Moon gravitational perturbations and solar radiation pressure with eclipses for the GTO. Also in this case, as for the MEO test case (Sect. 6.6), results in Fig. 20 show an error in the equinoctial parameter Q 2 at the end of the propagation. This is due to the absence of the modelling of nutation and precession, and the variation of SRP with distance from the Sun, as explained in Sects. 6.3 and 6.6. The apparent shift of the semi-major axis of the averaged semi-analytical propagator with respect to the average of the GMAT solution is caused by the behaviour of the short-periodic oscillation, as shown in more details in Fig. 21. It is evident that the short-periodic oscillations spend most of the time near their minimum; therefore, the average will be close to the minimum as well. Finally, Fig. 22 presents the results of the propagation with low-thrust acceleration following the inverse square law 1/r 2 (Sect. 5). Figure 23 shows the results of a propagation with perturbations due to J 2 , J 3 , J 4 , J 5 , Sun and Moon and solar radiation pressure with eclipses for the GEO. Also in this case, as for the MEO and GTO test cases (Sects. 6.6 and 6.7), results in Fig. 23 show an error in the equinoctial parameter Q 2 at the end of the propagation. Figure 24 presents the results of the propagation with low-thrust acceleration following the inverse square law 1/r 2 (Sect. 5). Figure 25 shows the results of a propagation with perturbations due to J 2 , J 3 , J 4 , J 5 , Sun and Moon and solar radiation pressure with eclipses for the orbit HEO. Also in this case, as   for the MEO, GTO and GEO test cases, results in Fig. 25 show an error in the equinoctial parameter Q 2 at the end of the propagation. Figure 26 presents the results of the propagation with low-thrust acceleration changing as 1/r 2 (Sect. 5).

Conclusion
This paper has presented a collection of analytical formulae for the propagation of the motion of a spacecraft under the effect of a number of disturbing accelerations. The analytical formulae were derived from a first-order expansion in the magnitude of the perturbing acceleration. It was shown how the analytical formulae can be used to compute the average variation of the orbital elements. The average variation was then propagated numerically forward in time. The propagation of the average variation of the orbital elements was validated against highprecision numerical orbit propagators (including the NASA open-source software GMAT). It was shown that all types of orbits could be propagated over moderate lengths of time faster than with GMAT. All comparisons show a good agreement with the numerical propagation for all the perturbations and accelerations in LEO and SSO, while in MEO, GTO, GEO and HEO a small error is evident in the analytical formulae due to the absence of a model for nutation and precession and the assumption of a constant distance of the Sun. This will be addressed in future works, by periodically updating the orientation between reference frames during the propagation and adjusting for the actual distance from the Sun. The analytical formulae proposed in this paper allow for the direct solution of the averaging integral in closed form over a complete revolution. In doing so, it was assumed that the motion of the third body was contained and could be considered constant. This assumption limits the use of the averaged solution to cases in which the orbital motion is faster than the motion of the third body in the sky. Furthermore, effects induced by the motion of the third body are not captured. It was, however, suggested that one could partition the orbital period in segments and apply the analytical formulae on each segment with a different position of the third body. This procedure would allow one to compute the full averaging integral accounting for the motion of the third body, although at a higher computational expense.
Future works will also consider a higher-order expansion of the variation of the orbital parameters, the effect of the perturbation on the true longitude, the motion of the third-body during averaging and possibly the effect of the variation of the SRP force with the distance from the Sun. In addition, the implementation of the code will be optimised for increased efficiency.

Compliance with ethical standards
Conflict of Interest The authors declare that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Integrals
This appendix gives the expressions for the integrals used to compute the analytical equations for the motion of the satellite subject to the Earth's gravitational potential perturbations.

B.2 Integrals
The expressions for the average integrals, I c mn and I s mn , which were obtained using Mathematica, have a secular term and a short-periodic term. The short-periodic term is too large to be presented here, but the secular term is much more manageable. In Tables 4 and 5, we have the averaged integralsĪ c mn andĪ s mn defined in Eq. (93). The secular components, S c m,n and S s m,n , can be expressed using the averaged integrals according to Eq. (94). The complete results are collected in a Mathematica notebook that can be found at: https://doi. org/10.15129/0d41a34d-2558-4eaf-96a7-4bfd5942f42c Table 4 Values ofĪ c m,n . An auxiliary variable is used, V = P 2 1 − P 2 2 . These expressions and their LaTeX representation were obtained using Mathematica m/n 1 2 3 4 5 0 e 2 +2 2B 5 3e 2 +2 2B 7 3e 4 +24e 2 +8 8B 9 15e 4 +40e 2 +8 8B 11 Table 5 Values ofĪ s m,n . An auxiliary variable is used, V = P 2 1 − P 2 2 . These expressions and their LaTeX representation were obtained using Mathematica m/n 1 2 3 4 5 In this section, and in the notebook, P 1 , P 2 , e and B are used as short hand for P 10 , P 20 , e 0 and B 0 .Ī Finally, note that the average variation for the semi-major axis a is zero for each Legendre polynomial, that is,