我要计算没有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都增加了。在

它没有显示椭圆轨道。我如何修正我的代码或写一个新的代码来描述地球和月球引力系统的完美椭圆轨道。在

谢谢你的帮助！在