|
|
Python代码模拟的太阳系，包括了水星(Mercury), 金星(Venus)，地球(Earth)，月球(Moon)，火星(Mars)
上面的动画是我用86行Python代码模拟的一个比较真实的太阳系，用到的基本原理就是N体问题，这是计算天体物理里面的一个重要问题和方法。计算天体物理是用计算机来模拟宇宙中的天体和天文现象，从而对它们进行科学研究。真正用于科研的模拟程序大多是非常复杂的，有着上万行甚至更多代码。这里我写了一个非常简短的代码模拟太阳系行星的运动轨迹，目的只是为了好玩以及展现一下计算天体物理的奇妙之处。
我们知道，行星沿着主要由太阳引力决定的轨道运动。要计算行星的轨道，我们需要计算来自太阳的引力并将其积分以获得速度和位置。我们在高中的时候学过牛顿万有引力定律和牛顿运动定律：
\vec a = \frac{GM}{|\vec{r}|^2} \frac{\vec r}{|\vec{r}|} \\
然后通过下面的公式来演化行星的位置和速度：
\vec{r}_{\rm new} = \vec{r}_{\rm old} + \vec{v}_{\rm old}~{\rm d}t\\ \vec{v}_{\rm new} = \vec{v}_{\rm old} + \vec{a}_{\rm new}~{\rm d}t
这里用的是欧拉积分法的一个改进方法--反向欧拉法。这个方法不同于欧拉法的是，它是稳定的，所以行星的轨道是闭合的。这个方法既简单又比较精确性。
上面的三个向量方程可以写成分量形式：
\begin{align} \begin{bmatrix} x \\ y \\ z \end{bmatrix} &= \begin{bmatrix} x \\ y \\ z \end{bmatrix} + \begin{bmatrix} v_x \\ v_y \\ v_z \end{bmatrix} \cdot {\rm d}t \\ \begin{bmatrix} a_x \\ a_y \\ a_z \end{bmatrix} &= \begin{bmatrix} x \\ y \\ z \end{bmatrix} \cdot \frac{GM}{(x^2+y^2+z^2)^{3/2}} \\ \begin{bmatrix} v_x \\ v_y \\ v_z \end{bmatrix} &= \begin{bmatrix} v_x \\ v_y \\ v_z \end{bmatrix} + \begin{bmatrix} a_x \\ a_y \\ a_z \end{bmatrix} \cdot {\rm d}t \end{align}\\
用Python编程的一个好处就是Python可以很优雅地处理向量，不同于C/C++。上面的九个方程式（三个向量方程式，每个具有三个分量）可以很容易地用三行代码实现：
p.r += p.v * dt
acc = -2.959e-4 * p.r / np.sum(p.r**2)**(3./2) # in units of AU/day^2
p.v += acc * dt现在我们可以演化行星的轨道了，我们还缺的是初始条件，也就是某个时间所有行星的位置和速度。我没有用随意假设的初始条件去做一个艺术性的动画，而是从NASA的JPL实验室的Horizons程序里调用太阳系的数据[1]，获取一个时刻的行星的精确位置和速度作为模拟的初始条件
from astroquery.jplhorizons import Horizons
obj = Horizons(id=1, location="@sun", epochs=Time("2017-01-01").jd, id_type='id').vectors()其中 id=1 代表水星（2代表金星，依此类推）。然后 obj['x'] 就是水星在2017年1月1日的x坐标。
程序计算了4个内行星（水星、金星、地球和火星）加上月球在2019年和2020年的轨道。对应的现实世界的时间显示在左上角。这个时间段可以在代码的前两行随意设置。
将上面这些内容放在一起，再加上一些 matplotlib 的 animation 动画包和其他一些设置，就完成了全部的代码。完整的代码在文章的最后。
这个模拟哪些是真实的，哪些是不真实的？
程序计算的所有星体的位置和速度是准确的，当然会存在一些模型和算法上的误差，比如没有考虑行星之间的引力。星体的大小不是真实的比例，不然就什么也看不见了，但这它们的相对大小是对的。也是这个原因，月球（小黄点）才会出现在地球的上部。另外，程序的计算是3维的，动画显示的只是在地球轨道平面的投影。如果你感兴趣的话，画面地球轨道的最左边对应的是春分，也就是3月21日，然后最下边是夏至。一般太阳系的坐标系是这样定义的。
程序是实时计算行星轨迹，还是从互联网调用数据然后做成动画？
是实时计算的。初始条件是用的网上的数据。
为什么水星的轨迹这么丑？
呃，它就是这么丑，我有什么办法。水星跟太阳系中其他所有行星不同，它的轨道的偏心率很高，高达0.2，意思就是它的轨道不是很圆。
可以做成3D动画吗？
2D的就足够了吧，因为行星基本都在同一个平面上。我估计写成3D的也不是非常复杂，大概可以把代码控制在100行以内。如果很多人有兴趣的话，我以后可能会更新一次。
为什么只有4个行星？
因为从第5颗行星木星开始，轨道半径变得非常大（大于地球的5倍）。加上它们的话，就几乎看不见地球了。
最后附上代码。下载请去GitHub：https://github.com/chongchonghe/Python-solar-system
#!/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from astropy.time import Time
from astroquery.jplhorizons import Horizons
sim_start_date = "2019-01-01" # simulating a solar system starting from this date
sim_duration = 2 * 365 # (int) simulation duration in days
m_earth = 5.9722e24 / 1.98847e30 # Mass of Earth relative to mass of the sun
m_moon = 7.3477e22 / 1.98847e30
class Object:
def __init__(self, name, rad, color, r, v):
self.name = name
self.r = np.array(r, dtype=np.float)
self.v = np.array(v, dtype=np.float)
self.xs = []
self.ys = []
self.plot = ax.scatter(r[0], r[1], color=color, s=rad**2, edgecolors=None, zorder=10)
self.line, = ax.plot([], [], color=color, linewidth=1.4)
class SolarSystem:
def __init__(self, thesun):
self.thesun = thesun
self.planets = []
self.time = None
self.timestamp = ax.text(.03, .94, 'Date: ', color='w', transform=ax.transAxes, fontsize='x-large')
def add_planet(self, planet):
self.planets.append(planet)
def evolve(self):
dt = 1.0
self.time += dt
plots = []
lines = []
for i2, p in enumerate(self.planets):
p.r += p.v * dt
acc = -2.959e-4 * p.r / np.sum(p.r**2)**(3./2) # in units of AU/day^2
if p.name == 399: # Force from the Moon to Earth
dr = p.r - self.planets[3].r # trick here, assuming moon is the 4th object
acc += -2.959e-4 * m_moon * dr / np.sum(dr**2)**(3./2)
if p.name == 301: # Force from earth to the Moon
dr = p.r - self.planets[2].r
acc += -2.959e-4 * m_earth * dr / np.sum(dr**2)**(3./2)
p.v += acc * dt
p.xs.append(p.r[0])
p.ys.append(p.r[1])
p.plot.set_offsets(p.r[:2])
plots.append(p.plot)
if i2 != 3: # ignore trajectory lines of the Moon
p.line.set_xdata(p.xs)
p.line.set_ydata(p.ys)
lines.append(p.line)
if len(p.xs) > 10000:
raise SystemExit("Stopping after a long run to prevent memory overflow")
self.timestamp.set_text('Date: ' + Time(self.time, format='jd', out_subfmt='date').iso)
return plots + lines + [self.timestamp]
plt.style.use('dark_background')
fig = plt.figure(figsize=[6, 6])
ax = plt.axes([0., 0., 1., 1.], xlim=(-1.8, 1.8), ylim=(-1.8, 1.8))
ax.set_aspect('equal')
ax.axis('off')
ss = SolarSystem(Object("Sun", 28, 'red', [0, 0, 0], [0, 0, 0]))
ss.time = Time(sim_start_date).jd
colors = ['gray', 'orange', 'blue', 'yellow', 'chocolate']
sizes = [0.38, 0.95, 1., 0.27, 0.53]
names = ['Mercury', 'Venus', 'Earth-Moon', 'Earth-Moon', 'Mars']
for i, nasaid in enumerate([1, 2, 399, 301, 4]):
obj = Horizons(id=nasaid, location="@sun", epochs=ss.time, id_type='id').vectors()
ss.add_planet(Object(nasaid, 20 * sizes, colors,
[np.double(obj[xi]) for xi in ['x', 'y', 'z']],
[np.double(obj[vxi]) for vxi in ['vx', 'vy', 'vz']]))
ax.text(0, - ([.47, .73, 1, 1.02, 1.5] + 0.1),
names[2] if i in [2, 3] else names, color=colors, zorder=1000,
ha='center', fontsize='large')
def animate(i):
return ss.evolve()
ani = animation.FuncAnimation(fig, animate, repeat=False, # init_func=init,
frames=sim_duration, blit=True, interval=20,)
plt.show()
# ani.save('solar_system_6in_200dpi.mp4', fps=40, dpi=200)长一点的视频：
Python代码模拟的太阳系
https://www.zhihu.com/video/1199877038029737984
视频：Python代码模拟的太阳系，包括了水星(Mercury), 金星(Venus)，地球(Earth)，月球(Moon)，火星(Mars) 参考
- ^https://ssd.jpl.nasa.gov/?horizons
|
|
窥视卡
雷达卡
发表于 2023-1-1 20:33:46
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
显身卡
