Extension of the King-Hele orbit contraction method for accurate, semi-analytical propagation of non-circular orbits

TL;DR

This paper extends the King-Hele semi-analytical orbit contraction method to better handle Earth's complex atmospheric density profiles, improving accuracy for non-circular orbit propagation while maintaining low computational costs.

Contribution

The authors introduce a novel atmosphere model derivative and superimposed exponential atmospheres to enhance the KH method's accuracy for realistic density profiles.

Findings

Extended KH method achieves accuracy comparable to numerical quadrature.

Method effectively propagates non-circular orbits with low computational load.

Refined series expansion improves the handling of complex atmospheric conditions.

Abstract

Numerical integration of orbit trajectories for a large number of initial conditions and for long time spans is computationally expensive. Semi-analytical methods were developed to reduce the computational burden. An elegant and widely used method of semi-analytically integrating trajectories of objects subject to atmospheric drag was proposed by King-Hele (KH). However, the analytical KH contraction method relies on the assumption that the atmosphere density decays strictly exponentially with altitude. If the actual density profile does not satisfy the assumption of a fixed scale height, as is the case for Earth's atmosphere, the KH method introduces potentially large errors for non-circular orbit configurations. In this work, the KH method is extended to account for such errors by using a newly introduced atmosphere model derivative. By superimposing exponentially decaying partial atmospheres, the superimposed KH method can be applied accurately while considering more complex density profiles. The KH method is further refined by deriving higher order terms during the series expansion. A variable boundary condition to choose the appropriate eccentricity regime, based on the series truncation errors, is introduced. The accuracy of the extended analytical contraction method is shown to be comparable to numerical Gauss-Legendre quadrature. Propagation using the proposed method compares well against non-averaged integration of the dynamics, while the computational load remains very low.

Tables

Table 1. Table 2: Relative density fitting errors ∀ h ∈ [ 100 , 2500 ] for-all ℎ 100 2500 \forall h\in[100,2500]  km.

${\eta }_{\rho }$	$T_{\infty }=750$ K	$T_{\infty }=1000$ K	$T_{\infty }=1250$ K
	$\forall h>239$ km	$\forall h>308$ km	$\forall h>306$ km
	$\forall h>134$ km	$\forall h>153$ km	$\forall h>154$ km
	$\forall h>119$ km	$\forall h>119$ km	$\forall h>130$ km
${\eta }_{\rho ,max}$	$1.6\%$ ($h=115$ km)	$1.8\%$ ($h=115$ km)	$1.9\%$ ($h=115$ km)

Table 2. Table 3: Smooth atmosphere model parameters resulting from a fit to the Jacchia-77 model, valid for altitudes h ∈ [ 100 , 2500 ] ℎ 100 2500 h\in[100,2500]  km.

	$T_{\infty }=750$ K		$T_{\infty }=1000$ K		$T_{\infty }=1250$ K
$p$	$H_p$	${\hat{\rho }}_p$	$H_p$	${\hat{\rho }}_p$	$H_p$	${\hat{\rho }}_p$
	[km]	[kg/m3]	[km]	[kg/m3]	[km]	[kg/m3]
$1$	$4.9948$	$2.4955\text{e}+02$	$4.9363$	$3.1632\text{e}+02$	$4.9027$	$3.6396\text{e}+02$
$2$	$10.471$	$8.4647\text{e}-04$	$11.046$	$5.2697\text{e}-04$	$11.437$	$3.8184\text{e}-04$
$3$	$21.613$	$9.1882\text{e}-07$	$24.850$	$3.7354\text{e}-07$	$25.567$	$2.8928\text{e}-07$
$4$	$37.805$	$1.2530\text{e}-08$	$46.462$	$1.0839\text{e}-08$	$44.916$	$1.2459\text{e}-08$
$5$	$49.967$	$1.3746\text{e}-09$	$64.435$	$1.0880\text{e}-09$	$76.080$	$9.2530\text{e}-10$
$6$	$174.23$	$1.5930\text{e}-13$	$147.46$	$3.8122\text{e}-13$	$111.09$	$1.6667\text{e}-11$
$7$	$315.15$	$1.1290\text{e}-14$	$314.53$	$4.8431\text{e}-14$	$354.23$	$5.9225\text{e}-14$
$8$	$1318.1$	$3.8065\text{e}-16$	$1214.6$	$4.2334\text{e}-16$	$892.19$	$1.7378\text{e}-15$

Table 3. Table 4: Comparison of 1081 1081 1081 propagations being subject to ρ 𝜌 \rho J or ρ 𝜌 \rho S , in total number of function evaluations, N f t ​ o ​ t superscript subscript 𝑁 𝑓 𝑡 𝑜 𝑡 N_{f}^{tot} , total integration evaluation time, t C ​ P ​ U t ​ o ​ t superscript subscript 𝑡 𝐶 𝑃 𝑈 𝑡 𝑜 𝑡 t_{CPU}^{tot} , and the minimum and maximum lifetime estimation error, η t L , m ​ i ​ n / m ​ a ​ x subscript 𝜂 subscript 𝑡 𝐿 𝑚 𝑖 𝑛 𝑚 𝑎 𝑥 \eta_{t_{L},min/max} .

$T_{\infty }$ [K]	$\rho$	$N_f^{tot}$ [$-$]	$t_{CPU}^{tot}$ [s]	${\eta }_{t_L,min/max}$ [$\%$]
$750$	${\rho }_J$	$593255$	$1086.9$
	${\rho }_S$	$592124$	$169.5$	$-0.060/0.051$
$1000$	${\rho }_J$	$568140$	$986.3$
	${\rho }_S$	$568140$	$153.6$	$-0.077/0.056$
$1250$	${\rho }_J$	$550021$	$789.9$
	${\rho }_S$	$549063$	$149.3$	$-0.074/0.048$

Table 4. Table 5: Performance of the different propagation and contraction methods, for a) t L = 30 subscript 𝑡 𝐿 30 t_{L}=30  days and b) t L = 360 subscript 𝑡 𝐿 360 t_{L}=360  days and various relative integration tolerances, γ r ​ e ​ l subscript 𝛾 𝑟 𝑒 𝑙 \gamma_{rel} . All figures are unit less.

	${ℳ}_1$	${ℳ}_2$	${\eta }_{t_L,50\%}$	${\eta }_{t_L,100\%}$	${ϵ}_{N_f}$	${ϵ}_{t_{CPU}}$
a)	SI-KH/${10}^{-6}$	SI-KH/${10}^{-12}$	$3.2\text{e}-6$	$8.4\text{e}-5$	$3.0\text{e}-1$	$2.9\text{e}-1$
	GL/${10}^{-6}$	GL/${10}^{-12}$	$3.3\text{e}-6$	$7.0\text{e}-5$	$3.7\text{e}-1$	$3.7\text{e}-1$
	NA/${10}^{-6}$	NA/${10}^{-12}$	$1.3\text{e}-3$	$2.5\text{e}-2$	$3.4\text{e}-1$	$3.4\text{e}-1$
	NA/${10}^{-9}$	NA/${10}^{-12}$	$1.6\text{e}-6$	$3.1\text{e}-5$	$6.0\text{e}-1$	$6.2\text{e}-1$
	SI-KH/${10}^{-6}$	NA/${10}^{-12}$	$8.7\text{e}-4$	$1.8\text{e}-3$	$1.1\text{e}-2$	$2.2\text{e}-2$
	GL/${10}^{-6}$	NA/${10}^{-12}$	$8.7\text{e}-4$	$1.7\text{e}-3$	$1.0\text{e}-2$	$3.6\text{e}-2$
b)	SI-KH/${10}^{-6}$	SI-KH/${10}^{-12}$	$3.2\text{e}-6$	$6.9\text{e}-5$	$3.1\text{e}-1$	$3.2\text{e}-1$
	GL/${10}^{-6}$	GL/${10}^{-12}$	$3.7\text{e}-6$	$6.9\text{e}-5$	$3.9\text{e}-1$	$4.0\text{e}-1$
	NA/${10}^{-6}$	NA/${10}^{-12}$	$1.6\text{e}-2$	$2.6\text{e}-1$	$3.4\text{e}-1$	$3.7\text{e}-1$
	NA/${10}^{-9}$	NA/${10}^{-12}$	$1.9\text{e}-5$	$4.1\text{e}-4$	$6.1\text{e}-1$	$6.4\text{e}-1$
	SI-KH/${10}^{-6}$	NA/${10}^{-12}$	$7.0\text{e}-5$	$3.2\text{e}-4$	$5.8\text{e}-4$	$1.1\text{e}-3$
	GL/${10}^{-6}$	NA/${10}^{-12}$	$7.2\text{e}-5$	$4.9\text{e}-4$	$5.8\text{e}-4$	$2.1\text{e}-3$

$E$	Eccentric anomaly [rad or deg]
$F$	$10.7$ cm solar flux [sfu]
$\overline{F}$	Smoothed $10.7$ cm solar flux [sfu]
$H$	Atmosphere density scale height [m or km]
$I_n$	Modified Bessel function of the first kind of order $n$ [$-$]
$P$	Orbit period [s]
$R$	Mean Earth radius [m or km]
$T$	Temperature [K]
$T_{\infty }$	Exospheric temperature [K]
$\Delta a$	Contraction over a full orbit period in a [m or km]
$\Delta e$	Contraction over a full orbit period in e [$-$]
$a$	Semi-major axis [m or km]
$e$	Eccentricity [$-$]
$e_b$	Boundary in e for selection of integral approximation method [$-$]
$h$	Height above Earth surface [m or km]
$h_a$	Apogee altitude [m or km]
$h_p$	Perigee altitude [m or km]
$n_p$	Number of partial atmospheres [$-$]
$r$	Radial distance from Earth’s center [m or km]
$t$	Time [seconds, days or years]
$t_L$	Lifetime [seconds, days or years]
$v$	Intertial velocity [m/s or km/s]
$z$	Auxilary variable for integration of the decay rate of highly eccentric orbits [$-$]
$\delta$	Inverse ballistic coefficient [m2/kg]
${\eta }_{\rho }$	Relative atmospheric density error [$-$ or $\%$]
${\eta }_{\Delta x}$	Relative integral approximation error [$-$ or $\%$]
${\eta }_{t_L}$	Relative lifetime error [$-$ or $\%$]
${\gamma }_{rel}$	Relative integration tolerance [$-$ or $\%$]
$\mu$	Earth gravitational parameter [m3/s2 or km3/s2]
$\rho$	Atmosphere density [kg/m3 or kg/km3]
$\hat{\rho }$	Atmosphere base density [kg/m3 or kg/km3]
$𝒞$	Contraction method
$ℳ$	Contraction and integration method with given ${\gamma }_{rel}$
$𝒪$	Order of series truncation error
J	Index corresponding to the Jacchia-77 atmosphere model
S	Index corresponding to the smooth atmosphere model
p	Index corresponding to the partial smooth atmosphere model
NS	Index corresponding to the non-smooth atmosphere model

Table 6. Table 7: Parameters to calculate 𝒂 𝒂 \boldsymbol{a} as a function of T ~ ∞ subscript ~ 𝑇 \tilde{T}_{\infty} . The factors are of unit [km -1 ].

$p$	$a_{p0}$	$a_{p1}$	$a_{p2}$
$1$	$-1.98541\text{e}-1$	$-1.40701\text{e}-2$	$1.87647\text{e}-2$
$2$	$-9.71648\text{e}-2$	$7.16062\text{e}-3$	$4.77822\text{e}-2$
$3$	$-5.05069\text{e}-2$	$3.33725\text{e}-2$	$-1.85987\text{e}-2$
$4$	$-2.83356\text{e}-2$	$1.64584\text{e}-2$	$-3.32683\text{e}-2$
$5$	$-2.18893\text{e}-2$	$8.84693\text{e}-3$	$5.46460\text{e}-2$
$6$	$-6.24488\text{e}-3$	$4.90041\text{e}-3$	$-6.03999\text{e}-3$
$7$	$-2.82771\text{e}-3$	$-3.17505\text{e}-3$	$1.93697\text{e}-3$
$8$	$-8.53512\text{e}-4$	$7.92640\text{e}-4$	$-1.24063\text{e}-3$
$p$	$a_{p3}$	$a_{p4}$	$a_{p5}$
$1$	$-1.72925\text{e}-2$	$2.77798\text{e}-2$	$-9.95750\text{e}-2$
$2$	$-1.51184\text{e}-1$	$3.51432\text{e}-1$	$-7.02642\text{e}-1$
$3$	$-1.03728\text{e}-1$	$5.51289\text{e}-1$	$-1.41638\text{e}+0$
$4$	$8.69501\text{e}-2$	$-6.20406\text{e}-2$	$-3.36952\text{e}-1$
$5$	$-2.34999\text{e}-1$	$5.47095\text{e}-1$	$-8.27779\text{e}-1$
$6$	$-7.24190\text{e}-2$	$5.32824\text{e}-1$	$-1.79828\text{e}+0$
$7$	$4.29619\text{e}-2$	$-1.78919\text{e}-1$	$3.53528\text{e}-1$
$8$	$4.65874\text{e}-3$	$-1.87465\text{e}-2$	$8.70408\text{e}-3$
$p$	$a_{p6}$	$a_{p7}$	$a_{p8}$
$1$	$1.76679\text{e}-1$	$-1.37542\text{e}-1$	$3.94618\text{e}-2$
$2$	$9.01640\text{e}-1$	$-6.03103\text{e}-1$	$1.59691\text{e}-1$
$3$	$1.87770\text{e}+0$	$-1.22379\text{e}+0$	$3.11852\text{e}-1$
$4$	$8.28293\text{e}-1$	$-6.99209\text{e}-1$	$2.06734\text{e}-1$
$5$	$7.76841\text{e}-1$	$-4.02671\text{e}-1$	$8.74533\text{e}-2$
$6$	$2.85818\text{e}+0$	$-2.11311\text{e}+0$	$5.91400\text{e}-1$
$7$	$-3.82857\text{e}-1$	$2.16923\text{e}-1$	$-5.02721\text{e}-2$
$8$	$3.62357\text{e}-2$	$-4.73838\text{e}-2$	$1.66805\text{e}-2$

Table 7. Table 8: Parameters to calculate 𝒃 𝒃 \boldsymbol{b} as a function of T ~ ∞ subscript ~ 𝑇 \tilde{T}_{\infty} . The factors are of unit [ln(kg/m 3 )].

$p$	$b_{p0}$	$b_{p1}$	$b_{p2}$
$1$	$5.35674\text{e}+0$	$1.36142\text{e}+0$	$-1.71993\text{e}+0$
$2$	$-6.96022\text{e}+0$	$-1.71534\text{e}-1$	$-6.26282\text{e}+0$
$3$	$-1.33334\text{e}+1$	$-4.29240\text{e}+0$	$1.12545\text{e}+0$
$4$	$-1.78792\text{e}+1$	$-2.89047\text{e}+0$	$3.93500\text{e}+0$
$5$	$-2.09320\text{e}+1$	$8.52674\text{e}+0$	$-5.08863\text{e}+1$
$6$	$-2.93700\text{e}+1$	$5.68339\text{e}-2$	$-2.61029\text{e}+1$
$7$	$-3.29807\text{e}+1$	$4.90080\text{e}+0$	$1.78391\text{e}+1$
$8$	$-3.51561\text{e}+1$	$-2.66659\text{e}+0$	$1.73783\text{e}+0$
$p$	$b_{p3}$	$b_{p4}$	$b_{p5}$
$1$	$1.48408\text{e}+0$	$-2.43815\text{e}+0$	$9.19988\text{e}+0$
$2$	$1.70218\text{e}+1$	$-3.66333\text{e}+1$	$7.26606\text{e}+1$
$3$	$1.41418\text{e}+1$	$-6.27283\text{e}+1$	$1.53398\text{e}+2$
$4$	$1.67754\text{e}+1$	$-1.15289\text{e}+2$	$3.24667\text{e}+2$
$5$	$1.56893\text{e}+2$	$-3.21951\text{e}+2$	$4.61948\text{e}+2$
$6$	$2.90804\text{e}+2$	$-1.47321\text{e}+3$	$3.87334\text{e}+3$
$7$	$-9.35850\text{e}+1$	$2.24591\text{e}+2$	$-3.60868\text{e}+2$
$8$	$-4.98942\text{e}+0$	$2.71676\text{e}+1$	$4.15537\text{e}+1$
$p$	$b_{p6}$	$b_{p7}$	$b_{p8}$
$1$	$-1.64492\text{e}+1$	$1.28147\text{e}+1$	$-3.67526\text{e}+0$
$2$	$-9.47544\text{e}+1$	$6.43396\text{e}+1$	$-1.72245\text{e}+1$
$3$	$-2.00134\text{e}+2$	$1.29740\text{e}+2$	$-3.30267\text{e}+1$
$4$	$-4.59063\text{e}+2$	$3.15704\text{e}+2$	$-8.42405\text{e}+1$
$5$	$-4.34126\text{e}+2$	$2.32404\text{e}+2$	$-5.27733\text{e}+1$
$6$	$-5.21125\text{e}+3$	$3.43718\text{e}+3$	$-8.85649\text{e}+2$
$7$	$3.73065\text{e}+2$	$-2.15221\text{e}+2$	$5.18052\text{e}+1$
$8$	$-1.88208\text{e}+2$	$1.86631\text{e}+2$	$-5.96266\text{e}+1$

Equations

$$H=-\frac{\rho }{\mathrm{d}\rho /\mathrm{d}h}$$

$$H_J(h)=-\frac{\rho (h)\Delta h}{\rho (h+\Delta h)-\rho (h)}\Delta h=1\text{ m}$$

$$T_{\infty }=5.48{\overline{F}}^{\frac{4}{5}}+101.8F^{\frac{2}{5}}$$

$${\rho }_{NS}(h)={\hat{\rho }}_i\exp \frac{h-h_i}{H_i}{h}_i<h<h_{i+1}$$

$$f_T=-\frac{1}{2}\rho v^2\delta$$

$$\text{derivative}at$$

$$\text{derivative}et$$

$$\text{derivative}Et$$

$$r$$

$$v$$

$$\Delta a$$

$$\Delta e$$

$$F_x=\text{derivative}xt\approx \frac{\Delta x}{P}x\in [a,e]$$

$$P=2\pi \sqrt{\frac{a^3}{\mu }}$$

$${\int }_0^{2\pi }f(E)dE\approx \pi \sum _i{w}_{i}f(E_i),E_i=(x_i+1)\pi$$

$$w_i=\frac{2}{(1-x_i^2)[P_n^{\prime }(x_i){]}^2}$$

$$\frac{1}{z}=\frac{H}{ae}$$

$$I_n(x)=\frac{1}{\pi }{\int }_0^{\pi }\exp (x\cos \theta )\cos (n\theta )d\theta$$

$${\rho }_S(h)=\sum _{p=1}^{n_p}{\rho }_p(h)=\sum _{p=1}^{n_p}{\hat{\rho }}_p{e}^{-h/H_p}$$

$$H_S(h)=-\frac{{\rho }_S(h)}{\mathrm{d}{\rho }_S/\mathrm{d}h}=\frac{{\sum }_{p=1}^{n_p}{\rho }_p(h)}{{\sum }_{p=1}^{n_p}{\rho }_p(h)/H_p}$$

$$h_i=\frac{h_0+h_1}{2}+\frac{h_1-h_0}{2}\cos \left(\frac{2i-1}{2N}\pi \right)i=1,\dots ,N$$

$$C=\sqrt{\frac{1}{N}\sum _{i=1}^N\ln {\left(\frac{{\rho }_S(h_i)}{{\rho }_J(h_i)}\right)}^2}$$

$${\eta }_{\rho }(h)=\frac{|{\rho }_S(h)-{\rho }_J(h)|}{{\rho }_J(h)}$$

$$a_p^i$$

$$b_p^i$$

$$a_p({\tilde{T}}_{\infty })$$

$$b_p({\tilde{T}}_{\infty })$$

$${\tilde{T}}_{\infty }=\frac{T_{\infty }-T_0}{T_1-T_0}$$

$$\mathbf{a}$$

$$\mathbf{b}$$

$$T_i=\frac{T_0+T_1}{2}+\frac{T_1-T_0}{2}\cos \left(\frac{2i-1}{2N}\pi \right)i=1,\dots ,N$$

$$H_p(T_{\infty })$$

$${\hat{\rho }}_p(T_{\infty })$$

$$\Delta a$$

Full text

In this work, the KH method is extended to account for such errors by using a newly introduced atmosphere model derivative. By superimposing exponentially decaying partial atmospheres, the superimposed KH method can be applied accurately while considering more complex density profiles. The KH method is further refined by deriving higher order terms during the series expansion. A variable boundary condition to choose the appropriate eccentricity regime, based on the series truncation errors, is introduced. The accuracy of the extended analytical contraction method is shown to be comparable to numerical Gauss-Legendre quadrature. Propagation using the proposed method compares well against non-averaged integration of the dynamics, while the computational load remains very low.

††journal: Advances in Space Research

1 Introduction

Numerical integration of the full orbital dynamics, including short-periodic variations, can be demanding from a computational point of view. For this reason, Semi-Analytical (SA)111The abbreviations used herein are, in alphabetical order; CIRA: COSPAR International Reference Atmosphere, CNES: Centre National d’Études Spatiales, COSPAR: Committee on Space Research, DTM: Drag Temperature Model, GL: Gauss-Legendre, KH: King-Hele, NA: Non-Averaged, NRLMSISE: Naval Research Laboratory Mass Spectrometer, Incoherent Scatter Radar Extended, SA: Semi-Analytical, SI-KH: SuperImposed King-Hele methods were developed to perform this task in a less demanding manner (e.g. Liu, 1974). Such methods remove the short-term periodic effects by averaging the variational equations, thereby reducing the stiffness of the problem. This is especially desired when orbits are to be propagated for many initial conditions and over long lifetimes, e.g. for estimating the future space debris environment.

The calculation of the orbit contraction – i.e. the reduction in semi-major axis and eccentricity – induced by atmospheric drag requires the integration of the atmosphere density along the orbit. Half a century ago, King-Hele (KH) derived analytical approximations to these integrals (King-Hele, 1964). Depending on the eccentricity of the orbit, e.g. circular, near-circular, low eccentric and highly eccentric, different series expansions were derived. Recommendations are given, found empirically, on when to use which formulation. Vinh et al. (1979) improved the theory by removing the ambiguity arising from the regions of validity in eccentricity and by applying the more mathematically rigorous Poincaré method for integration. The classical theory was adapted to non-singular elements, mitigating the problems that theories formulated in Keplerian elements have with vanishing eccentricities (Sharma, 1999; Xavier James Raj and Sharma, 2006).

The advantage of these methods is that the averaged contraction can be computed analytically using only a single density evaluation at the perigee. However, the analytical methods assume exponential decay of the atmosphere density above the perigee height. This fixed scale height assumption potentially introduces large errors, especially for highly eccentric orbits, if compared to propagation using quadrature.

Averaging methods based on quadrature solve the integral numerically. No assumption on the shape of the density profile is required, however, the density needs to be evaluated at many nodes along the orbit, slowing down the integration of the trajectory.

This work proposes modelling the atmosphere density by superimposing exponential functions, each with a fixed scale height. The KH formulation is then used for the calculation of the contraction of each individual component. As the assumption of a fixed scale height is satisfied for each component, the resulting decay rate is estimated with great accuracy. Finally, each individual contribution is summed up, resulting in the global contraction of the overall not strictly exponentially decaying atmosphere density. This superimposed approach is not limited to the KH method and can also be applied to the other analytical methods described above.

The proposed method is applied during propagation of different initial conditions from circular to highly elliptical orbits and compared against propagations using numerical quadrature of the contraction as well as against Non-Averaged (NA) integration. The smooth atmosphere derivative introduced here is independent of the underlying atmosphere model and can be extended to include time-variations, as shown here for the case of solar activity.

2 Background on Atmospheric Models

The atmosphere models discussed here can be divided into reference models and the derivatives thereof. The reference models commonly give the temperature, $T$222The nomenclature of all the variables can be found in A., and – more importantly for calculating the drag force – the density, $\rho$, of Earth’s atmosphere as a function of the altitude, $h$, and other input parameters. Examples are, in increasing degree of complexity, the COSPAR International Reference Atmosphere (CIRA), the Jacchia atmosphere (Jacchia, 1977), the Drag Temperature Model (DTM) (Bruinsma, 2015) and the Naval Research Laboratory Mass Spectrometer, Incoherent Scatter Radar Extended (NRLMSISE) model (Picone et al., 2002), all of which are (semi-)empirical models.

Of these reference models, derivatives can be obtained through fitting for two purposes: appropriate simplification of the mathematical formulation can lead to significant speed increases for a density evaluation; and adequate reformulation of the model improves the accuracy of analytical SA contraction methods as will become apparent in Section 3.3.

Sections 2.1 and 2.2 briefly introduce the Jacchia-77 atmosphere model and a derived non-smooth exponential atmosphere model.

2.1 Jacchia-77 Reference Atmosphere Model

The Jacchia-77 reference atmosphere (Jacchia, 1977) estimates the temperature and density profiles of the relevant atmospheric constituents as a function of the exospheric temperature, $T_{\infty}$. The density profile, $\rho$J, is based on the barometric equation and an empirically derived temperature profile in order to comply with observations of satellite decay. The static model is valid for altitudes $90<h<2500$ km and exospheric temperatures $500<T_{\infty}<2500$ K.

The computation of $\rho$J cannot be performed analytically and requires numerical integration for each of the 4 constituents, nitrogen, oxygen, argon and helium, plus integration of atomic nitrogen and oxygen. A fast, closed-form approximation is available (De Lafontaine and Hughes, 1983), but it was not considered here, as its modelled atmosphere does not purely decay exponentially.

The scale height, $H$, is defined as

[TABLE]

and numerically approximated for the Jacchia model scale height, $H$J, as

[TABLE]

Several thermospheric variations can be taken into account, such as solar cycle, solar activity, seasonal or daily variations. Generally, the objects of interest for SA propagation dwell on-orbit for several months to hundreds of years. Thus, only the variation with the $11$-year solar cycle is of interest here. The Jacchia reference uses the solar radio flux at $10.7$ cm, $F$, as an index for the solar activity (see Figure 1, source for data: Goddard Space Flight Center, 2018). From $F$, $T_{\infty}$ can be inferred as (Jacchia, 1977)

[TABLE]

where $\overline{F}$ is a smoothed $F$, commonly centred over an interval of several solar rotations. Jacchia recommended to use a smooth Gaussian mean based on weights which decay exponentially with time. Figure 1 shows the solar flux, $F$, and the Gaussian mean with a standard deviation of $\sigma=3$ solar rotations, i.e. 81 days, considering a window of $\pm 3\sigma$. More recent models such as NRLMSISE or DTM require $\overline{F}$ to be a moving mean of 3 solar rotations (ISO, 2013).

2.2 Non-Smooth Exponential Atmosphere Model

One very simple representation of the atmosphere density is using a piece-wise exponentially decaying model, by dividing the altitude range into bins. Each bin is defined by a lower altitude (base) and an upper altitude (base of the next bin), $h$i and $h$i+1, respectively, the base density, $\hat{\rho}$i, at $h$i and a scale height, $H$i, chosen such that the density is continuous over the limits of each bin. Then, within each altitude bin, the density, $\rho$NS, can be evaluated at each altitude $h$ as follows

[TABLE]

Such a model can be derived from any atmospheric model. Herein, the values given in Vallado (2013, Chapter 8.6) – fitting the CIRA-72 model at T_{\infty}$$=1000 K – are used for a comparison of models.

A problem with the non-smooth atmosphere model is that it is non-physical, with discontinuities in $H$. At each change of altitude bin, $H$ jumps from $H$i to $H$i+1. This non-smooth behaviour poses a problem to the (variable-step size) integrator, as the step size needs to be reduced to accurately describe the sudden change in contraction rate of the orbit. Thus, the number of function evaluations and the total time to propagate the orbit increases. An example is given in Figure 2b, comparing the number of steps required for propagation of an object subject to the non-smooth $\rho$NS to one using the smooth $\rho$J as a function of altitude. Evidently, each change of bin forces the integrator to reduce the step size.

The equally simple parametric model introduced in Section 4.1 does not suffer from these discontinuities.

3 Background on Semi-Analytical Orbit Contraction Methods

During SA propagation of an object trajectory subject to air-drag forces, the integrated change in the orbital element space, i.e. the contraction of the orbit, over a full revolution is of interest. This requires the integration of the (weighted) density along the orbit, which can either be done numerically using quadrature, or analytically.

Many quadrature rules exist (e.g. see Abramowitz and Stegun, 1972, p. 885–895) and they are independent of the underlying function, making them versatile. However, they require the evaluation of the density at multiple nodes along the orbit, increasing the computational load of the function evaluations during integration.

Analytical formulations, such as the one derived by D. King-Hele more than half a century ago (King-Hele, 1964) require the density to be evaluated only once per iteration in correspondence of the perigee altitude. Other examples of analytical formulations are the ones derived by Vinh et al. (1979), Sharma (1999) and Xavier James Raj and Sharma (2006). While offering improvements to the classical formulation of KH, such as being mathematically more rigorous and non-singular, they still suffer from the same assumption of a fixed scale height. The method proposed in Section 4 addresses this problem for any of the analytical formulations. For the sake of brevity, it is only applied to the KH method.

Section 3.1 introduces the system dynamics used throughout this work and discusses its averaging. Sections 3.2 and 3.3 introduce two averaging methods; the numerical Gauss-Legendre (GL) quadrature and the analytical KH method.

3.1 Dynamical System and Averaging

The main focus of this work is on correcting the errors arising from the fixed scale height assumption. Important effects of an oblate Earth, such as a non-spherical atmosphere or gravitational coupling (e.g. see  Brower and Hori, 1961), are not considered here. The superimposed approach does not replace the averaging method, rather it transforms one of the inputs, i.e. the atmosphere density, to fit its assumptions. Hence, it is also applicable to more elaborate theories.

The dynamical system used here is based on Lagrange’s planetary equations, given in Keplerian elements, stating the changes in the elements as a function of the applied forces from any small perturbations (see King-Hele, 1964, for more information). Only the tangential force induced by the aerodynamic drag is considered, i.e.

[TABLE]

with the density, $\rho$, the inertial velocity, $v$, and the effective area-to-mass ratio (i.e. the inverse of the ballistic coefficient), $\delta$, defined as $\delta=c_{D}A/m$, where $c_{D}$ is the drag coefficient, $A$ is the surface normal to $v$, and $m$ is the mass. Atmospheric rotation is ignored here, but could be taken into account by multiplying the right hand side of Equation 5 with the appropriate factor.

The variations of the semi-major axis, $a$, the eccentricity, $e$, and the eccentric anomaly, $E$, with respect to time, $t$, are

[TABLE]

with Earth’s gravitational parameter, $\mu$, the radius, $r$, and $v$ given as

[TABLE]

In order to reduce the stiffness of the problem, Equation 6 is averaged over a full orbit revolution, under the assumption that $a$ and $e$ remain constant. The resulting contractions, $\Delta a$ and $\Delta e$, for $a$ and $e$ respectively are

[TABLE]

with the altitude, $h=r-R$, given the mean Earth radius, $R$.

For SA propagation of the orbit, the derivatives of the variables with respect to time are approximated by the change over one revolution divided by the time required to cover the revolution

[TABLE]

with the orbit period, $P$, defined as

[TABLE]

3.2 Numerical Approximation

The integrals in Equation 8 can be approximated numerically using quadrature, e.g. GL quadrature (Abramowitz and Stegun, 1972, p. 887)

[TABLE]

where the node $x_{i}$ is the $i^{th}$ root of the Legendre Polynomial $P_{n}(x)$. The weights $w_{i}$ are given as

[TABLE]

and $P^{\prime}_{n}$ is the derivative of $P_{n}(x)$ with respect to $x$. The nodes and weights remain constant during the propagation, so they are calculated (or read from a table) only once upon initialisation. Routines to calculate ($x_{i}$, $w_{i}$) are available for various scientific programming tools, such as Matlab (MathWorks, 2018) and NumPy (Oliphant, 2006).

Advantages of a numerical approximation of the integrals in Equation 8 is that it can be found for any atmospheric model and that no series expansions are required. Disadvantages are the need of multiple density evaluations and the loss of an analytic formulation. E.g. the Jacobian cannot be inferred analytically, but requires another quadrature.

3.3 Classical King-Hele Approximation

Here, only a brief summary of the formulation is given. The treatment of the full theory behind the KH formulation can be found in King-Hele (1964). The integrals in Equation 8 can be approximated analytically by expanding the integrands as a power series in $e$ for low eccentric orbits, and in the inverted auxiliary variable, $z$

[TABLE]

for highly eccentric orbits, and cutting off at the appropriate degree.

With the assumption that the density, $\rho$, decreases strictly exponentially with altitude, i.e. with a fixed $H$, each expanded integrand can be represented by the modified Bessel function of the first kind, $I_{n}$, which for $n\in\mathbb{N}_{0}$ is given as (Abramowitz and Stegun, 1972, p. 376)

[TABLE]

In B, the resulting equations are given up to $5$th order, higher than the $2$nd order given originally by KH.

The KH formulation is fast as it can be evaluated analytically and requires only a single density evaluation for each computation of the contraction. The main problem with the fixed $H$ assumption is the underestimation of $\rho$ at altitudes above the perigee altitude, $h_{p}$, which for eccentric orbits can induce large errors. Figure 2a shows the trajectories of an object in an initially eccentric orbit with perigee and apogee height of $h_{p}\times h_{a}=750\times 2000$ km. They were propagated with two different atmosphere models, $\rho$NS and $\rho$J, and using two different contraction methods, GL quadrature and the KH formulation. For both atmosphere models, the KH method overestimates the density decay above perigee along the orbit, leading to an overestimation of the lifetime of up to $40\%$, compared to the propagation with the GL method. This is true – albeit sometimes less pronounced – for any object in a non-circular orbit subject to a non-strictly exponentially decaying atmosphere.

It has to be noted here that KH was aware of this problem and suggested a way to calculate the contraction of an orbit with a varying scale height (see King-Hele, 1964, Chapter 6). To keep the equations analytically integrable, he approximates the varying $H$ linearly, with a constant slope parameter. Linear approximation of the true $H$ is valid only locally. For low eccentric orbit configurations this might be sufficient, but high eccentricities will re-introduce the errors. Using a constant slope parameter will thus lead to a new over- or underestimation of the drag depending on $e$.

Another issue of the KH formulation is that it relies on series expansion. As the eccentricity grows, the formulation to calculate the contraction needs to switch from low to high eccentric orbits. This introduces discontinuities, at a classically fixed boundary eccentricty, $e_{b}$.

4 Proposed new model for the semi-analytical computation of the orbit contraction due to atmospheric drag

The proposed method of taking into account atmospheric drag for SA integration of trajectories consists of two parts: an atmosphere model based on constant scale heights, introduced in Section 4.1; and the extension of the KH formulation to reduce the errors induced by an atmosphere which in its sum does not decay exponentially, described in Section 4.2.

Table 1 shows an overview of how the proposed extension fits into the existing scheme of atmosphere models and SA orbit contraction methods. As mentioned earlier, the technique presented here is not limited to the KH method, but could be applied to any averaging method which is based on the fixed scale height assumption.

4.1 Smooth Exponential Atmosphere Model

The smooth atmosphere model proposed here does not in any way attempt to replace existing atmosphere density models. Instead, it is a derivation of those models. Nor is the idea of modelling the atmosphere as a sum of exponentials new: the Jacchia-77 reference model reduces – for each atmospheric constituent – to such a mathematical formulation if the vertical flux terms are neglected (Bass, 1980). The novelty of this work is the combination of the atmosphere model with the extended, superimposed KH formulation. Sections 4.1.1 and 4.1.2 introduce the static and variable atmosphere model, respectively.

4.1.1 Static Model

The smooth exponential atmosphere model, $\rho$S, is modelled by superimposing exponentials functions as

[TABLE]

where the number of partial atmospheres, $n_{p}$, the partial base densities, $\hat{\rho}$p, and the partial scale heights, $H$p, are fitting parameters. Note that the subscript $p$ does not stand for altitude bins, but for one of the partial atmospheres, each of which is valid for the whole altitude range. While it potentially could stand for a single atmosphere constituent, it is not restricted as such. The superimposed scale height, $H$S, is

[TABLE]

The derivative of $H$S with respect to $h$ is monotonically increasing, as $H$p is enforced to be larger than 0 for all $p$. Hence, the smooth atmosphere model can only be fitted to atmosphere models in altitude ranges where $\frac{\mathrm{d}{H}}{\mathrm{d}{h}}>0$. Above $h=100$ km, this is the case for $\rho$J for a wide range of $T_{\infty}$. Even if the underlying model shows slightly negative $H$ at the lower boundary $h_{0}$, a partial atmosphere with a small positive $H$p can still be fitted accurately.

To find the parameters, $H$p and $\hat{\rho}$p, the model in Equation 15 is fitted to $\rho$J for three different $T_{\infty}$: in accordance to a low solar activity, $T_{\infty}=750$ K; mean solar activity, $T_{\infty}=1000$ K; and high solar activity, $T_{\infty}=1250$ K (see Figure 1). The fit is performed in the logarithmic space as not to neglect lower densities at higher altitudes, using least squares minimisation at heights between $h_{0}=100$ km and the upper boundary, $h_{1}=2500$ km. To put more weights on the edges of the fit interval, the densities are evaluated at $N=100$ heights, $h_{i}$, distributed as Chebyshev nodes (Abramowitz and Stegun, 1972, p. 889)

[TABLE]

The number of partial atmospheres, $n_{p}$, is chosen to be $8$, as the cost function

[TABLE]

which is the root mean square of the logarithmic density fit residuals, stops improving (see Figure 3). For $T_{\infty}\in[750,1000,1250]$ K, the relative error, $\eta_{\rho}$, calculated as

[TABLE]

always remains below $0.1\%$ and $1\%$ for all $h>308$ km and $h>130$ km, respectively, and the maximum relative error, $\eta_{\rho}$,max, does not exceed $2\%$, as can be seen in Table 2. Hence, the density fit accurately represents the underlying model. The model parameters can be found in Table 3. Figure 4 shows a comparison between the underlying and fitted model, for T_{\infty}$$=1000 K.

A speed test for $2401$ density and scale height evaluations over the range $100\leq$ $h$ $\leq 2500$ km shows a near $60$-fold decrease in evaluation time for $\rho$S compared to $\rho$J. The implementation of the Jacchia-77 model used herein is written in the coding language C (taken from Instituto Nacional De Pesquisas Espaciais, 2018), and called from within Matlab, while the routine to calculate $\rho$S is implemented and called directly in Matlab. Thus, a further decrease of computational time could be expected if also the latter was implemented in C. The speed tests were performed using the same processor architecture.

4.1.2 Variable Model

Possible extensions to the smooth exponential atmosphere model are the inclusion of a temporal dependence, such as the solar cycle, annual or daily variations. Here, the model is extended to incorporate the variability in the atmosphere density due to a variable $T_{\infty}$. To conserve the mathematical formulation of the static model, the temperature dependence is introduced in the fitting parameters, $\hat{\rho}_{p}=\hat{\rho}_{p}(T_{\infty})$ and $H_{p}=H_{p}(T_{\infty})$.

$T_{\infty}$ is a function of the solar proxy $F$ (see Equation 3), so the fitting range is defined by $F$. Generally, the long-term predictions for $F$ – based on various numbers of previous solar cycles – remain between $F\in[60,230]$ sfu (Vallado and Finkleman, 2014; Dolado-Perez et al., 2015; Radtke and Stoll, 2016). This translates into $T_{\infty}\in[669,1321]$ K, as $\overline{F}$ per definition remains in the same range as $F$. The parameters for the variable smooth exponential atmosphere model derived below, and listed in  D, are valid for any $T_{\infty}\in[T_{0}=650,T_{1}=1350]$ K. They should not be used for $T_{\infty}$ outside this range, as polynomial fits tend to oscillate strongly outside the fitting interval.

The dependence on $T_{\infty}$ is incorporated using a polynomial least squares fit. Each partial atmosphere is fitted separately. The static parameters, fitted to the $i=1,2,\dots,M$ static atmospheres with different $T_{\infty}$, are converted

[TABLE]

and each time-variable partial atmosphere is fitted to two independent polynomials of order $l$ and $m$ respectively

[TABLE]

using a normalised and unit-less $\tilde{T}_{\infty}$, defined as

[TABLE]

In vector notation, Equation 21 can be written as

[TABLE]

To prevent over-fitting, the order of the polynomials should remain well below the number of fitted static atmospheres. Here, the model in Equation 21 is fitted to $M=50$ statically fitted models, distributed again as Chebyshev nodes between $T_{0}$ and $T_{1}$

[TABLE]

The orders are chosen to be $l=m=8$ such that the error remains below $0.5\%$ for all $h>155$ km and $T_{\infty}\in[650,1350]$ K.

If more accuracy is needed, the polynomial order can be increased and/or spline polynomial interpolation applied. Finally, the time-dependent atmosphere is recovered by inverting Equation 20

[TABLE]

Figure 5 compares the accuracy of the $T_{\infty}$-variable smooth exponential atmosphere model against the original Jacchia-77 model. It shows the ratio between $\rho_{S}(T_{\infty})/\rho_{J}(T_{\infty})$ for $T_{\infty}$ in the range from $650$ K to $1350$ K (left), and the corresponding parameters, $\hat{\rho}$p and $H$p as a function of $T_{\infty}$, including the underlying parameters of the static fits (right). Towards the lower edge of the temperature range (i.e. $T_{\infty}\rightarrow 650$ K), the polynomial fits for components $p=5-7$ do not well represent the underlying data. This leads to increased but still tolerable errors in the altitude range between $500$ and $1500$ km.

The advantage of this approach is, that the original structure of the model is maintained, so it can be used with the contraction model introduced in the next section.

4.2 Superimposed King-Hele Approximation

The extension of the KH contraction formulation into the SuperImposed King-Hele (SI-KH) formulation with a superimposed atmosphere is straightforward. Replacing $\rho$ from Equation 8 with the one defined in Equation 15 leads to

[TABLE]

i.e. each partial contraction reduces to the classical KH formulation with the partial exponential atmosphere $\rho$p. The important difference is that now $H$p is constant over the whole altitude range. The classical KH approximations – extended up to 5th order – can be found in B (dropping the subscript $p$). Finally, the rate of change is

[TABLE]

KH introduced the simple fixed boundary condition e_{b}$$=0.2 to select between the approximation method for low eccentric and high eccentric orbits, given in B.2 and B.3, respectively. However, as $H$p can be large, this condition is not always sufficient. Recall from Equation 13 that

[TABLE]

For low $a$ and high $H$, $z$ can approach unity at $e=0.2$, making the series expansion in $1/z$ inaccurate. Instead, it is proposed to define $e_{b}$ based on the truncation errors found in the formulations for the low and high eccentric orbits. The series truncation errors for the low eccentric orbit approximation (Equation 44), using the order notation, $\mathcal{O}$, are of the order of

[TABLE]

If $z$ is large (see justification below), $I_{0/1}(z)\to\exp(z)/\sqrt{2\pi z}$ and Equation 29 becomes

[TABLE]

For the high eccentric orbit approximation (Equation 46), the truncation errors are in the order of

[TABLE]

Assuming that the terms $\frac{(1+e)^{\frac{3}{2}}}{(1-e)^{\frac{1}{2}}}\frac{1}{(1-e^{2})^{6}}$ and $\left(\frac{1+e}{1-e}\right)^{\frac{1}{2}}\frac{1}{(1-e^{2})^{5}}$ are dominated by $1/z^{6}$ (see again below for a justification), Equation 31 simplifies to

[TABLE]

Equating the truncation errors from Equations 30 and 32, using Equation 13 and solving for $e$ results in the following condition

[TABLE]

Note that this boundary is most exact if the series expansions in both the low and high eccentric regimes are of the same order.

The derivation of the boundary condition required the assumptions of $z$ to be large, such that $I_{0/1}(z)\to\exp(z)/\sqrt{2\pi z}$ and such that $1/z^{6}$ dominates the other $e$ terms in Equation 31. To validate the assumptions, replace $a$ in Equation 33 with $a=(h_{p}+R_{E})/(1-e_{b})$ and solve for $e_{b}$, neglecting the negative solution

[TABLE]

Given $H_{min/max}=4.9/1320$ km (see Table 3) and the valid range for $h_{p}\in h\in[100,2500]$ km, the extrema in $e_{b}$ and $z_{b}=1/e_{b}$, are found to be

[TABLE]

For any $z>z_{b,min}$, $I_{0/1}$ remains close to $\exp(z)/\sqrt{2\pi z}$, being off only $+6\%$ and $-16\%$, respectively, at $z_{b,min}$. At the same time, $1/z^{6}$ dominates the terms dependent on $e$ in Equation 31 by two to three orders of magnitude $\forall e<e_{b,max}$. Thus, the assumptions made to derive $e_{b}$ are valid.

An advantage of an analytical expression of the dynamics is that the Jacobian of the dynamics can be derived analytically too, which can be used for uncertainty propagation. For a comprehensive discussion of the SI-KH method, the partial derivatives of the dynamics as derived by KH, with respect to $a$ and $e$, are given in C (again dropping the subscript $p$). As the SI-KH method is simply a summation of the individual contributions of the partial atmosphere, the derivatives can equally be summed up as

[TABLE]

5 Validation

The validation section is split into two parts: Section 5.1 validates the smooth exponential atmosphere, $\rho$S, by comparing it to the Jacchia-77 model, $\rho$J, during SA propagation using the GL contraction method; Section 5.2 validates the proposed SI-KH approach by comparing the contraction approximation along a single orbit, i.e. $\Delta a$ and $\Delta e$, to numerical quadrature. For completeness, propagations of a grid of initial conditions are performed using the GL and SI-KH methods and NA integration. The latter does not resort to any averaging technique, integrating the full dynamics of Equation 6, including $E$.

5.1 Validation of the Smooth Exponential Atmosphere Model

To validate $\rho$S against $\rho$J for $T_{\infty}$ $=750$, $1000$ and $1250$ K and at the same time distinguish it from the effects introduced by the SI-KH method on the resulting lifetime, $t_{L}$, the following orbits are propagated using the GL method only for the computation of the orbit contraction. All physically feasible initial orbit configurations on a $46\times 46$ grid from $250\leq h_{p}\leq 2500$ km and $250\leq h_{a}<2500$ km are propagated, using $\delta=1$ m2/kg. The lower limit, $250$ km, is selected as an object with such a large $\delta$ on a circular orbit at this altitude survives for a fraction of a day only at which point SA propagation becomes inaccurate. The upper limit, $2500$ km, is being imposed by definition of $\rho$J, but can be overcome by fitting to another model. The chosen $\delta$ is large, but does not limit the validity of this validation, as inaccuracies from the SA approach affect the propagation equally for both atmosphere models.

The SA propagation is performed using Matlab’s ode113 – a variable-step, variable-order Adams-Bashforth-Moulton integrator (Shampine and Reichelt, 1997) – and a relative error tolerance, \gamma_{rel}$$=10^{-6}, which is shown to be sufficient for different orbital scenarios in Section 5.2.

Figure 6 shows $t_{L}$ for the initial orbit grid, for propagation subject to $\rho$S, and the relative error, $\eta_{t_{L}}$, defined as

[TABLE]

comparing the propagations for each grid point using $\rho$S and $\rho$J, respectively. Table 4 contains information about the maximum error and the workload. Over the whole specified domain and for all $T_{\infty}$ $\in[750,1000,1250]$ K, $\eta_{t_{L}}$ remains within $[-0.1\%,0.1\%]$, which considering the uncertainties in atmospheric density modelling is more than accurate enough (Sagnieres and Sharf, 2017). Towards low perigees ($h_{p}<500$ km), the fitted $\rho$S starts to wobble around the underlying model (see Figure 4a), which is also apparent for the propagated orbits. A $6$-fold speed improvement can be observed, as no numerical integration is required when calculating the density with $\rho$S.

The reduction in function evaluations and computational time observable with an increasing $T_{\infty}$ is a consequence of the different density profiles. Increasing $T_{\infty}$ leads to an increased $\rho$, which increases the drag force and thus decreases the lifetime. However, the variable-step size integration method can compensate this by increasing the step size. Two possible explanations are: as the integrator is initialised with the same properties for all three cases, the initially set (small) step size favours shorter lifetimes; and the shape of the density profiles with high $T_{\infty}$ are more smooth, decreasing the number of failed function evaluation attempts.

5.2 Validation of the superimposed King-Hele method

The SA propagation relies on an accurate approximation of $\Delta a$ and $\Delta e$. Figure 7 shows – for different orbital configurations – the relative integral approximation error, $\eta_{\Delta x}$, defined as

[TABLE]

where $\mathcal{C}$ is the selected contraction method: $\mathcal{C}_{1}$ is the numerical GL method computed using 65 nodes; and $\mathcal{C}_{2}$ describes the analytical formulation, KH or SI-KH, using series expansion up to $5$th order.

Figure 7a reveals why orbits are predicted to re-enter much later using the classical KH contraction method: the density is underestimated at altitudes above $h_{p}$. The largest errors occur around $h_{p}=125$ km and $800$ km, where the rate of change in $H$ with respect to $h$ is large. Around these two altitudes, the contraction rate in $a$ is underestimated by more than $10\%$ and $20\%$, respectively, if $e>0.03$. Using the SI-KH the relative error remains well below $0.1\%$ $\forall h_{p}\in[100,2500]$ km and $\forall h_{a}\in[100,100000]$ km (see Figure 7b), a range that includes the vast majority of all Earth orbiting objects. Discontinuities can be found whenever $e$ passes through $e_{b}=e_{b}(H_{p})$. The biggest step occurs for the largest $H$p. Those discontinuities slightly increase the number of steps required during the integration. However, given the averaged dynamics, $\gamma_{rel}$ can be chosen large enough during integration mitigating the effects of the discontinuities.

To see how the SI-KH compares against GL in SA propagation and against NA propagation in terms of accuracy and computational power, the results from different initial orbit conditions are compared, for two scenarios:

1. a)

Short-term re-entry duration: $t_{L}=30$ days 2. b)

Mid-term re-entry duration: $t_{L}=360$ days

The reasons why long-term re-entry cases are not discussed here are two-fold: First, for long time spans, the NA integration requires small relative tolerances. If they are not met, the result cannot be trusted; Secondly, the longer the time spans, i.e. the smaller $\delta$, the more accurate the assumptions made for the SA propagation.

The initial conditions are spaced in $h_{p}\in[250,2500]$ km and $h_{a}\in[250,$ $100000]$ km and consist of all the $1558$ feasible solutions on a $46\times 46$ grid, where the grid spacing in $h_{a}$ is chosen to be logarithmic, as opposed to the equidistant grid in $h_{p}$. Two preliminary runs were performed using the SI-KH method to calculate the lifetimes. This way, the $\delta$ required to re-enter within the given time-span can be estimated. Figure 8 shows the grids of the resulting $\delta$ for both scenarios. Note that $\delta$ varies by almost 13 orders of magnitude.

The accuracy is described again as the relative lifetime, $\eta_{t_{L}}$, this time defined as

[TABLE]

where $\mathcal{M}$ is the selected contraction and integration method, combined with a given relative integrator tolerance, $\gamma_{rel}$, during integration. To give a feeling for the accuracy across all the different initial conditions, the $50\%$- and $100\%$-quantiles, i.e. the median and maximum denoted as $\eta_{t_{L},50\%}$ and $\eta_{t_{L},100\%}$, respectively, over all the $|\eta_{t_{L}}^{ij}|$ are given. The computation effort is compared via the total number of function calls, $N_{f}^{tot}$, and time required for the integration itself, $t_{CPU}^{tot}$

[TABLE]

Table 5 contains these figures comparing the different integration methods against each other. For both SI-KH and GL, the absolute maximum error over the whole grid and over both scenarios remains below $0.01$%, when decreasing $\gamma_{rel}$ from $10^{-6}$ to $10^{-12}$. Given this force model, it is therefore sufficient to use $\gamma_{rel}=10^{-6}$. For NA integration, this is not the case. While the maximum error remains modest ($0.18\%$) in the short-term case, it becomes large when the re-entry span is increased to one year ($26\%$), when decreasing $\gamma_{rel}$. Decreasing $\gamma_{rel}=10^{-9}$ and comparing to integration with $\gamma_{rel}=10^{-12}$, reduces the maximum error for the NA propagation in the mid-term case to $0.032\%$.

For the comparison of the SA techniques against NA propagation, the tolerance of the latter is set to $\gamma_{rel}=10^{-12}$. Again, SI-KH and GL fare very similar. For the short-term case, the boundaries of the SA propagation can be recognised for very high $\delta$, leading to still small maximum errors of $0.18\%$ and $0.17\%$, respectively. Figure 9a shows the resulting lifetime comparison for SI-KH and $t_{L}=30$ days. As $\delta$ increases to values above $10^{4}$ m2/kg, the assumption of constant $a$ and $e$ over one orbit starts to break down and small errors are introduced. This might be an issue for small debris such as multi-layer insulation fragments and paint flakes. For the mid-term scenario, the maximum error reduces by one order of magnitude for both SA methods tested. For high $h_{a}>10000$ km, the series expansion applied in the SI-KH method introduces small errors (see Figure 9b).

6 Conclusion

The classical KH orbit contraction method allows to analytically calculate the effects of drag on the orbit evolution averaged over an orbital period. However, it inaccurately estimates the orbital decay for eccentric orbits subject to a non-exponentially decaying atmosphere model. To improve the accuracy, a smooth exponential atmosphere model was proposed to be used in tandem with the new SI-KH orbit contraction method.

The classical KH method was extended to the SI-KH contraction method, making use of a superimposed atmosphere model to satisfy the assumption of a strictly decaying density for each component of the model. This greatly reduces the errors in the estimated decay rates of objects in eccentric orbits and subject to atmospheric density profiles with variable scale height. The analytical method was validated against an averaging technique based on numerical quadrature. Further, the semi-analytical propagation of orbits using the SI-KH method was validated against full numerical integration of the dynamics. The approach is applicable to any averaging techniques considering drag and based on the fixed scale height assumption above perigee. Finally, the Jacobian of the dynamics governed by the SI-KH method is given to be used for future applications such as uncertainty propagation.

7 Acknowledgements

The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme as part of project COMPASS (Grant agreement No 679086). The authors acknowledge the use of the Milkyway High Performance Computing Facility, and associated support services at the Politecnico di Milano, in the completion of this work. The datasets generated for this study can be found in the repository at the link www.compass.polimi.it/ publications.

8 References

Appendix A Nomenclature

Appendix B King-Hele Formulation

All the formulas presented here are explained and derived in the work of King-Hele (1964). The analytical formulas describe, for different eccentricities, the change in the semi-major axis, $a$, and the eccentricity, $e$, over one orbit as an approximation of Equation 8. Please note that one of the four cases was dropped, as it was introduced only due to the Bessel functions becoming inaccurate for small arguments. Today, the relevant mathematical software packages are accurate and fast enough to overcome this limitation.

Adaptations to the original formulation were made to

• •

find the change directly in $a$ and $e$, rather than $a$ and $x=ae$, to calculate the change in the variables of interest;

• •

find a more appropriate boundary condition, $e_{b}$, for the selection of the phase (see Section 4.2);

• •

increase the accuracy within each phase by taking into account more terms after the series expansion.

The two functions, $k_{a}$ and $k_{e}$, are introduced here for later use when describing the rate of change in all the eccentricity regimes described below, as

[TABLE]

with the effective area-to-mass ratio, $\delta$, the gravitational parameter, $\mu$, the atmospheric density, $\rho$, evaluated at the perigee altitude, $h_{p}$.

B.1 Circular Orbit

For circular orbits, no integration needs to be approximated, as the integral can be solved analytically as

[TABLE]

where $h_{p}$ reduces to the circular altitude. Dividing by the orbital period, $P$, according to Equation 9 and using the functions defined in Equation 41, the rate of change for circular orbits is

[TABLE]

B.2 Low Eccentric Orbit

For small $e<e_{b}(a,H)$, a series expansion in $e$ is performed and then integrated using the modified Bessel function of the first kind, $I_{n}$($z$), as

[TABLE]

with the auxiliary variable $z=ae/H$, the scale height, $H$, a single evaluation of the density at the perigee height, $h_{p}$, and the order of the series truncation error, $\mathcal{O}$, of $e^{6}$. The constant matrices are given as

[TABLE]

Dividing by $P$ according to Equation 9 and using the functions defined in Equation 41, the rate of change for low eccentric orbits is

[TABLE]

B.3 High Eccentric Orbit

Instead of performing the series expansion in $e$, which is infeasible for large values of $e>e_{b}(a,H)$, the expansion is performed for the substitute variable, $\lambda^{2}/z=1-\cos E$. KH truncated the series already after two powers. Here, instead, as $H$ can be large and as the formulation should be readily available for any $h_{p}<2500$ km, it is extended up to 5th power. The contractions over one orbit period are

[TABLE]

with the constant matrices

[TABLE]

[TABLE]

Plugging Equation 46 into Equation 9, using the functions defined in Equation 41, and introducing the new functions

[TABLE]

the rate of change for highly eccentric orbits is

[TABLE]

Appendix C Jacobian of Dynamics in $a$ and $e$

The partial derivatives of the dynamics with respect to $a$ and $e$ are given here for the three different regimes discussed in B. The partial derivatives of a partial atmosphere defined in Equation 15 (dropping the subscript $p$), and given $h_{p}=a(1-e)-R$, can be found as

[TABLE]

Thus, the partial derivatives of $k_{a}$ and $k_{e}$ (see Equation 41) with respect to $a$ and $e$ are

[TABLE]

C.1 Circular Orbit

For circular orbits, the rate and derivative in $e$ vanishes and the partial derivative of $F_{a}$ with respect to $a$, combining Equations 43 and 50, is

[TABLE]

C.2 Low Eccentric Orbit

For low eccentric orbits, with $e\leq e_{b}$, the partial derivative of $F_{a}$ and $F_{e}$ with respect to $a$ and $e$, combining Equations 45 and 50 and using the product rule, are

[TABLE]

where

[TABLE]

and

[TABLE]

C.3 High Eccentric Orbit

Using the partial derivatives of $c_{a}$ and $c_{e}$ from Equation 47 with respect to $a$ and $e$

[TABLE]

it follows that

[TABLE]

and the partial derivatives of $F_{a}$ and $F_{e}$ for high eccentric orbits (see Equation 48), with $e\geq e_{b}$, with respect to $a$ and $e$ become

[TABLE]

where

[TABLE]

Appendix D Variable Atmosphere Model Parameters

Tables 7 and 8 list the parameters to calculate $\boldsymbol{a}$ and $\boldsymbol{b}$ according to Equation 23 as a function of the normalised $\tilde{T}_{\infty}$. The two vectors are needed to recover $\hat{\rho}$p and $H$p $\forall p$, according to Equation 25. Note that the model should only be used for $T_{\infty}\in[T_{0}=650,T_{1}=1350]$ K.