我一直在尝试模拟绕太阳运动的行星和小行星，我发现了这个链接： 如何在已绘制的椭圆上绘制行星轨道随时间的变化？
我决定研究并尝试其中的代码。但似乎我要么使用错误，要么代码错误，因为当我绘制火星、地球、木星、水星和金星的轨道时，它们似乎与美国宇航局的在线模拟不一致，例如这个： https://eyes.nasa.gov/apps/solar-system/#/home
根据 NASA 模拟的当前位置： 美国宇航局模拟：
根据我的情节，当前位置： 我的模拟
代码：
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse, Circle
#Set axes aspect to equal as orbits are almost circular; hence square is needed
ax = plt.figure(0).add_subplot(111, aspect='equal')
#Setting the title, axis labels, axis values and introducing a grid underlay
#Variable used so title can indicate user inputed date
plt.title('Inner Planetary Orbits at[user input date]')
plt.ylabel('x10^6 km')
plt.xlabel('x10^6 km')
ax.set_xlim(-300, 300)
ax.set_ylim(-300, 300)
plt.grid()
#Creating the point to represent the sun at the origin (not to scale),
ax.scatter(0,0,s=200,color='y')
plt.annotate('Sun', xy=(0,-30))
#Implementing ellipse equations to generate the values needed to plot an ellipse
#Using only the planet's min (m) and max (M) distances from the sun
#Equations return '2a' (the ellipses width) and '2b' (the ellipses height)
def OrbitLength(M, m):
a=(M+m)/2
c=a-m
e=c/a
b=a*(1-e**2)**0.5
print(a)
print(b)
return 2*a, 2*b
#This function uses the returned 2a and 2b for the ellipse function's variables
#Also generating the orbit offset (putting the sun at a focal point) using M and m
def PlanetOrbit(Name, M, m):
w, h = OrbitLength(M, m)
Xoffset= ((M+m)/2)-m
Name = Ellipse(xy=((Xoffset),0), width=w, height=h, angle=0, linewidth=1, fill=False)
ax.add_artist(Name)
from math import *
EPSILON = 1e-12
def solve_bisection(fn, xmin,xmax,epsilon = EPSILON):
while True:
xmid = (xmin + xmax) * 0.5
if (xmax-xmin < epsilon):
return xmid
fn_mid = fn(xmid)
fn_min = fn(xmin)
if fn_min*fn_mid < 0:
xmax = xmid
else:
xmin = xmid
'''
Found something similar at this gamedev question:
https://gamedev.stackexchange.com/questions/11116/kepler-orbit-get-position-on-the-orbit-over-time?newreg=e895c2a71651407d8e18915c38024d50
Equations taken from:
https://en.wikipedia.org/wiki/Kepler%27s_laws_of_planetary_motion#Position_as_a_function_of_time
'''
def SolveOrbit(rmax, rmin, t):
# calculation precision
epsilon = EPSILON
# mass of the sun [kg]
Msun = 1.9891e30
# Newton's gravitational constant [N*m**2/kg**2]
G = 6.6740831e-11
# standard gravitational parameter
mu = G*Msun
# eccentricity
eps = (rmax - rmin) / (rmax + rmin)
# semi-latus rectum
p = rmin * (1 + eps)
# semi/half major axis
a = p / (1 - eps**2)
# period
P = sqrt(a**3 / mu)
# mean anomaly
M = (t / P) % (2*pi)
# eccentric anomaly
def fn_E(E):
return M - (E-eps*sin(E))
E = solve_bisection(fn_E, 0, 2*pi)
# true anomaly
# TODO: what if E == pi?
theta = 2*atan(sqrt((((1+eps)*tan(E/2)**2)/(1-eps))))
# if we are at the second half of the orbit
if (E > pi):
theta = 2*pi - theta
# heliocentric distance
r = a * (1 - eps * cos(E))
return theta, r
def DrawPlanet(name, rmax, rmin, t):
SCALE = 1e9
theta, r = SolveOrbit(rmax * SCALE, rmin * SCALE, t)
x = -r * cos(theta) / SCALE
y = r * sin(theta) / SCALE
planet = Circle((x, y), 8)
ax.add_artist(planet)
#These are the arguments taken from hyperphysics.phy-astr.gsu.edu/hbase/solar/soldata2.html
#They are the planet names, max and min distances, and their longitudinal angle
#Also included is Halley's Comet, used to show different scale and eccentricity
PlanetOrbit('Mercury', 69.8, 46.0)
PlanetOrbit('Jupiter', 741, 817)
PlanetOrbit('Venus', 108.9, 107.5)
PlanetOrbit('Earth', 152.1, 147.1)
PlanetOrbit('Mars', 249.1, 206.7)
# t = time since perihelion in seconds
DrawPlanet('Earth', 152.1, 147.1, 206*60*60*24)
DrawPlanet('Mars', 249.1, 206.7, 77*60*60*24)
DrawPlanet('Mercury', 69.8, 46.0, 43*60*60*24)
DrawPlanet('Jupiter', 741, 817, 552*60*60*24)
DrawPlanet('Venus', 108.9, 107.5, 16*60*60*24)
print(60*60*24*365)
plt.show()
这两个图像似乎没有明显对齐，所以我尝试了很多方法，例如：仅使用地球的近日点时间，重写 SolveOrbit 函数，以及使用 JPL HORIZONS api 查找近日点时间，以防万一它是错误的，我是这样的找不到。
我开始认为问题出在我输入的时间“t”而不是程序，因为我从这个网站得到了它们： https://in-the-sky.org/news.php?id=20210612_11_100
也许它不准确，但我不确定还能在哪里找到近日点时间。
用眼睛很难判断轨道，但是对于轨道的2D视图，它们看起来很好，请参阅此处的链接（回溯机）。提供的链接是 BPhO 提出的一系列挑战。