文章目录
- 1 坐标系的选择
- 1.1 黄道坐标系
- 1.2 三维空间直角坐标系
- 2 使用JPL星历表计算轨道
- 2.1 日期时间
- 2.2 特定时刻天体的位置
- 2.3 天体运行轨道
- 3 太阳系模型
- 3. 1 全家福
- 3.2 时间、距离和半径的缩放
- 3.3 黄道坐标系模型
天何所沓？十二焉分？日月安属？列星安陈？—— 屈原
远古时期的人类就对神秘幽邃的星空充满了好奇与敬畏。仰望星空，宇宙浩瀚。比宇宙更壮阔的，是人类对宇宙的不懈追问和沉淀在基因中的探索精神。
本文尝试用WxGL来回答“日月安属、列星安陈”这个古老的问题。太阳系天体轨道数据来源于JPL（美国喷气实验室）星历表，天体自转周期和自转轴倾角来源于网路。本模型天体公转轨道和自转轴倾角的准确性可以和以下两个网站相互印证。
- Heavens Above：太阳系
- NASA中文网：太阳系行星的倾角和自转
1 坐标系的选择
仰望星空，有一个非常重要却最容易被忽略的前提：站在哪里仰望？古人自然而然地选择了站在地球上仰望星空，于是诞生了地心说和天球模型。本模型选择站在银河系的某一点上，和太阳保持相对固定的位置关系，以上帝的视角来观察整个太阳系。
不管选择站在哪里，都需要一个坐标系——通常称为天体坐标系，用以标定和描述太阳系中各个天体的相对关系。常用的天体坐标系有很多种，本模型使用黄道坐标系描述天体运行轨迹，最终转化为WxGL使用的三维空间直角坐标系绘制输出。
1.1 黄道坐标系
所谓黄道，就是太阳在天球上的运动轨迹，黄道所在的平面称为黄道面。换成上帝视角，黄道面就是地球公转的轨道面。正如倾斜的地球仪呈现出的姿态那样，地球的赤道面和黄道面并不重合，二者存在23°26’的夹角，这就是黄赤交角。
黄道坐标系（Ecliptic Coordinate System）是以黄道面为基本面并且以春分点为原点组成的天球坐标系，其坐标不会随着时间或者观察者的地理位置而变化。
将X设为一个在天球中的星体的坐标，并且 X’ 为X在黄道面上的投影。那么X在黄道坐标系中的坐标就可以用一下2个参数来定义：
- 黄经（Ecliptic Longitude），简称λ；也就是角VX’，从春分点开始向逆时针方向测量，取值范围为 0° ~ 360°。
- 黄纬（Ecliptic Latitude），简称β；也就是角 XX’，即天体和太阳的连线与黄道面的夹角（线面角） ，取值范围为+90° ~ -90°。
1.2 三维空间直角坐标系
和OpenGL一样，WxGL使用右手坐标系，默认y轴为高度轴。使用haxis参数可以设置z轴为高度轴，使用azim参数和elev参数可以改变初始的方位角和仰角，使用fovy参数可以设置初始的水平视野宽度。本模型以z轴为高度轴。
import wxgl
app = wxgl.App(haxis='z', azim=15, elev=10, fovy=50)
app.axes()
app.show()
下图左右分别以y轴和z轴为高度轴建立坐标系，初始方位角都是15°，初始仰角都是10°。
2 使用JPL星历表计算轨道
JPL星历表给出了太阳、月球、八大行星以及冥王星过去和将来一段时间内的位置信息，并且是开放可使用的。JPL星历表在20世纪60年代由喷气推进实验室建立，最初用作行星探测导航的目的，随着观测技术的不断提高，新的观测数据不断获得，JPL星历表仍在不断修正和完善。
本文使用DE405星历表，下载地址为https://github.com/skyfielders/python-skyfield/blob/master/ci/de405.bsp，该文件约65M，涵盖了公元1600年～公元2200年的太阳系数据，可满足一般适用要求。
Python有很多模块都可以读取星历表数据，我喜欢使用skyfield，安装也很简单。
pip install skyfield
2.1 日期时间
JPL星历表使用儒略日（儒略纪元），因此一个在儒略日和UTC时刻之间转换的工具是必要的，skyfield模块提供了这个工具。
>>> from datetime import datetime, timedelta
>>> from skyfield.api import load, utc
>>> ts = load.timescale() # 创建一个时间处理工具
>>> ts.now() # 当前时刻的儒略历对象
<Time tt=2460055.7605328355>
>>> ts.utc(2023, 4, 21, 8, 0, 0) # 指定日期时间生成儒略历对象
<Time tt=2460055.834134074>
>>> dt = datetime.fromisoformat('2023-04-21 08:00:00') # 生成一个datetime对象
>>> t = ts.from_datetime(dt.replace(tzinfo=utc)) # 从datetime对象转化为儒略日对象
>>> t.utc_iso() # 转为UTC字符串
'2023-04-21T08:00:00Z'
>>> ts.utc(2023, 4, 21, range(8)) # 生成一个长度为8时间序列，间隔1小时
<Time tt=[2460055.500800741 ... 2460055.7924674074] len=8>
2.2 特定时刻天体的位置
以地球为例，计算2023年5月1日零时（UTC时刻）的位置，代码如下。
from skyfield.api import load
ts = load.timescale()
t = ts.utc(2023, 5, 1, 0, 0 ,0)
planets = load('res/de405.bsp') # 读星历表
earth = planets['earth'] # 地球对象
print(earth.at(t).ecliptic_xyz().au) # t时刻地球的位置，使用天文单位（日地距离）
print(earth.at(t).ecliptic_xyz().km) # t时刻地球的位置，使用天文单位（日地距离）
运行结果如下。
xufive@xuxiangwudeMacBook-Pro solar % python3 orbit_demo.py
[-7.79775078e-01 -6.49301881e-01 2.49865838e-04]
[-1.16652691e+08 -9.71341788e+07 3.73793973e+04]
2.3 天体运行轨道
下面的代码给出了地球从2023年5月1日零时（UTC时刻）开始的1小时内的运行轨迹，时间间隔1秒。
from skyfield.api import load
ts = load.timescale()
t = ts.utc(2023, 5, 1, 0, 0, range(3600))
planets = load('res/de405.bsp') # 读星历表
earth = planets['earth'] # 地球对象
orbit = earth.at(t).ecliptic_xyz().km
print(orbit.shape)
计算结果orbit是一个ndarray类型的二维数组。
xufive@xuxiangwudeMacBook-Pro solar % python3 orbit_demo.py
(3, 3600)
3 太阳系模型
3. 1 全家福
设想的太阳系模型包括太阳、月球和八大行星，其中最大的天体自然非太阳莫属，半径约69.6万千亩，最小的是月球，半径只有1738千米，大约是太阳的1/400。然而，这样的差异相比天体间距离的差异只能算是小巫见大巫：海王星距离太阳约45亿千米，几乎是地球和月球距离（38万千米）的1.2万倍！
这样大的动态范围，使得按照等比例绘制出来的太阳系模型基本没有价值，因此大多数的模型都会对天体半径和距离做缩放处理。本模型也不例外，但缩放比例可以自由选择，如果选择1，就是一个等比例的模型。
为了建立直观的概念，先来一张八大行星排队站在太阳面前的合影，看看它们的身形在视觉上差距有多大。前排从左向右分别是水星（Mercury）、金星（Venus）、地球（Earth）、火星（Mars）、木星（Jupiter）、土星（Saturn）、天王星（Uranus）和海王星（Neptune），最小的月球没有资格拍这张合影。
生成这张图片的代码很简单，只有区区13行。
import wxgl
app = wxgl.App(fovy=18)
app.sphere((0,0,0), 696300, texture='res/sun.jpg', transform=lambda t:((0,1,0,(0.002*t)%360),))
app.sphere((-535304,0,594516), 2440, texture='res/mercury.jpg')
app.sphere((-400000,0,692820), 6052, texture='res/venus.jpg')
app.sphere((-247214,0,760845), 6371, texture='res/earth.jpg')
app.sphere((-83623,0,795618), 3398, texture='res/mars.jpg')
app.sphere((83623,0,795618), 71492, texture='res/jupiter.jpg')
app.sphere((247214,0,760845), 60268, texture='res/saturn.jpg')
app.sphere((400000,0,692820), 25559, texture='res/uranus.jpg')
app.sphere((535304,0,594516), 24718, texture='res/neptune.jpg')
app.show()
虽然简单，却体现的WxGL的一个特定：使用函数作为参数。绘制出来的太阳是自转的，因为以参数形式传递了一个模型变换函数给圆球绘制函数sphere，该变换函数以时间t为参数，返回t时刻太阳绕(0,1,0)轴旋转的角度。
运行这段代码需要用到各个天体的纹理图片，下面给出的这些素材可以直接下载到本地使用，也可以从https://github.com/xufive/wxgl下载。
3.2 时间、距离和半径的缩放
地球绕太阳公转一周需要一年的时间，而模型则要快速实现地球公转，因此需要设置一个时间加速因子。如果用模型的1秒表示实际时间的1天，这个加速因子就是86400，这样模型中的地球大约6分钟就可以绕太阳公转一周。不过这会让模型中的地球1秒钟自转一周，快到无法看清了。因此，时间因子选择2800（模型的1秒相当于实际时间的8小时）是一个折衷的方案。
海王星等行星距离太阳太过遥远，若按照实际比例绘制模型的话，恐怕连太阳都小到不可见，因此等比例缩减行星距离太阳的距离是非常必要的。这个缩放因子选择1/50可以保证火星以内的行星不会相互挤压。
为了能够观察水星、月球等较小的天体，需要将除太阳外的天体半径放大，但也不宜让木星、土星等行星的大小超过太阳，因此半径的缩放因子选择20较为适宜。
放大天体半径、缩小天体间距离可能会导致地球和月球部分重叠，甚至月球被地球完全吞。在确定了距离和半径的缩放因子后，需要评估地球和月球是否重叠，如有，则需要增加地球和月球之间的距离。
3.3 黄道坐标系模型
太阳系模型代码大约180行。完整代码及图片、星历表等资源文件，可从https://github.com/xufive/wxgl下载。特别提醒：这段代码用到了WxGL最新版新增的功能，要运行代码的话，请务必将WxGL更新到0.9.12版本。
import wxgl
import numpy as np
from skyfield.api import load, utc
from datetime import datetime, timedelta
class SolarSystemModel:
"""太阳系天体轨道计算类"""
# 天体常量：半径r（km）、公转周期revo（太阳日）、自转周期spin（小时）和自转轴倾角tilt（度，相对于黄道面）
CONST = {
'sun': {'r': 696300, 'revo': 0, 'spin': 24*24.47, 'tilt': 7},
'mercury': {'r': 2440, 'revo': 87.99, 'spin': 24*58.6, 'tilt': 0},
'venus': {'r': 6052, 'revo': 224.70, 'spin': 24*243, 'tilt': 177.3},
'earth': {'r': 6371, 'revo': 365.2564, 'spin': 23.934, 'tilt': 23.43},
'mars': {'r': 3398, 'revo': 686.97, 'spin': 24.617, 'tilt': 25.2},
'jupiter': {'r': 71492, 'revo': 4332.71, 'spin': 9.833, 'tilt': 3.1},
'saturn': {'r': 60268, 'revo': 10759.5, 'spin': 10.233, 'tilt': 26.7},
'uranus': {'r': 25559, 'revo': 30685, 'spin': 17.24, 'tilt': 97.8},
'neptune': {'r': 24718, 'revo': 60194.25, 'spin': 15.966, 'tilt': 28.3},
'moon': {'r': 1738, 'revo': 27.32, 'spin': 24*27.32, 'tilt': 1.5424}
}
def __init__(self, de_file, t_factor=28800, d_factor=1/50, r_factor=20, start_dt=None):
"""构造函数
de_file - JPL星历表文件
t_factor - 时间加速因子，默认模型中的1秒对应实际时间的28800秒（8小时）
d_factor - 距离缩放因子，默认以实际天体间距离的1/50绘制模型
r_factor - 除太阳外其他天体半径缩放因子，默认以实际半径的20倍绘制模型
start_dt - 开始日期时间字符串（YYYY-MM-DD hh:mm:ss），默认None，表示当前UTC时刻开始
"""
self.t_factor = t_factor
self.d_factor = d_factor
self.r_factor = r_factor
self.start_dt = datetime.utcnow().replace(tzinfo=utc) if start_dt is None else datetime.fromisoformat(start_dt).replace(tzinfo=utc)
self.ts = load.timescale() # 创建处理时间的对象
self.planets = load(de_file) # 加载星历文件
self.k = 20000 * r_factor / (380000 * d_factor) # 地月距离缩放系数
def get_ecliptic_xyz_at_dt(self, planet_name, dt):
"""根据日期时间计算天体在黄道坐标系中的坐标
planet_name - 天体名称
dt - datetime类型的日期时间
"""
t = self.ts.from_datetime(dt)
name = '%s_barycenter'%planet_name if planet_name in ('jupiter', 'saturn', 'uranus', 'neptune') else planet_name
x, y, z = self.planets[name].at(t).ecliptic_xyz().km
return x*self.d_factor, y*self.d_factor, z*self.d_factor
def get_ecliptic_xyz(self, planet_name, time_delta):
"""计算天体在黄道坐标系中的坐标
planet_name - 天体名称
time_delta - 距离开始时刻的时间偏移量，以为毫秒单位
"""
seconds = self.t_factor * time_delta / 1000 # 模型时间转换为实际时间偏移量
dt = self.start_dt + timedelta(seconds=seconds)
return self.get_ecliptic_xyz_at_dt(planet_name, dt)
def get_ecliptic_orbit(self, planet_name):
"""计算天体单个公转周期的运行轨道，planet_name为天体名称"""
days = round(self.CONST[planet_name]['revo'] - 366)
dt = self.start_dt - timedelta(days=days) if days > 0 else self.start_dt
hours = np.linspace(0, self.CONST[planet_name]['revo'], 1001) * 24
t = self.ts.utc(dt.year, dt.month, dt.day, hours)
name = '%s_barycenter'%planet_name if planet_name in ('jupiter', 'saturn', 'uranus', 'neptune') else planet_name
x, y, z = self.planets[name].at(t).ecliptic_xyz().km
return np.stack((x*self.d_factor, y*self.d_factor, z*self.d_factor), axis=1)
def get_sphere(self, planet_name):
"""返回天体球面网格的顶点坐标，planet_name为天体名称"""
r = self.CONST[planet_name]['r'] if planet_name == 'sun' else self.CONST[planet_name]['r'] * self.r_factor
gv, gu = np.mgrid[np.pi/2:-np.pi/2:37j, 0:2*np.pi:73j]
zs = r * np.sin(gv)
xs = r * np.cos(gv) * np.cos(gu)
ys = r * np.cos(gv) * np.sin(gu)
return xs, ys, zs
def get_axis(self, planet_name):
"""返回天体旋转轴的顶点坐标，planet_name为天体名称"""
r = self.CONST[planet_name]['r'] * self.r_factor
return [[0, 0, 1.5*r], [0, 0, -1.5*r]]
def dt_func(self, t):
"""格式化日期时间的函数，用于在UI的状态栏显示模型对应的UTC时间"""
return self.ts.from_datetime(self.start_dt + timedelta(seconds=self.t_factor*t/1000)).utc_iso()
def tf_sun(self, t):
"""太阳模型变换函数，实现自转"""
speed = 0.36 / (self.CONST['sun']['spin'] * 3600 / self.t_factor)
phi = (t * speed) % 360
return ((0, 0, 1, phi), )
def tf_moon(self, t):
"""月球模型变换函数，跟随地球运动的同时实现自转、旋转轴倾斜和公转"""
rotate = (0, 0, 1, (t * 0.36 / (self.CONST['moon']['spin'] * 3600 / self.t_factor)) % 360)
tile = (-1, 0, 0, self.CONST['moon']['tilt'])
xm, ym, zm = self.get_ecliptic_xyz('moon', t)
xe, ye, ze = self.get_ecliptic_xyz('earth', t)
shift = xe+(xm-xe)*self.k, ye+(ym-ye)*self.k, ze+(zm-ze)*self.k
return (rotate, tile, shift)
def tf_factory(self, planet_name):
"""天体模型变换函数生成器，返回实现天体自转、旋转轴倾斜、公转的变换函数"""
def tf(t):
rotate = (0, 0, 1, (t * 0.36 / (self.CONST[planet_name]['spin'] * 3600 / self.t_factor)) % 360)
tile = (-1, 0, 0, self.CONST[planet_name]['tilt'])
shift = self.get_ecliptic_xyz(planet_name, t)
return (rotate, tile, shift)
return tf
def show_ecs(self):
"""绘制黄道坐标系模型"""
app = wxgl.App(haxis='z', elev=15, fovy=35, backend='qt')
app.title('太阳系模型')
app.info(time_func=self.dt_func) # 在状态栏中显示日期时间信息
# 绘制太阳
xs, ys, zs = self.get_sphere('sun')
app.mesh(xs, ys, zs, texture='res/sun.jpg', light=wxgl.BaseLight(), transform=self.tf_sun)
# 绘制月球
xs, ys, zs = self.get_sphere('moon')
app.mesh(xs, ys, zs, texture='res/moon.jpg', light=wxgl.BaseLight(), transform=self.tf_moon)
# 绘制8个行星
names = ('mercury', 'venus', 'earth', 'mars', 'jupiter', 'saturn', 'uranus', 'neptune')
colors = ('dodgerblue', 'gold', 'cyan', 'firebrick', 'burlywood', 'darksalmon', 'lightgray', 'lightskyblue')
for name, color in zip(names, colors):
# 绘制行星
xs, ys, zs = self.get_sphere(name)
app.mesh(xs, ys, zs, texture='res/%s.jpg'%name, light=wxgl.BaseLight(), transform=self.tf_factory(name))
# 绘制行星自转轴
vs = self.get_axis(name)
app.line(vs, color=color, stipple='dash-dot', transform=self.tf_factory(name))
# 绘制行星公转轨道线
orbit = self.get_ecliptic_orbit(name)
app.line(orbit, color=color)
# 绘制春分、秋分、夏至和冬至标识
for dt_str, word in (('03-21','春分'), ('06-22','夏至'), ('09-22','秋分'), ('12-23','冬至')):
dt = datetime.fromisoformat('%d-%s'%(self.start_dt.year, dt_str)).replace(tzinfo=utc)
x, y, z = self.get_ecliptic_xyz_at_dt('earth', dt)
d = 4000 * self.r_factor
box = [[x-2*d, y, z-d], [x-2*d, y, z-2*d], [x+2*d, y, z-2*d], [x+2*d, y, z-d]]
app.line([[x, y, z+d], [x, y, z-d]], color='white', width=2)
app.text3d(word, box, align='center')
app.show()
if __name__ == '__main__':
ssm = SolarSystemModel('res/de405.bsp', d_factor=1/50, r_factor=20, start_dt='1962-02-05 01:00:00') # 该日期七星连珠，排列在9.3°范围内
ssm.show_ecs()
绘制模型时，除了星历表文件参数不可省略，其他参数都是可选的。代码中给出的日期，是历史上有名的七星连珠日。省略日期，则从当前时刻开始绘制模型。
这是太阳系的全景图，火星以内的行星和太阳位于中间，只有外围的木星、土星、天王星和海王星轨道可以看清楚。
滚动鼠标滚轮，终于看到内行星了。
继续滚动滚轮，月球也显现出来了。
继续放大，终于看到了太阳-地球-月球系统的运动。这算是“日月安属、列星安陈”这个问题的答案吧。
