★ FEATURED ARTICLE 1. 从“抬头看天”到“动手算天”这个项目到底在做什么“地球围绕太阳转”这句话大概每个人在小学阶段就听过。但真正要把这件事从一句常识变成一套能算、能验证、能可视化的东西中间隔着的可不只是几个公式。我做的这个项目核心目标就一个用代码把日地系统的运动关系完整地模拟出来从最粗糙的圆形轨道开始一步步逼近真实情况最后能算出任意时刻地球在轨道上的位置、速度甚至能画出轨道形状、验证开普勒定律。说白了这不是一个天文科普项目而是一个“用计算思维重新理解天体运动”的实操项目。它适合谁呢如果你正在学物理、对天体力学感兴趣、想练手数值计算、或者单纯好奇“科学家到底怎么知道地球现在在哪儿”这个项目都能给你一条从零到一的路径。不需要你有天文背景但需要你会一点编程基础Python或者JavaScript都行我用的是Python因为科学计算生态成熟numpy和matplotlib能省掉大量造轮子的时间。这个项目解决的问题很具体把“地球绕太阳转”这个定性描述变成定量计算。你知道地球绕太阳一圈要365.25天但你知道地球在1月3日离太阳最近、7月4日离太阳最远吗你知道地球公转速度在近日点比远日点快大约3.4%吗这些数字背后是一整套力学和几何的配合。我踩过的坑包括一开始用正圆轨道算结果发现算出来的季节变化对不上后来加了偏心率又发现轨道倾角不能忽略再后来想算精确位置发现还得考虑其他行星的引力扰动。每一步都是“为什么这么算”和“这么算够不够用”的权衡。提示这个项目的核心不是追求绝对精确而是理解“精度是怎么一步步提升的”。你不需要一上来就搞N体问题从二体问题开始把基础打牢后面加东西才不慌。2. 核心思路拆解为什么从圆轨道开始又为什么必须放弃它2.1 圆轨道假设最简模型的价值与局限任何计算项目的第一步都是“先让它跑起来”。我最初写的第一版代码就是把地球轨道当成一个完美的圆半径取1个天文单位约1.496亿公里周期取365.25天然后让地球以恒定角速度绕太阳转。代码不到20行画出来的图也挺好看一个漂亮的圆地球在上面匀速跑。这个模型的价值在于它让你快速建立坐标系概念。太阳在原点地球位置用极坐标表示就是(r, θ)r固定θ随时间线性增加。你立刻就能算出地球在任意时刻的坐标x r·cos(θ)y r·sin(θ)。这个思路简单到可以用计算器手算但它已经能回答一些基本问题比如“春分时地球在哪个方向”。但问题很快就来了。我用这个模型去算地球到太阳的距离发现全年都是1.496亿公里可实际数据告诉我1月初地球离太阳最近约1.471亿公里7月初最远约1.521亿公里。差了整整500万公里这可不是小数目。更关键的是如果轨道是正圆地球公转速度应该恒定但实际观测显示地球在近日点附近跑得快远日点附近跑得慢速度差能达到每秒1公里左右。所以圆轨道模型只能用来“入门”不能用来“做事”。它的局限根源在于真实轨道是椭圆而且太阳不在椭圆中心而是在一个焦点上。这个事实是开普勒第一定律告诉我们的但知道定律和用代码实现它中间还差着参数换算和坐标系变换。2.2 椭圆轨道与开普勒定律从几何到力学的跨越从圆到椭圆核心变化是引入偏心率e。地球轨道的偏心率大约是0.0167这个数字很小意味着地球轨道其实非常接近圆但就是这0.0167的偏差造成了近日点和远日点500万公里的距离差。椭圆的半长轴a就是平均距离半短轴b a·√(1-e²)。太阳位于椭圆的一个焦点上另一个焦点是空的。这里有个容易搞混的地方椭圆有两个焦点太阳在其中一个那另一个焦点有什么用答案是没什么用它只是几何定义的一部分。但地球到太阳的距离r并不是常数而是随角度变化的。用极坐标表示r a(1-e²)/(1e·cos(θ))其中θ是从近日点开始算起的真近点角。这个公式是开普勒轨道方程的基础我第一次推导的时候卡了很久因为θ和时间的对应关系不是线性的。开普勒第二定律说地球和太阳的连线在相等时间内扫过相等面积。这意味着角速度不是常数近日点附近角速度大远日点附近角速度小。用数学表达就是r²·(dθ/dt) 常数这个常数等于轨道角动量除以地球质量。这条定律直接解释了为什么地球在近日点跑得快——因为r小dθ/dt就必须大才能保持r²·(dθ/dt)不变。但真正让计算变得可行的是开普勒第三定律T² ∝ a³。对于太阳系内所有行星轨道周期的平方和半长轴的立方成正比。这个比例常数只和太阳质量有关具体是4π²/(G·M_sun)。地球的a1 AUT1年所以这个常数在“年-AU”单位制下就是1。这意味着你只要知道半长轴就能算出周期反过来也一样。注意开普勒定律描述的是“行星怎么动”但没有解释“为什么这么动”。牛顿后来用万有引力给出了力学解释但在计算层面开普勒定律已经足够精确除非你要考虑其他行星的引力扰动。2.3 数值积分 vs 解析解为什么我最终选了数值方法理论上二体问题有解析解你可以用开普勒方程直接算出任意时刻的位置。开普勒方程是M E - e·sin(E)其中M是平近点角和时间成正比E是偏近点角e是偏心率。这个方程没有初等函数解只能用迭代法求数值解。我试过牛顿迭代法收敛很快一般迭代5到6次就能达到1e-10的精度。但解析解有个前提只有两个天体且没有其他力。真实情况是地球还受到月球、木星、金星等其他天体的引力影响虽然这些力很小但长期累积下来会造成轨道进动。比如地球近日点的进动速率大约是每世纪1.7度其中大部分可以用牛顿力学解释但有一小部分约每世纪0.43角秒需要广义相对论才能解释。所以我的选择是先用解析解开普勒方程把二体问题算准作为基准。然后在这个基础上用数值积分加入其他行星的引力扰动。数值积分我用的是四阶龙格-库塔法RK4步长取1天。为什么是1天因为地球公转周期是365天1天的步长意味着每圈有365个采样点对于轨道形状的还原已经足够。如果步长取太大比如10天能量会漂移轨道会慢慢变形如果取太小比如1小时计算量增加24倍但精度提升有限。实测下来RK4在1天步长下100年内的能量误差小于1e-6轨道半长轴的变化小于1公里。这个精度对于大多数应用已经绰绰有余。如果你要算几千年尺度的轨道变化那就得用更高阶的方法比如辛积分器它能更好地保持能量守恒。3. 核心细节解析参数、坐标系与时间系统3.1 关键参数一览每个数字都有来头做这个项目你首先得把一堆参数搞清楚。我把核心参数整理成表方便你对照。参数名称符号数值单位说明半长轴a1.00000011AU地球轨道平均半径偏心率e0.0167086无量纲轨道扁平程度轨道倾角i0.00005度相对黄道面近日点幅角ω114.20783度近日点方向升交点黄经Ω-11.26064度轨道面与黄道面交线平近点角M358.617度历元时刻的角度轨道周期T365.256363天恒星年引力常数G6.67430e-11m³/(kg·s²)万有引力常数太阳质量M_sun1.98892e30kg太阳质量天文单位AU1.495978707e11m日地平均距离这些参数里a和e是最核心的。a决定了轨道大小e决定了轨道形状。i、Ω、ω决定了轨道在三维空间中的朝向。M决定了地球在轨道上的起始位置。这六个参数合起来叫“轨道六根数”是描述二体问题的标准参数集。我一开始只用了a和e结果发现算出来的位置和实际观测差了几千公里。后来加上i、Ω、ω误差降到几百公里。再后来发现这些参数本身也在随时间变化比如a每年增加约1.5厘米e每年变化约1e-9。对于短期计算这些变化可以忽略对于千年尺度就必须考虑。提示如果你只是做可视化演示a和e就够了。如果你要算具体日期地球的位置至少需要六个参数。如果你要算几百年后的位置还需要考虑参数的时间导数。3.2 坐标系选择黄道坐标系 vs 赤道坐标系坐标系的选择直接决定了你算出来的结果怎么用。我最初用的是黄道坐标系因为地球轨道平面就是黄道面在这个坐标系里地球的z坐标几乎为零计算简单。黄道坐标系的x轴指向春分点z轴垂直于黄道面y轴在黄道面内与x轴垂直。但问题来了如果你要算地球在天空中的位置比如从地面观测站看太阳的位置那就得用赤道坐标系。赤道坐标系的z轴指向北天极x轴也指向春分点但黄道面和赤道面之间有23.44度的夹角叫黄赤交角。这个夹角造成了四季变化也是为什么地球轨道倾角虽然只有0.00005度但季节变化却如此明显。我的做法是在黄道坐标系里算轨道然后把结果转换到赤道坐标系。转换矩阵涉及黄赤交角ε具体是绕x轴旋转ε。这个旋转矩阵是[1, 0, 0] [0, cos(ε), -sin(ε)] [0, sin(ε), cos(ε)]把黄道坐标(x, y, z)乘以这个矩阵就得到赤道坐标。实测下来这个转换在Python里用numpy的dot函数一行就能搞定但要注意角度单位是弧度还是度。3.3 时间系统儒略日与历元转换时间系统是天文计算里最容易出错的地方。我用的是儒略日Julian Day因为它是一个连续的日数从公元前4713年1月1日12时开始算起没有月份和年份的复杂性。比如2024年1月1日0时对应的儒略日是2460310.5。为什么用儒略日因为轨道计算里所有的时间导数都是连续的如果你用日历日期每个月天数不一样闰年还要特殊处理很容易出错。儒略日就是一个浮点数加减乘除都方便。转换公式也不复杂Python的datetime模块可以转但要注意时区。我一般用UTC时间避免夏令时干扰。历元Epoch是另一个关键概念。轨道六根数是在某个特定时刻测定的这个时刻就是历元。我用的是J2000.0即2000年1月1日12时UTC对应的儒略日是2451545.0。所有计算都从这个时刻开始算出来的平近点角M要加上历元时刻的M0再乘以时间差。注意儒略日和简化儒略日MJD差2400000.5天别搞混。我一开始用MJD算结果差了240万天轨道直接飞到太阳系外了。4. 实操过程从零写出一个能跑的轨道计算器4.1 环境准备与依赖安装我用的是Python 3.10核心依赖就三个numpy做数值计算matplotlib做可视化scipy做方程求解。安装命令很简单pip install numpy matplotlib scipy如果你不想装scipy开普勒方程可以用牛顿迭代法自己写代码不超过10行。但scipy的fsolve更省事而且有现成的优化算法。我建议还是装上后面算数值积分也用得上。开发环境我用的是Jupyter Notebook因为可以分块运行边算边看结果。如果你习惯用IDEPyCharm或者VS Code都行但记得把matplotlib的后端设成TkAgg或者Qt5Agg不然图弹不出来。4.2 开普勒方程求解牛顿迭代法的实现细节开普勒方程M E - e·sin(E)是超越方程没有解析解。我用牛顿迭代法求E迭代公式是E_{n1} E_n - (E_n - e·sin(E_n) - M) / (1 - e·cos(E_n))初始值取E_0 M对于e0.1的情况迭代5次就能收敛到1e-12。我写了个函数def solve_kepler(M, e, tol1e-12, max_iter100): E M for i in range(max_iter): f E - e * np.sin(E) - M fp 1 - e * np.cos(E) dE -f / fp E dE if abs(dE) tol: return E raise ValueError(Kepler equation did not converge)这个函数里M是平近点角单位是弧度。如果你用度记得先转弧度。e是偏心率tol是收敛容差max_iter是最大迭代次数。实测下来对于地球轨道e0.0167迭代3次就收敛了非常快。但有个坑当e接近1时比如彗星轨道e0.99牛顿迭代法可能不收敛因为函数在E0附近导数接近零。这时候需要用二分法或者Brent法。地球轨道e很小所以不用担心这个问题。4.3 从偏近点角到真近点角位置计算的完整链条算出E之后还需要转成真近点角θ才能算距离r。转换公式是tan(θ/2) √((1e)/(1-e)) · tan(E/2)这个公式在Eπ附近会有奇点因为tan(E/2)趋于无穷。实际编程时我用的是atan2函数避免除零def eccentric_to_true(E, e): return 2 * np.arctan2(np.sqrt(1e) * np.sin(E/2), np.sqrt(1-e) * np.cos(E/2))然后距离r a(1 - e·cos(E))这个公式比用θ算r更稳定因为不需要处理tan的奇点。最后位置坐标x r · cos(θ) y r · sin(θ)在轨道平面内这就是地球的位置。如果要转到黄道坐标系还需要考虑ω、Ω、i的旋转。旋转顺序是先绕z轴转ω再绕x轴转i再绕z轴转Ω。这个顺序不能乱因为旋转矩阵不满足交换律。我一开始把顺序搞反了结果算出来的轨道面倾斜方向完全不对。后来查了资料才知道标准顺序是3-1-3即绕z、绕x、再绕z。用scipy的Rotation模块可以省掉手写矩阵的麻烦from scipy.spatial.transform import Rotation as R rot R.from_euler(zxz, [ω, i, Ω], degreesTrue) pos_ecliptic rot.apply([x, y, 0])这段代码把轨道平面内的位置转到黄道坐标系非常简洁。实测下来和手写矩阵的结果一致但代码可读性高很多。4.4 数值积分加入扰动RK4的实现与步长选择二体问题算准了接下来加扰动。我用RK4积分地球的运动方程d²r/dt² -G·M_sun·r/|r|³ Σ G·M_planet·(r_planet - r)/|r_planet - r|³第一项是太阳引力第二项是其他行星的引力扰动。我加了木星、金星、火星的引力因为这三颗对地球轨道影响最大。木星质量大金星离得近火星轨道交叉。RK4的步长我取1天因为地球轨道周期365天1天步长意味着每圈365步。实测下来100年后地球位置误差约1000公里相对于1.5亿公里的轨道半径相对误差小于1e-5。这个精度对于大多数应用足够了。如果你要算更长时间比如10000年步长得取小一点或者用自适应步长。我试过用scipy的solve_ivp它自带自适应步长但速度比手写RK4慢因为每次都要估计误差。对于固定步长够用的情况手写RK4更快。提示数值积分最怕能量漂移。你可以每算一步就检查总能量动能势能如果能量变化超过1e-6说明步长太大了。我一般会画能量随时间的变化图如果是一条水平线说明积分稳定。5. 常见问题与排查技巧实录5.1 算出来的轨道是螺旋线而不是闭合椭圆这是新手最容易遇到的问题。原因通常是数值积分误差累积导致能量不守恒。地球每绕一圈半长轴就增加一点几百圈下来轨道就变成螺旋线了。排查思路先检查步长是不是太大。如果步长是10天改成1天试试。如果还不行检查积分方法是不是太低阶。欧拉法一阶精度差RK4四阶一般够用如果你要算几千年用辛积分器比如Verlet算法它能严格保持能量守恒。我踩过的坑一开始用欧拉法步长1天结果100圈后轨道半径增加了1%。换成RK4后1000圈半径变化小于0.001%。所以方法的选择比步长更关键。5.2 开普勒方程不收敛怎么办对于e0.8的情况牛顿迭代法基本都能收敛。如果e0.8比如彗星轨道牛顿法可能震荡。这时候可以用二分法先找到E的区间然后不断二分。二分法收敛慢但绝对稳定。另一个办法是用Danby的初始猜测E_0 M e·sin(M)/(1 - sin(Me) sin(M))。这个猜测比E_0M更接近真解能减少迭代次数。我实测下来对于地球轨道E_0M就够了。对于e0.9的轨道Danby猜测能把迭代次数从20次降到5次。5.3 坐标系转换后位置对不上这个问题通常出在旋转顺序或者角度单位上。先检查角度是度还是弧度。numpy的sin和cos默认用弧度如果你传度进去结果肯定错。我习惯在函数开头统一转弧度ω np.radians(ω_deg) i np.radians(i_deg) Ω np.radians(Ω_deg)然后检查旋转顺序。黄道坐标系到赤道坐标系的转换是绕x轴转ε但如果你先转了别的轴结果就不对。我建议用scipy的Rotation模块它明确指定旋转轴和顺序不容易出错。还有一个隐蔽的坑春分点的定义。黄道坐标系的x轴指向春分点但春分点本身在进动每年移动约50角秒。对于短期计算可以忽略对于千年尺度必须考虑。我用的是J2000.0历元的春分点如果你用当前历元结果会差一点。5.4 常见问题速查表问题现象可能原因排查方法解决方案轨道不闭合数值积分误差检查能量守恒换RK4或辛积分器开普勒方程不收敛e太大或初始值差打印迭代过程用二分法或Danby猜测位置偏差几千公里坐标系转换错误检查旋转顺序用scipy Rotation季节对不上黄赤交角没加检查z坐标加23.44度旋转时间差一天儒略日转换错误检查闰年处理用datetime转儒略日速度方向不对真近点角公式错检查atan2参数用eccentric_to_true函数注意排查问题时先算一个已知时刻的位置比如J2000.0历元的地球位置和NASA JPL的Horizons系统对比。如果差得不多说明大方向对了如果差得离谱肯定是某个参数或公式错了。6. 精度提升与扩展方向从二体到多体从牛顿到爱因斯坦6.1 加入月球地月系统的质心修正地球其实不是绕太阳转而是地月质心绕太阳转。地月质心离地心约4670公里在地球内部。这个修正对地球位置的影响大约是每秒0.01公里一年累积下来约300公里。对于高精度计算这个修正不能忽略。我的做法是把地球和月球当成一个双星系统先算地月质心的轨道然后在地月质心轨道上叠加月球绕质心的运动。月球轨道也是椭圆偏心率0.0549半长轴384400公里。这个计算量不大但能让位置精度提升一个量级。6.2 广义相对论修正水星近日点进动的启示牛顿力学算水星轨道近日点进动每世纪比观测少43角秒。这个差异用广义相对论才能解释。地球轨道也有类似效应但小得多大约每世纪0.0001角秒。对于当前精度要求这个修正可以忽略。但如果你要算几万年后的轨道相对论效应会累积到可观测的程度。广义相对论修正项是a_rel -G·M_sun·r/(c²·r³) · (4·G·M_sun/r - v²) · r。这个项在近日点附近最大因为r最小。我试过加上这个项计算量增加不多但代码复杂度上升。对于地球轨道加不加这个项100年内位置差小于1米。所以除非你做高精度天体测量否则没必要。6.3 可视化从静态图到动态模拟算出来的数据不画出来等于白算。我用matplotlib画了三种图轨道形状图、距离随时间变化图、速度随时间变化图。轨道形状图就是x-y平面上的椭圆太阳在焦点上。距离图能清楚看到近日点和远日点。速度图能验证开普勒第二定律。如果你想做动态模拟可以用matplotlib的animation模块每帧更新地球位置。帧率设成30fps一年365天每帧走12天30秒能模拟一年。这个速度肉眼能看清地球在轨道上的运动也能看到速度变化。提示画图时记得把坐标轴比例设成相等不然椭圆会被拉成圆。用ax.set_aspect(equal)就行。6.4 扩展方向从地球到其他行星这套代码稍加修改就能算其他行星。水星e0.2056金星e0.0068火星e0.0934。你只需要改轨道六根数其他逻辑不变。我试过算水星因为e大开普勒方程迭代次数多但牛顿法依然收敛。算出来的水星轨道近日点进动和观测值对比牛顿力学差43角秒这个差异正好是广义相对论修正的量级。如果你把八大行星都算出来就能画一个太阳系俯视图。注意行星轨道倾角不同水星倾角7度金星3.4度地球0度。在三维空间里画出来能看到轨道面并不重合。这个可视化效果很震撼也能帮你理解为什么行星在天空中的轨迹是复杂的。我个人在实际操作中的体会是这个项目最大的价值不是算得多准而是让你亲手把“地球绕太阳转”这句话拆解成参数、公式、代码和图表。每解决一个bug你对天体运动的理解就深一层。踩过几次坑之后你会发现那些看似高深的天文计算底层逻辑其实很朴素——就是牛顿第二定律加万有引力再加一点几何变换。
阅读完成 · 觉得有帮助？
