sim.py

#!/usr/bin/env python2

import numpy as np

import matplotlib.pyplot as plt
import matplotlib.patches as ptch
import matplotlib.lines as lne
import matplotlib.animation as ani


#Constants#
G = 6.6734810*10**(-11)   # Universal grav constant G=6.6734810x10^-11 m^3/kgs^2
g = 9.81                # Gravity, g 

#Mass
massS = 1.989*10**30 #kg
#massS = 1000
massE = 5.972*10**24 #kg
#massE=1
#Velocity

num_steps = 5000

# These data arrays are [t, x, y, x', y', x'', y'']
#                   or [t, x, y, vx, vy, ax, ay]
earth_data = np.zeros((num_steps, 7))
earth_data[0] = [0, 149.59787*10**9, 0, 0, 29800, 0, 0]
dt = 36000
earth_data[:, 0] = dt * np.arange(num_steps)

for i in range(1, num_steps):
    t, x, y, vx, vy, ax, ay = earth_data[i - 1]
    r = np.sqrt( x**2 + y**2 )   # distance from Sun (0,0) to Earth, this should be constant in this sim
    ax = -(G*massS / r**3)*x     #acceleration in x direction
    ay = -(G*massS / r**3)*y     #acceleration in y direction
    vy = vy + ay*dt                         #updated Velocity with acc. in the y direction
    vx = vx + ax*dt                         #updated Velocity with acc. in the x direction
    x = x + vx*dt                           #updated position with vel. in the x direction
    y = y + vy*dt
    earth_data[i] = [earth_data[i, 0], x, y, vx, vy, ax, ay]

def update(i, earth):
    earth.set_data(earth_data[i, 1], earth_data[i, 2])
    return earth,

fig = plt.figure()
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])

# Sun will not move, don't replot
sun = ax.plot(0, 0, marker='o', ms=100, color='y')
earth, = ax.plot(earth_data[0, 1], earth_data[0, 2], marker='o', ms=10, color='g')

# Set figure dimensions
ax.set_xlim(-earth_data[0, 1], earth_data[0, 1])
ax.set_ylim(-earth_data[0, 1], earth_data[0, 1])

anim = ani.FuncAnimation(fig,         # Figure to operate on
                         update,      # Update function to call
                         num_steps,   # Number of "Frames" to generate
                         fargs=(earth,), # Additional arguments to pass to update function
                         interval=10, # Interval (in ms) between frames
                         blit=True)   # interpolate between frames to keep to 30fps
plt.show()