Python算法繪制特洛伊小行星群實(shí)現(xiàn)示例
書接上文
用Python搓一個(gè)太陽系
你們要的3D太陽系
3體人真的存在嗎
太長不看版
最小勢能點(diǎn)
在由兩個(gè)大質(zhì)量物體構(gòu)成的重力系統(tǒng)中，有一些特殊的區(qū)域會(huì)在兩個(gè)天體的頂級拉扯之下達(dá)到平衡，這些點(diǎn)就是拉格朗日點(diǎn)。而所謂平衡并非受力平衡，而是要求這個(gè)區(qū)域的物體會(huì)跟著雙星系統(tǒng)以相同的角速度運(yùn)動(dòng)。
根據(jù)上帝是個(gè)胖子這個(gè)假定，狀態(tài)穩(wěn)定意味著低勢能。所以在解析求解拉格朗日點(diǎn)之前，我們可以試著畫出這個(gè)雙星系統(tǒng)的勢能分布。
接下來搞一下太陽和木星：
可見木星在太陽的引力場下根本無法自己，但若把坐標(biāo)系調(diào)整一下，會(huì)看到木星雖小，卻還是有自己地盤的，畢竟也是有諸多衛(wèi)星的，這就意味著木星和太陽之間必然存在一些相對平衡的位置。
為了看得更加仔細(xì)，取對數(shù)是個(gè)不錯(cuò)的方法
import numpy as np import matplotlib.pyplot as plt from matplotlib import cm,ticker R = 7.1492e7 M1,M2 = 1.9891e30, 1.8982e27 G = 6.67e-11 mu = M2/M1 R1,R2 = np.array([mu,1])/(1+mu)*R x,y = np.meshgrid( np.arange(-0.5,1.5,1e-3)*R2 ,np.arange(-1,1,1e-3)*R2) H = -G*M1/np.sqrt((x+R1)**2+y**2) H -= G*M2/np.sqrt((x-R2)**2+y**2) H -= G*(M1+M2)*(x**2+y**2)/2/R**3 absH = np.abs(H) absH[absH>1e14] = 1e14 #去除奇點(diǎn) absH -= (np.min(absH)-1) print(np.max(absH),np.min(absH)) plt.contourf(x,y,np.log(absH),50,alpha=0.75, cmap=cm.PuBu_r) plt.show()
拉格朗日點(diǎn)
公式預(yù)警→_→
根據(jù)剛剛的圖可以看出，一般天體都會(huì)有一個(gè)屬于自己的私密區(qū)域，在這個(gè)區(qū)域里，別的天體的引力作用甚微，此區(qū)域即希爾球，拉格朗日點(diǎn)則是兩個(gè)天體希爾球的分界處。
在極坐標(biāo)下，可得
對于木星來說，五個(gè)拉格朗日點(diǎn)一般默認(rèn)為
特洛伊小行星群
參考此前的太陽系行星位置，得到其三維圖
from os import cpu_count import numpy as np from numpy.random import rand import matplotlib.pyplot as plt from matplotlib import animation 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])*ME*G r = np.array([0,0.387,0.723,1,1.524,5.203])*RE v = np.array([0,47.89,35.03,29.79,24.13,13.06])*1000 theta = rand(len(m))*np.pi*2 theta[-1] = 0#木星初始角度為0 cTheta,sTheta = np.cos(theta), np.sin(theta) xyz = r*np.array([cTheta, sTheta, 0*r]) #位置三分量 uvw = v*np.array([-sTheta, cTheta, 0*v]) #速度三分量 N_ast = 100 x_ast = xyz[0][-1]/2*( 1+(np.random.rand(100)-0.5)*0.1) y_ast = xyz[0][-1]/2*np.sqrt(3)*( 1+(np.random.rand(100)-0.5)*0.1) y_flag = np.random.rand(100)>0.5 y_ast = y_ast*(2*y_flag-1) m_ast = rand(N_ast)*1e20 r_ast = np.sqrt(x_ast**2+y_ast**2) v_ast = np.sqrt(G*3.32e5*ME/r_ast) #小行星速度sqrt(GM/R) theta = rand(N_ast)*np.pi*2 phi = (rand(N_ast)-0.5)*0.3 #給一個(gè)隨機(jī)的小傾角 cTheta,sTheta = x_ast/r_ast, y_ast/r_ast cPhi,sPhi = np.cos(phi),np.sin(phi) xyza = np.array([x_ast, y_ast, sPhi]) uvwa = v_ast*np.array([-sTheta*cPhi, cTheta*cPhi, sPhi]) name = "solar1.gif" fig = plt.figure(figsize=(10,10)) ax = fig.add_subplot(projection='3d') ax.grid() ax.set_xlim3d([-5.5*RE,5.5*RE]) ax.set_ylim3d([-5.5*RE,5.5*RE]) ax.set_zlim3d([-5.5*RE,5.5*RE]) traces = [ax.plot([],[],[],'-', lw=0.5)[0] for _ in range(len(m))] pts = [ax.plot([],[],[],marker='o')[0] for _ in range(len(m))] pt_asts = [ax.plot([],[],[],marker='.',lw=0.2)[0] for _ in range(N_ast)] N = 10000 dt = 3600*50 ts = np.arange(0,N*dt,dt) xyzs,xyzas = [],[] for _ in ts: xyz_ij = (xyz.reshape(3,1,len(m))-xyz.reshape(3,len(m),1)) r_ij = np.sqrt(np.sum(xyz_ij**2,0)) xyza_ij = (xyz.reshape(3,1,len(m))-xyza.reshape(3,N_ast,1)) ra_ij = np.sqrt(np.sum(xyza_ij**2,0)) for j in range(len(m)): for i in range(len(m)): if i!=j : uvw[:,i] += m[j]*xyz_ij[:,i,j]*dt/r_ij[i,j]**3 for i in range(N_ast): uvwa[:,i] += m[j]*xyza_ij[:,i,j]*dt/ra_ij[i,j]**3 xyz += uvw*dt xyza += uvwa*dt xyzs.append(xyz.tolist()) xyzas.append(xyza.tolist()) xyzs = np.array(xyzs).transpose(2,1,0) xyzas = np.array(xyzas).transpose(2,1,0) def animate(n): for i in range(len(m)): xyz = xyzs[i] traces[i].set_data(xyz[0,:n],xyz[1,:n]) traces[i].set_3d_properties(xyz[2,:n]) pts[i].set_data(xyz[0,n],xyz[1,n]) pts[i].set_3d_properties(xyz[2,n]) for i in range(N_ast): pt_asts[i].set_data(xyzas[i,0,n],xyzas[i,1,n]) pt_asts[i].set_3d_properties(xyzas[i,2,n]) return traces+pts+pt_asts ani = animation.FuncAnimation(fig, animate, range(1,N,100), interval=10, blit=True) plt.show() ani.save(name)
以上就是Python算法繪制特洛伊小行星群實(shí)現(xiàn)示例的詳細(xì)內(nèi)容，更多關(guān)于Python算法繪制特洛伊小行星群的資料請關(guān)注本站其它相關(guān)文章！
版權(quán)聲明：本站文章來源標(biāo)注為YINGSOO的內(nèi)容版權(quán)均為本站所有，歡迎引用、轉(zhuǎn)載，請保持原文完整并注明來源及原文鏈接。禁止復(fù)制或仿造本網(wǎng)站，禁止在非www.sddonglingsh.com所屬的服務(wù)器上建立鏡像，否則將依法追究法律責(zé)任。本站部分內(nèi)容來源于網(wǎng)友推薦、互聯(lián)網(wǎng)收集整理而來，僅供學(xué)習(xí)參考，不代表本站立場，如有內(nèi)容涉嫌侵權(quán)，請聯(lián)系alex-e#qq.com處理。
關(guān)注官方微信