问题背景

你一定会好奇，月球绕着地球做圆周运动，地球又绕着太阳转，太阳绕着银河系…….那么以任意其中一个旋转中心为相对坐标中心，月球的运动轨迹是怎么样的？
假设太阳为中心，那么月球的运动轨迹又是怎么样的呢？
本文将从数学构建坐标系的方法，阐述下任意星体相对其它星体轨迹方程，并运用python matplotlib这个强大的绘图库进行日地月运动模型的模拟

模型构建

现在假设历太阳-地球-月亮这个系统，太阳是旋转中心，地球相对太阳做角速度为ω1ω1\omega _1的匀速圆周运动，月球相对地球做ω2ω2\omega_2的匀速圆周运动。

以太阳为原点构建三维坐标系，即太阳坐标点(x0,y0,z0)(x0,y0,z0)\left(x_0,y_0,z_0\right) = (0,0,0)(0,0,0)\left(0, 0, 0\right)
现假设地球公转的轨道平面（黄道面）在x-y轴平面上。

已知月球的轨道平面（白道面）与黄道面（地球的公转轨道平面）保持著5.145 396°的夹角，即与x-y轴平面呈5.145 396°的夹角

由于与x-y轴平面的夹角可以是任意方向的，为建模方便设月球公转轨道面沿y轴方向向上倾斜5.145 396°且初始状态时：太阳、地球、月球在y-z平面上（后文用φφ\varphi 表示月球公转轨道平面与地球公转轨道平面的固定夹角）
假设地球公转半径为r1r1r_1, 月球公转半径为r2r2r_2 ，
对于地球，假设其在t时刻时的坐标为(x1,y1,z1)(x1,y1,z1)\left( x_1,y_1,z_1 \right) ，则有

⎧⎩⎨⎪⎪x1=x0+r1⋅cos(ω1t)y1=y0+r1⋅sin(ω1t)z1=z0+0{x1=x0+r1⋅cos⁡(ω1t)y1=y0+r1⋅sin⁡(ω1t)z1=z0+0

\left\{\begin{array}{lr} x_1 = x_0 + r_1\cdot \cos \left( \omega_1 t \right) \\y_1 = y_0 + r_1\cdot \sin \left( \omega_1 t \right) \\z_1 = z_0 + 0 & \end{array} \right.
现在对月球进行的运动轨迹进行建模: 任意t时刻，月球运行轨迹满足如下方程

{z2=z1+(y2−y1)⋅tanφ(x2−x1)2+(y2−y1)2+(z2−z1)2=r22{z2=z1+(y2−y1)⋅tan⁡φ(x2−x1)2+(y2−y1)2+(z2−z1)2=r22

\left\{\begin{array}{lr} z_2=z_1 + \left( y_2-y_1 \right) \cdot \tan \varphi \\\left( x_2-x_1 \right) ^2+\left( y_2-y_1 \right) ^2+\left( z_2-z_1 \right) ^2 =r_2^2 &\end{array} \right.
月球初始方向向量为 (0,cosφ,sinφ)(0,cos⁡φ,sin⁡φ)\left(0, \cos\varphi , \sin \varphi\right)
在t时刻，月球转动的夹角 ω2tω2t\omega_2 t, 此时的方向向量为 (x2−x1,y2−y1,z2−z1)(x2−x1,y2−y1,z2−z1)\left(x_2 - x_1, y_2 - y_1, z_2 - z_1\right)
由向量夹角公式有，

ω2t=(y2−y1)⋅cosφ+(z2−z1)⋅sinφω2t=(y2−y1)⋅cos⁡φ+(z2−z1)⋅sin⁡φ

\omega_2 t = \left(y_2 - y_1\right)\cdot\cos\varphi + \left(z_2 - z_1\right)\cdot\sin\varphi
联立上式，可得

⎧⎩⎨⎪⎪⎪⎪⎪⎪x2=x1+r2⋅sin(ω2t)y2=y1+r2⋅cos(ω2t)cosφ(1+tan2φ)z2=z1+(y2−y1)⋅tanφ{x2=x1+r2⋅sin⁡(ω2t)y2=y1+r2⋅cos⁡(ω2t)cos⁡φ(1+tan2⁡φ)z2=z1+(y2−y1)⋅tan⁡φ

\left\{\begin{array}{lr} x_2 = x_1 + r_2\cdot \sin\left( \omega_2 t \right) \\y_2 = y_1 + \frac{r_2\cdot \cos\left( \omega_2 t \right)}{\cos\varphi\left(1+ \tan^2\varphi\right) }\\z_2 = z_1 + \left(y_2 - y_1\right)\cdot\tan\varphi & \end{array} \right.

python实现

查阅相关资料可以知道
日地距离：约1.5亿千米,即一个天文单位.
月地距离：约38.4万千米
因此 r1r2=390.625r1r2=390.625\frac{r_1}{r_2} = 390.625

出于3D绘图的直观性，这里假设r1=10r1=10r_1 = 10 ,r2=1r2=1r_2 = 1
另外由于地球绕太阳公转一周，刚好十二个月，即月球绕地球公转角速度是地球绕太阳公转角速度的12倍
因此可设ω1=2πω1=2π\omega_1 = 2\pi, ω2=24πω2=24π\omega_2 = 24\pi

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 animmationr1 = 10
r2 = 1
omega1 = 2 * np.pi
omega2 = 24 * np.pi
phi = 5.1454 * np.pi / 180 def update(data):global line1, line2 , line3line1.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, line3ti = 0t = t_drange[np.mod(ti, t_dlen)]xt1 = x0 + r1 * np.cos(omega1 * t)yt1 = y0 + r1 * np.sin(omega1 * t)zt1 = z0 + 0xt2 = 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, phiglobal 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 = z0xt2 = 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演示了太阳-地球-月球运动轨迹模型动画
（红色球点表示太阳，蓝色球点表示地球，橙色球点表示月球，红色轨迹表示地球公转轨道，紫色轨迹表示月球相对地球轨道，蓝色轨迹表示月球相对太阳轨迹):

可以发现：月球在地球的公转轨道上是螺旋前进的,这与实际上是一致的。
这里的模型假设了太阳是静止的坐标系中心，如果再引入银河系中心，让太阳动起来，月球的轨迹会是怎么样呢？这实际并不难实现，只要对(x0,y0,z0)(x0,y0,z0)\left(x_0,y_0,z_0\right) 进行方程赋值即可，感兴趣的可以试试。

另外，你也可以用上面的模型进行自转模型的模拟，因为月球的自转周期与其绕地球公转的周期是一样的，
将上面代码稍作修改，将太阳位置替换成地球、地球位置替换成月球、原月球的公转轨道看成其自转轨迹，黄色球表示其某个面（可以设置初始时面对地球）。通过设置角速度一样，你就可以直观地看到——为什么我们在地球上总是看不到月球的另一面。

当然，如果你觉得模型太简单了，想更接近实际一点，也可以将地球公转轨道设置成椭圆，这也是很容易实现的，只要你知道椭圆的参数方程即可。

最后说下傅里叶里级数。是的，你应该发现了：星体运动无限迭代下去，就是一个曲线的傅里叶级数展开，只要迭代次数是有限的，星体的运动轨迹肯定是具有周期性的。

运用Python 模拟太阳-地球-月亮运动模型相关推荐

Python模拟登录，matplotlib模块，Python模拟太阳-地球-月亮运动模型
前言利用python模拟太阳-地球-月亮运动模型. 让我们愉快地开始吧~ 开发工具 **Python版本:**3.6.4 相关模块: pygame模块: matplotlib模块: numpy模块: ...
Python模拟太阳-地球-月亮运动模型
作者 | Charles,cv方向在读研究生.[Charles 的皮卡丘]专注于分享有趣好玩的Python小项目(AI.爬虫等等). 来源 | Charles 的皮卡丘编辑 | Jane [导语]春 ...
openGl编程实现一个太阳地球月亮的一个简单运动系统
一. 项目目的使用openGl编程实现一个太阳地球月亮的一个简单运动系统,要求实现三维转动.点光源变化.纹理映射及阴影等效果二. 任务实现 \1. 满足三者实际大小/距离的比例关系: \2. 满足 ...
webgl 绘制太阳地球月亮
目录 1.开发环境 2.内容说明 1.计算球体的坐标和纹理 2.求顶点索引 3.数据加载到缓存中 4.绘制球体 5.其他 3.运行结果及代码下载 1.开发环境浏览器 : 火狐 firefox(配置参 ...
HGE-7 模拟太阳和月亮的运转
1 /* 2 ** Haaf's Game Engine 1.8 3 ** Copyright (C) 2003-2007, Relish Games 4 ** hge.relishgames.com ...
Canvas模拟太阳地球月球的运动过程
先看效果图代码 package com.test.paintdemo.pathrelate;import android.content.Context; import android.graphi ...
CSS太阳地球月亮转圈圈loading
一个css例子 css太阳月亮地球转转转知识点总结: justify-content: center; justify-content 用于设置或检索弹性盒子元素在主轴(横轴)方向上的对齐方式 ...
matlab地球公转,Unity模拟太阳地球月球公转自转
1. 打开Unity编辑器,创建三个sphere,依次重命名为Sun,Earth,Moon.将三个球体大小比例控制在5:3:1,并适当调整位置. 2. 分别为三个小球附上材质 3. 新建脚本文件,重命 ...
three.js旋转，材质，灯光使用 —— 太阳地球月亮运动
import * as THREE from 'three' //引入控制器 import { OrbitControls } from 'three/examples/jsm/controls/Or ...
OpenGl太阳地球月亮运动系统
在讲解这个运动系统,首先我们的来讲解OpenGl里有关几个图形变换的知识,这里就以球为例,我们需要知道将球平移,旋转的2个知识.因为系统必须用到平移和旋转. 1 平移变换: glTranslatef( ...

运用Python 模拟太阳-地球-月亮运动模型

问题背景

模型构建

python实现

运用Python 模拟太阳-地球-月亮运动模型相关推荐

最新文章

热门文章