运用Python 模拟太阳-地球-月亮运动模型
2018-01-01 17:54
691 查看
问题背景
你一定会好奇,月球绕着地球做圆周运动,地球又绕着太阳转,太阳绕着银河系…….那么以任意其中一个旋转中心为相对坐标中心,月球的运动轨迹是怎么样的?假设太阳为中心,那么月球的运动轨迹又是怎么样的呢?
本文将从数学构建坐标系的方法,阐述下任意星体相对其它星体轨迹方程,并运用python matplotlib这个强大的绘图库进行日地月运动模型的模拟
模型构建
现在假设历太阳-地球-月亮这个系统,太阳是旋转中心,地球相对太阳做角速度为ω1的匀速圆周运动,月球相对地球做ω2的匀速圆周运动。以太阳为原点构建三维坐标系, 即太阳坐标点(x0,y0,z0) = (0,0,0)
现假设地球公转的轨道平面(黄道面)在x-y轴平面上。
已知月球的轨道平面(白道面)与黄道面(地球的公转轨道平面)保持著5.145 396°的夹角,即与x-y轴平面呈5.145 396°的夹角
由于与x-y轴平面的夹角可以是任意方向的,为建模方便设月球公转轨道面沿y轴方向向上倾斜5.145 396°且初始状态时:太阳、地球、月球在y-z平面上(后文用φ表示月球公转轨道平面与地球公转轨道平面的固定夹角)
假设地球公转半径为r1, 月球公转半径为r2,
对于地球,假设其在t时刻时的坐标为(x1,y1,z1),则有
⎧⎩⎨⎪⎪x1=x0+r1⋅cos(ω1t)y1=y0+r1⋅sin(ω1t)z1=z0+0
现在对月球进行的运动轨迹进行建模: 任意t时刻,月球运行轨迹满足如下方程
{z2=z1+(y2−y1)⋅tanφ(x2−x1)2+(y2−y1)2+(z2−z1)2=r22
月球初始方向向量为(0,cosφ,sinφ)
在t时刻,月球转动的夹角ω2t, 此时的方向向量为(x2−x1,y2−y1,z2−z1)
由向量夹角公式有,
ω2t=(y2−y1)⋅cosφ+(z2−z1)⋅sinφ
联立上式,可得
⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪x2=x1+r2⋅sin(ω2t)y2=y1+r2⋅cos(ω2t)cosφ(1+tan2φ)z2=z1+(y2−y1)⋅tanφ
python实现
查阅相关资料可以知道日地距离:约1.5亿千米,即一个天文单位.
月地距离:约38.4万千米
因此 r1r2=390.625
出于3D绘图的直观性,这里假设r1=10,r2=1
另外由于地球绕太阳公转一周,刚好十二个月,即月球绕地球公转角速度是地球绕太阳公转角速度的12倍
因此可设ω1=2π, ω2=24π
t 从 0 ~ 1 表示星体公转一周, python 代码模型实现如下:
import numpy as np import matplotlib as mpl mpl.use("TkAgg") from matplotlib import pyplot as plt from mpl_toolkits.mplot3d import Axes3D import matplotlib.animation as animmation r1 = 10 r2 = 1 omega1 = 2 * np.pi omega2 = 24 * np.pi phi = 5.1454 * np.pi / 180 def update(data): global line1, line2 , line3 line1.set_data([data[0], data[1]]) line1.set_3d_properties(data[2]) line2.set_data([data[3], data[4]]) line2.set_3d_properties(data[5]) line3.set_data([data[6], data[7]]) line3.set_3d_properties(data[8]) return line1,line2,line3, def init(): global line1, line2, line3 ti = 0 t = t_drange[np.mod(ti, t_dlen)] xt1 = x0 + r1 * np.cos(omega1 * t) yt1 = y0 + r1 * np.sin(omega1 * t) zt1 = z0 + 0 xt2 = xt1 + r2 * np.sin(omega2 * t) yt2 = yt1 + r2 * np.cos(omega2 * t)/(np.cos(phi) * (1 + np.tan(phi) ** 2)) zt2 = zt1 + (yt2 - yt1) * np.tan(phi) xt21 = xt1 + r2 * np.sin(2 * np.pi * t_range) yt21 = yt1 + r2 * np.cos(2 * np.pi * t_range)/(np.cos(phi) * (1 + np.tan(phi) ** 2)) zt21 = zt1 + (yt21 - yt1) * np.tan(phi) line1, = ax.plot([xt1], [yt1], [zt1], marker='o', color='blue',markersize=8) line2, = ax.plot([xt2], [yt2], [zt2], marker='o', color='orange',markersize=4) line3, = ax.plot(xt21, yt21, zt21, color='purple') return line1,line2,line3 def data_gen(): #global x0,y0,z0,ti, t_drang, t_range, omega1, omega2, phi global x0,y0,z0,t_dlen #while true: data = [] for ti in range(1,t_dlen): t = t_drange[ti] xt1 = x0 + r1 * np.cos(omega1 * t) yt1 = y0 + r1 * np.sin(omega1 * t) zt1 = z0 xt2 = xt1 + r2 * np.sin(omega2 * t) yt2 = yt1 + r2 * np.cos(omega2 * t)/(np.cos(phi) * (1 + np.tan(phi) ** 2)) zt2 = zt1 + (yt2 - yt1) * np.tan(phi) xt21 = xt1 + r2 * np.sin(2 * np.pi * t_range) yt21 = yt1 + r2 * np.cos(2 * np.pi * t_range)/(np.cos(phi) * (1 + np.tan(phi) ** 2)) zt21 = zt1 + (yt21 - yt1) * np.tan(phi) data.append([xt1, yt1, zt1, xt2, yt2, zt2, xt21, yt21, zt21]) return data #yield (xt1, yt1, zt1, xt2, yt2, zt2, xt21, yt21, zt21) t_range = np.arange(0, 1 + 0.005, 0.005) t_drange = np.arange(0, 1, 0.005 ) t_len = len(t_range) t_dlen = len(t_drange) #sun's coordination x0 = 0 y0 = 0 z0 = 0 #earth's orbit x1 = x0 + r1 * np.cos(omega1 * t_range) y1 = y0 + r1 * np.sin(omega1 * t_range) z1 = z0 + np.zeros(t_len) #moon's orbit x2 = x1 + r2 * np.sin(omega2 * t_range) y2 = y1 + r2 * np.cos(omega2 * t_range)/(np.cos(phi) * (1 + np.tan(phi) ** 2)) z2 = z1 + (y2 - y1) * np.tan(phi) f = plt.figure(figsize=(6,6)) ax = f.add_subplot(111,projection='3d') #plt.rcParams['animation.ffmpeg_path'] = r"C:\Program Files\ffmpeg\bin\ffmpeg" #plt.rcParams['animation.convert_path'] = r"C:\Program Files\ImageMagick-7.0.7-Q16\magick.exe" ax.set_aspect('equal') ax.set_title("Sun-Earth-Moon Model") ax.plot([0], [0], [0], marker='o', color= 'red', markersize=16) ax.plot(x1, y1, z1, 'r') ax.plot(x2, y2, z2, 'b') ax.set_xlim([-(r1 + 2), (r1 + 2)]) ax.set_ylim([-(r1 + 2), (r1 + 2)]) ax.set_zlim([-5, 5]) # line1 update Earth's track dynamically # line2 update Moon's track dynamically # line3 update Moon's orbit to earth line1, = ax.plot([], [], [], marker='o', color='blue',markersize=8,animated = True) line2, = ax.plot([], [], [], marker='o', color='orange',markersize=4,animated = True) line3, = ax.plot([], [], [], color='purple',animated = True) #red sphere for Sun, blue sphere for Earth, orange sphere for Moon ani = animmation.FuncAnimation(f, update, frames = data_gen(), init_func = init,interval = 20) #ffwriter = animmation.ffmpegwriter(fps = 200) #ani.save('planet.gif', writer='imagemagick', fps=40) #ani.save('planet.gif', writer = ffwriter) plt.show()
(注:上述代码注释部分的动画保存需要安装ffmpeg或者imagemagick的支持)
上述代码,通过matplotlib animation演示了太阳-地球-月球运动轨迹模型动画
[这里本来上要上传gif 图像,但是csdn 老是没有上传成功,所以作罢]
截取动画的几个页面如下(红色球点表示太阳,蓝色球点表示地球,橙色球点表示月球,红色轨迹表示地球公转轨道,紫色轨迹表示月球相对地球轨道,蓝色轨迹表示月球相对太阳轨迹:
可以发现:月球在地球的公转轨道上是螺旋前进的,这与实际上是一致的。
这里的模型假设了太阳是静止的坐标系中心,如果再引入银河系中心,让太阳动起来,月球的轨迹会是怎么样呢?这实际并不难实现,只要对(x0,y0,z0) 进行方程赋值即可,感兴趣的可以试试。
另外,你也可以用上面的模型进行自转模型的模拟,因为月球的自转周期与其绕地球公转的周期是一样的,
将上面代码稍作修改,将太阳位置替换成地球、地球位置替换成月球、原月球的公转轨道看成其自转轨迹,黄色球表示其某个面(可以设置初始时面对地球)。通过设置角速度一样,你就可以直观地看到——为什么我们在地球上总是看不到月球的另一面。
当然,如果你觉得模型太简单了,想更接近实际一点,也可以将地球公转轨道设置成椭圆,这也是很容易实现的,只要你知道椭圆的参数方程即可。
最后说下傅里叶里级数。是的,你应该发现了:星体运动无限迭代下去,就是一个曲线的傅里叶级数展开,只要迭代次数是有限的,星体的运动轨迹肯定是具有周期性的。
相关文章推荐
- 基于C#的2D太阳、地球、月亮运动轨迹模拟实现
- 关于OpenGL入门中地球,月亮,太阳运动模拟无法显示的问题(只显示黑框)
- OpenGL “太阳、地球和月亮”天体运动动画 例子
- SVG 模拟太阳 地球 月亮旋转
- 模拟太阳、地球和月亮
- openGL+VS2010的例程--太阳地球月球运动模型(三维)
- WPF太阳、地球、月球运动轨迹模拟
- OpenGL:太阳、地球、月球模型模拟动画
- openGL+VS2010的例程--太阳地球月球运动模型增强版(三维)
- Canvas模拟太阳地球月球的运动过程
- opengl: 太阳地球和月亮
- 太阳、地球、月亮
- OpenGL入门记录--“太阳,地球和月亮”系统代码
- OpenGL模拟太阳、地球、月球公转自转
- 太阳、地球、月亮
- [HTML5-SVG]使用svg、CSS3、raphaeljs:模拟月球绕地球,地球绕太阳旋转
- python蒙特卡洛算法模拟赌博模型
- python udp通信模型,模拟TFTP客户端
- Transform类 RotateAround() 太阳、地球、月亮围绕旋转
- Python-yield模拟生产消费者模型