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

Abstract

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.

1 Introduction

This paper is an extension of the work presented in Zuiani and Vasile (2012) 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.

In this paper, the analytical formulae presented in Zuiani and Vasile (2012, 2015) and Di Carlo et al. (2017b) are extended in two ways: we derive first-order analytical expressions for the zonal harmonics $J_3, J_4, J_5$ and the third-body gravitational perturbation, and we derive analytical formulae for the effect of an acceleration profile following an inverse square law. This acceleration profile is typically provided by solar electric propulsion systems or solar sails and will be referred to, in the rest of this paper, as low-thrust acceleration. Analytical formulae for the atmospheric drag were previously presented in Di Carlo et al. (2017b) and included the coupling with $J_2$.

One of the first works, in the literature, on the definition of analytical solutions for the orbital motion subject to Earth’s asymmetric gravitational perturbations was presented in Brouwer (1959). Brouwer proposed analytical solutions in canonical Delaunay variables for all inclinations but the critical value of 63.4 degrees. Contributions of zonal harmonics up to $J_5$ were considered, and equations were derived for the secular terms, the long periodic terms and the short periodic terms. In Kozai (1959), solutions were proposed in classical orbital elements for the perturbations of $J_2$, $J_3$ and $J_4$, using a different set of equations for the specific cases of small eccentricity and small inclination. Deprit (1981) proposed the elimination of the parallax, a technique by which, when the perturbation affecting a Keplerian motion is proportional to $r^{-n}$ with $n \ge 3$, a canonical transformation of Lie type will convert the system into a quasi Keplerian one. In this system, the perturbation is proportional to $r^{-2}$. The theory was originally formulated in Hill variables, but it was later shown in Lara et al. (2014) that it can be achieved in Delaunay variables Lara et al. (2014). Kechichian (2008) presented analytical expressions for the accelerations due to $J_3$ and $J_4$ expressed in equinoctial elements; he then derived analytical derivatives of the accelerations with respect to the equinoctial elements. No closed-form solution for the variation of the orbital elements was derived starting from the analytical equations presented by Kechichian. On the contrary, in this paper, closed-form solutions for the variation of the equinoctial elements are presented, for zonal harmonics up to $J_5$.

An expression for the gravitational potential of a third-body perturbing object was presented in Kaula (1962) as an expansion in the classical orbital elements, using Hansen functions. Cefola et al. (1974) presented equations for the third-body potential in non-singular equinoctial elements. In particular, the potential was expressed as the sum of terms dependent on Legendre polynomials of $\cos \psi$, where $\psi$ is the angle between the spacecraft vector and the third-body vector. The angle $\psi$ was expressed in terms of the true longitude and of two direction cosines of the third-body position vector relative to the equinoctial frame. Cefola provided equations for different terms of the expansion of the potential and for their derivatives with respect to the equinoctial elements, but he did not provide closed-form solutions for the variation of the orbital elements due to the third-body gravitational perturbations. An expansion in non-singular elements using Legendre polynomials was also presented in Giacaglia (1975). The mathematical formulation was, however, different from the one presented by Cefola.

Recent work by Bombardelli et al. (2011) derived analytical asymptotic equations for the motion of satellites under the effect of constant tangential acceleration in regularised orbital parameters.

This brief survey of the existing literature on the derivation of analytical formulas for the perturbed motion of a spacecraft shows that a number of solutions already exist in different parameterisations. However, to the author’s knowledge, analytical equations in non-singular equinoctial elements for the motion of a spacecraft subject to the forces considered in this paper are not present in the literature.

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 high-precision numerical orbit propagators. Section 7 concludes the paper.

2 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

$$\begin{aligned} \begin{aligned}&a, \\&P_1= e \sin \left( \Omega + \omega \right) , \\&P_2= e \cos \left( \Omega + \omega \right) , \\&Q_1 = \tan \frac{i}{2} \sin \Omega , \\&Q_2 = \tan \frac{i}{2} \cos \Omega , \\&L, \end{aligned} \end{aligned}$$

(1)

where $L = \Omega + \omega + \theta$ 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:

$$\begin{aligned} \mathbf {f} = \begin{bmatrix} f_R \\ f_T \\ f_N \\ \end{bmatrix} \ . \end{aligned}$$

(2)

The Gauss’ planetary equations, expressed in terms of equinoctial elements, equinoctial elements, are, as presented in Battin (1987):

$$\begin{aligned} \begin{aligned} \frac{da}{dt}&= \frac{2}{B} \sqrt{\frac{a^3}{\mu }} \left[ \left( P_2 \sin L - P_1 \cos L\right) f_R + \Phi (L) f_T \right] ,\\ \frac{dP_1}{dt}&= B \sqrt{\frac{a}{\mu }} \left[ -f_R \cos L + \left( \frac{P_1 + \sin L}{\Phi (L)} + \sin L\right) f_T - P_2 \frac{Q_1 \cos L - Q_2 \sin L}{\Phi (L)} f_N \right] ,\\ \frac{dP_2}{dt}&= B \sqrt{\frac{a}{\mu }} \left[ f_R \sin L + \left( \frac{P_2 + \cos L}{\Phi (L)} + \cos L\right) f_T + P_1 \frac{Q_1 \cos L - Q_2 \sin L}{\Phi (L)} f_N \right] ,\\ \frac{d Q_1}{dt}&= \frac{B}{2} \sqrt{\frac{a}{\mu }} \left( 1 + Q_1^2 + Q_2^2 \right) \frac{\sin L}{\Phi (L)} f_N,\\ \frac{d Q_2}{dt}&= \frac{B}{2} \sqrt{\frac{a}{\mu }} \left( 1 + Q_1^2 + Q_2^2 \right) \frac{\cos L}{\Phi (L)} f_N \ . \\ \end{aligned} \end{aligned}$$

(3)

In Eq. (3), B and $\Phi (L)$ are defined as:

$$\begin{aligned} B= & {} \sqrt{1-e^2} = \sqrt{1 - P_1^2 - P_2^2} , \end{aligned}$$

(4)

$$\begin{aligned} \Phi (L)= & {} 1 + e \cos \theta = 1 + P_1 \sin L + P_2 \cos L \ . \end{aligned}$$

(5)

The sixth Gauss’ equation for $\textit{dL}/\textit{dt}$ is, as presented in Kechichian (1997):

$$\begin{aligned} \frac{dL}{dt} = \sqrt{\frac{\mu }{a}}\frac{\Phi ^2}{B^3}-\sqrt{\frac{a^3}{\mu }}\frac{B}{\Phi }f_N \left( Q_1 \cos L - Q_2 \sin L\right) \ . \end{aligned}$$

(6)

Under the assumption that the perturbing acceleration is small compared to the local central body’s Keplerian gravitational acceleration, Zuiani et al. (2012) obtained the following approximation of Eq. (6):

$$\begin{aligned} \frac{dt}{dL} \approx \sqrt{\frac{a^3}{\mu }} \frac{B^3}{\Phi ^2(L)}\ . \end{aligned}$$

(7)

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.

Fig. 1

Terms of Eq. (6) for different values of L: SSO orbit (Table 2) subject to $J_2$ and to low-thrust acceleration directed along N

Combining Eqs. (3) and (7), the variations of the equinoctial elements with the true longitude can be expressed as:

$$\begin{aligned} \begin{aligned} \frac{da}{dL} =&\frac{2 a^3 B^2}{\mu } \left[ \frac{\left( P_2 \sin L - P_1 \cos L\right) }{\Phi ^2(L)} f_R + \frac{1}{\Phi (L)} f_T \right] , \\ \frac{d P_1}{dL} =&\frac{a^2 B^4}{ \mu } \left\{ \left[ -\frac{\cos L}{\Phi ^2(L)} f_R + \left( \frac{P_1 + \sin L}{\Phi ^3(L)} + \frac{\sin L}{ \Phi ^2(L)} \right) f_T \right] - P_2 \frac{Q_1 \cos L - Q_2 \sin L}{\Phi ^3 (L)} f_N \right\} ,\\ \frac{d P_2}{dL} =&\frac{a^2 B^4}{ \mu } \left\{ \left[ \frac{\sin L}{\Phi ^2(L)} f_R + \left( \frac{P_2 + \cos L}{\Phi ^3(L)} + \frac{\cos L}{ \Phi ^2(L)} \right) f_T \right] + P_1 \frac{Q_1 \cos L - Q_2 \sin L}{\Phi ^3 (L)} f_N \right\} ,\\ \frac{dQ_1}{dL} =&\frac{a^2 B^4}{2 \mu } \left( 1+ Q_1^2 + Q_2^2 \right) \frac{\sin L}{\Phi ^3(L)} f_N ,\\ \frac{dQ_2}{dL} =&\frac{a^2 B^4}{2 \mu } \left( 1+ Q_1^2 + Q_2^2 \right) \frac{\cos L}{\Phi ^3(L)} f_N \ .\\ \end{aligned} \end{aligned}$$

(8)

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 $\mathbf {X}$ the state of the spacecraft expressed in terms of equinoctial elements, $\mathbf {X} = [a, P_1, P_2, Q_1, Q_2, L]^T$, and the reduced vector with no true longitude $\widetilde{\mathbf {X}} = \left[ a, P_1, P_2, Q_1, Q_2 \right] ^T$, Eqs. (8) can be rewritten in compact form as

$$\begin{aligned} \frac{d\widetilde{\mathbf {X}}}{dL} = \epsilon F \left( \widetilde{\mathbf {X}}, L \right) \ , \end{aligned}$$

(9)

where $\epsilon = \Vert \mathbf {f} \Vert$. One can then look for the first-order approximation of the solution to Eq. (8) in the form

$$\begin{aligned} \widetilde{\mathbf {X}} \approx \widetilde{\mathbf {X}}_0 + \epsilon \widetilde{\mathbf {X}}_1 \left( \widetilde{\mathbf {X}}_0, L_0, L, F \right) \ , \end{aligned}$$

(10)

where

$$\begin{aligned} \widetilde{\mathbf {X}}_1 = \int _{L_0}^L F \left( \widetilde{\mathbf {X}}_0, \mathcal {L}\right) d \mathcal {L} \ , \end{aligned}$$

(11)

and $\widetilde{\mathbf {X}}_0$ represents the vector of initial conditions $\widetilde{\mathbf {X}}_0 = \left[ a_0, P_{10}, P_{20}, Q_{10}, Q_{20} \right] ^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

$$\begin{aligned} \bar{\mathbf {X}}(t) = \bar{\mathbf {X}}_0 + \int _{t_0}^t \dot{\bar{\mathbf {X}}}\left( \tau , \bar{\mathbf {X}}\left( \tau \right) \right) d\tau , \end{aligned}$$

(12)

where

$$\begin{aligned} \dot{\bar{\mathbf {X}}} = \epsilon \bar{F}(\bar{\mathbf {X}}_0), \end{aligned}$$

(13)

and

$$\begin{aligned} \bar{F}(\bar{\mathbf {X}}_0)=\frac{1}{2\pi }\int _0^{2\pi } F(\bar{\mathbf {X}}_0,\mathcal {L})d\mathcal {L} \ . \end{aligned}$$

(14)

We note that if one takes $\widetilde{\mathbf {X}}_0=\bar{\mathbf {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(\epsilon ^2)$ and remains small on a time scale that is proportional to $1/\epsilon$.

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_{\odot }^2$, where $r_{\odot }$ 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 inverse-square 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):

$$\begin{aligned} U = -\frac{G m_{\oplus }}{r} \sum _{l=2}^{\infty } J_l \left( \frac{R_{\oplus }}{r} \right) ^l \mathcal {P}_l \left( \sin \delta \right) \ . \end{aligned}$$

(15)

In Eq. (15), r is the distance of the considered point from the centre of mass of the Earth, $\delta$ is its declination, G is the gravitational constant, $m_{\oplus }$ is the mass of the Earth, $R_{\oplus }$ its radius, and $\mathcal {P}_l$ are the Legendre polynomials of order l in $\sin \delta$ as reported in Abramowitz and Stegun (1972):

$$\begin{aligned} \mathcal {P}_l(x) = \frac{1}{2^l l!} \frac{d^l}{dx^l}\left[ \left( x^2 - 1\right) ^l\right] \ . \end{aligned}$$

(16)

The expression for the coordinate z, $z = r \sin \delta$, is used to substitute $\sin \delta = 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

$$\begin{aligned} \mathbf {f}_{J_l} = \nabla U_{J_l} = \frac{\partial U_{J_l}}{\partial r} \hat{\mathbf {i}}_R + \frac{\partial U_{J_l}}{\partial z} \hat{\mathbf {k}} \ , \end{aligned}$$

(17)

where $\hat{\mathbf {i}}_R$ is the versor of the RTN reference frame and $\hat{\mathbf {k}}$ 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:

$$\begin{aligned} \begin{aligned} f_{{J_l}_R}&= \mathbf {f}_{J_l} \cdot \hat{\mathbf {i}}_R = \frac{\partial U_{J_l}}{\partial r} + \frac{\partial U_{J_l}}{\partial z} \hat{\mathbf {k}} \cdot \hat{\mathbf {i}}_R ,\\ f_{{J_l}_T}&= \mathbf {f}_{J_l} \cdot \hat{\mathbf {i}}_T = \frac{\partial U_{J_l}}{\partial z} \hat{\mathbf {k}} \cdot \hat{\mathbf {i}}_T ,\\ f_{{J_l}_N}&= \mathbf {f}_{J_l} \cdot \hat{\mathbf {i}}_N = \frac{\partial U_{J_l}}{\partial z} \hat{\mathbf {k}} \cdot \hat{\mathbf {i}}_N \ .\\ \end{aligned} \end{aligned}$$

(18)

The scalar products in the previous equations are

$$\begin{aligned} \begin{aligned} \hat{\mathbf {k}} \cdot \hat{\mathbf {i}}_R&= \sin i \sin u ,\\ \hat{\mathbf {k}} \cdot \hat{\mathbf {i}}_T&= \sin i \cos u ,\\ \hat{\mathbf {k}} \cdot \hat{\mathbf {i}}_N&= \cos i \ ,\\ \end{aligned} \end{aligned}$$

(19)

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$.

Table 1 Non-normalised values of $J_2, J_3, J_4$ and $J_5$ for Earth, Moon, Mars and Venus (Moon gravity model: LP165P, Mars gravity model: GMM-1 and Venus gravity model: Magellan 180U)

3.1 Analytical formulae for the effect of $J_3$

The potential due to $J_3$ is

$$\begin{aligned} U_{J_3} = - \frac{\mu _{\oplus }}{r} J_3 \left( \frac{R_{\oplus }}{r}\right) ^3 \mathcal {P}_3 \left( \sin \delta \right) = - \frac{\mu _{\oplus }}{r} J_3 \left( \frac{R_{\oplus }}{r} \right) ^3 \left( \frac{5}{2} \frac{z^3}{r^3} - \frac{3}{2} \frac{z}{r} \right) , \end{aligned}$$

(20)

where $\mu _{\oplus } = G m_{\oplus }$ is the Earth’s gravitational parameter and $P_3(x)$ is:

$$\begin{aligned} \mathcal {P}_3\left( x \right) = \frac{1}{2}\left( 5 x^3 - 3 x\right) \ . \end{aligned}$$

(21)

The derivatives of the potential with respect to r and z are:

$$\begin{aligned} \frac{d U_{J_3}}{d r}&=\frac{1}{2} \mu _{\oplus } J_3 \frac{R_{\oplus }^3}{r^5}\left( \frac{z^3}{r^3} - 15 \frac{z}{r} \right) , \end{aligned}$$

(22)

$$\begin{aligned} \frac{d U_{J_3}}{d z}&= \frac{3}{2} \mu _{\oplus } J_3 \frac{R_{\oplus }^3}{r^5} \left( 1 - 5 \frac{z^2}{r^2}\right) \ . \end{aligned}$$

(23)

Using Eqs. (18), the components of the perturbing acceleration due to $J_3$ are:

$$\begin{aligned} \begin{aligned} f_{{J_3}_R}&= 2 \mu _{\oplus } J_3 \frac{R_{\oplus }^3}{r^5} \sin i \sin u \left( 5 \sin ^2 i \sin ^2 u - 3 \right) ,\\ f_{{J_3}_T}&= \frac{3}{2} \mu _{\oplus } J_3 \frac{R_{\oplus }^3}{r^5} \sin i \cos u \left( 1 - 5 \frac{z^2}{r^2}\right) ,\\ f_{{J_3}_N}&= \frac{3}{2} \mu _{\oplus } J_3 \frac{R_{\oplus }^3}{r^5} \cos i \left( 1 - 5 \frac{z^2}{r^2} \right) \ . \end{aligned} \end{aligned}$$

(24)

The terms in i and u in Eqs. (24) can be expressed in terms of the equinoctial elements. Moreover, using

$$\begin{aligned} r = \frac{a (1-e^2)}{1 + e \cos \theta } = \frac{a B^2}{\Phi (L)} , \end{aligned}$$

(25)

it is possible to obtain the three components of the perturbing acceleration expressed in terms of equinoctial elements:

$$\begin{aligned} \begin{aligned} f_{{J_3}_R}&= 4 \mu _{\oplus } J_3 R_{\oplus }^3\frac{ \Phi ^5(L)}{a^5 B^{10}} \left( \frac{Q_2 \sin L - Q_1 \cos L}{S} \right) \left[ \frac{20 \left( Q_2 \sin L - Q_1 \cos L \right) ^2}{S^2} -3\right] ,\\ f_{{J_3}_T}&= 3 \mu _{\oplus } J_3 R_{\oplus }^3\frac{\Phi ^5(L)}{a^5 B^{10}} \left( \frac{Q_2 \cos L + Q_1 \sin L}{S} \right) \left[ 1-\frac{20 \left( Q_2 \sin L - Q_1 \cos L \right) ^2}{S^2} \right] , \\ f_{{J_3}_N}&= \frac{3}{2} \mu _{\oplus } J_3 R_{\oplus }^3 \frac{\Phi ^5(L)}{a^5 B^{10}} \frac{\left( 1 - Q_1^2 - Q_2^2 \right) }{S} \left[ 1 - \frac{20 \left( Q_2 \sin L - Q_1 \cos L\right) ^2}{S^2} \right] , \end{aligned} \end{aligned}$$

(26)

where $S=1 + Q_1^2 + Q_2^2$. The accelerations in Eqs. (26) are substituted in Eq. (8). After integration, the following expression for $\widetilde{\mathbf {X}}_1^{J_3}$ (Eq. 10) is obtained:

$$\begin{aligned} \widetilde{\mathbf {X}}_1^{J_3}= \left[ \begin{array}{ll} {a_0} \left( {B_0^8 S_0^3}\right) ^{-1} \left( 8 \ I_{J_3,a,1} + 6 \ I_{J_3,a,2} \right) \\ \left( {B_0^6 S_0^3}\right) ^{-1} \left[ I_{J_3,P_1} - \left( {3}/{2}\right) P_{20} \left( 1 - Q_{10}^2 - Q_{20}^2 \right) I_{J_3,P_1,P_2}\right] \\ \left( {B_0^6 S_0^3}\right) ^{-1} \left[ I_{J_3,P_2} + \left( {3}/{2}\right) P_{10} \left( 1 - Q_{10}^2 - Q_{20}^2 \right) I_{J_3,P_1,P_2}\right] \\ \left( {3}/{4}\right) \left( B_0^6 S_0^2\right) ^{-1} \left( 1 - Q_{10}^2 - Q_{20}^2\right) I_{J_3,Q_1}\\ \left( {3}/{4}\right) \left( B_0^6 S_0^2\right) ^{-1} \left( 1 - Q_{10}^2 - Q_{20}^2\right) I_{J_3,Q_2} \end{array}\right] \ . \end{aligned}$$

(27)

The contribution of $J_3$ to Eq. (10) can be obtained using $\widetilde{\mathbf {X}}_1^{J_3}$ and $\epsilon _{J_3} = J_3 R_{\oplus }^3 a_0^{-3}$. 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¹

3.2 Analytical formulae for the effect of $J_4$

The potential due to $J_4$ is expressed as:

$$\begin{aligned} U_{J_4} = \frac{\mu _{\oplus } J_4}{8} \frac{R_{\oplus }^4}{r^5} \left[ 35\frac{z^4}{r^4} - 30 \frac{z^2}{r^2} + 3 \right] \ . \end{aligned}$$

(28)

Following the same approach used for $J_3$, the components of the perturbing accelerations are:

$$\begin{aligned} \begin{aligned} f_{{J_4}_R} =&\frac{\mu _{\oplus } J_4 R_{\oplus }^4}{8 B^{12} a^6} \Phi (L)^6 \left[ \frac{2800}{S^4} \left( Q_2 \sin L - Q_1 \cos L \right) ^4 - \frac{600}{S^2} \left( Q_2 \sin L - Q_1 \cos L \right) ^2 + 15 \right] ,\\ f_{{J_4}_T} =&- \frac{\mu _{\oplus } J_4 R_{\oplus }^4}{S^2 B^{12} a^6} \Phi ^6(L) \left[ \frac{280}{S^2} \left( Q_2 \sin L - Q_1 \cos L \right) ^2 -30 \right] \left( Q_2 \cos L + Q_1 \sin L\right) \\&\left( Q_2 \sin L - Q_1 \cos L\right) ,\\ f_{{J_4}_N} =&- \frac{\mu _{\oplus } J_4 R_{\oplus }^4}{S^2 B^{12} a^6} \left( 1- Q_1^2 - Q_2^2\right) \Phi (L)^6 \left[ \frac{140}{S^2} \left( Q_2 \sin L - Q_1 \cos L \right) ^2 -15 \right] \\&\left( Q_2 \sin L - Q_1 \cos L\right) \ . \end{aligned} \end{aligned}$$

(29)

The resulting expression for $\widetilde{\mathbf {X}}_1^{J_4}$ is:

$$\begin{aligned} \widetilde{\mathbf {X}}_1^{J_4}= \left[ \begin{array}{ll} { a_0}\left( {B_0^{10} S_0^4}\right) ^{-1} \left( {{(1/4)} \ I_{J_4,a,1}} - 2 I_{J_4,a,2} \right) \\ - \left( {B_0^8 S_0^4 }\right) ^{-1} \left[ I_{J_4, P_1} +P_{20} (1-Q_{10}^2 - Q_ {20}^2) I_{J_4,P_1,P_2}\right] \\ \left( {B_0^8 S_0^4}\right) ^{-1} \left[ I_{J_4,P_2} + P_{10}(1-Q_{10}^2 - Q_{20}^2) I_{J_4,P_1,P_2}\right] \\ {(1/2)} \left( {2 B_0^8 S_0^3}\right) ^{-1} \left( 1 - Q_{10}^2 - Q_{20}^2 \right) I_{J_4,Q1}\\ {(1/2)} \left( {2 B_0^8 S_0^3}\right) ^{-1} \left( 1 - Q_{10}^2 - Q_{20}^2 \right) I_{J_4,Q2} \end{array}\right] \ . \end{aligned}$$

(30)

The contribution of $J_4$ to Eq. (10) can be obtained using $\widetilde{\mathbf {X}}_1^{J_4}$ and $\epsilon _{J_4} = J_4 R_{\oplus }^4 a_0^{-4}$. The integrals in Eqs. (30) are reported in Appendix A.

3.3 Analytical formulae for the effect of $J_5$

The potential due to $J_5$ is:

$$\begin{aligned} U_{J_5} = -\frac{\mu _{\oplus } J_5 R_{\oplus }^5}{8} \frac{z}{r^7} \left( 63 \frac{z^4}{r^4} - 70 \frac{z^2}{r^2} + 15 \right) \ . \end{aligned}$$

(31)

The components of the corresponding perturbing acceleration are:

$$\begin{aligned} f_{{J_5}_R} =&\frac{\mu _{\oplus } J_5 R_{\oplus }^5}{2 S B^{14} a^7} \Phi ^7(L) \left( Q_2 \sin L - Q_1 \cos L \right) \left[ \frac{3024}{S^4} \left( Q_2 \sin L - Q_1 \cos L\right) ^4 + \right. \nonumber \\&\left. - \frac{840}{S^2} \left( Q_2 \sin L - Q_1 \cos L \right) ^2 +45 \right] ,\nonumber \\ f_{{J_5}_T} =&- \frac{\mu _{\oplus } J_5 R_{\oplus }^5}{4 S B^{14} a^7} \Phi ^7(L) \left( Q_2 \cos L + Q_1 \sin L \right) \left[ \frac{5040}{S^4} \left( Q_2 \sin L - Q_1 \cos L\right) ^4 + \right. \nonumber \\&\left. - \frac{840}{S^2} \left( Q_2 \sin L - Q_1 \cos L \right) ^2 +15 \right] ,\nonumber \\ f_{{J_5}_N} =&- \frac{\mu _{\oplus } J_5 R_{\oplus }^5}{8 S B^{14} a^7} \Phi ^7(L) \left( 1- Q_1^2 - Q_2^2 \right) \left[ \frac{5040}{S^4} \left( Q_2 \sin L - Q_1 \cos L\right) ^4 + \right. \nonumber \\&\left. - \frac{840}{S^2} \left( Q_2 \sin L - Q_1 \cos L \right) ^2 +15 \right] .\nonumber \\ \end{aligned}$$

(32)

The resulting expression for $\widetilde{\mathbf {X}}_1^{J_5}$ is:

$$\begin{aligned} \widetilde{\mathbf {X}}_1^{J_5}= \left[ \begin{array}{ll} { a_0}\left( {S_0^5 B_0^{12}}\right) ^{-1} I_{J_5,a}\\ - {(1/2)} \left( {S_0^5 B_{0}^{10}}\right) ^{-1} \left[ I_{J_5,P_1} - {(1/4)} {P_{20}} \left( 1-Q_{10}^2 - Q_{20}^2\right) I_{J_5,P_1,P_2} \right] \\ - {(1/2)} \left( {S_0^5 B_{0}^{10}}\right) ^{-1} \left[ I_{J_5,P_2} + {(1/4)} {P_{10}} \left( 1-Q_{10}^2 - Q_{20}^2\right) I_{J_5,P_1,P_2}\right] \\ - \left( 1/16\right) \left( {B^{10} S_0^4}\right) ^{-1} \left( 1 - Q_{10}^2 - Q_{20}^2\right) I_{J_5,Q_1}\\ - \left( 1/16\right) \left( {B^{10} S_0^4}\right) ^{-1} \left( 1 - Q_{10}^2 - Q_{20}^2\right) I_{J_5,Q_2} \\ \end{array}\right] \ . \end{aligned}$$

(33)

The contribution of $J_5$ to Eq. (10) can be obtained using $\widetilde{\mathbf {X}}_1^{J_5}$ and $\epsilon _{J_5} = {J_5 R_{\oplus }^5}{a_0^{-5}}$. The integrals in Eqs. (33) are reported in Appendix A.

4 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 $\mathcal {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 and $\hat{\mathbf{r}}$ 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

$$\begin{aligned} R_{3^\mathrm{rd}} = \sum _{n=2}^\infty R_{\mathcal {P}_n} \ , \end{aligned}$$

(34)

where $R_{\mathcal {P}_n}$ refers to the disturbing third-body potential of Legendre polynomial order n and is given by:

$$\begin{aligned} R_{\mathcal {P}_n} = - \frac{\mu _3}{R_3} \sum _n \left( \frac{r}{R_3} \right) ^n \mathcal {P}_n( \hat{\mathbf{r}} \cdot \hat{\mathbf{R}}_3) \ . \end{aligned}$$

(35)

In Eq. (35), ${R}_3$ and $\hat{\mathbf{R}}_3$ are the norm and normalised vector of the position of the third-body relative to Earth. We can now obtain ${f_{\mathcal {P}_n}}$, the acceleration vector caused by this perturbation, by taking its gradient:

$$\begin{aligned} {\mathbf {f}_{\mathcal {P}_n}} = \frac{\mu _3}{R_3} \left( \frac{r}{R_3} \right) ^n \frac{1}{r} \left( n \hat{\mathbf{r}} \mathcal {P}_n( \hat{\mathbf{r}} \cdot \hat{\mathbf{R}}_3) + (\hat{\mathbf{R}}_3 - \hat{\mathbf{r}}(\hat{\mathbf{R}}_3 \cdot \hat{\mathbf{r}})) \mathcal {P}_n'( \hat{\mathbf{r}} \cdot \hat{\mathbf{R}}_3) \right) \ , \end{aligned}$$

(36)

where $\mathcal {P}_n'$ is the derivative of the Legendre polynomial. The exact third-body acceleration is $\mathbf {f}_{3^\mathrm{rd}} = \sum _{n=2}^\infty \mathbf {f}_{\mathcal {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 vector $\hat{\mathbf{r}}$ 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={\pi }/{2}$. In this frame, the quantity $\hat{\mathbf{r}}$ in Eq. (36) is

$$\begin{aligned} \hat{\mathbf{r}} = \left[ \cos L, \sin L, 0 \right] ^T \ , \end{aligned}$$

(37)

and we use the same direction cosines as in Cefola et al. (1974):

$$\begin{aligned} \begin{aligned} \hat{\mathbf{R}}_3&= \left[ \alpha _{3rd}, \beta _{3rd}, \gamma _{3rd} \right] ^T \ ,\\ \hat{\mathbf{r}} \cdot \hat{\mathbf{R}}_3&= \alpha _{3rd} \cos L + \beta _{3rd} \sin L \ ,\\ (\hat{\mathbf{R}}_3 - \hat{\mathbf{r}}(\hat{\mathbf{R}}_3 \cdot \hat{\mathbf{r}}))&= \left[ \begin{array}{ll} -\beta _{3rd} \cos L \sin L + \alpha _{3rd} \sin ^2 L \\ -\alpha _{3rd} \cos L \sin L + \beta _{3rd} \cos ^2 L \\ \gamma _{3rd} \end{array}\right] \quad .\\ \end{aligned} \end{aligned}$$

(38)

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

$$\begin{aligned} \mathbf {f}_{3^\mathrm{rd}}^{RTN} = R_z[-L] \mathbf {f}_{3^\mathrm{rd}}^\mathrm{equin} \ , \end{aligned}$$

(39)

where $R_z[-L]$ is the rotation matrix that rotates about the z axis of an angle L. Defining the function F as

$$\begin{aligned} F^{(n)}_{l}(L; c_{mn}, s_{mn}) = \frac{1}{\Phi ^{n+l}(L)} \sum _m c_{mn} \cos mL + s_{mn} \sin mL \ , \end{aligned}$$

(40)

the acceleration can be written as:

$$\begin{aligned} \begin{aligned} f^R_{3^{rd}}&= \sum _{n=2}^5 \frac{\mu _3 (a B^2)^{n-1}}{R_3^{n+1}} F^{(n)}_{-1}(L; c^R_{mn}, s^R_{mn}) \ , \\ f^T_{3^{rd}}&= \sum _{n=2}^5 \frac{\mu _3 (a B^2)^{n-1}}{R_3^{n+1}} F^{(n)}_{-1}(L; c^T_{mn}, s^T_{mn}) \ , \\ f^N_{3^{rd}}&= \sum _{n=2}^5 \frac{\mu _3 (a B^2)^{n-1}}{R_3^{n+1}} F^{(n)}_{-1}(L; c^N_{mn}, s^N_{mn}) \ . \\ \end{aligned} \end{aligned}$$

(41)

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,

$$\begin{aligned} \begin{aligned} \Phi ^k F_{l}(L; c_{mn}, s_{mn})&= F_{l-k}(L; c_{mn}, s_{mn}) \ , \\ \cos L ~ F_{l}(L; c_{mn}, s_{mn})&= \frac{1}{2} F_{l}(L; c_{m-1; n}+c_{m+1; n}, s_{m-1; n}+s_{m+1; n}) \ , \\ \sin L ~ F_{l}(L; c_{mn}, s_{mn})&= \frac{1}{2} F_{l}(L; s_{m+1; n}-s_{m-1; n}, c_{m-1; n}-c_{m+1; n}) \quad , \end{aligned} \end{aligned}$$

(42)

we can insert the formula for the acceleration in the Gauss’ equations, multiply by dt/dL (Eq. (7)) and obtain the following equations:

$$\begin{aligned} \begin{aligned} \frac{da}{dL}&= 2 \epsilon _{3^{rd}} a \sum _{n=2}^5 B^{2n} \frac{a^{n+1}}{a_0^2 R_3^{n-1}} F^{(n)}_{1}(L; c^{(1)}_{mn}, s^{(1)}_{mn}) \ ,\\ \frac{dP_1}{dL}&= \epsilon _{3^{rd}} \sum _{n=2}^5 B^{2n+4} \frac{a^{n+1}}{a_0^2 R_3^{n-1}} \left( F^{(n)}_{2}(L; c^{(2)}_{mn}, s^{(2)}_{mn}) + F^{(n)}_{1}(L; c^{(3)}_{mn}, s^{(3)}_{mn}) \right) \ , \\ \frac{dP_2}{dL}&= \epsilon _{3^{rd}} \sum _{n=2}^5 B^{2n+4} \frac{a^{n+1}}{a_0^2 R_3^{n-1}} \left( F^{(n)}_{2}(L; c^{(4)}_{mn}, s^{(4)}_{mn}) + F^{(n)}_{1}(L; c^{(5)}_{mn}, s^{(5)}_{mn}) \right) \ ,\\ \frac{d Q_1}{dL}&= \frac{1}{4} \epsilon _{3^{rd}} S \sum _{n=2}^5 B^{2n+4} \frac{a^{n+1}}{a_0^2 R_3^{n-1}} F^{(n)}_2(L; s^N_{m+1,n}-s^N_{m-1,n}, c^N_{m+1,n}+c^N_{m-1,n})\ , \\ \frac{d Q_2}{dL}&= \frac{1}{4} \epsilon _{3^{rd}} S \sum _{n=2}^5 B^{2n+4} \frac{a^{n+1}}{a_0^2 R_3^{n-1}} F^{(n)}_2(L; c^N_{m+1,n}+c^N_{m-1,n}, -s^N_{m+1,n}+s^N_{m-1,n}) \ , \\ \end{aligned} \end{aligned}$$

(43)

where $\epsilon _{3^{rd}}$ is given by

$$\begin{aligned} \epsilon _{3^{rd}} = \frac{\mu _3}{\mu } \frac{a_0^{2}}{R_3^{2}} \ , \end{aligned}$$

(44)

and the coefficients are defined as, using $c_m$ as short hand for $c_{mn}$:

$$\begin{aligned} c^{(1)}_{m}&= \frac{1}{2}\left( P_1(-c^R_{m+1} - c^R_{m-1}-s^T_{m-1}+s^T_{m+1})\right. \nonumber \\&\quad \left. + P_2(s^R_{m+1} - s^R_{m-1} + c^T_{m-1} + c^T_{m+1})\right) + c^T_m \ ,\nonumber \\ s^{(1)}_{m}&= \frac{1}{2}\left( P_1(s^R_{m+1} - s^R_{m-1} + c^T_{m-1} -c^T_{m+1}) \right. \nonumber \\&\left. P_2(-c^R_{m+1} + c^R_{m-1} + s^T_{m-1} + s^T_{m+1})\right) + s^T_m \ , \nonumber \\ c^{(2)}_{m}&= P_1c^T_m + \frac{1}{2}\left( s^T_{m+1}-s^T_{m-1} -P_2( +(c^N_{m+1}+c^N_{m-1})Q_1 - (s^N_{m+1}-s^N_{m-1})Q_2) \right) \ , \nonumber \\ s^{(2)}_{m}&= P_1s^T_m + \frac{1}{2}\left( c^T_{m+1}+c^T_{m-1} -P_2( -(c^N_{m+1}+c^N_{m-1})Q_2 - (s^N_{m+1}-s^N_{m-1})Q_1)\right) \ ,\nonumber \\ c^{(3)}_{m}&= \frac{1}{2}\left( -c^R_{m+1}-c^R_{m-1} + s^T_{m+1}-s^T_{m-1} \right) \ ,\nonumber \\ s^{(3)}_{m}&= \frac{1}{2}\left( +s^R_{m+1}-s^R_{m-1} + c^T_{m+1}+c^T_{m-1} \right) \ ,\nonumber \\ c^{(4)}_{m}&= P_2c^T_m + \frac{1}{2}\left( c^T_{m+1}+c^T_{m-1} +P_1( +(c^N_{m+1}+c^N_{m-1})Q_1 - (s^N_{m+1}-s^N_{m-1})Q_2) \right) \ ,\nonumber \\ s^{(4)}_{m}&= P_2s^T_m + \frac{1}{2}\left( -s^T_{m+1}+s^T_{m-1} +P_1( -(c^N_{m+1}+c^N_{m-1})Q_2 - (s^N_{m+1}-s^N_{m-1})Q_1)\right) \ , \nonumber \\ c^{(5)}_{m}&= \frac{1}{2}\left( s^R_{m+1}-s^R_{m-1} + c^T_{m+1} + c^T_{m-1} \right) \ ,\nonumber \\ s^{(5)}_{m}&= \frac{1}{2}\left( c^R_{m+1}+c^R_{m-1} - s^T_{m+1} + s^T_{m-1} \right) . \end{aligned}$$

(45)

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 third-body 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

$$\begin{aligned} \widetilde{\mathbf {X}}_1^{3^{rd}}= \left[ \begin{array}{ll} 2 a_0 \sum _{n=2}^5 B_0^{2n} \left( \frac{a_0}{R_3}\right) ^{n-1} G^{(n)}_{1}(L; c^{(1)}_{mn}, s^{(1)}_{mn}) \\ \sum _{n=2}^5 B_0^{2n+4} \left( \frac{a_0}{R_3}\right) ^{n-1} \left( G^{(n)}_{2}(L; c^{(2)}_{mn}, s^{(2)}_{mn}) + G^{(n)}_{1}(L; c^{(3)}_{mn}, s^{(3)}_{mn}) \right) \\ \sum _{n=2}^5 B_0^{2n+4} \left( \frac{a_0}{R_3}\right) ^{n-1} \left( G^{(n)}_{2}(L; c^{(4)}_{mn}, s^{(4)}_{mn}) + G^{(n)}_{1}(L; c^{(5)}_{mn}, s^{(5)}_{mn}) \right) \\ \frac{1}{4} S_0 \sum _{n=2}^5 B_0^{2n+4} \left( \frac{a_0}{R_3}\right) ^{n-1} G^{(n)}_2(L; s^N_{m+1}-s^N_{m-1}, c^N_{m+1}+c^N_{m-1})\\ \frac{1}{4} S_0 \sum _{n=2}^5 B_0^{2n+4}\left( \frac{a_0}{R_3}\right) ^{n-1} G^{(n)}_2(L; c^N_{m+1}+c^N_{m-1}, -s^N_{m+1}+s^N_{m-1}) \\ \end{array}\right] \ , \end{aligned}$$

(46)

where

$$\begin{aligned} G^{(n)}_{l}(L; c_{mn}, s_{mn}) = \sum _m c_{mn} I^c_{m, n+l} + s_{mn} I^s_{m, n+l} \ . \end{aligned}$$

(47)

The integrals $I^c_{m, n+l}$ and $I^s_{m, n+l}$ are defined as

$$\begin{aligned} \begin{aligned} I^c_{m, n+l}&= \int _{L_0}^L \frac{\cos m\mathcal {L}}{\Phi ^{n+l}(\mathcal {L})} d\mathcal {L} \ ,\\ I^s_{m, n+l}&= \int _{L_0}^L \frac{\sin m\mathcal {L}}{\Phi ^{n+l}(\mathcal {L})} d\mathcal {L} \ , \end{aligned} \end{aligned}$$

(48)

and their formulas are presented in Appendix B.2. The third-body perturbation is obtained by using $\widetilde{\mathbf {X}}_1^{3^{rd}}$ and $\epsilon _{3^{rd}}$ in Eq. (10).

5 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

$$\begin{aligned} \epsilon = \tilde{\epsilon } / (r/\tilde{r})^2 \ , \end{aligned}$$

(49)

where $\tilde{\epsilon }$ is a reference acceleration, delivered by the propulsion system at a distance from the Sun equal to $\tilde{r}$. The distance of the spacecraft from the central body, r, can be expressed as in Eq. (25). The magnitude of the acceleration is, therefore,

$$\begin{aligned} \tilde{\epsilon } \tilde{r}^2 \frac{\Phi (L)^2}{a^2 B^4} \ , \end{aligned}$$

(50)

and the acceleration vector can be expressed as:

$$\begin{aligned} \mathbf {f}_{LT} = \tilde{\epsilon } \tilde{r}^2 \frac{\Phi (L)^2}{a^2 B^4} \begin{bmatrix} \cos \beta \cos \alpha ' \\ \cos \beta \sin \alpha ' \\ \sin \beta \end{bmatrix} \ . \end{aligned}$$

(51)

In Eq. (51), $\alpha '$ and $\beta$ 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 $\alpha '$, where $\alpha ' = \pi /2 - \alpha$ and $\alpha$ is the angle measured from the transverse direction (see Fig. 2).

Fig. 2

Perturbing acceleration in the RTN reference frame and representation of the angles $\alpha$, $\alpha '$ and $\beta$

Due to the introduction of the angles $\alpha '$ and $\beta$, in the case of low-thrust acceleration, Eqs. (9) and (11) are replaced, respectively, by

$$\begin{aligned} \frac{d\widetilde{\mathbf {X}}}{dL} = \epsilon F \left( \widetilde{\mathbf {X}}, L, \alpha ', \beta \right) \ \end{aligned}$$

(52)

and

$$\begin{aligned} \widetilde{\mathbf {X}}_1 = \int _{L_0}^L F \left( \widetilde{\mathbf {X}}, \mathcal {L}, \alpha ', \beta \right) d \mathcal {L} \ . \end{aligned}$$

(53)

Substituting Eq. (51) into Eqs. (8) and integrating from $L_0$ to L with constant azimuth angle $\alpha '$ and elevation angle $\beta$ results in the following components for $\widetilde{\mathbf {X}}_1$ due to the low-thrust acceleration:

$$\begin{aligned} \widetilde{\mathbf {X}}_1^{LT} = \left[ \begin{array}{ll} \begin{array}{c} {2 a_0 }{ B_0^{-1}} \cos \beta \left[ \left( P_{20} \cos \alpha ' + P_{10} \sin \alpha ' \right) \left( \cos L_0 - \cos L \right) - \right. \\ \left. \left( P_{10} \cos \alpha ' - P_{20} \sin \alpha ' \right) \left( \sin L - \sin L_0 \right) + \sin \alpha ' \left( L - L_0 \right) \right] \ \\ \end{array} \\ \\ \begin{array}{c} \left[ -\cos \alpha ' \cos \beta \left( \sin L - \sin L_0 \right) + \sin \alpha ' \cos \beta \left( -\cos L + \cos L_0 \right) + \right. \\ \left. \left( \sin \alpha ' \cos \beta + P_{20} Q_{20} \sin \beta \right) I_{s1} - Q_{10} P_{20} \sin \beta I_{c1} + P_{10} \sin \alpha ' \cos \beta I_{11} \right] \ \\ \end{array} \\ \\ \begin{array}{c} \left[ \cos \alpha ' \cos \beta \left( \cos L_0 - \cos L \right) + \sin \alpha ' \cos \beta \left( \sin L - \sin L_0 \right) + \right. \\ \left. \left( \sin \alpha ' \cos \beta + P_{10} Q_{10} \sin \beta \right) I_{c1} - Q_{20} P_{10} \sin \beta I_{s1} + P_{20} \sin \alpha ' \cos \beta I_{11} \right] \\ \end{array} \\ \\ {(1/2)} \sin \beta S_0 I_{s1} \ \\ \\ {(1/2)} \sin \beta S_0 I_{c1} \ \\ \end{array}\right] \ . \end{aligned}$$

(54)

The contribution of the low-thrust acceleration to Eq. (10) can be obtained using $\widetilde{\mathbf {X}}_1^{LT}$ and $\epsilon _{LT} = \tilde{\epsilon } {\tilde{r}^2}/{ \mu }$. The expression for $I_{11}$ is reported in Zuiani and Vasile (2015). The integrals $I_{s1}$ and $I_{c1}$ are

$$\begin{aligned} I_{s1} = \frac{P_{10}}{e_0^2}\left( L - L_0 \right) - \frac{P_{20}}{e_0^2} \log \frac{\Phi (L)}{\Phi (L_0)} - \frac{P_{10}}{e_0^2} I_{11} \ , \end{aligned}$$

(55)

and

$$\begin{aligned} I_{c1} = \frac{P_{20}}{e_0^2}\left( L - L_0 \right) + \frac{P_{10}}{e_0^2} \log \frac{\Phi (L)}{\Phi (L_0)} - \frac{P_{20}}{e_0^2} I_{11} \ . \end{aligned}$$

(56)

If $e_0 \approx 0$, the following non-singular expressions are used:

$$\begin{aligned} I_{s1} = \cos L_0 - \cos L \ , I_{c1} = \sin L - \sin L_0 \ . \end{aligned}$$

(57)

6 Validation against high-precision numerical propagators

The equations presented in the previous sections and in Zuiani and Vasile (2012) 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)². 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)³. 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.⁴ 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).

6.1 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 elements $\bar{\mathbf{X}}(t)$ to the osculating orbital elements $\mathbf{X}(t)$ is, as in Vallado (2007):

$$\begin{aligned} \mathbf{X}(t) = \bar{\mathbf{X}}(t) + \eta (\bar{\mathbf{X}}(t), t) ~, \end{aligned}$$

(58)

where $\eta$ contains the short-periodic variations and is defined such that its average over one orbital period is zero. We write it as

$$\begin{aligned} \eta (t) = \mathbf{X}(\mathbf{X}_0, t) - \bar{\mathbf{X}}(\bar{\mathbf{X}}_0, t) ~, \end{aligned}$$

(59)

where $\mathbf{X}_0$ and $\bar{\mathbf{X}}_0$ are the initial osculating and mean orbital elements. The quantities $\mathbf{X}(\mathbf{X}_0, t)$ and $\bar{\mathbf{X}}(\bar{\mathbf{X}}_0, t)$ are the osculating and mean elements, obtained from the integration of Eq. (9) and Eq. (13), respectively, using $\mathbf{X}_0$ and $\bar{\mathbf{X}}_0$ as initial conditions. An iterative process is used to estimate $\bar{\mathbf{X}}_0$. Let $\bar{\mathbf{X}}_0^{(k)}$ be the k-th estimate of the initial mean elements. Since the estimate $\bar{\mathbf{X}}_0^{(k)}$ is not exact, the difference $\mathbf{X}(\mathbf{X}_0, t) - \bar{\mathbf{X}}(\bar{\mathbf{X}}_0^{(k)}, t)$ does not have a zero average over an orbital period. Thus, the short-periodic variation $\eta$ is estimated such that its average over a period is zero by subtracting the average of this quantity, that is,

$$\begin{aligned} \eta (t) \approx \mathbf{X}(\mathbf{X}_0, t) - \bar{\mathbf{X}}(\bar{\mathbf{X}}_0^{(k)}, t) - \frac{1}{T}\int _{t_0}^{t_0+T} \left[ \mathbf{X}(\mathbf{X}_0, t) - \bar{\mathbf{X}}(\bar{\mathbf{X}}_0^{(k)}, t) \right] dt ~. \end{aligned}$$

(60)

The $(k+1)$-th estimate of the initial mean elements is then obtained by solving Eq. (58) with $t=t_0$ for $\bar{\mathbf{X}}(t_0)=\bar{\mathbf{X}}_0$, using Eq. (60):

$$\begin{aligned} \bar{\mathbf{X}}_0^{(k+1)} = \bar{\mathbf{X}}_0^{(k)} + \frac{1}{T}\int _{t_0}^{t_0+T}\left[ \mathbf{X}(\mathbf{X}_0, t) - \bar{\mathbf{X}}(\bar{\mathbf{X}}_0^{(k)}, t)\right] dt ~. \end{aligned}$$

(61)

The process starts by using the initial osculating elements $\mathbf{X}_0$ as first estimate for $\bar{\mathbf{X}}_0^{(0)}$. The iterative process stops when the difference between consecutive estimates is less than $10^{-6}$. This condition is always met for $k\le 2$. For the purpose of obtaining the initial conditions of CALYPSO, the integral in Eq. (61) and the osculating elements $\mathbf{X}(\mathbf{X}_0, t)$ are computed numerically.

6.2 Numerical test set-up

Seven different types of orbits are considered in this test set: two low earth orbits (LEO), sun-synchronous 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.

Table 2 Initial osculating orbital element of the different orbits used for the validation of the semi-analytical propagator

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 $\text {m}^2/\text {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 at $\tilde{r} = R_{\oplus }$. The azimuth and elevation angles are set to $\alpha '=90$ deg and $\beta = 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, 2017b) 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.

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 Table 2, considering different perturbations

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.

6.3 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.

Fig. 3

Validation against GMAT for orbit LEO 1: $J_2$

Fig. 4

Validation against GMAT for orbit LEO 1: $J_2$ and $J_3$

Fig. 5

Validation against GMAT for orbit LEO 1: $J_2$, $J_3$ and $J_4$

Fig. 6

Validation against GMAT for orbit LEO 1: $J_2$, $J_3, J_4$ and $J_5$

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 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.

Fig. 7

Validation against GMAT for orbit LEO 1: $J_2$ to $J_5$, Sun, Moon, SRP and eclipses

Fig. 8

Validation against numerical propagation for orbit LEO 1: low-thrust acceleration (Sect. 5)

6.4 Results for the LEO 2 Orbit type

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).

Fig. 9

Validation against GMAT for orbit LEO 2: $J_2$ to $J_5$, Sun, Moon, SRP and eclipses

Fig. 10

Validation against numerical propagation for orbit LEO 2: low-thrust acceleration (Sect. 5)

6.5 Results for the SSO orbit type

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).

Fig. 11

Validation against GMAT for SSO: $J_2$ to $J_5$, Sun, Moon, SRP and eclipses

Fig. 12

Validation against numerical propagation for orbit SSO: low-thrust acceleration (Sect. 5)

6.6 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 against those of STK using the option “$J_2$ propagator”⁵. 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_{\odot }$). 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_{\odot }$, 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).

Fig. 13

Validation against GMAT for MEO: $J_2$ to $J_5$, Sun, Moon, SRP and eclipses

Fig. 14

Validation against GMAT for MEO: $J_2$

Fig. 15

Validation against STK “$J_2$ propagator” for MEO

Fig. 16

Validation against GMAT for MEO: Sun and Moon gravitational perturbations

Fig. 17

Validation against GMAT for MEO: SRP and eclipses

Fig. 18

Validation against GMAT and a numerical integrator for MEO: SRP and eclipses. The numerical integrator neglects the variation of the SRP force with the distance to the Sun

Fig. 19

Validation against numerical propagation for orbit MEO: low-thrust acceleration (Sect. 5)

6.7 Results for the GTO orbit type

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).

Fig. 20

Validation against GMAT for GTO: $J_2$ to $J_5$, Sun, Moon, SRP and eclipses

Fig. 21

Validation against GMAT for GTO: $J_2$ to $J_5$, Sun, Moon, SRP and eclipses. Zoom-in of the semi-major axis

Fig. 22

Validation against numerical propagation for orbit GTO: low-thrust acceleration (Sect. 5)

6.8 Results for the GEO orbit type

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).

Fig. 23

Validation against GMAT for GEO: $J_2$ to $J_5$, Sun, Moon, SRP and eclipses

Fig. 24

Validation against numerical propagation for GEO: low-thrust acceleration (Sect. 5)

6.9 Results for the HEO orbit type

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).

Fig. 25

Validation against GMAT for HEO: $J_2$ to $J_5$, Sun, Moon, SRP and eclipses

Fig. 26

Validation against numerical propagation for orbit HEO: low-thrust acceleration (Sect. 5)

7 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 high-precision 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.

Appendices

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. The integrals reported in the following subsections have been solved analytically using the software Mathematica; the results, too cumbersome to be reported here, have been exported directly from Mathematica to MATLAB. The complete results can be found at: https://doi.org/10.15129/0d41a34d-2558-4eaf-96a7-4bfd5942f42c.

1.1 Integrals for $J_3$

$$\begin{aligned} I_{J_3,a,1}= & {} \int _{L_0}^L \left\{ \Phi _0^3\left( \mathcal {L} \right) (P_{20} \sin \mathcal {L} - P_{10} \cos \mathcal {L}) \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) \right. \nonumber \\&\times \left. \left[ 20 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) ^2-3 S_0^2\right] \right\} d \mathcal {L} \end{aligned}$$

(62)

$$\begin{aligned} I_{J_3,a,2}= & {} \int _{L_0}^L \left\{ \Phi _0^4 \left( \mathcal {L} \right) \left( Q_{20} \cos \mathcal {L} + Q_{10} \sin \mathcal {L} \right) \right. \nonumber \\&\left. \times \left[ S_0^2- 20 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) ^2 \right] \right\} d\mathcal {L} \end{aligned}$$

(63)

$$\begin{aligned} I_{J_3,P_1}= & {} \int _{L_0}^L \left\{ -4 \cos (\mathcal {L}) \Phi _0^3(\mathcal {L}) \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) \right. \nonumber \\&\times \left. \left[ 20 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2-3 S_0^2\right] +\right. \nonumber \\&\left. +3\Phi _0^2(\mathcal {L}) \left( Q_{20} \cos \mathcal {L} + Q_{10} \sin \mathcal {L} \right) \left( P_{10} + \sin \mathcal {L} + \Phi _0(\mathcal {L}) \sin \mathcal {L} \right) \right. \nonumber \\&\times \left. \left[ S_0^2- 20 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2\right] \right\} d \mathcal {L} \end{aligned}$$

(64)

$$\begin{aligned} I_{J_3,P_1,P_2}= & {} \int _{L_0}^L \left\{ \left( Q_{10} \cos \mathcal {L} - Q_{20} \sin \mathcal {L} \right) \Phi _0^2(\mathcal {L})\right. \nonumber \\&\times \left. \left[ S_0^2- 20 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2\right] \right\} d\mathcal {L} \end{aligned}$$

(65)

$$\begin{aligned} I_{J_3,P_2}= & {} \int _{L_0}^L \left\{ 4 \sin (\mathcal {L}) \Phi _0^3(\mathcal {L}) \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) \right. \nonumber \\&\left. \left[ 20 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2-3 S_0^2\right] + \right. \nonumber \\&\left. 3 \Phi _0^2(\mathcal {L}) \left( Q_{20} \cos \mathcal {L} + Q_{10} \sin \mathcal {L} \right) \left( P_{20} + \cos \mathcal {L} + \Phi _0(\mathcal {L}) \cos \mathcal {L} \right) \right. \nonumber \\&\left. \left[ S_0^2- 20 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) ^2\right] \right\} d\mathcal {L} \end{aligned}$$

(66)

$$\begin{aligned} I_{J_3,Q_1}= & {} \int _{L_0}^L \sin (\mathcal {L}) \Phi _0^2(\mathcal {L}) \left[ S_0^2- 20 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) ^2 \right] d\mathcal {L} \end{aligned}$$

(67)

$$\begin{aligned} I_{J_3,Q_2}= & {} \int _{L_0}^L \cos (\mathcal {L}) \Phi _0^2(\mathcal {L}) \left[ S_0^2- 20 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2\right] d\mathcal {L} \end{aligned}$$

(68)

1.2 Integrals for $J_4$

$$\begin{aligned} I_{J_4,a,1}= & {} \int _{L_0}^L \left\{ \Phi _0^4(\mathcal {L}) \left[ 2800 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^4\right. \right. \nonumber \\&\left. \left. - 600 S_0^2 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 + 15 \right] \right. \nonumber \\&\times \left. \left( P_{20} \sin \mathcal {L} - P_{10} \cos \mathcal {L} \right) \right\} d\mathcal {L} \end{aligned}$$

(69)

$$\begin{aligned} I_{J_4,a,2}= & {} \int _{L_0}^L \left\{ \Phi _0^5(\mathcal {L}) \left( Q_{20} \cos \mathcal {L} + Q_{10} \sin \mathcal {L} \right) \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) \right. \nonumber \\&\times \left. \left[ 280 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 -30 S_0^2 \right] \right\} d\mathcal {L} \end{aligned}$$

(70)

$$\begin{aligned} I_{J_4,P_1}= & {} \int _{L_0}^{L} \left\{ \cos (\mathcal {L}) \frac{\Phi _0^4(\mathcal {L})}{8} \left[ 2800 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) ^4\right. \right. \nonumber \\&\left. \left. - 600 S_0^2 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 + \right. \right. \nonumber \\&\times \left. \left. 15 G_0^4\right] + \Phi _0^3(\mathcal {L}) \left( P_{10} + \sin \mathcal {L} + \Phi _0(\mathcal {L}) \sin \mathcal {L} \right) \left( Q_{20} \cos \mathcal {L} + Q_{10} \sin \mathcal {L}\right) \right. \nonumber \\&\times \left. \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) \left[ 280 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) ^2 - 30 S_0^2\right] \right\} d\mathcal {L} \end{aligned}$$

(71)

$$\begin{aligned} I_{J_4,P_1,P_2}= & {} \int _{L0}^L \Phi _0^3(\mathcal {L}) \left( Q_{10} \cos \mathcal {L} - Q_{20} \sin \mathcal {L}\right) ^2 \left[ 140 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) ^2 - 15 S_0^2\right] d\mathcal {L} \end{aligned}$$

(72)

$$\begin{aligned} I_{J_4,P_2}= & {} \int _{L_0}^{L} \left\{ \sin (\mathcal {L}) \frac{\Phi _0^4(\mathcal {L})}{8}\left[ 2800 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^4 - 600 S_0^2 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 + \right. \right. \nonumber \\&\times \left. \left. 15 S_0^4 \right] -\Phi _0^3(\mathcal {L}) \left( P_{20} + \cos \mathcal {L} + \Phi _0(\mathcal {L}) \cos \mathcal {L} \right) \left( Q_{20} \cos \mathcal {L} + Q_{10} \sin \mathcal {L}\right) \right. \nonumber \\&\times \left. \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) \left[ 280 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) ^2 - 30 S_0^2\right] \right\} d\mathcal {L} \end{aligned}$$

(73)

$$\begin{aligned} I_{J_4,Q_1}= & {} \int _{L_0}^L \Phi _0^3(\mathcal {L}) \sin \mathcal {L} \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) \left[ 15 S_0^2 - 140 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 \right] d\mathcal {L} \end{aligned}$$

(74)

$$\begin{aligned} I_{J_4,Q2}= & {} \int _{L_0}^L \Phi _0^3(\mathcal {L}) \cos \mathcal {L} \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) \left[ 15 S_0^2 -140 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 \right] d\mathcal {L} \end{aligned}$$

(75)

1.3 Integrals for $J_5$

$$\begin{aligned} I_{J_5,a}= & {} \int _{L_0}^L \left\{ \Phi _0^5(\mathcal {L}) \left[ 3024 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^4 - 840 S_0^2 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 +45 S_0^4 \right] \right. \nonumber \\&\left. \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) (P_{20} \sin \mathcal {L} - P_{10} \cos \mathcal {L})- \frac{1}{2}\Phi _0^6(\mathcal {L}) \left( Q_{20} \cos \mathcal {L} + Q_{10} \sin \mathcal {L} \right) ) \right. \nonumber \\&\left. \left[ 5040 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) ^4 - 840 S_0^2 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 +15 S_0^4 \right] \right\} d\mathcal {L} \end{aligned}$$

(76)

$$\begin{aligned} I_{J_5,P_1}= & {} \int _{L_0}^L \left\{ \Phi _0^5(\mathcal {L}) \cos L \left[ 3024 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^4 840 G_0^2 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 +45 G_0^4 \right] \right. \nonumber \\&\left. \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) + \frac{1}{2}\Phi _0^4(\mathcal {L}) \left( Q_{20} \cos \mathcal {L} + Q_{10} \sin \mathcal {L} \right) \left( P_{10} + \sin \mathcal {L} + \Phi _0(\mathcal {L}) \sin \mathcal {L}\right) \right. \nonumber \\&\left. \left[ 5040 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^4 + - 840 G_0^2 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 +15 G_0^4\right] \right\} d\mathcal {L} \end{aligned}$$

(77)

$$\begin{aligned} I_{J_5,P_1,P_2}= & {} \int _{L_0}^L \left\{ \Phi _0^4(\mathcal {L}) \left[ 5040 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^4 - 840 S_0^2 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 +15 S_0^4\right] \right. \nonumber \\&\left. \left( Q_{10} \cos \mathcal {L} - Q_{20} \sin \mathcal {L} \right) \right\} d\mathcal {L} \end{aligned}$$

(78)

$$\begin{aligned} I_{J_5,P_2}= & {} \int _{L_0}^L \left\{ - \Phi _0^5(\mathcal {L}) \sin \mathcal {L} \left[ 3024 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) ^4 + 840 S_0^2 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 +45 S_0^4 \right] \right. \nonumber \\&\left. \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) + \frac{1}{2}\Phi _0^4(\mathcal {L}) \left( Q_{20} \cos \mathcal {L} + Q_{10} \sin \mathcal {L} \right) \left( P_{20} + \cos \mathcal {L} + \Phi _0(\mathcal {L}) \cos \mathcal {L}\right) \right. \nonumber \\&\left. \left[ 5040 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) ^4 - 840 S_0^2 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 +15 S_0^4 \right] \right\} d\mathcal {L} \end{aligned}$$

(79)

$$\begin{aligned} I_{J_5,Q_1}= & {} \int _{L_0}^L \Phi _0^4(\mathcal {L}) \sin \mathcal {L} \left[ 5040 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) ^4 - 840 S_0^2 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 +15 S_0^4 \right] d \mathcal {L} \nonumber \\\end{aligned}$$

(80)

$$\begin{aligned} I_{J_5,Q_2}= & {} \int _{L_0}^L \Phi _0^4(\mathcal {L}) \cos \mathcal {L} \left[ 5040 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L}\right) ^4 - 840 S_0^2 \left( Q_{20} \sin \mathcal {L} - Q_{10} \cos \mathcal {L} \right) ^2 +15 S_0^4 \right] d \mathcal {L}\nonumber \\ \end{aligned}$$

(81)

Formulas for Third-Body Perturbation

1.1 Coefficients for acceleration

Equations (82-90) contain the values of $c_{nm}^{RTN}$ and $s_{nm}^{RTN}$. These coefficients are presented as matrices, where each matrix presents all the nonzero (and some zero) values of $c_{nm}^{RTN}$ or $s_{nm}^{RTN}$ for a fixed n. Each row corresponds to a fixed value of m, starting from $m=0$ for the top row. The columns are the radial (R), transverse (T) and normal (N) components, in this order. Each expression and its LaTeX representation were obtained using Mathematica. Missing coefficients are equal to zero. The direction cosines $\alpha _{3rd}, \beta _{3rd}, \gamma _{3rd}$ will be abbreviated $\alpha , \beta , \gamma$.

Notice that $c^R_{mn} = c^T_{mn} = 0$ for any values of m, n where $m>n$ or their sum $m+n$ is odd. Similarly, $c^N_{mn} = 0$ for any values of m, n where $m>n-1$ or their sum $m+n$ is even:

$$\begin{aligned} c^R_{mn}= & {} c^T_{mn} = 0 ~ \forall m,n\in \mathbb {N}: (m>n) \vee \left[ (m+n) ~\mathrm{odd}\right] \nonumber \\ c^N_{mn}= & {} 0 ~ \forall m,n\in \mathbb {N}: \left[ m>(n-1)\right] \vee \left[ (m+n) ~\mathrm{even}\right] \end{aligned}$$

(82)

$$\begin{aligned} c_2^{RTN}= & {} \left[ \begin{array}{ccc} \frac{1}{2} \left( 3 \alpha ^2+3 \beta ^2-2\right) &{} 0 &{} 0 \\ 0 &{} 0 &{} 3 \alpha \gamma \\ \frac{1}{2} \left( 3 \alpha ^2-3 \beta ^2\right) &{} 3 \alpha \beta &{} 0 \\ \end{array} \right] \end{aligned}$$

(83)

$$\begin{aligned} s_2^{RTN}= & {} \left[ \begin{array}{ccc} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 3 \beta \gamma \\ 3 \alpha \beta &{} -\frac{3}{2} \left( \alpha ^2-\beta ^2\right) &{} 0 \\ \end{array} \right] \end{aligned}$$

(84)

$$\begin{aligned} c_3^{RTN}= & {} \left[ \begin{array}{ccc} 0 &{} 0 &{} \frac{3}{4} \gamma \left( 5 \alpha ^2+5 \beta ^2-2\right) \\ \frac{9}{8} \alpha \left( 5 \alpha ^2+5 \beta ^2-4\right) &{} \frac{3}{8} \beta \left( 5 \alpha ^2+5 \beta ^2-4\right) &{} 0 \\ 0 &{} 0 &{} \frac{15}{4} \gamma (\alpha -\beta ) (\alpha +\beta ) \\ \frac{15}{8} \left( \alpha ^3-3 \alpha \beta ^2\right) &{} \frac{15}{8} \left( 3 \alpha ^2 \beta -\beta ^3\right) &{} 0 \\ \end{array} \right] \end{aligned}$$

(85)

$$\begin{aligned} s_3^{RTN}= & {} \left[ \begin{array}{ccc} 0 &{} 0 &{} 0 \\ \frac{9}{8} \beta \left( 5 \alpha ^2+5 \beta ^2-4\right) &{} -\frac{3}{8} \alpha \left( 5 \alpha ^2+5 \beta ^2-4\right) &{} 0 \\ 0 &{} 0 &{} \frac{15 \alpha \beta \gamma }{2} \\ \frac{15}{8} \left( 3 \alpha ^2 \beta -\beta ^3\right) &{} -\frac{15}{8} \left( \alpha ^3-3 \alpha \beta ^2\right) &{} 0 \\ \end{array} \right] \end{aligned}$$

(86)

$$\begin{aligned} c_4^{RTN}= & {} \left[ \begin{array}{ccc} \frac{3}{16} \left( 35 \alpha ^4+10 \alpha ^2 \left( 7 \beta ^2-4\right) +35 \beta ^4-40 \beta ^2+8\right) &{} 0 &{} 0 \\ 0 &{} 0 &{} \frac{15}{8} \alpha \gamma \left( 7 \alpha ^2+7 \beta ^2-4\right) \\ \frac{5}{4} (\alpha -\beta ) (\alpha +\beta ) \left( 7 \alpha ^2+7 \beta ^2-6\right) &{} \frac{5}{4} \alpha \beta \left( 7 \alpha ^2+7 \beta ^2-6\right) &{} 0 \\ 0 &{} 0 &{} \frac{35}{8} \alpha \gamma \left( \alpha ^2-3 \beta ^2\right) \\ \frac{35}{16} \left( \alpha ^4-6 \alpha ^2 \beta ^2+\beta ^4\right) &{} \frac{35}{4} \alpha \beta (\alpha -\beta ) (\alpha +\beta ) &{} 0 \\ \end{array} \right] \end{aligned}$$

(87)

$$\begin{aligned} s_4^{RTN}= & {} \left[ \begin{array}{ccc} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} \frac{15}{8} \beta \gamma \left( 7 \alpha ^2+7 \beta ^2-4\right) \\ \frac{5}{2} \alpha \beta \left( 7 \alpha ^2+7 \beta ^2-6\right) &{} -\frac{5}{8} (\alpha -\beta ) (\alpha +\beta ) \left( 7 \alpha ^2+7 \beta ^2-6\right) &{} 0 \\ 0 &{} 0 &{} -\frac{35}{8} \beta \gamma \left( \beta ^2-3 \alpha ^2\right) \\ \frac{35}{4} \alpha \beta (\alpha -\beta ) (\alpha +\beta ) &{} -\frac{35}{16} \left( \alpha ^4-6 \alpha ^2 \beta ^2+\beta ^4\right) &{} 0 \\ \end{array} \right] \end{aligned}$$

(88)

$$\begin{aligned} c_5^{RTN}= & {} \frac{15}{128} \left[ \begin{array}{ccc} 0 &{} 0 &{} 2 \gamma A_2 \\ 10 \alpha A_1 &{} 2 \beta A_1 &{} 0 \\ 0 &{} 0 &{} 56 \gamma (\alpha ^2 -\beta ^2 ) \left( 3 \alpha ^2+3 \beta ^2-2\right) \\ \frac{35}{3} \alpha \left( \alpha ^2-3 \beta ^2\right) \left( 9 \alpha ^2+9 \beta ^2-8\right) &{} -7 \beta \left( \beta ^2-3 \alpha ^2\right) \left( 9 \alpha ^2+9 \beta ^2-8\right) &{} 0 \\ 0 &{} 0 &{} 42 \gamma \left( \alpha ^4-6 \alpha ^2 \beta ^2+\beta ^4\right) \\ 21 \left( \alpha ^5-10 \alpha ^3 \beta ^2+5 \alpha \beta ^4\right) &{} 21 \left( 5 \alpha ^4 \beta -10 \alpha ^2 \beta ^3+\beta ^5\right) &{} 0 \\ \end{array} \right] \nonumber \\\end{aligned}$$

(89)

$$\begin{aligned} s_5^{RTN}= & {} \frac{15}{128} \left[ \begin{array}{ccc} 0 &{} 0 &{} 0 \\ 10 \beta A_1 &{} -2 \alpha A_1 &{} 0 \\ 0 &{} 0 &{} 112 \alpha \beta \gamma \left( 3 \alpha ^2+3 \beta ^2-2\right) \\ -\frac{35}{3} \beta \left( \beta ^2-3 \alpha ^2\right) \left( 9 \alpha ^2+9 \beta ^2-8\right) &{} -7 \alpha \left( \alpha ^2-3 \beta ^2\right) \left( 9 \alpha ^2+9 \beta ^2-8\right) &{} 0 \\ 0 &{} 0 &{} 168 \alpha \beta \gamma (\alpha -\beta ) (\alpha +\beta ) \\ 21 \left( 5 \alpha ^4 \beta -10 \alpha ^2 \beta ^3+\beta ^5\right) &{} -21 \left( \alpha ^5-10 \alpha ^3 \beta ^2+5 \alpha \beta ^4\right) &{} 0 \\ \end{array} \right] \nonumber \\ \end{aligned}$$

(90)

The following replacements were used to shorten the expressions:

$$\begin{aligned} A_1&= \left( 7 \left( 3 \alpha ^4+\alpha ^2 \left( 6 \beta ^2-4\right) +3 \beta ^4\right) -28 \beta ^2+8\right) \end{aligned}$$

(91)

$$\begin{aligned} A_2&= \left( 63 \alpha ^4+14 \alpha ^2 \left( 9 \beta ^2-4\right) +63 \beta ^4-56 \beta ^2+8\right) \end{aligned}$$

(92)

Table 4 Values of $\bar{I}^c_{m, n}$. An auxiliary variable is used, $V=P_1^2-P_2^2$. These expressions and their LaTeX representation were obtained using Mathematica

1.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 $\bar{I}^c_{mn}$ and $\bar{I}^s_{mn}$ defined in Eq. (93). The secular components, $\mathcal {S}^c_{m, n}$ and $\mathcal {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

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$.

$$\begin{aligned} \bar{I}^c_{m, n}= & {} \frac{1}{2\pi }\int _{L_0}^{L_0+2\pi } \frac{\cos m\mathcal {L}}{\Phi ^{n}(\mathcal {L})} d\mathcal {L} \nonumber \\ \bar{I}^s_{m, n}= & {} \frac{1}{2\pi }\int _{L_0}^{L_0+2\pi } \frac{\sin m\mathcal {L}}{\Phi ^{n}(\mathcal {L})} d\mathcal {L} \end{aligned}$$

(93)

$$\begin{aligned} \mathcal {S}^c_{m, n}= & {} 2 \bar{I}^c_{m, n}\left[ \Lambda (\mathcal {L}) \right] _{L_0}^{L} \nonumber \\ \mathcal {S}^s_{m, n}= & {} 2 \bar{I}^s_{m, n}\left[ \Lambda (\mathcal {L}) \right] _{L_0}^{L} \end{aligned}$$

(94)

where $\Lambda (\mathcal {L})$ is defined as

$$\begin{aligned} \Lambda (\mathcal {L}) = \arctan \left( \frac{P_1 + \left( 1-P_2\right) \tan \left( \frac{\mathcal {L}}{2}\right) }{B} \right) \end{aligned}$$

(95)

Table 5 Values of $\bar{I}^s_{m, n}$. An auxiliary variable is used, $V=P_1^2-P_2^2$. These expressions and their LaTeX representation were obtained using Mathematica

Finally, note that the average variation for the semi-major axis a is zero for each Legendre polynomial, that is,

$$\begin{aligned} G^{(n)}_{1}(L_0 + 2\pi ; c^{(1)}_{mn}, s^{(1)}_{mn}) - G^{(n)}_{1}(L_0; c^{(1)}_{mn}, s^{(1)}_{mn}) = 0 ~\forall n\in \mathbb {N} \end{aligned}$$

(96)