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.

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?