It turns out Newton's laws aren't quite right (they're superseded by Einstein's General Theory of Relativity). But it took a long time to realize it, since they're very close to accurately describing the behaviour of things that aren't too big or too fast. I spent a lot of secondary school applying Newtonian mechanics. So I thought it would be interesting to use them to build simulation of the solar system and see how accurate the simulation is.
The first step was to obtain exact locations and velocities for the various bodies of the solar system. I obtained these using the HORIZONS system by NASA's Jet Propulsion Laboratory, at http://ssd.jpl.nasa.gov/horizons.cgi.
I recorded the location, velocity and mass for the 16 largest objects in the solar system: earth, moon, sun, mercury, venus, mars, jupiter and it's four largest moons (io, ganymede, europa, callisto), saturn, uranus, neptune and its moon triton, and pluto.
So let's import everything we'll need, and then get the initial state of our bodies at midnight on January 1st, 2018
import math, time, pickle, requests
from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML
from functools import reduce
import seaborn as sns
G = 6.67408*10**-11 # (in m^3 kg^-1 s^-2)
AU = 1.4960e11 # (in meters)
d = 86400 # in seconds
planet_masses = [
1.98854e30, #sun
3.302e23, #mercury
4.8685e24, #venus
5.97219e24, #earth
734.9e20, #moon
6.4185e23, #mars
1.89813e27, #jupiter
8.933e22, #io
479.7e22, #europa
1.482e23, #ganymede
1.076e23, #calisto
5.68319e26, #saturn
8.68103e25, #uranus
1.024e26, #Neptune
13455.3e22, #triton
1.307e22] #pluto. Ok it's not a planet and small
#NASA uses numbers to identify planets and their moons.
planet_codes = [10,199,299,399,301,499,599,501,502,503,504,699,799,899,801,999]
planet_names = ['sun', 'mercury', 'venus','earth', 'moon','mars', 'jupiter', 'io',
'europa','ganymede', 'calisto','saturn', 'uranus', 'neptune', 'triton', 'pluto']
def get_vectors(body, start, end, step):
# given a planetary body code and a time span, returns cartesian position
# and velocity coordinates.
base_url = 'https://ssd.jpl.nasa.gov/horizons_batch.cgi?batch=l&TABLE_TYPE=%27VECTORS%27&CSV_FORMAT=%27YES%27&CENTER=%27500@0%27&OUT_UNITS=%27AU-D%27&REF_PLANE=%27ECLIPTIC%27REF_SYSTEM=%27J2000%27&TP_TYPE=%27ABSOLUTE%27&ELEM_LABELS=%27YES%27CSV_FORMAT=%27YES%27&OBJ_DATA=%27YES%27&COMMAND=%27'
url = base_url + str(body) + '%27&START_TIME=%27' + start + '%2000%3A00%27&STOP_TIME=%27' + end + '%2000%3A00%27&STEP_SIZE=%27' + step + '%27'
r = requests.get(url)
return parse_text(r.text)
def parse_text(ntxt):
## parses the text given by HORIZONS and returns a list of position/velocity coordinates.
steps = ntxt.split('$$SOE')[1].strip().split('$$EOE')[0].strip().split('\n')
steps = [step.split(' ') for step in steps]
steps = [[s for s in step if s][4:10] for step in steps]
steps = [[float(s[:-1]) for s in step] for step in steps]
return steps
## Since I'm curious how close we can get to a whole year in advance:
start = '2018-01-01'
end = '2019-01-01'
step_size ='1h'
t0 = time.time()
initial_state = []
total_state = []
## the horizons data is in AU and AU/d, so we'll convert to SI units. We'll also
## add the planet masses, and combine the entire initial state into one list
def to_SI (data):
# NASA data comes in AU/d
ndata = data[:]
ndata[0] = data[0]*AU
ndata[1] = data[1]*AU
ndata[2] = data[2]*AU
ndata[3] = (data[3]*AU)/d
ndata[4] = (data[4]*AU)/d
ndata[5] = (data[5]*AU)/d
if len(data)>6:
ndata[6] = data[6]
return ndata
for i, body in enumerate(planet_codes):
print ('Working on %s.' %planet_names[i])
body_data = get_vectors(body,start,end,step_size) #the whole thing
total_state.append(body_data)
body_initial = body_data[0]
body_initial.append(planet_masses[i])
initial_state.extend(to_SI(body_initial))
t1 = time.time()
delta = t1 - t0
print ('Took %s seconds' %delta)
# Saving both
with open('1yr_2018_01_01.p', 'wb') as f:
pickle.dump(total_state,f)
with open('initial_2018_01_01.p', 'wb') as f:
pickle.dump(initial_state,f)
with open('1yr_2018_01_01.p', 'rb') as f:
total_state = pickle.load(f)
with open('initial_2018_01_01.p', 'rb') as f:
initial_state = pickle.load(f)
nasa_data = []
for planet in total_state:
planetx = []
planety = []
for timepoint in planet:
planetx.append(timepoint[0]*AU)
planety.append(timepoint[1]*AU)
nasa_data.append([planetx,planety])
The key to simulating the motions of bodies in solar system is the `scipy.integrate.odeint()` method. This function takes three arguments: (a) an initial state, (b) an array of variable values to evaluate the system at (in this case, these are instants of time) and (c) a function `ode()` that, given the values that characterize the state at one point returns the derivatives for those values.
Since we need the derivatives for all the values that characterize the solar system, initial_state encodes the state of the solar system in one long list of 16*7 = 112 items. (16 objects, 3 position and 3 velocity vectors + mass for each object.) The `planet_seperator()` function seperates this long list into individual lists for each body.
How do we write `ode()`? Well,the derivatives for the position vectors (x,y,z) are just the velocity vectors (vx, vy, vz). The acceleration vectors are calculated according to Newton's Law of universal gravitation. For each pair of bodies, o1 and o2, the component acceleration of o1 due to o2 is G (the gravitational constant) times the mass of o2 divided by the square of the distance between o1 and o2. The resultant acceleration is the sum of all the component accelerations.
def distance(x1,y1,z1,x2,y2,z2):
# returns the distance between two points
return math.sqrt((x2-x1)**2 + (y2-y1)**2 + (z2-z1)**2)
def component_acceleration(x1,y1,z1,x2,y2,z2,m2):
# given location of two planets and mass of planet 2 and returns
# component acceration of planet 1
dist = distance(x1,y1,z1, x2, y2, z2)
a = G*m2/dist**2 #magnitude of acceleration
# normalized direction vectors:
dx = (x2-x1)/dist
dy = (y2-y1)/dist
dz = (z2-z1)/dist
ax = dx*a
ay = dy*a
az = dz*a
return [ax,ay,az]
def planet_separator(X):
# converts long list of params into list of lists for each body
numplanets = int(len(X)/7)
ans = [X[i*7:i*7+7] for i in range(numplanets)]
return ans
def ode(params, t):
# given the position and velocity of a planet, return list of same length with derivatives
deltas = []
planets = planet_separator(params)
for i, planet in enumerate(planets):
#for each body, calculate force due to gravity on it due to each other body
x,y,z = planet[:3]
comp_accelerations = [component_acceleration(x,y,z, other_planet[0], other_planet[1], other_planet[2], other_planet[6]) for j, other_planet in enumerate(planets) if j != i]
resultant_a = reduce(lambda a1,a2: [a1[0]+a2[0], a1[1]+a2[1], a1[2]+a2[2]], comp_accelerations)
deltas += [planet[3], planet[4], planet[5], resultant_a[0], resultant_a[1],resultant_a[2], 0]
return deltas
def integrator (initial, timespan):
start = time.time()
solution = integrate.odeint(ode,initial,timespan)
end = time.time()
delta = end - start
print('Integrating took %s seconds' %delta)
return solution
Let's test this out over a one year period. Ode() works by numerical integration: instead of a continuous process, it breaks the time into chunks and makes an estimation of the derivates during that chunk. Obviously the smaller the chunks, the more accurate this will be. But the smaller the chunks, the more of them and the more processing it'll take. I'll break up the time into 8670 one hour chunks over a year.
# Now let's integrate a full year
one_year = 60*60*24*365
hours_per_year = int(one_year/(60*60))
one_hour = 60*60
timespan = np.linspace(0, one_year, hours_per_year+1)
solution_1yr = integrator(initial_state, timespan)
with open('solar_system_integrated_one_year_n1.p', 'wb') as f:
pickle.dump(solution_1yr, f)
with open('solar_system_integrated_one_year_n1.p', 'rb') as f:
solution_1yr = pickle.load(f)
Ok, let's check that against the actual values obtained from NASA's HORIZONS. For clarity we'll just look at Earth's position:
calculated_x = solution_1yr[-1][3*7]
calculated_y = solution_1yr[-1][3*7+1]
nasa_x = nasa_data[3][0][-1]
nasa_y = nasa_data[3][1][-1]
x_error = ((calculated_x-nasa_x)/AU)*100 #errors as a % of AU
y_error = ((calculated_y-nasa_y)/AU)*100
print (x_error, y_error)
So the simulation is off by around 0.01%. That's not awful. To illustrate, let's plot both the JPL and the simulated Earth together:
fig, ax = plt.subplots(figsize=(10,10))
maxx = max(nasa_data[3][0])*1.1
minx = min(nasa_data[3][0])*1.1
maxy = max(nasa_data[3][1])*1.1
miny = min(nasa_data[3][1])*1.1
ax.set_xlim(minx,maxx)
ax.set_xlabel('Meters from Sun')
ax.set_ylabel('Meters from Sun')
ax.set_ylim(miny,maxy)
ax.plot(solution_1yr[:,[3*7]], solution_1yr[:,3*7+1], 'b', label="Earth Calculated Position")
ax.plot(nasa_data[3][0], nasa_data[3][1], 'g', label = "Actual position of Earth from NASA")
ax.legend()
earth_calc, = ax.plot([], [], 'bo', ms=1)
earth_nasa, = ax.plot([], [], 'go', ms=1)
earths = [earth_calc, earth_nasa]
def init():
earth_calc.set_data([], [])
earth_nasa.set_data([], [])
return earth_calc, earth_nasa
def animate(i):
global ax, fig
ms = 6
earth_calc.set_data(solution_1yr[i*40][3*7], solution_1yr[i*40][3*7+1])
earth_nasa.set_data(nasa_data[3][0][i*40],nasa_data[3][1][i*40])
earth_calc.set_markersize(ms)
earth_nasa.set_markersize(ms)
return earths
ani = animation.FuncAnimation(fig, animate, frames=219,
interval=20, blit=True, init_func=init)
HTML(ani.to_html5_video())
They're close enough that you can't tell there are two plots. Just for fun, let's go ahead and plot all the planets:
sol = solution_1yr
num_planets = 16
def find_plot_bounds(sol):
#Find x,y max/min loc for each planet
xs,ys,zs = [], [], []
numplanets = int(len(sol[0])/7)
for step in sol:
for i in range(16):
xs.append(step[i*7])
ys.append(step[i*7+1])
zs.append(step[i*7+2])
return {'minx': min(xs)*1.1, 'maxx': max(xs)*1.1, 'miny': min(ys)*1.1, 'maxy': max(ys)*1.1}
bounds = find_plot_bounds(sol)
minx = bounds['minx']*1.3
maxx = bounds['maxx']*1.3
miny = bounds['miny']
maxy = bounds['maxy']*2
fig, ax = plt.subplots(figsize=(10,10))
ax.set_xlim(minx,maxx)
ax.set_ylim(miny,maxy)
planets = []
annotations = []
lines = []
os = [(-2,2), (2,0),(1,2),(-2,0),(0,1),(1,-1),(0,-1),
(-1,-1),(-1,-2),(-0.75,3),(0,2),(1,1),(1,-1),(-1,1),(-1,-1),(1,1)]
for i in range(16):
# ax.plot(nasa_data[i][0], nasa_data[i][1], linewidth=0.5)
line, = ax.plot(sol[:,i*7],sol[:,i*7+1],label=planet_names[i], linewidth=0.5)
lines.append(line)
if i not in [4,7,8,9,10,14]:
name = planet_names[i]
xy = (sol[3000,i*7],sol[3000,i*7+1])
x_os = os[i][0]*0.4e12
y_os = os[i][1]*0.4e12
xytext = (sol[3000,i*7]+x_os,sol[3000,i*7+1]+y_os)
ax.annotate(name,xy=xy, xytext=xytext,
arrowprops={'arrowstyle': '-'}, va='center')
for i in range(16):
planet, = ax.plot([],[],'go', ms=2, label=planet_names[i])
planets.append(planet)
def init():
return planets
def animate(i):
global ax, fig
ms = 6
for j,p in enumerate(planets):
thisx = solution_1yr[i*40][j*7]
thisy = solution_1yr[i*40][j*7+1]
p.set_data(thisx,thisy)
p.set_markersize(2)
return planets
ani = animation.FuncAnimation(fig, animate, frames=219,
interval=20, blit=True, init_func=init)
HTML(ani.to_html5_video())