# Orbital Trajectory simulator

I have written a simple program to do trajectory simulation in the Earth-Moon system, it still has a long way to go I am working on making it more class oriented and am looking into implementing a better system than just the straight time step integration method.

I want to focus on having it more class-oriented. I am unsure where to draw a good line as to how much stuff should be class/function driven and how much should be in the main block. I am also looking for feedback on how to make it faster, as right now when I have the timestamp low it takes quite a while to run. In addition I am trying to find a good way to implement planets, as the way I have them now is not it, I could make a planet class, and initialize all the ones I am going to use, but I would rather just import/use them.

Here is some example output:

Console: (old output)

Days remaining: 28.84, (96.14% left)
Tics per second: 49286 est time remaining: 00:00:50
Days remaining: 27.69, (92.28% left)
Tics per second: 49555 est time remaining: 00:00:48
Days remaining: 26.53, (88.43% left)
Tics per second: 49653 est time remaining: 00:00:46
Days remaining: 25.37, (84.57% left)
Tics per second: 49770 est time remaining: 00:00:44
Days remaining: 24.21, (80.71% left)
Tics per second: 49777 est time remaining: 00:00:42
Days remaining: 23.06, (76.85% left)
Tics per second: 49789 est time remaining: 00:00:40
Days remaining: 21.90, (72.99% left)
Tics per second: 49844 est time remaining: 00:00:37
Days remaining: 20.74, (69.14% left)
Tics per second: 49901 est time remaining: 00:00:35
Days remaining: 19.58, (65.28% left)
Tics per second: 49796 est time remaining: 00:00:33
Days remaining: 18.43, (61.42% left)
Tics per second: 49725 est time remaining: 00:00:32
Days remaining: 17.27, (57.56% left)
Tics per second: 49682 est time remaining: 00:00:30
Days remaining: 16.11, (53.70% left)
Tics per second: 49632 est time remaining: 00:00:28
Days remaining: 14.95, (49.85% left)
Tics per second: 49645 est time remaining: 00:00:26
Days remaining: 13.80, (45.99% left)
Tics per second: 49553 est time remaining: 00:00:24
Days remaining: 12.64, (42.13% left)
Tics per second: 49541 est time remaining: 00:00:22
Days remaining: 11.48, (38.27% left)
Tics per second: 49519 est time remaining: 00:00:20
Days remaining: 10.32, (34.41% left)
Tics per second: 49523 est time remaining: 00:00:18
Days remaining: 9.17, (30.56% left)
Tics per second: 49538 est time remaining: 00:00:15
Days remaining: 8.01, (26.70% left)
Tics per second: 49512 est time remaining: 00:00:13
Days remaining: 6.85, (22.84% left)
Tics per second: 49532 est time remaining: 00:00:11
Days remaining: 5.69, (18.98% left)
Tics per second: 49501 est time remaining: 00:00:09
Days remaining: 4.54, (15.12% left)
Tics per second: 49518 est time remaining: 00:00:07
Days remaining: 3.38, (11.27% left)
Tics per second: 49516 est time remaining: 00:00:05
Days remaining: 2.22, (7.41% left)
Tics per second: 49533 est time remaining: 00:00:03
Days remaining: 1.06, (3.55% left)
Tics per second: 49526 est time remaining: 00:00:01
Total Tics per second: 49509.11
Total time: 00:00:52


Plot:

The blue line is the craft (whose orbit has been altered by the moon), and the Green line is the moons position over time.

(Also here in a git repo, please excuse the fluff I am trying to get started with readthedocs and sphinx, but it's not going so well).

### Orbital.py (main bit)

import time as t
import math
import sys

from astropy.time import Time
import numpy as np

#Custum objects live in sub directory
sys.path.append('./Objects/')
# Import Objects
from Objects import Craft
import moon
import earth

# Generate output for plotting a sphere
def drawSphere(xCenter, yCenter, zCenter, r):
# Draw sphere
u, v = np.mgrid[0:2 * np.pi:20j, 0:np.pi:10j]
x = np.cos(u) * np.sin(v)
y = np.sin(u) * np.sin(v)
z = np.cos(v)
# Shift and scale sphere
x = r * x + xCenter
y = r * y + yCenter
z = r * z + zCenter
return (x, y, z)

def plot(ship,planets):
"""3d plots earth/moon/ship interaction"""
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Initilize plot
fig = plt.figure()
ax.set_xlabel('X Km')
ax.set_ylabel('Y Km')
ax.set_zlabel('Z Km')
ax.set_xlim3d(-500000, 500000)
ax.set_ylim3d(-500000, 500000)
ax.set_zlim3d(-500000, 500000)

ax.plot(xs=ship.hist[0][0::10],
ys=ship.hist[1][0::10],
zs=ship.hist[2][0::10],
zdir='z', label='ys=0, zdir=z')

# Plot planet trajectory
for planet in planets:
ax.plot(xs=moon.hist[0],
ys=moon.hist[1],
zs=moon.hist[2],
zdir='z', label='ys=0, zdir=z')

# Plot Earth (plot is moon position relative to earth)
# also plotting to scale
(xs, ys, zs) = drawSphere(0, 0, 0, 6367.4447)
ax.plot_wireframe(xs, ys, zs, color="r")

plt.show()

def sim(startTime, endTime, step, ship, planets):
"""Runs orbital simulation given ship and planet objects as well as start/stop times"""

# Caculate moon planet update rate (1/10th as often as the craft)
plan_step = int(math.ceil(((endTime - startTime) / step)/100))

# Initilize positions of planets
for planet in planets:
planet.getRelPos(startTime)
planet.log()

start = t.time()
for i, time in enumerate(np.arange(startTime, endTime, step)):
# Every plan_step update the position estmation
if (i % plan_step == 0):
for planet in planets:
planet.getRelPos(time)
planet.log()

# Update craft_vol
for planet in planets:
ship.force_g(planet.mass, planet.pos[0], planet.pos[1], planet.pos[2])
ship.update()

# Log the position of the ship every 1000 steps
if (i % 1000) == 0:
# Append the position to the lists
ship.log()

# Print status update every 100,000 steps
if (i % 100000) == 0:
if t.time()-start > 1:
print "Days remaining: {0:.2f}, ({1:.2f}% left)".format((endTime - time), ((1-((time - startTime) / (endTime - startTime)))*100))
# Caculate estmated time remaining
# Total tics
tot_tics = int((endTime - startTime) / step)
# Tics per second till now
tics_s = int(math.ceil(i/(t.time()-start)))
# Estmated remaining seconds
sec_rem = (tot_tics-i)/tics_s
m, s = divmod(sec_rem, 60)
h, m = divmod(m, 60)
print "Tics per second: {0} est time remaining: {1:0>2}:{2:0>2}:{3:0>2}".format(tics_s, h, m, s)
print "Total Tics per second: {0:.2f}".format(((endTime - startTime) / step)/(t.time()-start))
tot_s = t.time()-start
m, s = divmod(tot_s, 60)
h, m = divmod(m, 60)
print "Total time: {0:0>2.0f}:{1:0>2.0f}:{2:0>2.0f}".format(h, m, s)

if __name__ == "__main__":

# Delta time for simulations (in days)
del_t = np.longdouble(1.0) / (24.0 * 60.0 * 60.0)

ship = Craft(del_t, x=35786, y=1, z=1, v_x=0, v_y=4.5, v_z=0, mass=12)
planets = [earth,moon]

# Initilize simulation time
start_time = Time('2015-09-10T00:00:00')
end_time = Time('2015-10-10T00:00:00')

# Initilize start date/time (Julian)
time = start_time.jd

sim(start_time.jd, end_time.jd, del_t, ship, planets)
plot(ship,[moon])


### ./Orbital/Objects.py

import math
import numpy as np
from jplephem.spk import SPK
import os

class Craft(object):

def __init__(self, delt_t, x=0.0, y=0.0, z=0.0,
v_x=0.0, v_y=0.0, v_z=0.0, mass=0):
# Pos is in km
self.pos = np.array([np.longdouble(x),
np.longdouble(y),
np.longdouble(z)])
# Vol is in km/s
self.vol = np.array([np.longdouble(v_x),
np.longdouble(v_y),
np.longdouble(v_z)])
# Mass is in kg
self.mass = mass
# Delta_t is in days
self.del_t = np.longdouble(delt_t) * 86400

# Initilize non input vars
# Initilize history of positions
self.hist = [[np.longdouble(x)],
[np.longdouble(y)],
[np.longdouble(z)]]
# Initilize force array
self.f = np.array([np.longdouble(0),
np.longdouble(0),
np.longdouble(0)])
# Initilize constants
# Gravitational constant (6.674*10-11 N*(m/kg)2)
self.G_const = np.longdouble(0.00000000006674)

def force_g(self, body_mass, body_x=0.0, body_y=0.0, body_z=0.0):
# Caculate x, y, z distances
xdist = (self.pos[0] - body_x)
ydist = (self.pos[1] - body_y)
zdist = (self.pos[2] - body_z)
# Caculate vector distance
d = math.sqrt((xdist**2) + (ydist**2) + (zdist**2)) * 1000
# Caculate comman force of gravity
g_com = ((self.G_const * body_mass * self.mass) / (d**3))

# Update forces array
self.f -= np.array([g_com * (xdist * 1000),
g_com * (ydist * 1000),
g_com * (zdist * 1000)])

# Function to step simulation based on force caculations
def update(self):
# Step velocity
self.vol += np.array([(((self.f[0] / self.mass) * self.del_t)) / 1000,
(((self.f[1] / self.mass) * self.del_t)) / 1000,
(((self.f[2] / self.mass) * self.del_t)) / 1000])
# Step position
self.pos += np.array([(self.del_t * self.vol[0]),
(self.del_t * self.vol[1]),
(self.del_t * self.vol[2])])
# Reset force profile
self.f = np.array([np.longdouble(0),
np.longdouble(0),
np.longdouble(0)])

def VV_update(self):
"""Parameters:
r is a numpy array giving the current position vector
v is a numpy array giving the current velocity vector
dt is a float value giving the length of the integration time step
a is a function which takes x as a parameter and returns the acceleration vector as an array"""
r_new = r + v*dt + a(r)*dt**2/2
v_new = v + (a(r) + a(r_new))/2 * dt

def dist(self, body_x=0.0, body_y=0.0, body_z=0.0):
return math.sqrt(((self.pos[0] - body_x)**2) +
((self.pos[1] - body_y)**2) +
((self.pos[2] - body_z)**2))

def log(self):
self.hist[0].append(self.pos[0])
self.hist[1].append(self.pos[1])
self.hist[2].append(self.pos[2])


### ./Objects/earth.py

import numpy as np
from jplephem.spk import SPK
import os

if not os.path.isfile('de430.bsp'):
kernel = SPK.open('de430.bsp')
# Initilize mass
mass = np.longdouble(5972198600000000000000000)
# Initilize history of positions
hist = [[],[],[]]
pos = np.array([np.longdouble(0),
np.longdouble(0),
np.longdouble(0)])

def getPos(time):
"""Returns the earths position relative to the solar system barycentere
"""
return kernel[3, 399].compute(time)
def getRelPos(time):
global pos
"""Returns relitive position of the earth (which is the origin
in an earth moon system)"""
pos = np.array([np.longdouble(0),
np.longdouble(0),
np.longdouble(0)])
return pos
def log():
global pos
"""log current position"""
hist[0].append(pos[0])
hist[1].append(pos[1])
hist[2].append(pos[2])


### ./Objects/moon.py

import numpy as np
from jplephem.spk import SPK
import os

if not os.path.isfile('de430.bsp'):
kernel = SPK.open('de430.bsp')
# Initilize mass
mass = np.longdouble(7.34767309 * 10**22)
# Initilize history of positions
hist = [[],[],[]]
pos = []

def getPos(time):
"""Returns the moons position relative to the solar system barycentere"""
global pos
pos = kernel[3, 301].compute(time)
return pos
def getRelPos(time):
"""Returns relitive position of the moon (relative to the earth)"""
global pos
pos = kernel[3, 301].compute(time) - kernel[3, 399].compute(time)
return np.array([np.longdouble(pos[0]),
np.longdouble(pos[1]),
np.longdouble(pos[2])])
def log():
"""log current position"""
global pos
global hist
hist[0].append(pos[0])
hist[1].append(pos[1])
hist[2].append(pos[2])


(Note: if you would like to run this code you will need the de430.bsp file you can get that here)

My updated code can be found here.

# The Earth and the Moon

You have a very strange pattern here where you're treating moon.py as if it's a class. You're using global variables to make the scoping work but it just looks funny. Instead, inside moon.py you should have a Moon class. This should contain relevant attributes and methods.

Everything in your moon file is relevant to a moon 'object', so let's make it a class. You have your imports as normal, but now the mass and kernel initialisation should be in __init__:

class CelestialBody(object):

def __init__(self, mass, position=None):
if not os.path.isfile('de430.bsp'):
self.kernel = SPK.open('de430.bsp')

self.mass = np.longdouble(mass)
self.hist = [[],[],[]]
if position is not None:
self.pos = position
else:
self.pos = []


And no I'm not happy with the name CelestialBody. Seems long and silly, if you have a better name for Planet/Moon then it'd be good to use there instead. Even though you shouldn't need it anywhere else.

Now this can actually apply to both the Moon and Earth, so it's a much more useful class. You have very similar functions too. log is identical in both, so that's not an issue. getPos is practically the same except that you have different values passed to the kernel. Since you also have different code for getRelPos, it seems best to have separate Earth and Moon classes inherit from CelestialBody:

class Moon(CelestialBody):

KERNAL_CONSTANT = 301

def __init__(self, mass):
super(Moon, self).__init__(mass)


KERNAL_CONSTANT is a bad name but I'm not familiar with what its relevance is. The point is, this value can be passed to the getPos method from CelestialBody:

def getPos(time):
"""Returns the position relative to the solar system barycentere"""
return kernel[3, self.KERNAL_CONSTANT].compute(time)


And of course the Earth class can have it's own constant defined. Importantly you can now set them to have entirely different getRelPos methods in their own classes.

And of course the Earth class can have it's own constant defined. Importantly you can now set them to have entirely i getRelPos methods in their own classes. Here's how I'd rewrite the two classes into one file:

import numpy as np
import os

from jplephem.spk import SPK

class CelestialBody(object):

def __init__(self, mass, position=None):
if not os.path.isfile('de430.bsp'):
self.kernel = SPK.open('de430.bsp')

self.mass = np.longdouble(mass)
self.hist = [[],[],[]]
if position is not None:
self.pos = position
else:
self.pos = []

def getPos(time):
"""Returns the position relative to the solar system barycentere"""
return kernel[3, self.KERNAL_CONSTANT].compute(time)

class Earth(CelestialBody):

KERNAL_CONSTANT = 399

def __init__(self, mass):
super(Earth, self).__init__(mass)

def getRelPos(time):
"""Returns relitive position of the earth (which is the origin
in an earth moon system)"""

location = np.longdouble(0)
self.pos = np.array([location, location, location])
return self.pos

class Moon(CelestialBody):

KERNAL_CONSTANT = 301

def __init__(self, mass):
super(Moon, self).__init__(mass)

def getRelPos(time):
"""Returns relitive position of the moon (relative to the earth)"""
self.pos = self.getPos(time) - kernel[3, 399].compute(time)
return np.array([np.longdouble(self.pos[0]),
np.longdouble(self.pos[1]),
np.longdouble(self.pos[2])])


Also shouldn't Moon's getRelPos actually call Earth's getPos rather than just reuse the code from there? It seems much cleaner to do it that way. And if Moon is so dependant on Earth then it should probably have Earth as an attribute to be easily accessed there.

And of course now you can get your classes like this:

import bodies

earth = bodies.Earth(5972198600000000000000000)
moon = bodies.Moon(7.34767309 * 10**22)


# Craft

In Craft you have a lot of noise inside __init__. Your pos is in km comments are littered through the function, but would be much more readable and useful implemented as a docstring:

def __init__(self, delt_t, x=0.0, y=0.0, z=0.0,
v_x=0.0, v_y=0.0, v_z=0.0, mass=0):
"""Units:

Position: km
Volume: km/s
Mass: kg
Delta_t: days
"""

self.pos = np.array([np.longdouble(x),
np.longdouble(y),
np.longdouble(z)])
self.vol = np.array([np.longdouble(v_x),
np.longdouble(v_y),
np.longdouble(v_z)])
self.mass = mass
self.del_t = np.longdouble(delt_t) * 86400


You also don't need to note that you're initialising values, that's what the whole point of __init__ is. The only other comment I'd keep is a note about what the gravitational constant is, but as a constant it should be at the class level, not inside __init__:

class Craft(object):

# Gravitational constant (6.674*10-11 N*(m/kg)2)
GRAVITY = np.longdouble(0.00000000006674)


This can still be referred to with self.GRAVITY as every instance will have this attribute. But it will also allow you to use Craft.GRAVITY which may at times be useful. You also no longer need to recreate this value for each craft.

Also in Python, names in UPPER_SNAKE_CASE denote that they refer to a constant value, so instead of g_const, GRAVITY denoting a constant value seems clearer to me. GRAVITATIONAL_CONSTANT would be extra clear, if also extra letters.

I have a question here too, you instantiate the history as a list of 3 lists for each co-ordinate. Is that necessary for the plotting you do later? It would make more sense to me if you were storing a list of three item tuples. Here's an example (with arbitrary data):

[
(231, 23, 12),
(21, 121, 89),
(451, 41, 21),
(23, 1, 7),
]


Instead yours would look like this:

[
[231, 21, 451, 23],
[23, 121, 41, 1],
[12, 89, 21, 7],
]


Which is harder to read. If it's possible to switch to my suggestion, I think you should as then you can access a single co-ordinate, and from that take out the three dimension values. I believe this would also be faster as you only need to index one long list, and then just have a short 3 item tuple.

In force_g you have more unnecessary comments. We can see that you're calculating x, y and z. And we could easily recognise that you're calculating the vector if you changed the name from d to vector_distance. That makes your equation below easier too. It is probably easier to just use a comment rather than try make a manageable name for common force of gravity, especially since it's so localised.

Lastly Update forces array should really be told in a docstring, so that it's easier to understand what the function is for, and seems to indicate that the function name should be changed:

def update_forces(self, body_mass, body_x=0.0, body_y=0.0, body_z=0.0):
"""Update forces array

Calculates x, y, z and the common force of gravity.
Using those calculations the forces array is updated."""

xdist = (self.pos[0] - body_x)
ydist = (self.pos[1] - body_y)
zdist = (self.pos[2] - body_z)
vector_distance = math.sqrt((xdist ** 2) + (ydist ** 2) + (zdist ** 2)) * 1000

# Caculate common force of gravity
g_com = (self.GRAVITY * body_mass * self.mass) / (d ** 3)

self.f -= np.array([g_com * (xdist * 1000),
g_com * (ydist * 1000),
g_com * (zdist * 1000)])


VV_update has some serious readability problems. I am not following. What's the point of r_new and v_new? As far as I can tell you just create them and then neither return them nor set them as attributes. Effectly doing nothing. Perhaps there's some quirk of numpy that means values are being modified in place here, but I don't see it. Aside from that, you don't tell me what VV is or what updates it. There's explanations for the very confusing names, but those explanations don't actually indicate what the use of this function is. Having clarification on how a function works is good for a docstring, but if there isn't an explanation of the function's use then you've got problems.

Also, you could think about putting the Craft class in the same file as the others. They don't have any overlap, so there's no particular benefit aside from less files to import from and being slightly neater as long as the file doesn't get too long.

# Main

In your imports you import time as t. I'm not fond of that alias as single letter variables are usually used for throwaway values in brief contexts (like a for loop). You may have done this to avoid clashing with Time. Though you technically don't need to (since capitalisation means they refer to two separate values) it's a good idea to distinguish them. I'd recommend instead doing something that clarifies the specific use of time here, and from time import time as current_time. Now you've distinguished them while simultaneously showing exactly what you need time for.

You also put some imports statements inside of functions. This is a bad idea, as it's common style to put all of them at the start so it's easier for people to read and know what you're using in the script.

Odd as it may seem, you don't need parentheses for Python's multiple assignment and it's actually confusing to use them. Just remove the parentheses here:

xs, ys, zs = drawSphere(0, 0, 0, 6367.4447)


That's how people usually write it, so it's clearer to follow.

• I'll probably right more since I didn't touch the craft in your objects file. If you have particular concerns or questions about what I said here, let me know and I can expand on them when I add more. – SuperBiasedMan Dec 2 '15 at 17:55
• That is terrific, I was struggling as to how to implement the differences properly, I would love a general review of both the object's class and my main code in regards to style, format, odd things as you noted earlier, also I am looking for a good way to move the update function away from the Object class, and into it's on class? I am unsure of the best way to do that, the idea is to have different ways of updating, some slower and more accurate and some faster and less accurate, but again thanks for the amazing feedback! – Mark Omo Dec 2 '15 at 18:23
• I went ahead and added CelestialBody, Earth, and Moon to Objects.py, and when I tried to use them I got TypeError: must be type, not classobj on the super function, I am trying to resolve it but I am having trouble – Mark Omo Dec 2 '15 at 21:19
• @MarkOmo: Glad to help! I'll probably have time to review more tomorrow and answer about the update function. And ah, that error is due to using Python 2, I should have noted that if you're using Python 2 you need to get CelestialBody to inherit from object for everything to work properly. I fixed my example code if you want to take a look at the changes. – SuperBiasedMan Dec 2 '15 at 21:30
• @MarkOmo I added more thoughts, I found your Craft object to be pretty solid so I had more general thoughts. Did you mean you wanted to move the update away from Craft? Because I don't think you should, it fits with OOP where it is as a self reliant method. It's easy to say craft.update() and let it adjust its own values without needing any input. – SuperBiasedMan Dec 4 '15 at 11:49

Let us start with some code review, before improving the performance, and finish off with a comparision against the original code.

## Code review

As a gesture to us reviewing your code, it would have been nice to know which modules where external to a default installation, as I had to plunder a little to install the astropy and jplephem module. However I got thing up and running, and then I looked at code and found the following:

• Naming is still not entirely snake_case – According to PEP8 the recommendation is to use snake_case, and not abbreviations. Your code breaks both of these, with names like pos, f, endTime, getRelPos(), and so on
• Function should be named accordingly to what they do – Your function getRelPos() does not get anything, but it does calculate a new position. In other words, update_position(time) would be a better name, as it clearly conveys what is happening.
• Avoid import's in functions – It's not considered good style to import modules within functions, like you do in plot(). It hides the dependencies, and can be time consuming if the function is called often.
• Add spaces around operators – Instead of self.f-=np.array(...), open it up using self.f -= np.array(...)
• Use enough parenthesis in formulas – Use enough, but not too many. Some places you can rely on operator precedence to avoid a lot of parenthesis clogging up, like math.sqrt(xdist**2 + ydist**2 + zdist**2), there is no need to have the power expression in parenthesis.
• Try to avoid costly functions in top-level of modules – In both earth.py and moon.py you do the SPK.open('de430.bsp') which is >100 MB file. This should be moved into main code, and the top level code of the modules should be within a initialization functions.

## Performance issues

There are various issues related to performance which we can address:

• Ideally simplify the de430.bsp – I don't know how this large file is handled, but besides having it open in just one place (at main), you should look into whether it can be trimmed down to something more relevant for you, and of smaller size. You only seem to use a small part of it, which could possibly be extracted somehow.
• Do you really need numpy.arrays and numpy.longdouble? – I kind of question if you really need to use this, or if you could have benefitted from using ordinary arrays and standard ints. Especially when I see code where you use the gravity constant with only 3 decimals, 6.674, then using numpy.longdouble's seems inadequate and imprecise...

However, I've left them as is in my refactored code below, but it could be worthwhile to look into simplifications of this.

• Simpler initialization of numpy arrays – A lot of places you reset your arrays to a triplet of numpy.longdouble zeros. This can easier, and faster, be done using np.zeros(3, dtype=np.longdouble)
• Simpler multiplication with numpy arrays – Instead of using distinct variables, you multiply the entire vector immediately using self.pos += np.dot(self.vol, self.del_t)
• Precalculate constant factors in heavy formulas – For example in update(self) you do (((self.f[0] / self.mass) * self.del_t)) / 1000 for each of the f's, when it suffices doing np.dot(self.f, self.VELOCITY_FACTOR), where self.VELOCITY_FACTOR = self.del_t / self.mass / 1000. This removes 2 of 3 multiplications, and that is a time saver.

I used kernprof to profile your code, and found that almost 90% of the time was divided on doing update() and force(), so therefore I spent some time simplifying these amongst others to reduce number of calculations executed.

• If possibly, avoid math.sqrt – Taking the square root is rather expensive, and if not needed try avoiding it. (You might also here get a precision issue as the is uses normal precision here?). In your code you do d = math.sqrt(...) followed by g_com = (...) / d**3. This can be simplifed to one operation of g_com = (...) / d**1.5. Still expensive, but not as expensive as the combination of sqrt and **3.
• Calculate and store only what is needed – In your code you execute the simulation loop 2592000 times with full calculation on the ship, and for each 10th you update the planets position. You store all planets positions, but you only store every 1000th position of the ship, and in addition you only plot every 10th of these again. This concludes to a lot of extraneous movement of the ship not to be used again.

This has an effect on the memory usage, and on execution time to calculate the positions. And in the end, you don't need it to plot the graph you want to plot. Similarily you update a non-moving earth, and store the positions of it, but in the end, you don't use the trajectory of the earth anyways. Waste of memory and time.

## Refactored code

I started of simplify the earths.py and moon.py into a single file, planets.py. I then continued updating iteratively the two other files. This resulted in the following files, with a few comments on each of them.

### File: planets.py

import numpy as np

class Planet(object):

def __init__(self, kernel, planet_id, relative_id, mass):
if kernel is None:
raise ValueError("Need a base kernel, like de430.bsp")

self._planet_id = planet_id
self._relative_id = relative_id
self._kernel = kernel
self.mass = mass
self.position = np.zeros(3, dtype=np.longdouble)

# Initialize position history
self.trajectory = [[], [], []]

def log_position(self):
"""Append current position to trajectory."""

if self._planet_id != self._relative_id:
for i in xrange(3):
self.trajectory[i].append(self.position[i])

def update_position(self, time):
"""Update the position according to a related position and self movement."""

if self._planet_id != self._relative_id:
self.position = (  self._kernel[3, self._planet_id].compute(time)
- self._kernel[3, self._relative_id].compute(time))


I moved the kernel into main code, so this is passed as a parameter. From reading your code I decided to use planet_id and relative_id to indicate the id of the objects. I used Planet althought that is slightly inadequate.

I started out with making the Earth and Moon to be distint classes, but found that the differences are so small, that there was no need for them to be separate classes. To allow for changing reference point, away from earth, the classes defaults to not storing/updating the positions if the planet_id and relative_id are the same. This since they would zero out, anyway.

### File: Objects.py

import math
import numpy as np

def profile(func):  return func

class Craft(object):

def __init__(self, delta_t, x=0.0, y=0.0, z=0.0,
v_x=0.0, v_y=0.0, v_z=0.0, mass=0):

# Pos is in km
self.position = np.array([np.longdouble(x),
np.longdouble(y),
np.longdouble(z)])

# Velocity is in km/s
self.velocity = np.array([np.longdouble(v_x),
np.longdouble(v_y),
np.longdouble(v_z)])
# Mass is in kg
self.mass = mass

# Delta_t is in days
self.delta_t = np.longdouble(delta_t) * 86400

# Initialize trajectory history of posisitons
self.trajectory = [[np.longdouble(x)],
[np.longdouble(y)],
[np.longdouble(z)]]

# Initialize force array
self.force = np.zeros(3, dtype=np.longdouble)

# Define some class constants, relative to self.mass
g_const = np.longdouble(0.00000000006674)
self.FORCE_FACTOR  = g_const * self.mass / 1000000
self.VELOCITY_FACTOR = self.delta_t / self.mass / 1000

@profile
def force_g(self, body):
"""Updates the force the body has on this craft."""

dist = self.position - body.position
gravity_component = ( self.FORCE_FACTOR * body.mass /
((dist[0]**2 + dist[1]**2 + dist[2]**2)**1.5) )

# Update forces array
self.force -= np.dot(dist, gravity_component)

@profile
def update(self):
"""Steps up simulation based on force calculations."""
# Step up velocity and position
self.velocity += np.dot(self.force, self.VELOCITY_FACTOR)
self.position += np.dot(self.velocity, self.delta_t)

# Reset force profile
self.force = np.zeros(3, dtype=np.longdouble)

def log(self):
for i in range(3):
self.trajectory[i].append(self.position[i])


The main change in this file is consolidation of factors into the precomputed VELOCITY_FACTOR and FORCE_FACTOR, in addition to using simpler and faster numpy alternatives.

The somewhat strange def profile(funct): return func is a little hack related to running the code with or without kernprof. The code runs as it stands, but to profile it you need to comment the def profile(). If you don't intend to profile it, you can remove the def profile() and the @profile decorators.

### File: Orbital.py

import time as t
import math
import sys
from jplephem.spk import SPK
import os
from astropy.time import Time
import numpy as np

#Custum objects live in sub directory
sys.path.append('./Objects/')

# Import Objects
from Objects import Craft
from planets import Planet

import matplotlib.pyplot as plt
# The following is needed for the projection='3d' to work in plot()
from mpl_toolkits.mplot3d import Axes3D

def profile(func): return func

# Generate output for plotting a sphere
def drawSphere(xCenter, yCenter, zCenter, r):
# Draw sphere
u, v = np.mgrid[0:2 * np.pi:20j, 0:np.pi:10j]
x = np.cos(u) * np.sin(v)
y = np.sin(u) * np.sin(v)
z = np.cos(v)
# Shift and scale sphere
x = r * x + xCenter
y = r * y + yCenter
z = r * z + zCenter
return (x, y, z)

def plot(ship,planets):
"""3d plots earth/moon/ship interaction"""

# Initialize plot
fig = plt.figure()
ax.set_xlabel('X Km')
ax.set_ylabel('Y Km')
ax.set_zlabel('Z Km')
ax.set_xlim3d(-500000, 500000)
ax.set_ylim3d(-500000, 500000)
ax.set_zlim3d(-500000, 500000)

ax.plot(xs=ship.trajectory[0][0::10],
ys=ship.trajectory[1][0::10],
zs=ship.trajectory[2][0::10],
zdir='z', label='ys=0, zdir=z')

# Plot planet trajectory
for planet in planets:
ax.plot(xs=planet.trajectory[0],
ys=planet.trajectory[1],
zs=planet.trajectory[2],
zdir='z', label='ys=0, zdir=z')

# Plot Earth (plot is moon position relative to earth)
# also plotting to scale
(xs, ys, zs) = drawSphere(0, 0, 0, 6367.4447)
ax.plot_wireframe(xs, ys, zs, color="r")

plt.show()

@profile
def run_simulation(start_time, end_time, step, ship, planets):
"""Runs orbital simulation given ship and planet objects as well as start/stop times"""

# Calculate moon & planet update rate (1/10th as often as the craft)
planet_step_rate = int(math.ceil(((end_time - start_time) / step)/100))

# Initialize positions of planets
for planet in planets:
planet.update_position(start_time)
planet.log_position()

start = t.time()
total_time = end_time - start_time
# Total tics
total_steps = int((end_time - start_time) / step)

print("Days in simulation: {:6.2f}, number of steps: {}".format(
total_time, total_steps))

for i, time in enumerate(np.arange(start_time, end_time, step)):

# Every planet_step_rate update the position estmation
if (i % planet_step_rate == 0):
for planet in planets[1:]:
planet.update_position(time)
planet.log_position()

# Update craft velocity and force vectors related to planets
for planet in planets:
ship.force_g(planet)

# Update ship position according to sum of affecting forces
ship.update()

# Log the position of the ship every 1000 steps
if (i % 1000) == 0:
# Append the position to the lists
ship.log()

# Print status update every 100,000 steps
if (i % 100000) == 0:
if t.time()-start > 1:
time_since_start = time - start_time

# Calculate estmated time remaining
remaining_days_of_simulation = end_time - time
remaining_percentage_of_simulation = (1 - time_since_start / total_time) * 100

# Tics per second till now
steps_pr_second = int(math.ceil( i / (t.time()-start)))
# Estmated remaining seconds
remaining_simulation_time_in_secs = (total_steps - i) / steps_pr_second

# Split remaining_simulation_time_in_secs into hours, minutes and second
minutes, seconds = divmod(remaining_simulation_time_in_secs, 60)
hours, minutes = divmod(minutes, 60)
print("  Simulation remaining: {:6.2f} days, {:5.2f}%.  Script statistics - steps/s: {:6}, est. time left: {:0>2}:{:0>2}:{:0>2}".format(
remaining_days_of_simulation, remaining_percentage_of_simulation,
steps_pr_second, hours, minutes, seconds))

end = t.time()
elapsed = end - start
print("Final steps per second: {0:.2f}".format((total_time / step)/elapsed))

minutes, seconds = divmod(elapsed, 60)
hours, minutes = divmod(minutes, 60)
print("Total running time: {:0>2.0f}:{:0>2.0f}:{:0>2.0f}".format(
hours, minutes, seconds))

def main():

# Delta time for simulations (in days)
delta_t = np.longdouble(1.0) / (24.0 * 60.0 * 60.0)

ship = Craft(delta_t, x=35786, y=1, z=1, v_x=0, v_y=4.5, v_z=0, mass=12)

# Initialize simulation time
simulation_start = Time('2015-09-10T00:00:00')
simulation_end = Time('2015-10-10T00:00:00')

if not os.path.isfile('de430.bsp'):
kernel = SPK.open('de430.bsp')

planets = [Planet(kernel, 399, 399, np.longdouble(5972198600000000000000000)),
Planet(kernel, 301, 399, np.longdouble(7.34767309 * 10**22))]

run_simulation(simulation_start.jd, simulation_end.jd, delta_t, ship, planets)
plot(ship, [planets[1]])

if __name__ == "__main__":
main()


I've renamed some of the variables, combined the log output, and introduced a main() to tidy it a little more. But more or less this code is as is, and could most likely benefit a lot from reducing number of steps and what gets stored in the arrays. I skipped updating of the position of earth, as that is the base of your system.

In total I've reduced the number of lines in your code, even though I've added a lot of blank lines here and there to enhance readability.

## Performance comparision

Your original code when I ran it on my computer gave this output:

Days in simulation:  30.00, number of steps: 2592000
Days remaining: 28.84, (96.14% left)
Tics per second: 15317 est time remaining: 00:02:42
Days remaining: 27.69, (92.28% left)
Tics per second: 15348 est time remaining: 00:02:35
Days remaining: 26.53, (88.43% left)
Tics per second: 15353 est time remaining: 00:02:29
Days remaining: 25.37, (84.57% left)
Tics per second: 15364 est time remaining: 00:02:22
...
Days remaining: 4.54, (15.12% left)
Tics per second: 15371 est time remaining: 00:00:25
Days remaining: 3.38, (11.27% left)
Tics per second: 15376 est time remaining: 00:00:18
Days remaining: 2.22, (7.41% left)
Tics per second: 15380 est time remaining: 00:00:12
Days remaining: 1.06, (3.55% left)
Tics per second: 15386 est time remaining: 00:00:05
Total Tics per second: 15386.26
Total time: 00:02:48


And my code gave the following:

Days in simulation:  30.00, number of steps: 2592000
Simulation remaining:  28.84 days, 96.14%.  Script statistics - steps/s:  21574, est. time left: 00:01:55
Simulation remaining:  27.69 days, 92.28%.  Script statistics - steps/s:  21610, est. time left: 00:01:50
Simulation remaining:  26.53 days, 88.43%.  Script statistics - steps/s:  21648, est. time left: 00:01:45
Simulation remaining:  25.37 days, 84.57%.  Script statistics - steps/s:  21640, est. time left: 00:01:41
...
Simulation remaining:   4.54 days, 15.12%.  Script statistics - steps/s:  21727, est. time left: 00:00:18
Simulation remaining:   3.38 days, 11.27%.  Script statistics - steps/s:  21729, est. time left: 00:00:13
Simulation remaining:   2.22 days,  7.41%.  Script statistics - steps/s:  21734, est. time left: 00:00:08
Simulation remaining:   1.06 days,  3.55%.  Script statistics - steps/s:  21739, est. time left: 00:00:04
Final steps per second: 21741.83
Total running time: 00:01:59


In short even without reducing the high number of steps to execute, I still shaved of almost a minute, and increased the steps pr second from 15400 steps/seconds to 21700 steps/second. Hopefully you should see the same percentage improvement on your computer, even though my number initially are somewhat slower than on your computer.

• I appreciate your feedback, I much prefer SuperBiasedMan's method for planets, but I love the performance improvements at first glance around 30%, I was going to make a case for camelCase but I realised I was inconstant in my naming regardless and switched over to underscores, I am working on incorporating everyone's suggestions, again thank you for your feedback! – Mark Omo Dec 3 '15 at 5:30
• @MarkOmo, You fail to realise that in structure both mine and SuperBiasedMan's solution are quite close. I've left out the getPos() as you didn't use them, but they could easily be readded. You could also make explicit classes for the earth and the moon, but what have you gained? – holroy Dec 3 '15 at 14:47
• down the road I plan on having the ability to simulate more planets and their associated bodies, as well as providing other functionality, that makes defining bodies in their own classes attractive thus allowing me to just add planets without having to reference there associated constants and to use different ephemeris files for better accuracy/different time spans. – Mark Omo Dec 3 '15 at 20:02
• @MarkOmo, I'm not sure if making each and every new planetary object a new class is in the true spirit of OOP, but it can be done in both solutions. More importantly, you can achieve having different ephemeris files and/or constants whilst still using one file as well. It's just a matter of instantiation. – holroy Dec 3 '15 at 20:04
• what would you suggest? my hope is to just be able to call/initlilize Objects.earth/mars/luna/ect as oposed to having to provide ephemeris referances or mass or diamater is there a better way thay you would recommend? – Mark Omo Dec 3 '15 at 20:08

Here is the code modified for running a given number of times, and as it turns out only around 3000 times is adequate to get a similar figure. How many times you need to run it does depend on the velocity of the different object and turning rate (like when ship sling shots around the moon or earth).

The timings gets quite a lot better, of course, when reducing from 2.5 millions steps to 3000 steps...

The planets.py and Objects.py files are as in my other answer.

### File: Orbital.py

import time as t
import math
import sys
from jplephem.spk import SPK
import os
from astropy.time import Time
import numpy as np

#Custum objects live in sub directory
sys.path.append('./Objects/')

# Import Objects
from Objects import Craft
from planets import Planet

import matplotlib.pyplot as plt
# The following is needed for the projection='3d' to work in plot()
from mpl_toolkits.mplot3d import Axes3D

def profile(func): return func

# Generate output for plotting a sphere
def drawSphere(xCenter, yCenter, zCenter, r):
# Draw sphere
u, v = np.mgrid[0:2 * np.pi:20j, 0:np.pi:10j]
x = np.cos(u) * np.sin(v)
y = np.sin(u) * np.sin(v)
z = np.cos(v)
# Shift and scale sphere
x = r * x + xCenter
y = r * y + yCenter
z = r * z + zCenter
return (x, y, z)

def plot(ship,planets):
"""3d plots earth/moon/ship interaction"""

# Initialize plot
fig = plt.figure()
ax.set_xlabel('X Km')
ax.set_ylabel('Y Km')
ax.set_zlabel('Z Km')
ax.set_xlim3d(-500000, 500000)
ax.set_ylim3d(-500000, 500000)
ax.set_zlim3d(-500000, 500000)

ax.plot(xs=ship.trajectory[0][0::10],
ys=ship.trajectory[1][0::10],
zs=ship.trajectory[2][0::10],
zdir='z', label='ys=0, zdir=z')

# Plot planet trajectory
for planet in planets:
ax.plot(xs=planet.trajectory[0],
ys=planet.trajectory[1],
zs=planet.trajectory[2],
zdir='z', label='ys=0, zdir=z')

# Plot Earth (plot is moon position relative to earth)
# also plotting to scale
(xs, ys, zs) = drawSphere(0, 0, 0, 6367.4447)
ax.plot_wireframe(xs, ys, zs, color="r")

plt.show()

@profile
def run_simulation(start_time, end_time, step, ship, planets):
"""Runs orbital simulation given ship and planet objects as well as start/stop times"""

# Calculate moon & planet update rate (1/10th as often as the craft)
planet_step_rate = int(math.ceil(((end_time - start_time) / step)/100))

# Initialize positions of planets
for planet in planets:
planet.update_position(start_time)
planet.log_position()

start = t.time()
total_time = end_time - start_time
# Total tics
total_steps = int((end_time - start_time) / step)

print("Days in simulation: {:6.2f}, number of steps: {}".format(
total_time, total_steps))

for i, time in enumerate(np.arange(start_time, end_time, step)):

# Every planet_step_rate update the position estmation
if (i % planet_step_rate == 0):
for planet in planets[1:]:
planet.update_position(time)
planet.log_position()

# Update craft velocity and force vectors related to planets
for planet in planets:
ship.force_g(planet)

# Update ship position according to sum of affecting forces
ship.update()

# Log the position of the ship every 1000 steps
if (i % 1000) == 0:
# Append the position to the lists
ship.log()

# Print status update every 100,000 steps
if (i % 100000) == 0:
if t.time()-start > 1:
time_since_start = time - start_time

# Calculate estmated time remaining
remaining_days_of_simulation = end_time - time
remaining_percentage_of_simulation = (1 - time_since_start / total_time) * 100

# Tics per second till now
steps_pr_second = int(math.ceil( i / (t.time()-start)))
# Estmated remaining seconds
remaining_simulation_time_in_secs = (total_steps - i) / steps_pr_second

# Split remaining_simulation_time_in_secs into hours, minutes and second
minutes, seconds = divmod(remaining_simulation_time_in_secs, 60)
hours, minutes = divmod(minutes, 60)
print("  Simulation remaining: {:6.2f} days, {:5.2f}%.  Script statistics - steps/s: {:6}, est. time left: {:0>2}:{:0>2}:{:0>2}".format(
remaining_days_of_simulation, remaining_percentage_of_simulation,
steps_pr_second, hours, minutes, seconds))

end = t.time()
elapsed = end - start
print("Final steps per second: {0:.2f}".format((total_time / step)/elapsed))

minutes, seconds = divmod(elapsed, 60)
hours, minutes = divmod(minutes, 60)
print("Total running time: {:0>2.0f}:{:0>2.0f}:{:0>2.0f}".format(
hours, minutes, seconds))

def main():

# Delta time for simulations (in days)
delta_t = np.longdouble(1.0) / (24.0 * 60.0 * 60.0)

ship = Craft(delta_t, x=35786, y=1, z=1, v_x=0, v_y=4.5, v_z=0, mass=12)

# Initialize simulation time
simulation_start = Time('2015-09-10T00:00:00')
simulation_end = Time('2015-10-10T00:00:00')

if not os.path.isfile('de430.bsp'):
kernel = SPK.open('de430.bsp')

planets = [Planet(kernel, 399, 399, np.longdouble(5972198600000000000000000)),
Planet(kernel, 301, 399, np.longdouble(7.34767309 * 10**22))]

run_simulation(simulation_start.jd, simulation_end.jd, delta_t, ship, planets)
plot(ship, [planets[1]])

if __name__ == "__main__":
main()


This updates the planets and the ship at every step, and prints all positions in history. One could add a for loop around the ship update to make it update slightly more often, both there is almost no gain doing that as it is the ships calculation which are time consuming.

Timed output:

\$ time python Orbital.py
Days in simulation:  30.00, number of steps: 5000
Simulation remaining:  30.00 days, 100.00%.  Script statistics - remaining rounds: 5000
Simulation remaining:  27.00 days,  90.00%.  Script statistics - steps/s:   2225, est. time left: 00:00:02
Simulation remaining:  24.00 days,  80.00%.  Script statistics - steps/s:   2237, est. time left: 00:00:01
Simulation remaining:  21.00 days,  70.00%.  Script statistics - steps/s:   2231, est. time left: 00:00:01
Simulation remaining:  18.00 days,  60.00%.  Script statistics - steps/s:   2236, est. time left: 00:00:01
Simulation remaining:  15.00 days,  50.00%.  Script statistics - steps/s:   2230, est. time left: 00:00:01
Simulation remaining:  12.00 days,  40.00%.  Script statistics - steps/s:   2236, est. time left: 00:00:00
Simulation remaining:   9.00 days,  30.00%.  Script statistics - steps/s:   2232, est. time left: 00:00:00
Simulation remaining:   6.00 days,  20.00%.  Script statistics - steps/s:   2230, est. time left: 00:00:00
Simulation remaining:   3.00 days,  10.00%.  Script statistics - steps/s:   2226, est. time left: 00:00:00
Final steps per second: 2225.10
Total running time: 00:00:02

real    0m6.939s
user    0m3.143s
sys 0m0.352s


That is 7 seconds including compiling the Python script, and accessing that large de430.bsp file. And the image rotated to almost the same as OP's rotation gives this image:

In a comment you state that you want accuracy, and as such want to calculate more. There is a balance between doing 2.5 millions steps, and 3000 steps, which is related to the speed of the various objects and distance. The further away the less impact they have on each other.

Another point regarding accuracy is that you claim to use the earth as your reference point, whilst your code doesn't actually do that, and here is why:

• Your ship starts at a given point, (35786, 1, 1), which presumably is the Earth's position at simulation start. But you don't correct the position of the ship against the moved Earth position
• When calculating the position of the Moon position you used the [3, 399] vs [3, 301], but you don't counteract that this is the positioning according to the Earth´s barycenter vs Solar System barycenter (aka [0, 3]). See https://pypi.python.org/pypi/jplephem/, the paragraph on "learning the position of Mars with respect to the Earth".
• The same applies to some extent to the Earth position, if I'm not mistaken

So you if aiming for accuracy you might want to reevaluate the mathematics, and what your base reference is supposed to be, aka if you should use the Earths barycenter or the actual position of the Earth, or possibly the Solar System barycenter. In the code I've provided so far, I assumed that your mathematics where correct, so I've not counteracted for these inaccuracies.

• I initialize the position of the ship at (35786, 1, 1) arbitrarily, just for testing, not because it was the position of earth, the positions provided of the earth and the moon are relative the the earth-moon barycentre, so in order to get the relative position of the earth to the moon I only need to subtract the positions, or I would have the earth wobbling around" and the moon in a less circular (appearing) orbit but because I am only focusing on the earth moon system I only care about the position of the moon in the earth's inertial frame. – Mark Omo Dec 3 '15 at 20:14
• Still you disregard the movement of the ship according to general movement of earth and/or moon with respect to the Solar System barycenter which the de430.bsp uses as base... – holroy Dec 3 '15 at 20:16
• that's because at this point, I only care about the movement of the ship relative to the earth/moon not the solar system. Only the earth/moon barycenter is referenced to the solar system barycenter, along with all the other planets, and the sun. See here:pypi.python.org/pypi/jplephem#getting-started-with-de430 it shows the relations including what is referenced to what. – Mark Omo Dec 3 '15 at 20:20