I'm trying to make a computer simulation for Earth, Moon and Asteroid as they orbit around the Sun.
The data r_earth, r_moon and the r_asteroid are already given (they are an array). The xlim and ylim are so big because I have to use real values from the universe, otherwise I wouldn't see anything on the animation.
import numpy as np
def RK4_step(y, t, f, dt): #for numerical part I decided to use RK4 method
k1 = dt * f(y, t)
k2 = dt * f(y + 0.5 * k1, t + 0.5 * dt)
k3 = dt * f(y + 0.5 * k2, t + 0.5 * dt)
k4 = dt * f(y + k3, t + dt)
return (t + dt, y + (k1 + 2. * (k2 + k3) + k4) / 6.)
def integrate_ode(y0, derivative_fn, step_fn, a, b, n):
dt = (b - a) / (n - 1)
t = a
y = np.array(y0)
res = [y]
for _ in range(1, n):
t, y = step_fn(y, t, derivative_fn, dt)
res.append(y)
return (np.linspace(a, b, n), np.array(res))
# np.array(res) is 2D array of dimensions (<time steps>, <um variables>) that holds the values for each time step.
# Initial parameters for body masses(Sun, Earth, Moon, Asteroid)
G = 6.674e-11
m_s = 1.989e30
m_e = 5.972e24
m_m = 7.346e22
m_a = 1e25
def derivative_fn(y, t): # vectors necessary for 2.Newton's Law
d1 = y[1] - y[0]
vec01 = d1 / np.linalg.norm(d1) ** 3
d2 = y[2] - y[0]
vec02 = d2 / np.linalg.norm(d2) ** 3
d3 = y[2] - y[1]
vec12 = d3 / np.linalg.norm(d3) ** 3
der = np.array([
y[3],
y[4],
y[5],
- G * m_s * y[0] / np.linalg.norm(y[0]) ** 3 + G * m_m * vec01 + G * m_a * vec02,
- G * m_s * y[1] / np.linalg.norm(y[1]) ** 3 - G * m_e * vec01 + G * m_a * vec12,
- G * m_s * y[2] / np.linalg.norm(y[2]) ** 3 - G * m_e * vec02 - G * m_m * vec12])
return der
# numerically computed data in the array form
# y = [r_e, r_m, r_a, v_e, v_m, v_a] -> 2D array of size (6, 2)
# the initial distances between the bodies,
# for Earth and Moon we compute it from 2nd Newton's Law (orbital velocity)
r_e0 = 1.519e11
v_e0 = np.sqrt(G * m_s / r_e0)
r_m0 = 1.521e11
v_m0 = np.sqrt(G * m_e / abs(r_m0 - r_e0))
r_a0 = 1e11
v_a0 = 31293 # m/s
y0 = np.array([[r_e0, 0.], [r_m0, 0.], [r_a0, 0.], [0., v_e0], [0., v_e0 + v_m0], [0., v_e0 + v_a0]])
# all the parameters are now in the array y0,
# the last two are velocities relative to the Earth
n = 10000
ts, ys = integrate_ode(y0, derivative_fn, RK4_step, 0., 5.154e7, n)
r_earth = ys[:, 0, :]
r_moon = ys[:, 1, :]
r_asteroid = ys[:, 2, :]
# Plotting
import matplotlib.pyplot as plt
plt.xlabel("x")
plt.ylabel("y")
plt.plot(r_earth[:, 0], r_earth[:, 1], label="Earth")
plt.plot(r_moon[:, 0], r_moon[:, 1], label="Moon")
plt.plot(r_asteroid[:, 0], r_asteroid[:, 1], label="Asteroid")
plt.show()
from matplotlib.animation import FuncAnimation
fig, ax = plt.subplots()
ax.set_xlim([-0.3e12, 0.3e12])
ax.set_ylim([-0.3e12, 0.3e12])
# xlim and ylim are so big because we wouldn't see the animation.
sun, = ax.plot(0., 0., 'oy', ms=20)
earth, = ax.plot(r_earth[0, 0], r_earth[0, 1], 'og', ms=10)
moon, = ax.plot(r_moon[0, 0], r_moon[0, 0], 'ob', ms=3)
asteroid, = ax.plot(r_asteroid[0, 0], r_asteroid[0, 0], 'ok', ms=1)
num_frames = 10000
frame_duration = n / num_frames
def animation_frame(frame):
index = int(frame_duration * frame)
earth.set_data(r_earth[index, 0], r_earth[index, 1])
moon.set_data(r_moon[index, 0] / r_earth[index, 0], r_moon[index, 1] / r_earth[index, 1])
asteroid.set_data(r_asteroid[index, 0] / r_earth[index, 0], r_asteroid[index, 1] / r_earth[index, 1])
return earth, moon, asteroid
animation = FuncAnimation(fig, func=animation_frame, frames=range(num_frames), interval=50)
plt.show()
When I start my code, the Moon is in the same position as Sun, both are static. The Earth is moving around the Sun very slowly and there is no Asteroid on the animation.
How to fix this? I still want to use real data, but wish that animation would go 100x faster(but still smooth animation). How do I do that?