Orbital Trajectory simulator

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

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

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