一篇文章教你用Python繪畫一個太陽系
更新時間：2021年10月14日 12:05:09 作者：微小冷
這篇文章主要給大家介紹了關于如何利用Python繪畫一個太陽系，文中通過示例代碼介紹的非常詳細，對大家的學習或者工作具有一定的參考學習價值，需要的朋友們下面隨著小編來一起學習學習吧
圖片上傳之后不知為何幀率降低了許多。。。
日地月三體
所謂三體，就是三個物體在重力作用下的運動。由于三點共面，所以三個質點僅在重力作用下的運動軌跡也必然無法逃離平面。
三體運動所遵循的規(guī)律就是古老而經(jīng)典的萬有引力
則對于 m i 而言，
且
將其寫為差分形式
由于我們希望觀察三體運動的復雜形式，而不關系其隨對應的宇宙星體，所以不必考慮單位制，將其在二維平面坐標系中拆分，則
#后續(xù)代碼主要更改這里的參數(shù) m = [1.33e20,3.98e14,4.9e12] x = np.array([0,1.5e11,1.5e11+3.8e8]) y = np.array([0,0,0]) u = np.array([0,0,0]) v = np.array([0,2.88e4,1.02e3])
由于地月之間的距離相對于日地距離太近，所以在畫圖的時候將其擴大100倍，得到圖像
盡管存在誤差，但最起碼看到了地球圍繞太陽轉，月球圍繞地球轉。。。代碼為
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
m = [1.33e20,3.98e14,4.9e12]
x = np.array([0,1.5e11,1.5e11+3.8e8])
y = np.array([0.0,0,0])
u = np.array([0.0,0,0])
v = np.array([0,2.88e4,2.88e4+1.02e3])
fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(xlim=(-2e11,2e11),ylim=(-2e11,2e11))
ax.grid()
trace0, = ax.plot([],[],'-', lw=0.5)
trace1, = ax.plot([],[],'-', lw=0.5)
trace2, = ax.plot([],[],'-', lw=0.5)
pt0, = ax.plot([x[0]],[y[0]] ,marker='o')
pt1, = ax.plot([x[0]],[y[0]] ,marker='o')
pt2, = ax.plot([x[0]],[y[0]] ,marker='o')
k_text = ax.text(0.05,0.85,'',transform=ax.transAxes)
textTemplate = 't = %.3f days\n'
N = 1000
dt = 36000
ts = np.arange(0,N*dt,dt)/3600/24
xs,ys = [],[]
for _ in ts:
x_ij = (x-x.reshape(3,1))
y_ij = (y-y.reshape(3,1))
r_ij = np.sqrt(x_ij**2+y_ij**2)
for i in range(3):
for j in range(3):
if i!=j :
u[i] += (m[j]*x_ij[i,j]*dt/r_ij[i,j]**3)
v[i] += (m[j]*y_ij[i,j]*dt/r_ij[i,j]**3)
x += u*dt
y += v*dt
xs.append(x.tolist())
ys.append(y.tolist())
xs = np.array(xs)
ys = np.array(ys)
def animate(n):
trace0.set_data(xs[:n,0],ys[:n,0])
trace1.set_data(xs[:n,1],ys[:n,1])
#繪圖時的地月距離擴大100倍，否則看不清
tempX2S = xs[:n,1]+100*(xs[:n,2]-xs[:n,1])
tempY2S = ys[:n,1]+100*(ys[:n,2]-ys[:n,1])
trace2.set_data(tempX2S,tempY2S)
pt0.set_data([xs[n,0]],[ys[n,0]])
pt1.set_data([xs[n,1]],[ys[n,1]])
tempX = xs[n,1]+100*(xs[n,2]-xs[n,1])
tempY = ys[n,1]+100*(ys[n,2]-ys[n,1])
pt2.set_data([tempX],[tempY])
k_text.set_text(textTemplate % ts[n])
return trace0, trace1, trace2, pt0, pt1, pt2, k_text
ani = animation.FuncAnimation(fig, animate,
range(N), interval=10, blit=True)
plt.show()
ani.save("3.gif")
日地火
m = [1.33e20,3.98e14,4.28e13]
x = np.array([0,1.5e11,2.28e11])
y = np.array([0.0,0,0])
u = np.array([0.0,0,0])
v = np.array([0,2.88e4,2.4e4])
### 由于火星離地球很遠，所以不必再改變尺度
def animate(n):
trace0.set_data(xs[:n,0],ys[:n,0])
trace1.set_data(xs[:n,1],ys[:n,1])
trace2.set_data(xs[:n,2],ys[:n,2])
pt0.set_data([xs[n,0]],[ys[n,0]])
pt1.set_data([xs[n,1]],[ys[n,1]])
pt2.set_data([xs[n,2]],[ys[n,2]])
k_text.set_text(textTemplate % ts[n])
return trace0, trace1, trace2, pt0, pt1, pt2, k_text
得到
這個運動要比月球的運動簡單得多——前提是開上帝視角，俯瞰太陽系。如果站在地球上觀測火星的運動，那么這個運動可能相當帶感
所以這都能找到規(guī)律，托勒密那幫人也真夠有才的。
太陽系
由于太陽和其他星體之間的質量相差懸殊，所以太陽系內(nèi)的多體運動，都將退化為二體問題，甚至如果把太陽當作不動點，那就成了單體問題了。
盡管如此，我們還是盡可能地模仿一下太陽系的運動情況
| 質量 | 半長軸(AU) | 平均速度(km/s) | |
|---|---|---|---|
| 水星 | 0.055 | 0.387 | 47.89 |
| 金星 | 0.815 | 0.723 | 35.03 |
| 地球 | 1 | 1 | 29.79 |
| 火星 | 0.107 | 1.524 | 24.13 |
| 木星 | 317.8 | 5.203 | 13.06 |
| 土星 | 95.16 | 9.537 | 9.64 |
| 天王星 | 14.54 | 19.19 | 6.81 |
| 海王星 | 17.14 | 30.07 | 5.43 |
| 冥王星 |
除了水星偏心率為0.2，對黃道面傾斜為7°之外，其余行星的偏心率皆小于0.1，且對黃道面傾斜普遍小于4°。由于水星的軌道太小，偏不偏心其實都不太看得出來，所以就當它是正圓也無所謂了，最后得圖
au,G,RE,ME = 1.48e11,6.67e-11,1.48e11,5.965e24
m = np.array([3.32e5,0.055,0.815,1,
0.107,317.8,95.16,14.54,17.14])*ME*6.67e-11
r = np.array([0,0.387,0.723,1,1.524,5.203,
9.537,19.19,30.7])*RE
theta = np.random.rand(9)*np.pi*2
x = r*np.cos(theta)
y = r*np.sin(theta)
v = np.array([0,47.89,35.03,29.79,
24.13,13.06,9.64,6.81,5.43])*1000
u = -v*np.sin(theta)
v = v*np.cos(theta)
name = "solar.gif"
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(xlim=(-31*RE,31*RE),ylim=(-31*RE,31*RE))
ax.grid()
traces = [ax.plot([],[],'-', lw=0.5)[0] for _ in range(9)]
pts = [ax.plot([],[],marker='o')[0] for _ in range(9)]
k_text = ax.text(0.05,0.85,'',transform=ax.transAxes)
textTemplate = 't = %.3f days\n'
N = 500
dt = 3600*50
ts = np.arange(0,N*dt,dt)
xs,ys = [],[]
for _ in ts:
x_ij = (x-x.reshape(len(m),1))
y_ij = (y-y.reshape(len(m),1))
r_ij = np.sqrt(x_ij**2+y_ij**2)
for i in range(len(m)):
for j in range(len(m)):
if i!=j :
u[i] += (m[j]*x_ij[i,j]*dt/r_ij[i,j]**3)
v[i] += (m[j]*y_ij[i,j]*dt/r_ij[i,j]**3)
x += u*dt
y += v*dt
xs.append(x.tolist())
ys.append(y.tolist())
xs = np.array(xs)
ys = np.array(ys)
def animate(n):
for i in range(9):
traces[i].set_data(xs[:n,i],ys[:n,i])
pts[i].set_data(xs[n,i],ys[n,i])
k_text.set_text(textTemplate % (ts[n]/3600/24))
return traces+pts+[k_text]
ani = animation.FuncAnimation(fig, animate,
range(N), interval=10, blit=True)
plt.show()
ani.save(name)
由于外圈的行星軌道又長速度又慢，而內(nèi)層的剛好相反，所以這個圖很難兼顧，觀感上也不太好看。
如果只畫出木星之前的星體，順便加上小行星帶，可能會好一些。
通過這個圖就能看出來，有一顆小行星被木星彈了過來，直沖沖地向地球趕來，幸好又被太陽彈了出去，可見小行星還是挺危險的，好在這只是個假想圖。
相關文章
這篇文章主要介紹了django 簡單實現(xiàn)登錄驗證給你,文中通過示例代碼介紹的非常詳細，對大家的學習或者工作具有一定的參考學習價值,需要的朋友可以參考下2019-11-11
pygame.display.flip()和pygame.display.update()的區(qū)別及說明
這篇文章主要介紹了pygame.display.flip()和pygame.display.update()的區(qū)別及說明,具有很好的參考價值,希望對大家有所幫助,如有錯誤或未考慮完全的地方,望不吝賜教2024-03-03
在本篇文章里小編給大家整理了關于django3.02模板中的超鏈接配置實例代碼內(nèi)容，需要的朋友們可以學習參考下。2020-02-02
這篇文章主要介紹了CoordConv實現(xiàn)卷積加上坐標實例詳解，有需要的朋友可以借鑒參考下，希望能夠有所幫助，祝大家多多進步，早日升職加薪2023-03-03
python編程scrapy簡單代碼實現(xiàn)搜狗圖片下載器
這篇文章主要為大家介紹了使用python scrapy簡單代碼實現(xiàn)搜狗圖片下載器示例詳解，有需要的朋友可以借鑒參考下，希望能夠有所幫助2021-11-11
Python實現(xiàn)內(nèi)存泄露排查的示例詳解
一般在python代碼塊的調(diào)試過程中會使用memory-profiler、filprofiler、objgraph等三種方式進行輔助分析，今天這里主要介紹使用objgraph對象提供的函數(shù)接口來進行內(nèi)存泄露的分析，感興趣的可以了解一下2023-01-01
這篇文章主要介紹了Python3通過chmod修改目錄或文件權限的方法示例，文中通過示例代碼介紹的非常詳細，對大家的學習或者工作具有一定的參考學習價值，需要的朋友們下面隨著小編來一起學習學習吧2020-06-06
Jupyter?Notebook出現(xiàn)不是內(nèi)部或外部的命令解決方案
這篇文章主要介紹了Jupyter?Notebook出現(xiàn)不是內(nèi)部或外部的命令解決方案，具有很好的參考價值，希望對大家有所幫助。如有錯誤或未考慮完全的地方，望不吝賜教2023-06-06
主流被使用的地理坐標系并不統(tǒng)一，常用的有WGS84、GCJ02（火星坐標系）、BD09（百度坐標系）以及百度地圖中保存矢量信息的web墨卡托，本文利用Python編寫相關類以實現(xiàn)4種坐標系統(tǒng)之間的互相轉換，希望對大家有所幫助2023-02-02
這篇文章主要介紹了基于Python編寫一個計算器程序，實現(xiàn)簡單的加減乘除和取余二元運算，文中通過示例代碼介紹的非常詳細，對大家的學習或者工作具有一定的參考學習價值，需要的朋友們下面隨著小編來一起學習學習吧2021-09-09