How to perform a soft landing on Mars

So, let’s say that you’re a spacecraft. Or in one. And you just arrived at the Martian Sphere of Influence, with a 15000km circular orbit. How do you land safely? You have no parachute, and a carton of eggs that must not break (a friendly reminder that there are probably no hens on Mars). A good idea would perhaps be to lower your periapsis enough to touch the Martian atmosphere, and utilize the resistance of the air to lower your apoapsis without spending any fuel.
How much fuel will I need? How long will it take? Are you sure that my eggs won’t break?
All these questions can be answered by formulating a model for the system and performing a simulation!
We start with the Martian atmosphere. A quick duckduckgoing (equivalent to googling but less invasive) points us to this NASA model, which we happily implement, only to plot it and realize that there is an error in the model, and the temperature goes below -273.15°C at an altitude of around 110km (which is probably not a very good idea). Disclaimer: this is an outdated model based on some quite old curve fits, and NASA provides a multitude of other more realistic and accurate options but since (thankfully) nobody assigned me the task of safely landing some eggs on Mars this model will be deemed sufficient. The side effect of the below absolute zero temperature is that the density asymptotically spikes on both sides around this point:
Atmospheric model error
Luckily, another duckduckgoing indicates that the folks at the University of Colorado have already devised a solution that does not include implementing a more advanced model, a polynomial curve fitting (the document disappeared as I was writing this post, please use the wayback machine). After patching up our atmospheric function, we move on to the spacecraft.
There is no upper limit for the complexity of the spacecraft dynamics model, and the lower limit is dictated by one’s conscience. So let’s just consider that the spacecraft is an infinitely small point in 2 dimensional space with fixed aerodynamic characteristics! Just like real life.
By conveniently placing Mars in the origin of the Cartesian coordinate system (0,0), the distance of the spacecraft from the center of Mars (yes, it is a perfect sphere with uniformly distributed mass in case you were wondering) becomes:

\(r = \sqrt{x^2+y^2}\\\)

where x,y are the coordinates of the spacecraft

In similar fashion we calculate the absolute velocity:

\(v = \sqrt{v_x^2+v_y^2}\\\)

where \(v_x,v_y\) are the velocity components of the spacecraft

The magnitude of the force of gravity is calculated using Newton’s law of universal gravitation:

\(F_{grav} = \mu \frac{m}{r^2}\\\)

where \(\mu\) is the standard gravitational parameter for Mars (pick the right units!) and m the mass of the spacecraft

Since we have the speed of the spacecraft and its altitude (remember to subtract Mars’ radius from r) to plug it in the atmosphere model, we can calculate the magnitudes of the aerodynamic forces:

\(
F_D =\frac{1}{2} \rho v^2 C_D A \\
F_L =\frac{1}{2} \rho v^2 C_L A \\
C_L = \frac{L}{D} C_D \\
\)

where \(F_D\) is the magnitude of the drag force, \(F_L\) the magnitude of the lift force, \(\rho\) the density derived by the atmospheric model, \(C_D\) the drag coefficient, \(C_L\) the lift coefficient, \(A\) the aerodynamically active area, \(v\) the magnitude of the velocity and \(\frac{L}{D}\) the lift to drag ratio

And yes, all the aerodynamic characteristics are constant since we’re aiming for maximum realism, we’ve been over that already!
I assume that you are familiar with drag and lift, \(\frac{L}{D}\) is just a design parameter specifying the relationship between the two in a capsule. Together with the capsule’s ballistic coefficient, they are the two major factors that dictate the trajectory of the spacecraft when aerobraking, or in general whenever it is flying (or falling) through the atmosphere.
In order to not break our precious eggs, we of course need to strap a huge ass rocket motor to the capsule to decelerate.
Every time the engine is producing thrust, it does so by spewing massive quantities of exhaust gasses with furious anger from the nozzle. To do so, it has to expend significant amounts of fuel and oxidizer (in the case of liquid rocket motors), which of course means that the spacecraft’s mass is dropping. The rate of change of the mass can be calculated:

\(\dot{m} = \frac{F_{thrust}}{g_0 I_{sp}}\\\)

where \(\dot{m}\) the rate of change of the spacecraft mass, \(F_{thrust}\) the magnitude of the thrust, \(g_0\) is the standard gravity and \(I_{sp}\) the specific impulse of the motor

Specific impulse is just a measure of how efficiently the rocket uses up its propellant, which in reality directly relates to the velocity of the exhaust gasses. Higher exhaust gas velocities produce more momentum with the same mass (amount of propellant used) which means more furious anger with less weight loss.

Now that we have the magnitudes of everything, we need to add the directions. Do do so we can just normalize the vector components to get the unit vector that’s pointing in the right direction. For example with velocity:

\(
v_{1x} = \frac{v_x} {v} \\
v_{1y} = \frac{v_y} {v} \\\)

where the resulting \(v_{1x}, v_{1y}\) are the components of the normalized velocity vector

This vector will be used to set the direction of most of the forces, namely the drag, lift and motor thrust. The drag and motor thrust act in the opposite direction of this vector. Drag does that by definition, and we pick that for the thrust because we want to decelerate (the term for that is that we are burning in the retrograde direction). For the lift, there are two cases. Since we considered the aerodynamic characteristics constant, it can either point to the sky (lift-up) or to the ground (lift-down). A configuration with downward lift results in more abrupt deceleration and higher forces, as the spacecraft dives faster in the thicker layers of the lower atmosphere. On the other hand, when the lift vector is pointing upwards the vertical descent of the spacecraft is slower and it is subjected to less stress, with the disadvantage of a longer aerobraking duration. To minimize the stress on our zygotes, we pick the lift-up configuration and calculate the direction as a right angle to the velocity vector:

\(
l_{1x} = v_{1x}cos(\theta_{L})-v_{1y}sin(\theta_{L}) \\
l_{1y} = v_{1x}sin(\theta_{L})+v_{1y}cos(\theta_{L}) \\\)

where \(\theta_{L}\) is the angle of the lift vector relative to the velocity

Finally, Newton’s second law of motion can be used to calculate the acceleration of the spacecraft from the acting forces, leading to the complete differential equation system. As mentioned, the thrust and drag are always applied in the retrograde direction (hence the use of the velocity unit vector), and the gravitational attraction to the direction of the origin (position unit vector).

\(
a_x = \frac{\frac{x}{r} F_{grav} – v_{1x} F_D – v_{1x} F_{thrust} + l_{1x} F_L}{m} \\
a_y = \frac{\frac{y}{r} F_{grav} – v_{1y} F_D – v_{1y} F_{thrust} + l_{1y} F_L}{m} \\
\)

And the complete system:

\(
\dot{y} = \left[ \begin{array}{c} v_x \\ v_y \\ a_x \\ a_y \\ \dot{m} \end{array} \right], y_0 = \left[ \begin{array}{c} x_0 \\ y_0 \\ v_{x0} \\ v_{y0} \\ m_0 \end{array} \right], \quad
\dot{y}(t) = f(t,y(t),p), \quad y(t_0) = y_0 \\
\)

where \(y_0\) are the initial states for each variable and p is the parameters vector (coefficients and the such)

To calculate the initial states, we place the spacecraft at \(x_0 = 15000, y_0 = 0\) (just a point in the desired 15000km circular orbit, on the x axis for simplicity) and calculate the initial velocity using

\(
v_{apo0} = \sqrt{\mu_{Mars} \frac{2}{R_{A0}} – \frac{1}{a}} \\
\)

where \(R_{A0}\) is the apoapsis radius (15000km), and a the semi-major axis of the orbit (again 15000km because it is circular)

To control when the motor is active or not, a simple if checks if the current time is within the pre-specified burn intervals.

After implementing everything and putting them together, the system is integrated using MATLAB’s (I’m sorry open source ecosystem, I truly am) built in explicit Dormand-Prince 5-4 (ode45). While dynamic orbital systems are already stiff and the addition of thrust intervals with abrupt changes make the system even stiffer, the simulation duration is small enough to allow a very small tolerance (hence time step) to ensure sufficiently correct integration.

After educated guesses on the spacecraft’s characteristics based on a scaled SpaceX crew dragon capsule and sufficient fiddling around with the manipulated variables, we get the following output (of course some prettifying stuff are added such as highlight the powered parts of the trajectory and a gradient of the logarithm of the atmospheric density):

Mars landing simulation
Mars landing simulation

I played enough with the numbers to achieve a touchdown velocity of 0.295 m/s (THE EGGS DID NOT BREAK) with a fuel mass fraction of 0.718. This means that 72% of the initial mass is fuel that got burnt in the process.

Next step? Strap an optimizer (most likely fmincon) to find the optimal solution with the smallest mass fraction. Because there are discrete phases in the simulation (thus discontinuities) this is not as simple as I originally anticipated.

That’s it! Thank you for reading!

Hanging the NASA space travel posters in an orderly fashion

I assume that the reader is familiar with the NASA JPL Visions of the Future posters. Because all cool people know about them and only cool people read my blog.

Well after years of anticipation, I finally found an opportunity to actually put them up on my wall -as I recently relocated- and I was struck by an unforeseen problem.
How does one arrange fifteen colored posters with different dominant colors?
It really seems that the theme of this blog is answering unavailing (at least according to some) questions.

Initially we need the images to be sorted, in our case the thumbnails of the posters, which we shall retrieve with a rather crude but fully functional (as of 6/9/17) chain of commands:

curl -s https://www.jpl.nasa.gov/visions-of-the-future/ | grep "small.jpg' alt" | cut -d "'" -f 2 | xargs -I @ wget https://www.jpl.nasa.gov/visions-of-the-future/@

Ok now what?
As it turns out color sorting is not that trivial. After sufficient research, my decision was to sort by hue instead of playing with color distance (maybe this should be given a try as well).

Step one is to create a loop that iterates over jpgs and pngs in a given (by command line argument) folder using os.walk:

for (path, dirs, files) in os.walk(sys.argv[1]):
    for fname in files:
        if fname.endswith(('.jpg', '.png')):

Now using pillow (non-dead version of Python Imaging Library) we can open our images and get the RGB channels values of the pixels with getdata.
By dividing the sum with the len of the pixels, we have our median poster color, followed by a division with 255 since our conversion library prefers a range from 0-1. (thrice, once for each color channel)

meanRGB = [(float(sum(list(im.getdata(channel)))) / len(list(im.getdata(channel))))/255 for channel in range(3)]

Next, we convert this lonely pixel to HSV using colorsys, and fill a list with that and the filename of the image:

meanHSV.append([colorsys.rgb_to_hsv(meanRGB[0],meanRGB[1],meanRGB[2]),fname])

Finally, when the loops are done we can sort this list and print out the hue sorted filenames:

meanHSV.sort()
print '\n'.join([q[1] for q in meanHSV])

In case you were wondering:

~/py/visionssort$ ./sort.py .
kepler186f-small.jpg
mars-small.jpg
51pegasib-small.jpg
kepler16b-small.jpg
venus-small.jpg
titan-small.jpg
grand_tour-small.jpg
earth-small.jpg
superearth-small.jpg
ceres-small.jpg
europa-small.jpg
nightlife-small.jpg
jupiter-small.jpg
trappist-small.jpg
enceladus-small.jpg

Since the right arrangement for my wall is 3×5, this pattern was preferred (i.e. I couldn’t come up with a better solution):

And the final result is:



I can now hang my posters in peace, knowing that they are arranged.
Full code

Sunset – Sunrise monitor on a raspberrypi

Imagine that for some reason it is imperative that you know at any given time when the Sun will rise or set next. (Maybe you are a vampire.)

For whatever reasons – who am I to judge – you need that, you can easily achieve it with a raspberry pi, a 2×16 character LCD screen and python!

And as is almost always the case with every python mini project, we are going to utilize a powerful library to achieve our puny goals. In this case Ephem.

Ephem is a quote: Python package for performing high-precision astronomy computations.

I will not be covering how to interface the LCD module with the raspberrypi or how to use the Adafruit library as those topics have been exhausted already.

Let’s start by importing the libraries we are going to need (the number of which will inevitably increase):

import ephem
from datetime import datetime
import Adafruit_CharLCD as LCD

Some pin definitions required by the LCD library which I don’t know why but I didn’t pass them directly to the function:

lcd_rs = 25
lcd_en = 24
lcd_d4 = 23
lcd_d5 = 17
lcd_d6 = 27
lcd_d7 = 22
lcd_backlight = 4
lcd_columns = 16
lcd_rows = 2
 
lcd = LCD.Adafruit_CharLCD(lcd_rs, lcd_en, lcd_d4, lcd_d5, lcd_d6, lcd_d7,lcd_columns, lcd_rows, lcd_backlight)

In order to cram as much information as possible in such a small screen, a nice way to delimit our data is by using custom characters. The hd44780 controller is nice like that.

Using online tools that let you draw the pixels and acquire the byte sequence (such as this one) and some Bob Ross videos, we can create the characters we need.

I made one that is supposedly a Sun on the horizon (yes I know.) and another one for the temperature:

lcd.create_char(0, [0,0,14,27,17,31,31,31])
lcd.create_char(1, [0,4,10,10,10,17,31,14])

Which translate to those two characters (yes it does definitely not look like a Sun we’ve been over that):

And now we can proceed in computing the time difference from now to the next sunset / sunrise.

First, we create a Sun and an observer:

sun = ephem.Sun()
 
obs = ephem.Observer()
obs.lat = '-49.55'
obs.long= '69.83'
obs.date = datetime.utcnow()

And now by taking advantage of the wonder that is modern day computing, we can get what wanted in just two lines of code:

delta = min(obs.next_setting(sun).datetime(),obs.next_rising(sun).datetime()) - datetime.utcnow()
riseset = str(delta.seconds//3600) + ":" + str(delta.seconds%3600//60).zfill(2)

This works by calling two functions from the ephem library that return the next sunset and sunrise, finding out which one is closer to the present using min, and subtracting the present daytime in order to get a timedelta object (the time difference).

Then, with a series of mods and divs that should be self-explanatory to get the hours and minutes from the seconds in a pretty way, we have our value!

I also happen to have a dirty bash script that scrapes the current outside temperature from a website for my city, so I’ll throw that in as well (more imports):

temp = subprocess.check_output("./gettemp.sh").rstrip()

And finally send everything to the LCD:

lcd.clear()
lcd.message('\x00 ' + riseset + '\n' + '\x01 ' + temp)

‘\x00’ is the sunset character because in the function create_char we specified the first memory position for it to be stored, and ‘\x01’ the little thermometer.

Now throw everything after the coordinates assignment in a loop with a time delay and Bob’s your uncle:

Thanks for reading!
Code

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

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

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 Fg:

Orbit graph
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), m1,2 the masses of the two bodies and r the distance between them. G and m1,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 x1,2, y1,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

Note that the output of DistBetween is multiplied by 103 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:

direcion vector between two points

\( \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))

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))

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

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)

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)

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)

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()

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)

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.
Periapsis in exosphere

…Well as it turns out we need 5*1016 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):
Two body demonstration
or a binary system:
Two body example2

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.

Hello world?

This is a blog! Yay.

The purpose of this blog is.. I guess it only exists because I wanted to try out installing WP on a virtual machine. Debian wiki instructions unclear, dick stuck in /dev/null.

If I continue to publish content here, it will probably be about electronics, linux, random thoughts etc. How uncommon eh?