Python代碼模擬的太陽系，包括了水星(Mercury), 金星(Venus)，地球(Earth)，月球(Moon)，火星(Mars)

上面的動畫是我用86行Python代碼模擬的一個比較真實的太陽系，用到的基本原理就是N體問題，這是計算天體物理里面的一個重要問題和方法。計算天體物理是用計算機來模擬宇宙中的天體和天文現象，從而對它們進行科學研究。真正用于科研的模擬程序大多是非常復雜的，有著上萬行甚至更多代碼。這里我寫了一個非常簡短的代碼模擬太陽系行星的運動軌跡，目的只是為了好玩以及展現一下計算天體物理的奇妙之處。

我們知道，行星沿著主要由太陽引力決定的軌道運動。要計算行星的軌道，我們需要計算來自太陽的引力并將其積分以獲得速度和位置。我們在高中的時候學過牛頓萬有引力定律和牛頓運動定律：

然后通過下面的公式來演化行星的位置和速度：

這里用的是歐拉積分法的一個改進方法--反向歐拉法。這個方法不同于歐拉法的是，它是穩(wěn)定的，所以行星的軌道是閉合的。這個方法既簡單又比較精確性。

上面的三個向量方程可以寫成分量形式：

用Python編程的一個好處就是Python可以很優(yōu)雅地處理向量，不同于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日，然后最下邊是夏至。一般太陽系的坐標系是這樣定義的。

程序是實時計算行星軌跡，還是從互聯(lián)網調用數據然后做成動畫？

是實時計算的。初始條件是用的網上的數據。

為什么水星的軌跡這么丑？

呃，它就是這么丑，我有什么辦法。水星跟太陽系中其他所有行星不同，它的軌道的偏心率很高，高達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 Horizonssim_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.98847e30class Object:def __init__(self, name, rad, color, r, v):self.name = nameself.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 = thesunself.planets = []self.time = Noneself.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.0self.time += dtplots = []lines = []for i2, p in enumerate(self.planets):p.r += p.v * dtacc = -2.959e-4 * p.r / np.sum(p.r**2)**(3./2) # in units of AU/day^2if p.name == 399: # Force from the Moon to Earthdr = p.r - self.planets[3].r # trick here, assuming moon is the 4th objectacc += -2.959e-4 * m_moon * dr / np.sum(dr**2)**(3./2)if p.name == 301: # Force from earth to the Moondr = p.r - self.planets[2].racc += -2.959e-4 * m_earth * dr / np.sum(dr**2)**(3./2)p.v += acc * dtp.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 Moonp.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[i], colors[i], [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][i] + 0.1),names[2] if i in [2, 3] else names[i], color=colors[i], 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

總結

以上是生活随笔為你收集整理的python一百行代码多少钱_用86行Python代码模拟太阳系的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。