程序员炫酷溜娃!用代码画地球日月的动态轨道模型

Posted AI科技大本营

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了程序员炫酷溜娃!用代码画地球日月的动态轨道模型相关的知识,希望对你有一定的参考价值。

作者|天元浪子

来源|Python作业辅导员

无言相守45亿年,太阳、地球和月球这三个好基友究竟是怎样的关系呢?从孩提时代我就一直在想,要是能有一个可以直观演示太阳、地球和月球运行轨迹的模型就好了。今天,我终于实现了小时候的梦想:用WxGL画出了太阳、地球和月球的动态轨道模型。配上简单的解说,小朋友也可以秒懂四季更迭、日蚀月蚀、黄赤交角。

在开始绘制模型前,让我们先来了解一下太阳、地球和月球的起源,以及它们的大小、远近和行踪路线。

1、主流的大爆炸理论认为,宇宙大爆炸发生在138.17亿年前。

2、约66亿年前,一颗超新星爆炸后产生了一团星云,逐渐形成了太阳系的雏形。

3、约46亿年前,这团星云最中心的巨大氢球,开始产生聚变反应,太阳就此诞生。

4、约45亿年前,一颗火星大小的行星“忒伊亚”撞上了地球的雏形,体积增加了近一倍,地球诞生。

5、在这次撞击中,“忒伊亚”残骸的一部分形成了月球。

6、太阳半径696000km,地球半径6371km,月球半径1738km,三者之比约为400:3.67:1。

7、太阳绕银河系中心公转周期约2.5亿年。

8、太阳也在自转,不过因为是等离子体,所以不同纬度有不同的自转速度。赤道区域自转最快,周期为24.47天。

9、地球自转周期为1天,公转周期为365.2564天(恒星年)或365.2422天(回归年)。

10、地球公转轨道是一个近似圆的椭圆,长半轴为149600000km,短半轴为149580000km,曲率为0.016722。

11、地球自转轴并不垂直于地球公转轨道面,地球公转轨道面(黄道面)与地球赤道面夹角为23.43° 。

12、月球公转周期为27.32日(恒星月)或29.53天(朔望月),月球自转周期为27.32日(恒星月)。

13、月球轨道面对黄道面的倾角是5.145°,月球赤道面对黄道面的倾角是1.543°,月球轨道面与月球赤道面夹角为6.688° 。

14、月球轨道是一个椭圆,长半轴为385000km,离心率是0.0549。

15、地球公转和自转方向、月球公转和自转方向、太阳的自转方向,是一致的,都是自西向东。

了解了这些,就可以开始绘制模型了。关于WxGL模块的安装,请参考我近期的文章,这里不再赘述。下面是完整的代码,大约一百余行

# -*- coding: utf-8 -*-

import numpy as np
from scipy.spatial.transform import Rotation as sstr
from wxgl import wxplot as plt

R_SUN = 1                           # 以太阳半径696000km为1个单位
R_EARTH = 30 * 6371/696000          # 地球半径6371km,此处放大30倍
R_MONTH = 50 * 1738/696000          # 月球半径1738km,此处放大50倍
A_EARTH = (149600000/696000)/100    # 地球公转轨道长半轴149600000km,此处缩小100倍
E_EARTH = 0.016722                  # 地球公转轨道离心率
A_MONTH = 385000/696000             # 月球公转轨道长半轴385000km
E_MONTH = 0.0549                    # 月球公转轨道离心率

T = 36                              # 每天36个计数周期
T_SUN_SELF = int(24.47 * T)         # 太阳自转周期24.47天
T_EARTH_SELF = 1 * T                # 地球自转周期1天
T_MONTH_SELF = int(27.32 * T)       # 月球自转周期27.32天
T_EARTH = int(365.2564/3 * T)       # 地球公转周期365.2564天,此处取1/3
T_MONTH = int(27.32 * T)            # 月球公转周期27.32天

A_E_ORBIT = 23.43                   # 地球轨道面与地球赤道面夹角
A_M_ORBIT = 23.43 + 5.145           # 月球轨道面与地球赤道面夹角
A_M_EQUATOR = 23.43 - 1.543         # 月球赤道面与地球赤道面夹角

# 地球轨道倾角旋转器
ROTATOR_E_ORBIT = sstr.from_euler('xyz', (0, -A_E_ORBIT, 0), degrees=True)

# 月球轨道倾角旋转器
ROTATOR_M_ORBIT = sstr.from_euler('xyz', (0, -A_M_ORBIT, 0), degrees=True)

def rotate_merge(av1, av2):
    """将两个轴角旋转合并为一个

    av1         - 元组,首元素(浮点型)为旋转角度(逆时针为正,右手定则),尾元素(列表或元组)为旋转向量
    av2         - 元组,首元素(浮点型)为旋转角度(逆时针为正,右手定则),尾元素(列表或元组)为旋转向量
    """

    a1, v1 = av1
    a2, v2 = av2
    v1, v2 = np.array(v1), np.array(v2)

    r1 = sstr.from_rotvec(np.radians(a1)*v1/np.linalg.norm(v1))
    r2 = sstr.from_rotvec(np.radians(a2)*v2/np.linalg.norm(v2))

    m = np.dot(r1.as_matrix(), r2.as_matrix())
    r = sstr.from_matrix(m)
    vec = r.as_rotvec()
    phi = np.degrees(np.linalg.norm(vec))

    return phi, vec

def get_ellipse_orbit(a, e, n):
    """计算椭圆轨道,a为长半轴, e为离心率,n为轨道点数"""

    t = np.linspace(0, 2*np.pi, n)
    r = a*(1-e*e)/(1-e*np.cos(t))
    xs = r*np.cos(t)
    ys = r*np.sin(t)
    zs = np.zeros(r.shape)

    return xs, ys, zs

def rotate_s(i):
    """太阳自转函数:沿Z轴旋转,T_SUN_SELF个计数周期旋转一周"""

    return (i%T_SUN_SELF)*360/T_SUN_SELF, (0,0,1)

def rotate_e_orbit(i):
    """地球轨道旋转函数:沿y轴旋转23.43°(黄赤夹角)"""

    return -A_E_ORBIT, (0,1,0)

def rotate_e(i):
    """地球自转函数:沿z轴旋转,T_EARTH_SELF个计数周期旋转一周"""

    return (i%T_EARTH_SELF)*360/T_EARTH_SELF, (0,0,1)

def translate_e(i):
    """地球位移函数:T_EARTH个计数周期循环一次"""

    phi = (i%T_EARTH)*2*np.pi/T_EARTH
    r = A_EARTH*(1-E_EARTH*E_EARTH)/(1-E_EARTH*np.cos(phi))
    d = np.array((r*np.sin(phi), -r*np.cos(phi), 0))

    return ROTATOR_E_ORBIT.apply(d)

def rotate_m_orbit(i):
    """月球轨道旋转函数:沿y轴旋转28.575°后再跟随月球自转"""

    av1 = -A_M_ORBIT, (0,1,0)
    av2 = (i%T_MONTH_SELF)*360/T_MONTH_SELF, (0,0,1)

    return rotate_merge(av1, av2)

def translate_m_orbit(i):
    """月球轨道位移函数:T_EARTH个计数周期循环一次"""

    phi = (i%T_EARTH)*2*np.pi/T_EARTH
    r = A_EARTH*(1-E_EARTH*E_EARTH)/(1-E_EARTH*np.cos(phi))
    d = np.array((r*np.sin(phi), -r*np.cos(phi), 0))

    return ROTATOR_E_ORBIT.apply(d)

def translate_m(i):
    """月球位移函数:T_EARTH个计数周期循环一次"""

    phi = (i%T_EARTH)*2*np.pi/T_EARTH
    r = A_EARTH*(1-E_EARTH*E_EARTH)/(1-E_EARTH*np.cos(phi))
    d1 = np.array((r*np.sin(phi), -r*np.cos(phi), 0))
    d1 = ROTATOR_E_ORBIT.apply(d1)

    phi = (i%T_MONTH)*2*np.pi/T_MONTH
    r = A_MONTH*(1-E_MONTH*E_MONTH)/(1-E_MONTH*np.cos(phi))
    d2 = np.array((r*np.sin(phi), -r*np.cos(phi), 0))
    d2 = ROTATOR_M_ORBIT.apply(d2)

    return d1 + d2

def rotate_m(i):
    """月球自转函数:T_MONTH_SELF个计数周期旋转一周"""

    return (i%T_MONTH_SELF)*360/T_MONTH_SELF, (0,0,1)


# 初始化画布
plt.figure(elevation=5, azimuth=-25)
plt.axis(grid=False)

# 绘制太阳及其自转轴
plt.sphere((0,0,0), R_SUN, texture='res/sun.jpg', rotate=rotate_s, order='R', light=0)
plt.plot((0,0), (0,0), (1.5,-1.5), color='red', style='dash-dot')

# 绘制地球公转轨道
xs_e, ys_e, zs_e = get_ellipse_orbit(A_EARTH, E_EARTH, T_EARTH)
plt.plot(xs_e, ys_e, zs_e, color='cyan', width=1, rotate=rotate_e_orbit, order='R')

# 标注冬夏至和春秋分点
i_summer, i_winter = np.argmin(xs_e), np.argmax(xs_e)
i_autumn, i_spring = np.argmin(ys_e), np.argmax(ys_e)
plt.text('夏至', size=32, pos=ROTATOR_E_ORBIT.apply((xs_e[i_summer], ys_e[i_summer], zs_e[i_summer])))
plt.text('冬至', size=32, pos=ROTATOR_E_ORBIT.apply((xs_e[i_winter], ys_e[i_winter], zs_e[i_winter])))
plt.text('春分', size=32, pos=ROTATOR_E_ORBIT.apply((xs_e[i_spring], ys_e[i_spring], zs_e[i_spring])))
plt.text('秋分', size=32, pos=ROTATOR_E_ORBIT.apply((xs_e[i_autumn], ys_e[i_autumn], zs_e[i_autumn])))

# 绘制地球及其自转轴
plt.sphere((0,0,0), R_EARTH, texture='res/earth.jpg', rotate=rotate_e, translate=translate_e, order='TR', light=0)
plt.plot((0,0), (0,0), (0.5,-0.5), color='cyan', style='dash-dot', rotate=rotate_e, translate=translate_e, order='TR')

# 绘制月球公转轨道
xs_m, ys_m, zs_m = get_ellipse_orbit(A_MONTH, E_MONTH, T_MONTH)
plt.plot(xs_m, ys_m, zs_m, color='magenta', width=1, rotate=rotate_m_orbit, translate=translate_m_orbit, order='TR')

# 绘制月球及其自转轴
plt.sphere((0,0,0), 0.1, texture='res/month.jpg', rotate=rotate_m, translate=translate_m, order='TR', light=0)
plt.plot((0,0), (0,0), (0.3,-0.3), color='magenta', style='dash-dot', rotate=rotate_m, translate=translate_m, order='TR')

plt.show()

代码中用到了太阳、地球和月球的纹理图片,读者可自行下载或替换为自己喜欢的图片。


更多精彩推荐
GitHub 遭抵制!AI 代码生成神器竟成“抄袭工具”?

基于 Python 的 8 种常用抽样方法

太酷了,Python 制作足球可视化图表 | 代码干货

以上是关于程序员炫酷溜娃!用代码画地球日月的动态轨道模型的主要内容,如果未能解决你的问题,请参考以下文章

WebGL实现炫酷地球--数据可视化

请问各位高手,如何用matlab画地球模型,最好带有经纬线的那种!非常感谢!

MATLAB | 老版本也能用,默认设置让简单的代码画出炫酷的图像

MATLAB | 老版本也能用,默认设置让简单的代码画出炫酷的图像

月球轨道相对于太阳

基于jQuery带动画的超炫酷弹出对话框