我要计算没有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都增加了。在
它没有显示椭圆轨道。我如何修正我的代码或写一个新的代码来描述地球和月球引力系统的完美椭圆轨道。在
谢谢你的帮助！在