用Python搓一个太阳系
文章目录
- 日地月三体
- 日地火
- 太阳系
你们要的3D太阳系
3D太阳系+特洛伊小行星群
图片上传之后不知为何帧率降低了许多。。。
日地月三体
所谓三体,就是三个物体在重力作用下的运动。由于三点共面,所以三个质点仅在重力作用下的运动轨迹也必然无法逃离平面。
三体运动所遵循的规律就是古老而经典的万有引力
F⃗=Gmimjr2e⃗r\vec F=\frac{Gm_im_j}{r^2}\vec e_r F=r2Gmimjer
则对于mim_imi而言,
midv⃗idt=Gmimjrij3r⃗ijm_i\frac{\text d\vec v_i}{\text dt}=\frac{Gm_im_j}{r_{ij}^3}\vec r_{ij} midtdvi=rij3Gmimjrij
且
dr⃗idt=v⃗i\frac{\text d\vec r_i}{\text dt}=\vec v_i dtdri=vi
将其写为差分形式
v⃗i=∑j≠iGmjrij3r⃗ijdtr⃗i=v⃗idt\begin{aligned} \vec v_i&=\sum_{j\not=i}\frac{Gm_j}{r_{ij}^3}\vec r_{ij}\text dt\\ \vec r_i&= \vec v_i\text dt \end{aligned} viri=j=i∑rij3Gmjrijdt=vidt
由于我们希望观察三体运动的复杂形式,而不关系其随对应的宇宙星体,所以不必考虑单位制,将其在二维平面坐标系中拆分,令v⃗=(u,v)\vec v=(u,v)v=(u,v),则
ui+=∑j≠iGmj(xj−xi)dt(xi−xj)2+(yi−yj)23vi+=∑j≠iGmj(yj−yi)dt(xi−xj)2+(yi−yj)23xi+=u⃗idtyi+=v⃗idt\begin{aligned} u_i&+=\sum_{j\not=i}\frac{Gm_j(x_j-x_i)\text dt}{\sqrt{(x_i-x_j)^2+(y_i-y_j)^2}^3}\\ v_i&+=\sum_{j\not=i}\frac{Gm_j(y_j-y_i)\text dt}{\sqrt{(x_i-x_j)^2+(y_i-y_j)^2}^3}\\ x_i&+= \vec u_i\text dt\\ y_i&+= \vec v_i\text dt \end{aligned} uivixiyi+=j=i∑(xi−xj)2+(yi−yj)23Gmj(xj−xi)dt+=j=i∑(xi−xj)2+(yi−yj)23Gmj(yj−yi)dt+=uidt+=vidt
太阳、地球和月亮就是一个典型的三体系统,其中太阳质量为1.989×1030kg1.989×10^{30}kg1.989×1030kg,地球质量为5.965×1024kg5.965×10^{24}kg5.965×1024kg,月球质量为7.342✕1022kg7.342✕10^{22}kg7.342✕1022kg,万有引力常数为G=6.67×10−11N⋅m2/kg2G=6.67×10^{-11}N·m2/kg^2G=6.67×10−11N⋅m2/kg2。地月距离为3.8×108m3.8\times10^8m3.8×108m;日地距离为1.5×1011m1.5\times10^{11}m1.5×1011m;地球公转速度为28.8km/s28.8km/s28.8km/s;月球公转速度为1km/s1km/s1km/s,则各参数初始化为
#后续代码主要更改这里的参数
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 animationm = [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*dty += v*dtxs.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_textani = animation.FuncAnimation(fig, animate, range(N), interval=10, blit=True)plt.show()
ani.save("3.gif")
日地火
质量MMM | GMGMGM | 与太阳距离 | 公转速度 | |
---|---|---|---|---|
地球 | 5.965×1024kg5.965×10^{24}kg5.965×1024kg | 3.98×10143.98×10^{14}3.98×1014 | 1.5×1011m1.5\times10^{11}m1.5×1011m | 28.8km/s28.8km/s28.8km/s |
火星 | 6.4171✕1023kg6.4171✕10^{23}kg6.4171✕1023kg | 4.28×10134.28×10^{13}4.28×1013 | 1.52A.U.=2.28×10111.52 A.U.=2.28\times10^{11}1.52A.U.=2.28×1011 | 24km/s24km/s24km/s |
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
得到
这个运动要比月球的运动简单得多——前提是开上帝视角,俯瞰太阳系。如果站在地球上观测火星的运动,那么这个运动可能相当带感
所以这都能找到规律,托勒密那帮人也真够有才的。
太阳系
由于太阳和其他星体之间的质量相差悬殊,所以太阳系内的多体运动,都将退化为二体问题,甚至如果把太阳当作不动点,那就成了单体问题了。
尽管如此,我们还是尽可能地模仿一下太阳系的运动情况
质量 | 半长轴(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.965e24m = 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*2x = 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*dty += v*dtxs.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)
由于外圈的行星轨道又长速度又慢,而内层的刚好相反,所以这个图很难兼顾,观感上也不太好看。
如果只画出木星之前的星体,顺便加上小行星带,可能会好一些。
通过这个图就能看出来,有一颗小行星被木星弹了过来,直冲冲地向地球赶来,幸好又被太阳弹了出去,可见小行星还是挺危险的,好在这只是个假想图。
用Python搓一个太阳系相关推荐
- 用Python搓一个黑洞
文章目录 简介 单位制 观测绘图 简介 黑洞图像大家都知道,毕竟前几年刚发布的时候曾火遍全网,甚至都做成表情包了. 问题在于,凭什么认为这就是黑洞的照片,而不是一个甜甜圈啥的给整模糊了得到的呢?有什么 ...
- Python 用Ursina 3D引擎做一个太阳系行星模拟器
这次,我们再来用Ursina引擎来做一个太阳系行星模拟器吧! 想要了解Ursina 3D引擎的基本使用方法的话,查看我的另一篇文章: 手把手教你用Python编一个<我的世界> 1. 认识 ...
- 用Python画一个3D太阳系
用Python画一个平面的太阳系得到一些朋友的欣赏,然后有同学提出了绘制三维太阳系的要求. 从Python画图的角度来说,三维太阳系其实并不难,问题在于八大行星对黄道面的倾斜太小,所以尽管画个三维的图 ...
- python绘制太阳系模型_画一个太阳系的模型
1 .画一个太阳系的模型? 2 .地球公转对地球有什么影响? 3 .摆在摆钟里是怎样工作的?摆在摆钟里起到了什么作用?人们为什么会选择摆作为摆钟 的控制核心?摆的快慢与那些因素有关? 4 . Xx 家 ...
- python绘制动态模拟图-如何利用Python动态模拟太阳系运转
前言 提到太阳系,大家可能会想到哥白尼和他的日心说,或是捍卫.发展日心说的斗士布鲁诺,他们像一缕光一样照亮了那个时代的夜空,对历史感兴趣的小伙伴可以深入了解一下,这里就不多说了. 太阳以巨大的引力使周 ...
- 用html+ccs3就能 做出一个太阳系行星
2019独角兽企业重金招聘Python工程师标准>>> 做一个太阳系八大行星的运转动画,不包括行星的卫星,所有行星围绕太阳公转,行星采用纯色,暂时没有自转. 效果静态图: 动画中包括 ...
- 如何利用Python动态模拟太阳系运转
更多编程教程请到:菜鸟教程 https://www.piaodoo.com/ 友情链接:好看站 http://www.nrso.net/ 高州阳光论坛https://www.hnthzk.com/ 前 ...
- python绘制太阳系_如何利用Python动态模拟太阳系运转
前言 提到太阳系,大家可能会想到哥白尼和他的日心说,或是捍卫.发展日心说的斗士布鲁诺,他们像一缕光一样照亮了那个时代的夜空,对历史感兴趣的小伙伴可以深入了解一下,这里就不多说了. 太阳以巨大的引力使周 ...
- python怎么画地球绕太阳转_如何利用Python动态模拟太阳系运转
前言 提到太阳系,大家可能会想到哥白尼和他的日心说,或是捍卫.发展日心说的斗士布鲁诺,他们像一缕光一样照亮了那个时代的夜空,对历史感兴趣的小伙伴可以深入了解一下,这里就不多说了. 太阳以巨大的引力使周 ...
最新文章
- iOS 疑难杂症— — 收到推送显示后自动消失的问题
- Oracle buffer状态深入剖析
- maven shade
- 竟然被尤雨溪点赞了:我给Vue生态贡献代码的这一年
- 操作系统文件编程知识
- 利用xlwt写excel并进行单元格的合并
- linux 和服务通讯,Android 的Activity和Service之间的通信
- easyexcel和poi是否有版本冲突_easyexcel--解决poi大文件发生OOM问题
- python是什么时候出现的_python诞生于什么时候
- tensorflow之variables_to_restore
- C# 判断输入的字符是不是数字
- VS2015 卸载与安装相关问题(包丢失以及没有WIN控制台)
- openpyxl进行excel的整行复制
- 计算机图形学(六)-光栅化、采样、走样与反走样、滤波与卷积
- 三相全桥整流电路_三相桥式全控整流电路的工作原理
- 《关键对话》教你如何摆脱沟通困境
- cherry键盘win键锁定的问题
- pandas删除行删除列,增加行增加列
- 人生=亲情+爱情+金钱+理想+友情?
- STM32单片机实现连接USB摄像头