The aim of this paper is to explain why it is inherently dicult to predict the point of
reentry for a de-orbiting satellite. Ultimately, your explanation may be as simple as that
satellites travels an enormous distance over a short period of time, leaving room for big errors,
or it may involve a more subtle reason. Whatever the case, create a dierential equations
model of a satellite’s position that includes reasonable physical assumptions. Solve the model
equations numerically and use numerical experimentation (or even better, analysis) to
support your conclusions. Two to three single-spaced pages is an appropriate length.
Where to start? The equations of motion in vacuum are
m1x00 = ?
Gm1m2x
jxj3 ;
where G is the gravitational constant, m1 is the body’s mass and m2 is the mass of the
earth. We know that all orbits are conic sections and so you may as well assume that
x(t) = hx1(t); x2(t)i is a vector in the plane. Consider that the earth radius is around 6000
km and that the International Space Station circles the earth in 92 minutes. This suggests
expressing length in thousands of km (megameters, Mm) and time in minutes. Put in all
the units and make the appropriate cancelations to end up with the relevant dimensionless
equations. They take the form
x00
1 = f1(t; x1; x2; x0
1; x0
2); x00
2 = f2(t; x1; x2; x0
1; x0
2);
which you then convert in the usual fashion to a system of rst order equations by introducing
the new variables v1(t) = x0
1(t) and v2(t) = x0
2(t):
The atmosphere high above the earth’s surface is not completely vacuum and so satellites
eventually come out of orbit due to atmospheric drag, a process called orbital decay. (Solar
storms are another factor.) To account for atmospheric drag, add the drag force1
FD = ?1
2CDAjvjv; v = hx0
1; x0
2i;
to the force due to gravity in the equations of motion. Here CD = 2:2 is a good value
for the drag coecient, A is the cross-sectional area (hundreds of meters squared) and is
atmospheric density. For density ; we can use the approximation
(x1; x2) = 0 exp
R ? r
8:5 103m
; 0 = 1:3 kg m?3: (SEE NEXT PAGE!)
Here r =
p
x21
+ x22
is the distance of the satellite to earth’s center and R = 6:4 Mm
(megameters). This way, atmospheric density is equal to sea-level density 0 whenever
r = R and decays exponentially as r ! 1:
1https://www.wired.com/2013/04/air-resistance-in-orbit/
Corrigendum. The form of the prior expression for is adequate for altitudes up to
about 100 km but experimental models show that the decay rate actually changes at higher
altitudes. In fact, the rate of decay for higher altitudes is less than for lower altitudes. The
following modied expression accounts for the two dierent rates and is a more accurate
representation of the experimental values. Let a = r ? R be altitude (still a function of x1
and x2) and take
(a) = 0 exp(f(a))
where the piecewise function
f(a) =
(
? 137:71a; 0 a < 0:14
? 19:28 ? 17:63(a ? 0:14); 0:14 a < 1
accounts for the dierent rates.
A Matlab/Octave code that executes this is
z = rho0*exp(-137.71*a.*(a<0.14) + (-19.28-17.63*(a-0.14)).*(a>=0.14));
