Two weeks ago, the USAF and DDS put up a Hack-A-Sat space hacking challenge. I only had time to check out a handful of the problems, one of which was Where’s the Sat that asked people to identify a satellite by location and time and then predict where it would be at another time. It turns out that with the right tools, such as Skyfield, the problem is straightforward (check, e.g., here). But the solution is not very satisfying unless one understands what goes on inside the black box… so if you are curious, read on.
To first approximation, Earth is spherically symmetric and satellites are pointlike, in which case orbits are ellipses with the center of Earth at one of the foci (recall Keplerian orbits). E.g., if you picture the “standard” ellipse in the x-y plane with semiaxes a and b along the x and y axes, and shift Earth to the origin, then the center of the ellipse will be at (-a e, 0), while the point of closest approach (perigee) at P = (a (1-e), 0), where e is the eccentricity.
Sats in the challenge were given as a list of NORAD two-line elements (TLEs). These contain the parameters of the orbit and where on the orbit the sat was at the given time, specifically, the
The list in the challenge had sats with ballpark 1.5-hour orbital periods, zipping around at roughly 6700 km from the center of Earth - such as the International Space Station (ISS).
With the help of some geometry (if you are rusty, check Wikipedia), you can reconstruct the ellipse and the position on it corresponding to the TLE. The result is not terrible but still misses by a lot; e.g., for the ISS it was off by 30 km already at epoch and by more than 150 km within ~5 orbits after epoch. Orbits are not ellipses because Earth’s gravitational field is not perfectly spherical plus there is atmospheric drag, so Kepler’s solution no longer applies (other bodies such as the Moon and the Sun also exert forces but for orbits near Earth those are less relevant). Short portions of orbits can always be considered elliptical but parameters change along the path. Moreover, one cannot use TLE data directly even near epoch because TLEs contain averaged quantities, not instantaneous values (thanks to Dr. Kelso at CelesTrack for pointing that out to me).
The most important correction is that Earth’s gravitational potential U (r,ϑ) varies with latitude. People usually expand the potential over Legendre polynomials in sin ϑ, and keep the first couple terms, say, up to the fourth-order term - the coefficients of which (J2, J3, J4) are determined by measurements (J1 is zero by construction). Gradients of the terms beyond zeroth order then give additional force contributions that one can plug into Newton’s a = F /m and find sat trajectories via integrating the differential equations numerically (about freshman physics at most places, though a high-order integrator does not hurt). However, this approach still requires initial conditions, so one cannot sidestep extracting instantaneous positions and velocities from TLEs… somehow.
TLEs are designed to work with a certain class of orbital propagation models. A popular one of these, applicable to near-Earth trajectories, is Simplified General Perturbations version 4 (SGP4). The model stems from approximate integration of Newton’s equations in convenient coordinates for Keplerian motion, which yields tractable math for small J2, J3, J4 corrections to gravity and a weak drag force. The simplified math really mattered in the old days when computers were slow. The algorithm with the resulting equations is publicly available (in Python, too, by now) but for the actual physics papers
behind it you will likely need access to a research library. TLEs contain averaged orbit parameters that are convenient input quantities for SGP codes, and the codes then give position and velocity at any time, both before and after epoch (with errors that grow with the time difference relative to epoch).
It seems like we are done: we can either use SGP4 for everything, or take initial conditions at epoch from SGP4 and then solve Newton’s equation numerically. However, if we try the latter for the ISS example above, we get a discrepancy of about 4 km per day from epoch compared to SGP4 predictions. This is worrisome because SGP4 is known to work very well in practice, and the task in Where’s the Sat is to predict positions at about 23 days prior to epoch, at which point our real physics calculation is off by nearly 100 km.
Earth’s atmosphere dilutes with altitude but even a small air resistance can matter for sat motion.
SGP4 uses quadratic
F = -const ρ|v|v drag, written as the acceleration a = -(B* /R) (ρ/ρ0) |v| v,
where B* is specified in the TLE, while R is the radius of Earth.
The density ratio is of a power-law form [(q0 - s)/(r - s)]**tau (see Lane’s paper),
where q0, s, tau are fit parameters. Specific values can be obtained by inspecting the code:
(for r < R + 158 km, SGP4 uses a different set). With air resistance included in Newton’s equations, the deviation from SGP4 reduces by half for the ISS to 2 km/day, which is good but not good enough.
Is 2 km/day simply the inherent accuracy of SGP4, or does the model have some other physics? Actually, no and no - the reason for the discrepancy is more subtle.
SGP4 is a good approximation for an extended time around epoch. But its answers are only approximate, and therein lies the problem when we want to use SGP4 positions and velocities as initial conditions for Newton’s equations. In fact, with just about one part per million adjustments to initial conditions from SGP4 (couple tens of meters in position and few millimeters per second in velocity), one can obtain Newtonian trajectories that stay close to the SGP4 solutions for a long time.
A reasonably intelligent way to find these adjustments is to systematically alter each component of the initial position and velocity by a small amount (six degrees of freedom), and then find the linear combination of those changes that brings the position predicted by Newton’s equations in agreement with SGP4 at two points in time, say, about one period before and one period after epoch. An even better way would be to minimize the sum of deviations at multiple points along one orbit (I have not explored that). Either way involves some straightforward linear algebra. From the adjusted initial conditions, the Newtonian trajectory agrees with SGP4 to about 100 meters/10 days from epoch for the ISS when both nonspherical gravity corrections and drag are included.
This translates to roughly 200-meter accuracy at twenty-some days (~300 orbits) before epoch, which would have been plenty for the Where’s the Sat challenge (as far as I can tell, the challenge was coded to tolerate up to 1% error, i.e., couple tens of kilometers).
For accurate satellite orbit predictions you must go beyond Keplerian orbits and take into account the nonspherical gravitational field of Earth and also atmospheric drag. Straightforward physics, in principle, but getting precise-enough initial conditions for Newton’s equations from TLEs is tricky. Thus, solving Where’s the Sat in the available timeframe was basically impossible without applying SGP4 or a similar propagation model.