一篇文章教你用Python繪畫一個(gè)太陽系
圖片上傳之后不知為何幀率降低了許多。。。
日地月三體
所謂三體,就是三個(gè)物體在重力作用下的運(yùn)動(dòng)。由于三點(diǎn)共面,所以三個(gè)質(zhì)點(diǎn)僅在重力作用下的運(yùn)動(dòng)軌跡也必然無法逃離平面。
三體運(yùn)動(dòng)所遵循的規(guī)律就是古老而經(jīng)典的萬有引力
則對于 m i 而言,
且
將其寫為差分形式
由于我們希望觀察三體運(yùn)動(dòng)的復(fù)雜形式,而不關(guān)系其隨對應(yīng)的宇宙星體,所以不必考慮單位制,將其在二維平面坐標(biāo)系中拆分,則
#后續(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])
由于地月之間的距離相對于日地距離太近,所以在畫圖的時(shí)候?qū)⑵鋽U(kuò)大100倍,得到圖像
盡管存在誤差,但最起碼看到了地球圍繞太陽轉(zhuǎn),月球圍繞地球轉(zhuǎn)。。。代碼為
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]) #繪圖時(shí)的地月距離擴(kuò)大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]) ### 由于火星離地球很遠(yuǎn),所以不必再改變尺度 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
得到
這個(gè)運(yùn)動(dòng)要比月球的運(yùn)動(dòng)簡單得多——前提是開上帝視角,俯瞰太陽系。如果站在地球上觀測火星的運(yùn)動(dòng),那么這個(gè)運(yùn)動(dòng)可能相當(dāng)帶感
所以這都能找到規(guī)律,托勒密那幫人也真夠有才的。
太陽系
由于太陽和其他星體之間的質(zhì)量相差懸殊,所以太陽系內(nèi)的多體運(yùn)動(dòng),都將退化為二體問題,甚至如果把太陽當(dāng)作不動(dòng)點(diǎn),那就成了單體問題了。
盡管如此,我們還是盡可能地模仿一下太陽系的運(yùn)動(dòng)情況
質(zhì)量 | 半長軸(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°。由于水星的軌道太小,偏不偏心其實(shí)都不太看得出來,所以就當(dāng)它是正圓也無所謂了,最后得圖
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)層的剛好相反,所以這個(gè)圖很難兼顧,觀感上也不太好看。
如果只畫出木星之前的星體,順便加上小行星帶,可能會(huì)好一些。
通過這個(gè)圖就能看出來,有一顆小行星被木星彈了過來,直沖沖地向地球趕來,幸好又被太陽彈了出去,可見小行星還是挺危險(xiǎn)的,好在這只是個(gè)假想圖。
相關(guān)文章
Python監(jiān)測屏幕界面內(nèi)容變化并發(fā)送通知方法詳解
這篇文章主要為大家介紹了Python監(jiān)測屏幕界面內(nèi)容變化并發(fā)送通知,有需要的朋友可以借鑒參考下,希望能夠有所幫助,祝大家多多進(jìn)步,早日升職加薪2023-06-06python數(shù)據(jù)處理之如何選取csv文件中某幾行的數(shù)據(jù)
這篇文章主要給大家介紹了關(guān)于python數(shù)據(jù)處理之如何選取csv文件中某幾行的數(shù)據(jù)的相關(guān)資料,文中通過示例代碼介紹的非常詳細(xì),對大家學(xué)習(xí)或者使用python具有一定的參考學(xué)習(xí)價(jià)值,需要的朋友們下面來一起學(xué)習(xí)學(xué)習(xí)吧2019-09-09python爬蟲模擬瀏覽器訪問-User-Agent過程解析
這篇文章主要介紹了python爬蟲模擬瀏覽器訪問-User-Agent過程解析,文中通過示例代碼介紹的非常詳細(xì),對大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值,需要的朋友可以參考下2019-12-12在python中l(wèi)ist作函數(shù)形參,防止被實(shí)參修改的實(shí)現(xiàn)方法
這篇文章主要介紹了在python中l(wèi)ist作函數(shù)形參,防止被實(shí)參修改的實(shí)現(xiàn)方法,具有很好的參考價(jià)值,希望對大家有所幫助。一起跟隨小編過來看看吧2020-06-06python3 pathlib庫Path類方法總結(jié)
這篇文章主要介紹了python3 pathlib庫Path類方法總結(jié),文中通過示例代碼介紹的非常詳細(xì),對大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值,需要的朋友可以參考下2019-12-12Python使用Pygame實(shí)現(xiàn)時(shí)鐘效果
這篇文章主要為大家詳細(xì)介紹了Python使用Pygame實(shí)現(xiàn)時(shí)鐘效果,文中示例代碼介紹的非常詳細(xì),具有一定的參考價(jià)值,感興趣的小伙伴們可以參考一下2021-08-08