moon.py

import numpy as np
import matplotlib.pyplot as plt
import math

num_steps = 150000

G = 6.67e-11 # N m2 / kg2
m_e = 5.97e24 # kg
m_m = 7.35e22 # kg
total_time = 27.32 * 24 * 3600 # s
radius = ((total_time**2 * G * m_e) / (4 * math.pi**2))**(1/3) # m
speed = math.sqrt((G * m_e) / radius) # m/s

def acceleration(moon_pos, ship_pos):
    d_em = np.linalg.norm(moon_pos)
    d_es = np.linalg.norm(ship_pos)
    d_sm = np.linalg.norm(ship_pos - moon_pos)
    acc_m = ((-moon_pos * G * m_e) / (d_em**3))
    acc_s = G * (((-ship_pos * m_e)/(d_es**3)) + ((m_m * (moon_pos - ship_pos))/(d_sm**3)))
    return acc_m, acc_s

def trajectory():
    h = total_time / num_steps
    x_m = np.zeros([num_steps + 1, 2])
    v_m = np.zeros([num_steps + 1, 2])
    x_s = np.zeros([num_steps + 1, 2])
    v_s = np.zeros([num_steps + 1, 2])

    x_m[0,0] = radius
    v_m[0,1] = speed
    x_s[0,1] = 82.e6
    v_s[0,0] = -2.e3

    for step in range(num_steps):
        acc_m, acc_s = acceleration(x_m[step], x_s[step])
        Ex_m = x_m[step] + (h * v_m[step])
        Ev_m = v_m[step] + (h * acc_m)
        Ex_s = x_s[step] + (h * v_s[step])
        Ev_s = v_s[step] + (h * acc_s)
        Hx_m = x_m[step] + (0.5 * h * (v_m[step] + Ev_m))
        Hx_s = x_s[step] + (0.5 * h * (v_s[step] + Ev_s))
        Hacc_m, Hacc_s = acceleration(Hx_m, Hx_s)
        Hv_m = v_m[step] + (0.5 * h * (acc_m + Hacc_m))
        Hv_s = v_s[step] + (0.5 * h * (acc_s + Hacc_s))       
        x_m[step + 1] = Hx_m
        x_s[step + 1] = Hx_s
        v_m[step + 1] = Hv_m
        v_s[step + 1] = Hv_s

    error = np.linalg.norm(x_m[-1] - x_m[0])

    return x_s, x_m, error


x_s, x_m, error = trajectory()
print("GTE = " + str(error))

def plot():
    plt.plot(x_s[:,0], x_s[:,1])
    plt.plot(x_m[:,0], x_m[:,1])
    plt.scatter(0,0)
    plt.show()


plot()