我一直在尝试模拟绕太阳运动的行星和小行星，我发现了这个链接： 如何在已经绘制的椭圆上绘制行星轨道作为时间的函数？
并且我决定研究并尝试其中的代码。但似乎我要么使用错误，要么代码错误，因为当我绘制火星、地球、木星、水星和金星的轨道时，它们似乎与美国宇航局的在线模拟不一致，例如这个： https://eyes.nasa.gov/apps/solar-system/#/home
根据美国航空航天局模拟的当前位置： 美国航空航天局模拟：
根据我的绘图的当前位置： 我的模拟
代码：
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
也许它不准确，但我不确定还能在哪里找到近日点时间。
在代码中发现了一个难题！虽然它看起来很复杂，但问题在于模拟中行星的起始位置。是对的，问题出在
t
（自近日点后的时间）上，但不是找到的值不准确。
假设所有行星的近日点同时发生。 这并非事实。 每个行星都有自己独立的近日点时间，该时间与其轨道周期不同。
以下是更正代码的方法：
-
找到每个行星的正确近日点时间： 可以从像 NASA HORIZONS 这样的可靠来源找到这些信息。
-
计算自近日点后的时间：
- 获取当前日期和时间。
-
为每个行星找到自其最近一次近日点以来的持续时间。这将是代码中的
t值。 -
将计算出的
t应用于每个行星： 在调用DrawPlanet时，将每个行星特有的t传入。
其他技巧：
- 验证的单位： 确保的所有单位（时间、距离等）在所有计算中都是一致的（例如，秒、米或天文单位）。
- 考虑时间跨度： 如果模拟较长的时间跨度，则可能需要考虑岁差和其他长期轨道变化等因素以获得更高的准确性。
- 逐步调试： 从一个行星开始，验证其在多个时间点的计算位置，然后添加更多行星。
通过更正自近日点以来的时间计算，模拟的准确性应该会提高，并且行星位置应该更接近于 NASA Eyes 等参考点。 祝好运，编码愉快！
标签：python,ellipse,orbital-mechanics From： 78800541