A few days ago, whilst gazing upon the midnight sky, a totally justified question crawled into my mind.

*How many SpaceX rockets does it take to deorbit the Moon?*

It went straight into one of those places where you can’t really ignore it or shake it off, you just *have to know*.

I could have probably just calculated the required *Δv* to reduce the *periapsis* (“lowest” point in orbit or closest approach) of the Moon enough for it to enter the atmosphere, but motivated by my need of a kinematics – gravity framework for a long standing project I decided to simulate the orbit of the Moon around the Earth in Python and brute-force the required amount of thrust (or number of rockets).

So here I am, registering domain names, configuring webservers and editing css stylesheets with the hope that my work will be helpful or educational to some, in an attempt to pay my long due obligation of giving back to the internet.

But enough of that, lets get to the point.

In order to model orbital motion, we need to be able to compute the position of a given body at each time step. This depends on its velocity which in turn depends on its acceleration, which is computed from the forces acting on it using *Newton’s second law of motion*:

\( \sum F=ma \\\)

Where *F* are the acting forces, *m* the mass of the body and *a* the acceleration.

Our first task is to create a class of an orbiting (or not) body, since this is meant to be reusable and not just a one night stand.

class Body(object):
def __init__(self, Pos0, Mass0, V0, A0, R0, F0=[0,0]):
self.Pos = Pos0
self.Mass = Mass0
self.V = V0
self.A = A0
self.R = R0
self.F = F0
def __call__(self):
print self.Pos, self.V, self.A, self.F |

class Body(object):
def __init__(self, Pos0, Mass0, V0, A0, R0, F0=[0,0]):
self.Pos = Pos0
self.Mass = Mass0
self.V = V0
self.A = A0
self.R = R0
self.F = F0
def __call__(self):
print self.Pos, self.V, self.A, self.F

The attributes of the object consist of its position vector, mass, velocity vector, acceleration vector, radius and force vector. All the vectors are two dimensional since we *assume that the Sun does not exist* for the scope of this problem, so we don’t have to bother with orbital inclinations and another dimension (who needs the Sun and a third dimension anyway?).

When the object is called, some information about its state gets printed.

In order for a mass *m* to orbit around a central point, it has to have a sufficient velocity *v* orthogonal to the centripetal force *F*_{g}:

At first, we need the gravity force vector, \(\vec{F_{g}}\), the magnitude of which is given to us by *Newton’s law of universal gravitation*:

\( |\vec{F_{g}}|=G\frac{m_1 m_2}{r^{2}} \\\)

Where *G* is the universal gravitational constant (not to be confused with *g*), *m*_{1,2} the masses of the two bodies and *r* the distance between them. *G* and *m*_{1,2} are known so we need to compute *r*. To do so, we shall utilize the *Pythagorean theorem*:

\( r= \sqrt{(x_1 -x_2)^2 + (y_1- y_2)^2} \\\)

Where *x*_{1,2}, *y*_{1,2} are the Cartesian coordinates of the two points (bodies). Now we can construct the functions *DistBetween* and *FgravMag:*

def DistBetween(body1, body2):
return math.sqrt((body1.Pos[0]-body2.Pos[0])**2+(body1.Pos[1]-body2.Pos[1])**2)
def FgravMag(body1, body2):
return G*body1.Mass*body2.Mass/(DistBetween(body1, body2)*10**3)**2 |

def DistBetween(body1, body2):
return math.sqrt((body1.Pos[0]-body2.Pos[0])**2+(body1.Pos[1]-body2.Pos[1])**2)
def FgravMag(body1, body2):
return G*body1.Mass*body2.Mass/(DistBetween(body1, body2)*10**3)**2

Note that the output of *DistBetween* is multiplied by 10^{3} in *FgravMag* because the distance unit of choice is Km but we need the force in Newtons.

In order to turn our gravity force magnitude to a force vector, we need a unit vector (one that has a magnitude of 1) \(\vec{u}\) that points from body1 to body2. This is exactly what we are going to get if we grab a body1 centered vector that points to body2 by subtracting their coordinates and divide that by the distance between them:

\( \vec{D}=(x_2 -x_1,y_2- y_1) \\ \vec{u}=\frac{\vec{D}}{|\vec{D}|} \\\)

Which is handled by the prosaically named function *PosUnitDirVct*:

def PosUnitDirVct(body1, body2):
return numpy.divide([body2.Pos[0]-body1.Pos[0], body2.Pos[1]-body1.Pos[1]], DistBetween(body1, body2)) |

def PosUnitDirVct(body1, body2):
return numpy.divide([body2.Pos[0]-body1.Pos[0], body2.Pos[1]-body1.Pos[1]], DistBetween(body1, body2))

By multiplying this direction vector with the magnitude of the gravitational pull we get the gravity force vector we sought for:

\(\vec{F_{g}} = \vec{u}|\vec{F_{g}}| \\\)

And the equivalent code:

def FgravVct(body1, body2):
return numpy.multiply(PosUnitDirVct(body1, body2), FgravMag(body1, body2)) |

def FgravVct(body1, body2):
return numpy.multiply(PosUnitDirVct(body1, body2), FgravMag(body1, body2))

We are half way there!

To reduce the periapsis of the Moon enough, in order to reach the atmosphere or crash into the earth, we need to apply a force in the *retrograde* direction, which means opposite to the direction of its velocity. This task is made a hell of a lot easier by the fact that the Moon is *tidally locked* to the Earth (orbital frequency same as rotational frequency, or why there is an album called *The Dark Side of the Moon [sic]*), so we can just park our boosters nose down on the “leading” point of the Moon that is always facing to the direction of its velocity *(prograde)* and fire them up at the apoapsis. To do so, we first employ a velocity direction vector (prograde vector) formulated in an orderly fashion, which we will later invert, resulting in a retrograde direction vector. By multiplying this new vector with the thrust of a Falcon9 v1.1 first stage booster (in Newtons!) and the count N of the boosters, we get our deorbiting force vector \( \vec{F_r} \):

def VelUnitDirVct(body):
return numpy.divide(body.V, math.sqrt(body.V[0]**2+body.V[1]**2)) # math.sqrt(V[0]**2+V[1]**2) is equivalent to numpy.linalg.norm(V)
...
Fr = numpy.multiply(VelUnitDirVct(Moon),-1) * F9thrust * N |

def VelUnitDirVct(body):
return numpy.divide(body.V, math.sqrt(body.V[0]**2+body.V[1]**2)) # math.sqrt(V[0]**2+V[1]**2) is equivalent to numpy.linalg.norm(V)
...
Fr = numpy.multiply(VelUnitDirVct(Moon),-1) * F9thrust * N

Now we have both of our acting forces and we can proceed with the computation of the position at each step. Since we will be working with discrete time progression in the form of a modest *for* loop, this is only a matter of integrating the acceleration given to us by Newton’s second law of motion twice, without forgetting to add its previous velocity/position:

\(\vec{a}_t=\frac{\sum{}\vec{F}_t}{m}

\\\\

\vec{v}_{t+n} = \int_{t}^{t+n} \vec{a_t}dt + \vec{v}_{t}

\\\\

\vec{x}_{t+n} = \int_{t}^{t+n} \vec{v_t}dt + \vec{x}_{t}\\\)

Where *x* the position, *t* is a given point in time and *n* our time step.

We can implement that directly as methods of our objects:

def updateAcc(self):
self.A = numpy.divide(self.F, self.Mass)/10**3
def updateVel(self, time):
self.V = [self.V[0] + self.A[0]*time, self.V[1] + self.A[1]*time]
def updatePos(self, time):
self.Pos = [self.Pos[0] + self.V[0]*time, self.Pos[1] + self.V[1]*time]
def Move(self, time):
self.updateAcc()
self.updateVel(time)
self.updatePos(time) |

def updateAcc(self):
self.A = numpy.divide(self.F, self.Mass)/10**3
def updateVel(self, time):
self.V = [self.V[0] + self.A[0]*time, self.V[1] + self.A[1]*time]
def updatePos(self, time):
self.Pos = [self.Pos[0] + self.V[0]*time, self.Pos[1] + self.V[1]*time]
def Move(self, time):
self.updateAcc()
self.updateVel(time)
self.updatePos(time)

I should also note at this point that this method of time discretization (Solving the differential equations using *Euler method*) has an inherent error when working with time varying quantities, directly related to the size of our step.

*All models are wrong, but some are useful.* – George Box

Finally it’s time to glue everything together. Okay but how do we do that?

First we should define the two bodies:

Earth = Body([0,0], 5.97237*10**24, [0,0], [0,0], 6371)
Moon = Body([0,405400], 7.341*10**22, [0,0], [0,0], 1737) |

Earth = Body([0,0], 5.97237*10**24, [0,0], [0,0], 6371)
Moon = Body([0,405400], 7.341*10**22, [0,0], [0,0], 1737)

We place the Earth at 0,0 (how arrogant of us, casually placing the Earth) and the apoapsis of the Moon at the top of our coordinate system. Mass is in Kg and distance in Km. Every position of the Moon during our simulation will be placed in a list named *Moonposplot (yeah I know).*

The first stage of a Falcon9 has a burn duration of 180 seconds, so after that the only acting force will be Earth’s gravity:

for t in numpy.linspace(0,endtime,endtime/step):
if t >= 180:
Moon.F = FgravVct(Moon, Earth)
else:
Fr = numpy.multiply(VelUnitDirVct(Moon),-1) * F9thrust * N
Moon.F = numpy.add(FgravVct(Moon, Earth), Fr)
Moon.Move(step)
Moonposplot.append(Moon.Pos) |

for t in numpy.linspace(0,endtime,endtime/step):
if t >= 180:
Moon.F = FgravVct(Moon, Earth)
else:
Fr = numpy.multiply(VelUnitDirVct(Moon),-1) * F9thrust * N
Moon.F = numpy.add(FgravVct(Moon, Earth), Fr)
Moon.Move(step)
Moonposplot.append(Moon.Pos)

And using matplotlib to output an animated plot:

Ecircle = plt.Circle((Earth.Pos[0], Earth.Pos[1]), Earth.R, color='g')
Mcircle = plt.Circle((Moon.Pos[0], Moon.Pos[1]), Moon.R, color='grey')
Esphere = plt.Circle((Earth.Pos[0], Earth.Pos[1]), Earth.R+10000, color=(0,0,0.5,0.1))
Thsphere = plt.Circle((Earth.Pos[0], Earth.Pos[1]), Earth.R+700, color=(0,0,0.6,0.15))
Msphere = plt.Circle((Earth.Pos[0], Earth.Pos[1]), Earth.R+80, color=(0,0,0.7,0.2))
Ssphere = plt.Circle((Earth.Pos[0], Earth.Pos[1]), Earth.R+50, color=(0,0,0.8,0.25))
Trsphere = plt.Circle((Earth.Pos[0], Earth.Pos[1]), Earth.R+12, color=(0,0,0.9,0.3))
fig = plt.figure()
ax = plt.axes(xlim=(-415.4*10**3, 415.4*10**3),ylim=(-415.4*10**3, 415.4*10**3))
line, = ax.plot(map(list, zip(*Moonposplot))[0],map(list, zip(*Moonposplot))[1], lw=0.5)
def init():
ax.add_artist(Esphere)
ax.add_artist(Thsphere)
ax.add_artist(Msphere)
ax.add_artist(Ssphere)
ax.add_artist(Trsphere)
ax.add_artist(Ecircle)
ax.add_artist(Mcircle)
return Mcircle, Esphere, Thsphere, Msphere, Ssphere, Trsphere, Ecircle, line
def anim(i):
Mcircle.center = (Moonposplot[i][0], Moonposplot[i][1])
return Mcircle, Esphere, Thsphere, Msphere, Ssphere, Trsphere, Ecircle, line
ani = animation.FuncAnimation(fig, anim, init_func=init, frames=int(endtime/step), interval=30, blit=True)
plt.show() |

Ecircle = plt.Circle((Earth.Pos[0], Earth.Pos[1]), Earth.R, color='g')
Mcircle = plt.Circle((Moon.Pos[0], Moon.Pos[1]), Moon.R, color='grey')
Esphere = plt.Circle((Earth.Pos[0], Earth.Pos[1]), Earth.R+10000, color=(0,0,0.5,0.1))
Thsphere = plt.Circle((Earth.Pos[0], Earth.Pos[1]), Earth.R+700, color=(0,0,0.6,0.15))
Msphere = plt.Circle((Earth.Pos[0], Earth.Pos[1]), Earth.R+80, color=(0,0,0.7,0.2))
Ssphere = plt.Circle((Earth.Pos[0], Earth.Pos[1]), Earth.R+50, color=(0,0,0.8,0.25))
Trsphere = plt.Circle((Earth.Pos[0], Earth.Pos[1]), Earth.R+12, color=(0,0,0.9,0.3))
fig = plt.figure()
ax = plt.axes(xlim=(-415.4*10**3, 415.4*10**3),ylim=(-415.4*10**3, 415.4*10**3))
line, = ax.plot(map(list, zip(*Moonposplot))[0],map(list, zip(*Moonposplot))[1], lw=0.5)
def init():
ax.add_artist(Esphere)
ax.add_artist(Thsphere)
ax.add_artist(Msphere)
ax.add_artist(Ssphere)
ax.add_artist(Trsphere)
ax.add_artist(Ecircle)
ax.add_artist(Mcircle)
return Mcircle, Esphere, Thsphere, Msphere, Ssphere, Trsphere, Ecircle, line
def anim(i):
Mcircle.center = (Moonposplot[i][0], Moonposplot[i][1])
return Mcircle, Esphere, Thsphere, Msphere, Ssphere, Trsphere, Ecircle, line
ani = animation.FuncAnimation(fig, anim, init_func=init, frames=int(endtime/step), interval=30, blit=True)
plt.show()

And let’s test that for a second without any boosters…

Oh no! A premature crash!

That happened because we did not specify an initial velocity *V0*, and it is as if a stationary Moon just *appeared* and fell straight into the Earth.

The power of modern search engines immediately gives us 0.964 Km/s as the answer to the question *“Moon minimum orbital velocity”* (Minimum because at t=0 we are at the apoapsis). Isn’t that great? This power is readily available to everyone, yet some people choose to use the internet just to post photos of their meals.

Anyway, let’s initialize the Moon in a meaningful way now:

Moon = Body([0,405400], 7.341*10**22, [0.964,0], [0,0], 1737) |

Moon = Body([0,405400], 7.341*10**22, [0.964,0], [0,0], 1737)

And we have achieved orbit!

A *well respected and accurate* verification can be performed using the cursor and the output of matplotlib’s figure to check out if the periapsis is where it should, which surprisingly enough it is! Another test that can be performed is playing with the duration of the simulation until exactly one revolution around the Earth is completed and then compare that to the orbital period of the Moon. Using a quick and dirty *for* loop to accomplish that, we verify that the found value is indeed really close to a *sidereal month*.

Getting the answer to the initial question is only a matter of applying a force that is a multiple of 5885 kN (Falcon9 1.1 stage 1 thrust) until the Moon reaches the exosphere! (Then orbital decay will take over.) This can be achieved with another *for* loop that modifies the multiplier *N* until the distance between the two bodies is sufficiently small.

…Well as it turns out we need *5*10*^{16} Falcon9 stage 1 boosters so that’s not happening any time soon.

Such a significant amount of boosters adds enough mass to our body to play an important role in our calculations. It is trivial to account for that, but if we decide to be *that* precise (let me remind you that *we ignored the Sun*) we would also need to take in account the reduction of the boosters’ mass, something we could do because SpaceX was kind enough to provide us with the *specific impulse* of Falcon9’s stage 1. So we shall deem our previous answer satisfactory.

By developing the program in such way, it is now trivial to demonstrate multi-body situations such as the two body problem (because *nothing is really stationary*):

or a *binary system*:

Well, I can now consider my question answered!

Thank you for reading this, and if you want to check out the complete code, it is on Github.

Todo: Implement a higher order symplectic & time-reversible integrator such as JANUS.