Created
October 23, 2016 01:26
-
-
Save TheDataLeek/210025af1e49d530124e7116f4680abc to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/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() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment