Как я могу решить орбиту Земли и Луны в питоне в численном моделировании без vpython?

Question

Я должен рассчитать 2 гравитационные системы тела (Земля и Луна) без vpython.

Целью этого кода является практика численного расчета и есть мой код.

import numpy as np
from math import *

from astropy.constants import *

import matplotlib.pyplot as plt

import time

start_time = time.time()

"""
G = Gravitational constant
g0 = Standard acceleration of gravity ( 9.8 m/s2)
M_sun = Solar mass
M_earth = Earth mass
R_sun = Solar darius
R_earth = Earth equatorial radius
au = Astronomical unit


"""

M_moon = 7.342E22
R_moon = 1.737E6
# Mean radius of moon.


M_earth = M_earth.value
R_earth = R_earth.value
G = G.value

perigee, apogee = 3.626E8, 4.054E8

position_E = np.array([0,0])
position_M = np.array([(perigee+apogee)/2.,0])
position_com = (M_earth*position_E+M_moon*position_M)/(M_earth+M_moon)

rel_pE = position_E - position_com
rel_pM = position_M - position_com


p_E = {"x":rel_pE[0], "y":rel_pE[1],"v_x":0, "v_y":10}
p_M = {"x":rel_pM[0], "y":rel_pM[1],"v_x":0, "v_y":-100}

t = range(0,365)

data_E , data_M = [0]*len(t), [0]*len(t)

def s(initial_velocity, acceleration, time):
    result = initial_velocity*time + 0.5*acceleration*time**2
    return result

def v(initial_velocity, acceleration, time):
    result = initial_velocity + acceleration*time
    return result

dist = float(sqrt((p_E["x"]-p_M['x'])**2 + (p_E["y"]-p_M["y"])**2))

# position data of Earth and Moon. make new list to make easy to draw plot

xE=[]
yE=[]
xM=[]
yM=[]

for i in t:
    dist = float(sqrt((p_E["x"]-p_M["x"])**2 + (p_E["y"]-p_M["y"])**2))


    a_Ex = -G*M_moon*p_E["x"]/(dist**2)
    a_Ey = -G*M_moon*p_E["y"]/(dist**2)

    data_E[i] = p_E

    p_E["x"] += s(p_E['v_x'], a_Ex, 24*3600)
    p_E["v_x"] += v(p_E['v_x'], a_Ex, 24*3600)
    p_E["y"] += s(p_E['v_y'], a_Ey, 24*3600)
    p_E["v_y"] += v(p_E['v_y'], a_Ey, 24*3600)

    xE += [p_E["x"]]
    yE += [p_E["y"]]

    a_Mx = -G*M_earth*p_M["x"]/(dist**2)
    a_My = -G*M_earth*p_M["y"]/(dist**2)

    data_M[i] = p_M

    p_M["x"] += s(p_M['v_x'], a_Mx, 24*3600)
    p_M["v_x"] += v(p_M['v_x'], a_Mx, 24*3600)
    p_M["y"] += s(p_M['v_y'], a_My, 24*3600)
    p_M["v_y"] += v(p_M['v_y'], a_My, 24*3600)

    xM += [p_M["x"]]
    yM += [p_M["y"]]

print("\n Run time \n --- %d seconds ---" %(time.time()-start_time))

Но этот код увеличивает как x, так и y.

Это не показывает эллиптическую орбиту. Как я могу исправить свой код или написать новый код, описывающий идеальную эллиптическую орбиту относительно гравитационной системы Земля и Луна.

Спасибо за вашу помощь!