python - 如何使用python制作行星轨道

Question

我已经为 'r'(Radius) 编写了一般轨道方程的代码。值，并尝试为 x=r*cos(theta) 绘制图形和 y=r*sin(theta)但是我的代码为不同的“r”值提供了相同的“x”值。所以我要结束enter code here用直线代替椭圆。

from numpy import *
import matplotlib.pyplot as plt
import pprint

L = 9.11*10**38     #L = angular momentum
m = 3.28*10**23     #m = mass of mercury
M = 1.99*10**30     #M = mass of sun
a = 5.8*10**7       #a = semi-major axis
G = 6.674*10**-11   #G = Gravitationl constant
k = G*M*m        
E = -k/(2*a)        #E = energy
p = L**2/(m*k)   
c = 1 + (2*E*L**2)/m*k**2
e = sqrt(-c)        #e = eccentricity

def fx(x):
    r = p/(1 + e*cos(x))
    return r

n = 1000
phi =linspace(0,2*pi,n)
radius = zeros([n])
theta = zeros([n])
x = zeros([n])
y = zeros([n])

for i in range(0,n):
    radius[i] = fx(phi[i])
    theta[i] = 180*phi[i]/pi

for i in range(0,n):
    x[i] = radius[i]*cos(phi[i])

for i in range(0,n):
    y[i] = radius[i]*sin(phi[i])

print('r =',radius)
print('x =',x)
print('y =',y)
plt.plot(x,y)
plt.show()

这是我的数学书的截图:

Answer 1

您的代码(至少)有两个问题:

1) a = 5.8*10**7 应该是 a = 5.8*10**10 因为(给定你的 G)您需要以米而不是公里为单位的距离。

2) 正如@JohnO 所指出的，您需要 1 + (2*E*L**2)/(m*k**2) 而不是 1 + ( 2*E*L**2)/m*k**2

如果你做了这两个改变，你会得到一个合理的椭圆: