0

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?

thunder
  • 113
  • 4
  • I do not understand the problem. This runs smoothly on my machine (close to 20 FPS as specified). The earth move slowly but this is due to the computed values. If you want to be faster, then you can skip frames (by multiplying `index` by 2 or 3). You can also reduce the interval to 30 ms rather than 50 ms to get more FPS. – Jérôme Richard Aug 08 '21 at 12:09
  • Yes, it works!!! – thunder Aug 08 '21 at 14:26

0 Answers0