Orbital Trajectory simulator

Question

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 = fig.add_subplot(111, projection='3d')
    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'):
    raise ValueError('de430.bsp Was not found!')
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'):
    raise ValueError('de430.bsp Was not found!')
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.

Answer 1

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'):
            raise ValueError('de430.bsp Was not found!')
        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'):
            raise ValueError('de430.bsp Was not found!')
        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.

Answer 2

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 vertical space – Both between functions, and within a function, you could add blank lines to help your code read better.
• 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 = fig.add_subplot(111, projection='3d')
    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'):
        raise ValueError('de430.bsp Was not found!')
    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.

Answer 3

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 = fig.add_subplot(111, projection='3d')
    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'):
        raise ValueError('de430.bsp Was not found!')
    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:

Addendum on accuracy

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.