Python算法绘制特洛伊小行星群实现示例 Python算法绘制特洛伊小行星群

最小势能点
拉格朗日点
特洛伊小行星群

书接上文
用Python搓一个太阳系
你们要的3D太阳系
3体人真的存在吗
太长不看版

文章图片

文章图片

最小势能点在由两个大质量物体构成的重力系统中，有一些特殊的区域会在两个天体的顶级拉扯之下达到平衡，这些点就是拉格朗日点。而所谓平衡并非受力平衡，而是要求这个区域的物体会跟着双星系统以相同的角速度运动。

文章图片

根据上帝是个胖子这个假定，状态稳定意味着低势能。所以在解析求解拉格朗日点之前，我们可以试着画出这个双星系统的势能分布。

文章图片

接下来搞一下太阳和木星：

文章图片

文章图片

可见木星在太阳的引力场下根本无法自己，但若把坐标系调整一下，会看到木星虽小，却还是有自己地盘的，毕竟也是有诸多卫星的，这就意味着木星和太阳之间必然存在一些相对平衡的位置。
为了看得更加仔细，取对数是个不错的方法

文章图片

import numpy as npimport matplotlib.pyplot as pltfrom matplotlib import cm,tickerR = 7.1492e7M1,M2 = 1.9891e30, 1.8982e27G = 6.67e-11mu = M2/M1R1,R2 = np.array([mu,1])/(1+mu)*Rx,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**3absH = np.abs(H)absH[absH>1e14] = 1e14#去除奇点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()

拉格朗日点 【Python算法绘制特洛伊小行星群实现示例】公式预警→_→
根据刚刚的图可以看出，一般天体都会有一个属于自己的私密区域，在这个区域里，别的天体的引力作用甚微，此区域即希尔球，拉格朗日点则是两个天体希尔球的分界处。

文章图片

在极坐标下，可得

文章图片

文章图片

对于木星来说，五个拉格朗日点一般默认为

文章图片

特洛伊小行星群

文章图片

参考此前的太阳系行星位置，得到其三维图

文章图片

from os import cpu_countimport numpy as npfrom numpy.random import randimport matplotlib.pyplot as pltfrom matplotlib import animationau,G,RE,ME = 1.48e11,6.67e-11,1.48e11,5.965e24m = np.array([3.32e5,0.055,0.815,1,0.107,317.8])*ME*Gr = np.array([0,0.387,0.723,1,1.524,5.203])*REv = np.array([0,47.89,35.03,29.79,24.13,13.06])*1000theta = rand(len(m))*np.pi*2theta[-1] = 0#木星初始角度为0cTheta,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 = 100x_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.5y_ast = y_ast*(2*y_flag-1)m_ast = rand(N_ast)*1e20r_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*2phi = (rand(N_ast)-0.5)*0.3#给一个随机的小倾角cTheta,sTheta = x_ast/r_ast, y_ast/r_astcPhi,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 = 10000dt = 3600*50ts =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]**3for i in range(N_ast):uvwa[:,i] += m[j]*xyza_ij[:,i,j]*dt/ra_ij[i,j]**3xyz += uvw*dtxyza += uvwa*dtxyzs.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_astsani = animation.FuncAnimation(fig, animate, range(1,N,100), interval=10, blit=True)plt.show()ani.save(name)

以上就是Python算法绘制特洛伊小行星群实现示例的详细内容，更多关于Python算法绘制特洛伊小行星群的资料请关注脚本之家其它相关文章！