SDN-AP: Mean-Element Propagation from Low Earth Orbit to Interplanetary Cruise

Executive summary

SDN-AP, the Space Data Network advanced propagator, turns element sets into states for a whole catalog at once. One WebAssembly module carries three theories (section 1): SGP4, the theory of the public catalog; the SGP4XP mode, which implements the extended theory the US Space Force publishes; and SDNAP, the Space Data Network's own theory. The module returns the same bytes under its browser host and under native WasmEdge, at any thread count (section 12.2).

The SDNAP theory contributes:

Held out against precise orbits, sets fitted to 24 h of data have median day-3 position errors of 0.16 to 0.47 km in low Earth orbit (laser-ranged satellites), 0.13 to 0.18 km for navigation satellites in medium Earth orbit and 0.10 to 0.14 km in geostationary orbit (section 12.3). The SGP4XP mode reproduces the reference implementation of its theory to a median of 0.37 µm (section 12.1).

1 The theory family

1.1 SGP4

The SGP4 base is written from Spacetrack Report No. 3¹ and from the account of the theory by Hoots, Schumacher and Glover².

1.2 The SGP4XP mode

The mode implements SGP4-XP as the US Space Force describes it³:

1.3 The SDNAP theory

SDNAP keeps the state, integrator, grid, checkpoints and output chain of section 2 and changes the forces:

2 Mean elements and averaging

2.1 State and rates

The state is $y = (n, a_f, a_g, \chi, \psi, \lambda)$ with

regular at zero eccentricity and zero inclination⁴. An acceleration with radial, transverse and normal components $(R, T, N)$ changes the elements by Gauss's equations in equinoctial form¹⁴. With $L$ the true longitude, $w = 1 + a_f\cos L + a_g\sin L$, $s^2 = 1+\chi^2+\psi^2$, $p$ the semilatus rectum and $h = \sqrt{\mu p}$:

A force enters the theory as its average over the mean anomaly,

and gravity enters through Lagrange's equations for the averaged disturbing function (section 3). A mean-motion error grows into an along-track error $s(t) = a\int\delta n\,dt$. A fractional error $\delta$ in the drag, held for a time $t$, gives

with $\dot a_D = -\rho B\sqrt{\mu a}$ the drag decay of the semi-major axis: 3.4 km after three days for a 10 % error at 550 km, where $\dot a_D = -53$ m/day. Most of the error of a low orbit lies along track, so each term below is judged by what it does to the mean motion.

2.2 Integration and output

2.3 Orbit averages of drag

The density at each node is taken at the osculating radius: the mean orbit's plus its first-order J2 short-period radial terms. Their constant part moves the mean height by −0.43 km at 300 km and 53° inclination, about 1 % of the density.

2.4 The inverse map

The mean elements whose output is a given osculating state $x$ at $t$ solve

where $q$ gives the osculating equinoctial elements of a state: the osculating-to-mean iteration¹⁷. The short-periodic corrections are of order J2, so the iteration converges to $10^{-9}$ km and $10^{-12}$ km/s in about four steps. Maneuvers, detection and launch sets use it (sections 8 to 10).

3 Gravity

Kaula's disturbing function¹² is

with $S_{lmpq}$ the cosine and sine of $\psi_{lmpq}$ weighted by the field coefficients $\bar C_{lm}$ and $\bar S_{lm}$. The terms split by $j = l-2p+q$, the multiple of the mean anomaly in the angle:

The field is EGM96 to degree 70¹⁸, with EGM2008¹⁹ as an option. Each truncation is the degree beyond which the median day-3 error against full-force Starlink truths stops falling.

3.1 Zonal

J2 to J5 enter in closed-form averaged rates (section 1.2). J6 to J20 add first-order rates: the secular terms of even degree, and the long-period terms of odd degree with angle $\pm\omega$, which are of order $e$ in $G$ but of order one in the eccentricity vector. Fitted to one full-force truth with zonal gravity only, the day-3 error is 0.35 km with J2 to J5, 0.14 km to J10, 0.047 km to J20 and 0.042 km to J40.

3.2 m-daily terms

For $j = 0$, $q = 0$ and even $l$ ($p = l/2$) the angle is $\psi = m(\Omega - \theta_G)$, with a period of about one $m$-th of a day. With the slow elements frozen over that period and $\psi' = m(\dot\Omega_{J2} - \omega_E)$, Lagrange's equations integrate in closed form. Since $\partial R/\partial M = 0$, the semi-major axis, the mean motion and the eccentricity do not oscillate:

with $A = \mu R_E^{l}/a^{l+1}$, $F = \bar F_{l,m,l/2}(i)$ at the epoch inclination, $F' = dF/di$, $G = \langle(a/r)^{l+1}\rangle$, $\eta = \sqrt{1-e^2}$, $S$ the term's combination of $\bar C_{lm}$ and $\bar S_{lm}$ and $S_1$ its antiderivative in $\psi$: the equinoctial form of Danielson et al.²⁰. The corrections are added to the elements at output time. The normalized inclination functions are Fourier coefficients. On a circular orbit, $\bar P_{lm}(\sin i\sin u)\,e^{\mathrm i m\Delta\lambda(u)} = \sum_p \bar F_{lmp}\,e^{\mathrm i(l-2p)u}$, with $\Delta\lambda(u)$ the longitude from the node, is a trigonometric polynomial of degree $l$; so $2l + 1$ samples give every coefficient exactly, and the normalized Legendre recurrence²¹ keeps the values stable to degree 70.

The m-daily set is a few hundred meters along track in every low orbit (EGM2008, even $l \le 8$, exact inclination functions):

3.3 Eccentricity-vector terms

Odd-degree terms with $q = \pm1$ and $j = 0$ have the angle $\psi = w\,\omega + m(\Omega-\theta_G)$, $w = -q$. They are of order $e$ in $G$ but of order one in the eccentricity vector, through the $1/e$ of Lagrange's equations, and they oscillate at the m-daily frequencies, amplified by $n/\dot\psi$. They are added at output time, to first order and to degree 7. Fitted to one full-force truth with tesserals to degree 3, they take the fit residual from 0.274 to 0.051 km and the day-3 error from 2.6 to 0.37 km.

3.4 Shallow resonance

The angle $\psi = (M+\omega) + m(\Omega-\theta_G)$ of order $m$ drifts slowly when about $m$ revolutions fill a nodal day. For each order whose angle moves slower than $\varepsilon n$, with $\varepsilon = 0.045$, the terms of odd degree to 31 enter the averaged rates, with $F$ and $dF/di$ at the epoch inclination and $G$ at the epoch eccentricity. The grid resolves their periods, $1/\varepsilon \approx 22$ revolutions and longer. With $x = \lvert\dot\psi\rvert/(\varepsilon n)$ from the current elements, each term carries the weight

so a term that crosses the band edge as drag raises $n$ fades instead of switching; terms are selected at the epoch out to $x = 1.3$. The 530 to 570 km Starlink shells lie within about 0.2 rev/day of the 15:1 commensurability. On the 53° shell $\dot\psi \approx -0.164$ rev/day: a 6.1-day period and, in EGM2008, an along-track amplitude near 0.8 km. The 97.6° shell sits on the resonance, where the term acts as a nearly constant transverse force of up to $2.7\times10^{-8}$ m/s².

3.5 Deep resonance

The one-day, half-day, 16 h and 8 h classes of section 1.2 take $G_{lpq}(e)$ as the Hansen coefficients themselves^22,12, valid at every eccentricity. Power series in $e$ diverge beyond the Laplace limit $e = 0.6627$²³, and Molniya orbits sit near $e = 0.72$. The derivatives of the inclination functions are exact. The functions are held at their epoch values.

3.6 How an omitted term enters the drag

Take an along-track term $A\cos(\omega_m t + \varphi)$ that a theory omits, on a fit arc of length $T$ that holds $C = \omega_m T/2\pi$ cycles, with the arc mapped to $x \in [-1, 1]$. It projects onto the quadratic signature of drag with the Legendre coefficient

from $e^{\mathrm i kx} = \sum_l (2l+1)\,\mathrm i^l j_l(k) P_l(x)$, with $j_l$ the spherical Bessel functions. By equation (1) the drag itself has

which is 627 m at 24 h, 2,507 m at 48 h and 5,641 m at 72 h for $\dot a_D = -53$ m/day at 550 km. The fit biases the drag by $\delta B/B = c_2/c_{2,D}$ and carries the bias forward as $c_2 P_2(1 + 2L/T)$ at a lead $L$ past the arc; at the end of day 3 after a 24 h fit, $P_2(7) = 73$. For the J22 term of the 53° shell (period 11.8 h, along-track amplitude 454 m) the bias is 22 %, 2.1 % and 0.64 % at 24, 48 and 72 h (RMS over the phase). The $m = 1$ terms, with periods near a day, alias harder: $\lvert 5 j_2(\pi C)\rvert = 1.53$ against 0.436 at 24 h. The tesseral terms close this path in the theory (sections 3.2 to 3.4) and longer spans close it in the fit (section 7).

4 Atmosphere and density calibration

4.1 Density and environment

Beyond the last row or node the values in force are held. The atmosphere co-rotates with the Earth.

4.2 The catalog-shared density kernel

The atmosphere at a given time is the same for every object. NRLMSISE-00 is therefore evaluated once per kernel time on a grid, every object's quadrature nodes interpolate it, and $e^{s(t)}$ is applied at the query time:

Ninety-six grids (10.7 days) share one allocation of 6.3 MiB, so a catalog propagated over a week builds each grid once. A grid depends only on its time and the environment, never on which object first needs it. Against direct evaluation along 53° orbits, $\ln(\rho_{\mathrm{kernel}}/\rho_{\mathrm{direct}})$ has an RMS of 0.21 %, 0.33 % and 0.47 % at 350, 450 and 550 km, and at most 1.7 %.

4.3 Effective ballistic coefficient and thrust

A drag-makeup thrust proportional and opposite to the drag, $\mathbf a_T = -\kappa\,\mathbf a_D$, leaves the rates of drag with

so the averaged theory carries an operator's thrust, storms included, through one fitted number. $B_{\mathrm{eff}}$ is signed: net thrust, $\kappa > 1$, is a negative $B_{\mathrm{eff}}$. It acts in the drag regime, perigee below 2,500 km. A ballistic arc observes only $B_{\mathrm{eff}}$. Given a prior $B \sim N(\hat B, \sigma_B^2)$, the set reports $\kappa = 1 - B_{\mathrm{eff}}/\hat B$ with

4.4 Empirical transverse terms

with $u$ the argument of latitude, $t_e$ the end of the fit arc and $\tau_e = 24$ h. To first order in $e$ the averaged rates are $\dot a = 2\tau_0 g/n$ and $\tfrac{d}{dt}(a_f + \mathrm i a_g) = e^{\mathrm i\Omega}(\tau_c + \mathrm i\tau_s)\,g/(na)$. The constant term is an along-track quadratic $-\tfrac32\tau_0 t^2$; the pair turns the eccentricity vector in a fixed direction of the orbit frame, the once-per-revolution empirical acceleration of precise orbit determination^29,30. $\tau_0$ and $B_{\mathrm{eff}}$ share the $t^2$ signature, correlated as $\langle\rho\rangle/\sqrt{\langle\rho^2\rangle}$, above 0.99 for a 10 % density variation, so one of the two is free per object: $B_{\mathrm{eff}}$, or $\tau_0$ with $B$ at its prior where the configuration selects sustained thrust (section 11). Above the drag regime, $\tau_c$ and $\tau_s$ carry the once-per-revolution part of radiation pressure and are held after the arc. On 54 GNSS windows (ESA final orbits) the day-3 median is 0.176 km with them held, 0.243 km with them decaying and 0.219 km without them.

5 Radiation pressure and third bodies

The lightness number $\beta = P_{1\,\mathrm{AU}}\,(1\,\mathrm{AU})^2\,\mathrm{AGOM}/\mu_\odot = 7.69\times10^{-4}\,\mathrm{AGOM}$, with AGOM in m²/kg, measures radial radiation pressure against the Sun's gravity. For high lunar orbits the Earth's quadrupole drives the Kozai–Lidov exchange of eccentricity and inclination^32,33, which the averaged rates carry. The DE440 positions are refitted to Chebyshev records: the geocentric Moon to 1.6 m, the Sun to 3.3 km and the planets to between 2.6 m (Venus) and 0.25 km (Jupiter and Saturn).

6 Moon-centered, cislunar and heliocentric sets

6.1 Moon-centered mean elements

The Moon turns once in 27.3 days, so the tesseral angles $m(\Omega-\theta_M)$ are slow for every lunar orbit, and the lunar tesserals belong in the averaged rates rather than in the output:

The averaged degree sets the lunar error. On 8 windows of each orbiter, day 3 is 6.5, 3.7, 1.95 and 1.83 km for LRO and 22.6, 6.3, 5.9 and 4.0 km for Danuri at degrees 10, 30, 50 and 70.

6.2 Cislunar and heliocentric segments

Where averaging about one body does not hold, a segment carries a fitted epoch state that is integrated numerically:

6.3 Reach and segmentation

A segment fit tries the whole arc and its latest half and quarter, each at least the shortest segment, and keeps the longest whose residual is within $\max(2r,\ r + 0.05\ \mathrm{km})$ of the smallest residual $r$, so a burn or a change of thrust ends the earlier segment. Each segment is an SDS $OMM^35,36 with CENTER_NAME, REFERENCE_FRAME ICRF, the fitted epoch state as osculating elements with SEMI_MAJOR_AXIS and GM (an open orbit has $e > 1$ and $a < 0$), MEAN_ELEMENT_THEORY SDNAP, the regime (CISLUNAR or HELIO), the force parameters, and the time by which the epoch was carried past a maneuver or a sphere-of-influence passage.

7 Fitting and spans

The fit adjusts the epoch elements and the parameters $(B_{\mathrm{eff}}, \mathrm{AGOM}, \tau_c, \tau_s)$ to the positions of an SDS $OEM by Levenberg–Marquardt^37,38, with analytic Jacobians: the exact derivatives of two-body motion with the J2 secular drifts, and the drag column through the theory's $\dot n$ per unit $B$. A prior on $B_{\mathrm{eff}}$ enters as one more row: a previous set's value and sigma (the sigma floored at 10 %), or the configuration's. Each fit starts on two revolutions and grows the arc by factors of three to the full interval, each stage starting from the last solution. The span follows the data:

Longer spans shrink the aliasing of section 3.6 by $T^2$ through equation (7). On laser-ranged spheres in low orbit (Ajisai, LARES, LARES 2, Larets, Starlette, Stella and WESTPAC; 21 windows of ILRS orbits³⁹; fits ending at the same time; median 3D RMS, km):

A 6 h tail scores the fit's own short-period error rather than its drag estimate, so the tail is a day.

8 Maneuvers

8.1 Impulse and finite-burn maps

For an impulse $\Delta\mathbf v = (\Delta v_U, \Delta v_V, \Delta v_W)$ in the radial, in-track and cross-track axes of the pre-burn state at $t_b$, the position is continuous and the post-burn mean elements follow from the inverse map, equation (2):

To first order on a near-circular orbit, $\Delta a = 2\Delta v_V/n$, $\Delta(a_f + \mathrm i a_g) = e^{\mathrm i L}(2\Delta v_V - \mathrm i\Delta v_U)/(na)$, $\Delta(\psi + \mathrm i\chi) = s^2 e^{\mathrm i L}\Delta v_W/(2na)$ and $\Delta\lambda = -2\Delta v_U/(na)$: a radial impulse leaves the semi-major axis unchanged, and a cross-track one moves only $\chi$ and $\psi$. A burn of constant acceleration $(a_U, a_V, a_W)$ lasting $\Delta t_b$ sweeps $\Delta L = n\Delta t_b$ about its centroid $L_c$. Integrating the same rates over it gives

with $\operatorname{sinc}x = \sin x/x$. A finite burn is therefore the map (12) at its centroid with the eccentricity-vector and plane changes scaled by $\operatorname{sinc}(\Delta L/2)$, the arc loss of the impulsive approximation⁴⁰. Each burn yields a post-burn $OMM with its epoch at the centroid, and an SDS $MNV with the pre- and post-burn states and the characterization of $\Delta\mathbf v$.

8.2 Detection

The statistic

is the generalized likelihood ratio of a jump⁴², distributed as $\chi^2_3$ without a burn. A burn is declared when $\max_{t_b}\Lambda$ exceeds 20: about 5,500 candidate epochs a year in low orbit give $\alpha \approx 1.8\times10^{-4}$ for one false alarm per object-year. The epoch and $\Delta\mathbf v$ then come from the likelihood profiled over the nuisance $\theta$, the epoch elements and $B_{\mathrm{eff}}$, re-estimated with and without the burn:

with $A = \sum_k J_k^{\top}WJ_k$ and $\mathbf b = \sum_k J_k^{\top}We_k$ over the arc, $A_+$ and $\mathbf b_+$ the same sums over $t_k > t_b$ (suffix sums, so one pass serves every candidate epoch), $J_k = \partial e_k/\partial\theta$ and $\tilde N$ the burn's image in $\theta$. The $MNV gives the epoch where $\Lambda_p$ peaks and the interval within 3.84 of the peak, $\chi^2_1$ at 95 %⁴³; a burn-free arc reports the floor $\sqrt{20/F_{cc}}$ per axis at the best epoch. After a detection, a joint segment fit estimates $(y_0, \Delta\mathbf v, t_b)$ over the whole arc, the post-burn set being the map (12) of $y_0$ at $t_b$, so the early segment still informs $B_{\mathrm{eff}}$ and the post-burn state.

8.3 Measured

9 Reentry

9.1 Where averaging ends

Averaging needs the elements to change little over one integration step of two periods, against the scales the rates depend on. At nodes every 30 min along the mean trajectory, within 120 days:

$H$ is the scale height of the density model at the perigee height, and the per-revolution changes are the averaged rates over one period. Drag on an eccentric orbit acts mostly at perigee⁴⁶, so the osculating semi-major axis steps by $\Delta a_{\mathrm{rev}}$ at each pass while the mean theory ramps; the along-track difference is a sawtooth of $\tfrac{3\pi}{8}\lvert\Delta a_{\mathrm{rev}}\rvert$ peak to peak.

9.2 Hand-off and window

At the first node where a criterion trips, the set hands off, one way, to osculating integration: EGM96 to degree and order 20, the Sun and Moon, NRLMSISE-00 times $e^{s(t)}$ under the same $SPW and $DCL records with the set's $B_{\mathrm{eff}}$, and radial radiation pressure, by fourth-order Runge–Kutta at 10 s down to 80 km. The density model is the same on both sides of the boundary. A set flagged as lowering decays by its sustained transverse thrust $\tau_0$, with $B_{\mathrm{eff}}$ at its prior.

A fractional ballistic error $\delta B$ and a density error $\delta(t)$ move a lifetime $L$ by $\delta L = -L\,\delta B/B - \int_0^L\delta(t)\,dt$. With $\delta$ first-order Gauss–Markov, of variance $\sigma_\delta^2$ and correlation time $\tau$⁴⁷, the 1-sigma window is

with $\sigma_\delta = 0.15$, $\tau = 2$ days, $\sigma_B/B$ from the set's $B_{\mathrm{eff}}$ sigma or the configuration, and $L$ the time to go from the end of the fitted data. With $\sigma_B/B = 0.05$ the window is 9.9 % of the time to go at 10 days and 14.7 % at 1 day. Each issue is one SDS $RDM with the keywords of the CCSDS Re-entry Data Message⁴⁸: the nominal epoch and its window, the orbit lifetime from the set's epoch, BALLISTIC_COEFF $= 1/B_{\mathrm{eff}}$ in kg/m², and the impact point.

9.3 Measured

Twelve decays integrated by FFM-1 to 100 km (three ballistic coefficients at three levels of solar activity, and three eccentric orbits; 134 nodes), each fitted at successive perigee heights and handed to reentry:

10 Launch sets

From the $LAM uncertainties of perigee and apogee height, argument of latitude and node, the first set's sigmas for a near-circular insertion are

A deployment releases objects from a parent state $x_0$ at $t_0$, object $k$ with a separation impulse $\Delta\mathbf v_k$ and its own $B_k$. Its set is the map (12) of the parent set with $\Delta\mathbf v_k$, plus element offsets under a tight prior (about 100 m of position). The fit estimates the parent's elements, the shared $\tau_c$ and $\tau_s$ and the batch mean $\bar B$, and per object $\Delta\mathbf v_k$ and $B_k$, under

the hierarchical prior of Stein shrinkage^49,50, with $\sigma_{\mathrm{sep}} = 1$ m/s. The normal matrix is an arrowhead, the parent block coupled to one block per object, and its Schur complement solves it in time linear in the number of objects. Each object's along-track offset from the parent follows Clohessy–Wiltshire⁵¹ plus differential drag,

and the objects are associated with later catalog sets by sorting both by along-track position.

11 Configuration through $PCF

Each method takes one SDS $PCF³⁵ on its required configuration port. THEORY selects SGP4, SGP4XP or SDNAP, and each method reads its own tables:

The truncations of section 3, the kernel of section 4.2, the empirical terms of section 4.4, the segment rules of section 6, the span rule of section 7, the detection threshold of section 8, the criteria and window of section 9 and the launch priors of section 10 are fields of these tables, so the configuration of a run names every choice that shaped its result. A field set to a value the module cannot apply is refused by name, as "TABLE.FIELD: reason", never silently replaced.

12 Verification

12.1 Conformance of the SGP4XP mode

A seeded corpus of 10,200 synthetic element sets, 1,700 per regime, is propagated to eight times over 7 days by the module and by the US Space Force's Sgp4Prop library, version 9.9, run as a black box³. Differences are 3D positions at the times both answer:

Where both refuse a decayed set, the two refusal epochs agree within 1 % of the elapsed lifetime (164 of 164). Every module output is digested with SHA-256, so a rebuilt module is checked bit for bit against a recorded run.

12.2 Two-runtime parity and determinism

The module⁵³ is compiled once, with clang for wasm32-wasip1-threads and without exception handling, and the same module.wasm runs in the browser host and in native WasmEdge 0.16.4⁵⁴. propagate returns identical $OEM bytes in both. Batch propagation returns the same bytes at any thread count, and no state depends on the order of requests: the checkpoints of section 2.2 and the kernel grids of section 4.2 are keyed by time and environment alone.

12.3 Held-out accuracy

Windows are registered in a frozen plan before scoring⁵⁵. Each fits the first 24 h of its truth (30 days for cruise, 14 days for the halo orbit) under the space weather of GFZ and NOAA SWPC and IERS Earth orientation, with $s = 0$ and no ballistic prior, and is scored on the days that follow. A cell is the median over windows of each day's 3D RMS, in km:

The Starlink truth shares the theory's density model and weather, so that row measures the theory's own approximations against full-force integration; the other rows measure the real atmosphere and the real spacecraft. Further out, the halo orbit's segmented error is 17.2, 23.1, 36.4 and 48.4 km at days 4 to 7, and at days 7 and 15 of cruise it is 1.4 and 5.6 km for Psyche, 0.6 and 4.3 km for Lucy and 0.9 and 5.8 km for Juno (720 h fits).

12.4 Cost

Times are one thread of the WebAssembly module on a desktop workstation.

References

R1

Hoots, F. R., Roehrich, R. L. (1980). Models for propagation of NORAD element sets. Spacetrack Report No. 3, Aerospace Defense Command.

R2

Hoots, F. R., Schumacher, P. W., Glover, R. A. (2004). History of analytical orbit modeling in the U.S. space surveillance system. Journal of Guidance, Control, and Dynamics 27(2), 174–185.

R3

US Space Force. Astrodynamics Standards: the Sgp4Prop library and its release notes, version 8.0 and later (the public description of SGP4-XP).

R4

Broucke, R. A., Cefola, P. J. (1972). On the equinoctial orbit elements. Celestial Mechanics 5, 303–310.

R5

Kutta, W. (1901). Beitrag zur näherungsweisen Integration totaler Differentialgleichungen. Zeitschrift für Mathematik und Physik 46, 435–453.

R6

Brouwer, D. (1959). Solution of the problem of artificial satellite theory without drag. Astronomical Journal 64, 378–397.

R7

Jacchia, L. G. (1970). New static models of the thermosphere and exosphere with empirical temperature profiles. Smithsonian Astrophysical Observatory Special Report 313.

R8

NOAA, NASA, US Air Force (1976). U.S. Standard Atmosphere, 1976. U.S. Government Printing Office, Washington, D.C.

R9

Abramowitz, M., Stegun, I. A. (1964). Handbook of Mathematical Functions. National Bureau of Standards, Applied Mathematics Series 55.

R10

Escobal, P. R. (1965). Methods of Orbit Determination. John Wiley & Sons.

R11

Kaula, W. M. (1962). Development of the lunar and solar disturbing functions for a close satellite. Astronomical Journal 67, 300–303.

R12

Kaula, W. M. (1966). Theory of Satellite Geodesy. Blaisdell.

R13

Allan, R. R. (1967). Resonance effects due to the longitude dependence of the gravitational field of a rotating primary. Planetary and Space Science 15, 53–76.

R14

Walker, M. J. H., Ireland, B., Owens, J. (1985). A set of modified equinoctial orbit elements. Celestial Mechanics 36, 409–419.

R15

IAU SOFA Board. IAU SOFA Software Collection, and ERFA: Essential Routines for Fundamental Astronomy. ERFA

R16

Trefethen, L. N., Weideman, J. A. C. (2014). The exponentially convergent trapezoidal rule. SIAM Review 56(3), 385–458.

R17

Walter, H. G. (1967). Conversion of osculating orbital elements into mean elements. Astronomical Journal 72, 994–997.

R18

Lemoine, F. G., et al. (1998). The development of the joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) geopotential model EGM96. NASA/TP-1998-206861, Goddard Space Flight Center.

R19

Pavlis, N. K., Holmes, S. A., Kenyon, S. C., Factor, J. K. (2012). The development and evaluation of the Earth Gravitational Model 2008 (EGM2008). Journal of Geophysical Research 117, B04406.

R20

Danielson, D. A., Sagovac, C. P., Neta, B., Early, L. W. (1995). Semianalytic satellite theory. Naval Postgraduate School, NPS-MA-95-002.

R21

Holmes, S. A., Featherstone, W. E. (2002). A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions. Journal of Geodesy 76, 279–299.

R22

Hansen, P. A. (1855). Entwickelung des Products einer Potenz des Radius Vectors mit dem Sinus oder Cosinus eines Vielfachen der wahren Anomalie in Reihen. Abhandlungen der Königlich Sächsischen Gesellschaft der Wissenschaften 4, 182–281.

R23

Murray, C. D., Dermott, S. F. (1999). Solar System Dynamics. Cambridge University Press.

R24

Picone, J. M., Hedin, A. E., Drob, D. P., Aikin, A. C. (2002). NRLMSISE-00 empirical model of the atmosphere: statistical comparisons and scientific issues. Journal of Geophysical Research 107(A12), 1468.

R25

Tapping, K. F. (2013). The 10.7 cm solar radio flux (F10.7). Space Weather 11, 394–406.

R26

Matzka, J., Stolle, C., Yamazaki, Y., Bronkalla, O., Morschhauser, A. (2021). The geomagnetic Kp index and derived indices of geomagnetic activity. Space Weather 19, e2020SW002641.

R27

Marcos, F. A., Kendra, M. J., Griffin, J. M., Bass, J. N., Larson, D. R., Liu, J. J. F. (1998). Precision low Earth orbit determination using atmospheric density calibration. Journal of the Astronautical Sciences 46(4).

R28

Casali, S. J., Barker, W. N. (2002). Dynamic calibration atmosphere (DCA) tool for the high accuracy satellite drag model (HASDM). AIAA 2002-4888.

R29

Colombo, O. L. (1989). The dynamics of Global Positioning System orbits and the determination of precise ephemerides. Journal of Geophysical Research 94(B7), 9167–9182.

R30

Wu, S. C., Yunck, T. P., Thornton, C. L. (1991). Reduced-dynamic technique for precise orbit determination of low earth satellites. Journal of Guidance, Control, and Dynamics 14(1), 24–30.

R31

Park, R. S., Folkner, W. M., Williams, J. G., Boggs, D. H. (2021). The JPL planetary and lunar ephemerides DE440 and DE441. Astronomical Journal 161, 105.

R32

Kozai, Y. (1962). Secular perturbations of asteroids with high inclination and eccentricity. Astronomical Journal 67, 591–598.

R33

Lidov, M. L. (1962). The evolution of orbits of artificial satellites of planets under the action of gravitational perturbations of external bodies. Planetary and Space Science 9, 719–759.

R34

Lemoine, F. G., et al. (2013). High-degree gravity models from GRAIL primary mission data. Journal of Geophysical Research: Planets 118, 1676–1698.

R35

Space Data Standards, release 1.264.0: the $OMM, $OEM, $MNV, $RDM, $LAM, $BOV, $SPW, $PRW, $DCL and $PCF schemas. Space Data Standards

R36

CCSDS (2023). Orbit Data Messages. Recommended Standard CCSDS 502.0-B-3, Blue Book.

R37

Levenberg, K. (1944). A method for the solution of certain non-linear problems in least squares. Quarterly of Applied Mathematics 2, 164–168.

R38

Marquardt, D. W. (1963). An algorithm for least-squares estimation of nonlinear parameters. Journal of the Society for Industrial and Applied Mathematics 11(2), 431–441.

R39

Pearlman, M. R., Degnan, J. J., Bosworth, J. M. (2002). The International Laser Ranging Service. Advances in Space Research 30(2), 135–143.

R40

Robbins, H. M. (1966). An analytical study of the impulsive approximation. AIAA Journal 4(8), 1417–1423.

R41

Hampel, F. R. (1974). The influence curve and its role in robust estimation. Journal of the American Statistical Association 69, 383–393.

R42

Willsky, A. S., Jones, H. L. (1976). A generalized likelihood ratio approach to the detection and estimation of jumps in linear systems. IEEE Transactions on Automatic Control 21(1), 108–112.

R43

Wilks, S. S. (1938). The large-sample distribution of the likelihood ratio for testing composite hypotheses. Annals of Mathematical Statistics 9(1), 60–62.

R44

Space Data Network. FFM-1, the full-force numerical propagator: propagator/ffm. Source

R45

Gumbel, E. J. (1958). Statistics of Extremes. Columbia University Press.

R46

King-Hele, D. G. (1987). Satellite Orbits in an Atmosphere: Theory and Applications. Blackie.

R47

Emmert, J. T., Warren, H. P., Segerman, A. M., Byers, J. M., Picone, J. M. (2017). Propagation of atmospheric density errors to satellite orbits. Advances in Space Research 59(1), 147–165.

R48

CCSDS (2019). Re-entry Data Message. Recommended Standard CCSDS 508.1-B-1, Blue Book.

R49

James, W., Stein, C. (1961). Estimation with quadratic loss. Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability 1, 361–379.

R50

Efron, B., Morris, C. (1975). Data analysis using Stein's estimator and its generalizations. Journal of the American Statistical Association 70, 311–319.

R51

Clohessy, W. H., Wiltshire, R. S. (1960). Terminal guidance system for satellite rendezvous. Journal of the Aerospace Sciences 27(9), 653–658.

R52

SpaceX. Starlink satellite ephemerides: operator predictions, reissued about every 8 h.

R53

Space Data Network. The SDN-AP propagator: propagator/sdn-ap, Apache-2.0. Source

R54

WasmEdge runtime, version 0.16.4. WasmEdge

R55

Space Data Network. Held-out accuracy of SDN-AP: the frozen plan, windows, truth and results of the experiment sdn-gp-starlink-ffm-match. Repository

R56

NASA Johnson Space Center, Flight Operations Directorate. ISS trajectory data, CCSDS orbit ephemeris messages.

R57

Noll, C. E. (2010). The Crustal Dynamics Data Information System: a resource to support scientific analysis using space geodesy. Advances in Space Research 45(12), 1421–1440.

R58

Johnston, G., Riddell, A., Hausler, G. (2017). The International GNSS Service. In Teunissen, P. J. G., Montenbruck, O. (eds.), Springer Handbook of Global Navigation Satellite Systems, 967–982. Springer.

R59

Montenbruck, O., et al. (2017). The Multi-GNSS Experiment (MGEX) of the International GNSS Service (IGS): achievements, prospects and challenges. Advances in Space Research 59(7), 1671–1697.

R60

Acton, C. H. (1996). Ancillary data services of NASA's Navigation and Ancillary Information Facility. Planetary and Space Science 44(1), 65–70.

R61

Ricker, G. R., et al. (2015). Transiting Exoplanet Survey Satellite. Journal of Astronomical Telescopes, Instruments, and Systems 1(1), 014003.

R62

Giorgini, J. D., et al. (1996). JPL's on-line solar system data service. Bulletin of the American Astronomical Society 28(3), 1158.