1. 项目缘起与整体设计思路

1.1 为什么选择用Python模拟太阳系

第一次看到"150行代码模拟太阳系"这个目标时，我的直觉是：这要么是个玩具，要么是个陷阱。玩具是因为网上太多"几行代码画个圆就叫太阳系"的示例，陷阱是因为一旦涉及真实轨道力学，代码量很容易失控。但实际动手之后我发现，150行是一个相当精妙的目标——它逼着你砍掉所有不必要的抽象，只保留最核心的物理逻辑和渲染逻辑。

这个项目的本质是用数值积分方法求解N体问题的一个简化版本。严格意义上的N体问题没有解析解，但太阳系有一个天然优势：太阳的质量占整个系统质量的99.8%以上，所以我们可以把太阳当作固定不动的中心天体，其他行星只受太阳引力影响，忽略行星之间的相互作用。这就是所谓的"受中心力场约束的二体问题近似"，在这个近似下，每颗行星的运动方程可以独立求解，计算量从O(N²)降到O(N)。

选择Python而不是C++或JavaScript，理由很实际：numpy负责向量化计算，matplotlib负责动画渲染，pygame负责音频播放，三件套加起来让代码量压缩到极致。如果换成纯C++写OpenGL，光是窗口初始化和着色器编译就不止150行。Python在这个场景下的定位是"快速验证物理直觉的工具"，而不是"生产级仿真软件"。

1.2 150行代码的分配策略

我给自己定的分配方案是这样的：物理计算约40行，渲染与动画约50行，音频处理约20行，参数配置与主循环约40行。这个分配比例不是拍脑袋来的，而是根据"改动频率"决定的——物理公式一旦写对基本不用动，渲染代码需要反复调参，音频是锦上添花的部分，参数配置则是你调试时改得最多的区域。

提示：不要一上来就追求150行。先写一个能跑通的300行版本，确认物理正确、动画流畅之后，再逐步合并重复逻辑、删除调试代码，自然就能压到150行左右。强行压缩代码只会让你在调试时痛不欲生。

1.3 最终效果预期

跑起来之后你会看到：黑色背景上，太阳固定在中心，八大行星按照各自的轨道半径和公转周期运动，轨道用淡色椭圆线标出，行星大小按对数比例缩放（否则木星和地球的尺寸差异会让地球变成一个像素点）。背景播放一段循环的环境音乐，音量控制在30%左右，不抢视觉焦点。

这个效果适合三类人：想入门科学计算可视化的初学者、想给学生演示开普勒定律的教师、以及想找一个周末小项目练手的开发者。代码不依赖任何付费库，一台普通笔记本就能流畅运行。

2. 核心物理原理与参数计算

2.1 万有引力定律的数值实现

牛顿的万有引力公式是 F = G·M·m/r²，方向沿两物体连线。在我们的简化模型中，太阳固定，行星只受太阳引力，所以加速度 a = G·M/r²，方向指向太阳。这里M是太阳质量，r是行星到太阳的距离。

关键问题是：G和M的数值怎么选？真实值G=6.674×10⁻¹¹ N·m²/kg²，M=1.989×10³⁰ kg，两者相乘GM≈1.327×10²⁰ m³/s²。如果直接用国际单位制，行星轨道半径是10¹¹米量级，周期是10⁷秒量级，数值跨度太大，浮点精度容易出问题。

我的做法是归一化单位：令GM=4π²，距离单位取天文单位（AU），时间单位取年。这样地球的轨道半径r=1，周期T=1年，代入a=4π²/r²正好满足圆轨道条件a=v²/r=4π²/r，自洽。这个技巧在数值模拟里叫"无量纲化"，能显著提升计算稳定性。

GM = 4 * np.pi ** 2 # 归一化后的引力常数乘太阳质量

2.2 初始速度的设定逻辑

行星要绕太阳做圆周运动，需要满足向心力等于引力：v²/r = GM/r²，所以v = sqrt(GM/r)。对于椭圆轨道，在近日点速度略大于这个值，远日点略小。为了简化，我全部按圆轨道处理，初始位置放在x轴上，初始速度沿y轴方向。

planets = [ {"name": "水星", "r": 0.387, "color": "#b0b0b0", "size": 2}, {"name": "金星", "r": 0.723, "color": "#e8cda2", "size": 3}, {"name": "地球", "r": 1.000, "color": "#4a90d9", "size": 3}, {"name": "火星", "r": 1.524, "color": "#c1440e", "size": 2.5}, {"name": "木星", "r": 5.203, "color": "#d8ca9d", "size": 6}, {"name": "土星", "r": 9.537, "color": "#ead6b8", "size": 5.5}, {"name": "天王星", "r": 19.19, "color": "#b5e3e3", "size": 4}, {"name": "海王星", "r": 30.07, "color": "#3f54ba", "size": 4}, ]

注意海王星的轨道半径是30AU，如果按线性比例画图，水星到太阳的距离只有海王星的1.3%，几乎贴在太阳上。所以渲染时必须用对数坐标或者分段缩放，否则内行星全挤成一团。

2.3 数值积分方法的选择

最简单的积分方法是欧拉法：x_new = x + v·dt，v_new = v + a·dt。但欧拉法能量不守恒，跑几百步之后行星就会螺旋飞出去。我试过，地球大约在模拟30年之后轨道半径扩大了15%，完全不能用。

改进方案有两个：速度Verlet法和四阶Runge-Kutta法。Verlet法的优势是辛结构，能量长期守恒，而且每步只需要计算一次加速度，代码量极小：

def verlet_step(pos, vel, dt): acc = -GM * pos / np.linalg.norm(pos) ** 3 vel_half = vel + 0.5 * acc * dt pos_new = pos + vel_half * dt acc_new = -GM * pos_new / np.linalg.norm(pos_new) ** 3 vel_new = vel_half + 0.5 * acc_new * dt return pos_new, vel_new

RK4精度更高但每步要算四次加速度，代码量翻倍。对于这个可视化项目，Verlet法完全够用，而且它的能量守恒特性让长时间模拟不会出现轨道漂移。

2.4 时间步长的确定

dt取多大？这取决于你想让动画跑多快。如果dt=0.01年，地球公转一圈需要100步，按30fps渲染就是3.3秒一圈，视觉上比较舒服。但水星周期只有0.24年，24步就转完一圈，会显得飞快。这是真实物理的体现——水星确实比地球快得多。

注意：dt不能太大，否则Verlet法也会失稳。经验值是dt < 0.05×最短周期。水星周期0.24年，所以dt上限约0.012年。我取0.005年，留足安全余量。

3. 渲染与动画实现细节

3.1 matplotlib动画框架的搭建

matplotlib的FuncAnimation是标准方案，核心是定义一个update函数，每帧调用一次，返回需要重绘的艺术家对象。这里有个性能陷阱：如果每帧都调用ax.clear()再重新画所有元素，帧率会掉到个位数。正确做法是只更新已有对象的数据。

fig, ax = plt.subplots(figsize=(8, 8), facecolor="black") ax.set_facecolor("black") ax.set_xlim(-32, 32) ax.set_ylim(-32, 32) ax.set_aspect("equal") ax.axis("off") scatters = [] for p in planets: sc, = ax.plot([], [], "o", color=p["color"], markersize=p["size"]) scatters.append(sc)

每帧只需要调用sc.set_data([x], [y])，而不是重新创建对象。这个优化能让帧率从5fps提升到30fps以上。

3.2 轨道拖尾效果的实现

纯点状行星看起来太单调，加一条轨道拖尾会好看很多。实现方式是维护一个固定长度的历史位置队列，每帧把当前位置append进去，超出长度的pop掉。

from collections import deque trails = [deque(maxlen=200) for _ in planets]

200帧的拖尾大约对应6-7秒的动画，刚好能看出轨道弧度但不会糊成一团。拖尾用ax.plot画线，透明度设0.3，颜色和行星一致。

3.3 对数坐标的视觉处理

前面提到海王星30AU的问题。我的处理方式是：物理计算用真实距离，渲染时对坐标做压缩。具体做法是对位置向量取平方根再乘以符号：

def compress(pos): return np.sign(pos) * np.sqrt(np.abs(pos))

这样1AU变成1，30AU变成5.48，视觉比例从30:1压缩到5.5:1，内行星和外行星都能看清。缺点是轨道形状会略微变形，但对于科普演示来说完全可以接受。

3.4 太阳的渲染与光晕效果

太阳用一个大一点的黄色圆点表示，再加几层半透明的同心圆模拟光晕。matplotlib没有内置的glow效果，但可以用多个Circle对象叠加，透明度递减：

for radius, alpha in [(1.5, 0.3), (2.5, 0.15), (4.0, 0.05)]: glow = plt.Circle((0, 0), radius, color="yellow", alpha=alpha) ax.add_patch(glow)

这个技巧在暗色背景上效果很好，成本几乎为零。

4. 背景音乐的集成方案

4.1 音频库的选型对比

Python播放音频有几个选择：pygame.mixer、pydub+simpleaudio、playsound。我最终选pygame.mixer，理由是它支持循环播放、音量控制、异步播放，而且不需要额外安装ffmpeg之类的系统依赖。

4.2 音频加载与循环播放

import pygame pygame.mixer.init() pygame.mixer.music.load("space_ambient.mp3") pygame.mixer.music.set_volume(0.3) pygame.mixer.music.play(-1) # -1表示无限循环

play(-1)的参数是循环次数，-1代表无限循环。音量设0.3是因为环境音乐只需要烘托氛围，太大声会干扰思考。

4.3 音频与动画的同步问题

matplotlib的动画是在主线程运行的，pygame的音频在后台线程播放，两者天然异步。这其实是好事——你不需要精确同步，音乐只是背景氛围。但有一个坑：程序退出时音频不会自动停止，需要注册退出回调：

import atexit atexit.register(pygame.mixer.music.stop)

否则关掉窗口后音乐还在响，得手动去任务管理器杀进程。

提示：如果你没有合适的背景音乐文件，可以用numpy生成一段简单的正弦波环境音，频率在200-400Hz之间缓慢变化，叠加一点白噪声。这样就不依赖外部文件了，代价是多20行代码。

5. 完整代码结构与关键片段

5.1 整体代码骨架

整个程序分为五个部分：导入与配置、物理初始化、渲染初始化、动画更新函数、主循环启动。下面按顺序拆解关键片段。

import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation from collections import deque import pygame, atexit # ========== 配置区 ========== GM = 4 * np.pi ** 2 DT = 0.005 TRAIL_LEN = 200

5.2 物理状态初始化

每颗行星的状态用位置向量和速度向量表示，初始位置在x轴上，速度沿y轴：

states = [] for p in planets: r = p["r"] pos = np.array([r, 0.0]) vel = np.array([0.0, np.sqrt(GM / r)]) states.append({"pos": pos, "vel": vel, "trail": deque(maxlen=TRAIL_LEN)})

5.3 动画更新函数

这是整个程序的核心，每帧做三件事：推进物理、更新散点位置、更新拖尾：

def update(frame): for i, s in enumerate(states): pos, vel = verlet_step(s["pos"], s["vel"], DT) s["pos"], s["vel"] = pos, vel xy = compress(pos) scatters[i].set_data([xy[0]], [xy[1]]) s["trail"].append(xy) trail_lines[i].set_data( [p[0] for p in s["trail"]], [p[1] for p in s["trail"]] ) return scatters + trail_lines

5.4 启动与退出处理

ani = FuncAnimation(fig, update, interval=33, blit=True) atexit.register(pygame.mixer.music.stop) plt.show()

interval=33对应约30fps，blit=True只重绘变化的部分，能显著提升性能。

6. 常见问题与排查技巧实录

6.1 行星飞出去或者螺旋坠落

这是最常见的问题，根本原因是积分方法不稳定或者dt太大。排查顺序：先确认用的是Verlet而不是欧拉，再把dt减半试试。如果减半后正常，说明是步长问题；如果还是飞，检查加速度公式里的norm(pos)**3是不是写成了norm(pos)**2——这是最经典的错误，分母少了一次方，力的方向对了但大小错了。

6.2 动画卡顿、帧率低

三个可能原因：一是用了ax.clear()，改成只更新数据；二是拖尾太长，把maxlen从200降到100；三是blit=False，改成True。如果还卡，把行星数量从8颗减到4颗测试，确认是不是计算量问题。

6.3 音乐不播放或报错

pygame.mixer.init()在某些系统上需要指定采样率，加参数pygame.mixer.init(frequency=44100, size=-16, channels=2)试试。如果报"file not found"，检查音频文件路径是不是相对路径——matplotlib的动画可能改变工作目录，最好用绝对路径。

6.4 内行星挤在一起看不清

这是尺度问题，不是bug。解决方案前面提过：对数压缩坐标。如果不想改物理坐标，也可以给每颗行星单独设置缩放系数，但那样轨道就不是真实比例了，科普价值会打折扣。

6.5 一个容易被忽略的坑：浮点累积误差

即使Verlet法能量守恒，跑上几万步之后位置精度还是会下降。如果你打算让程序连续跑几个小时，建议每隔一万步把行星位置重新归一化到初始轨道半径上。这不是物理上正确的做法，但作为可视化演示，它能保证画面永远好看。

if frame % 10000 == 0: for s, p in zip(states, planets): s["pos"] = s["pos"] / np.linalg.norm(s["pos"]) * p["r"]

这个技巧是我在实际长时间运行中摸索出来的，常规教程里不会提，但对于做展示类项目非常实用。

7. 可扩展方向与个人实践体会

7.1 加入月球或卫星系统

在现有框架下加一颗绕地球转的月球，只需要把月球的加速度改成"太阳引力+地球引力"的叠加。代码增加约15行，但视觉效果会丰富很多。注意月球的dt需要更小，因为它的周期只有0.0748年，建议单独给它用dt/5的步长。

7.2 用真实星历数据初始化

目前用的是圆轨道近似，如果想更真实，可以从公开的星历数据里读取某一天的行星位置和速度作为初始条件。这样轨道会变成椭圆，近日点和远日点的速度差异也能看出来。代价是需要额外处理数据格式，代码量会增加30行左右。

7.3 交互式控制

matplotlib支持鼠标和键盘事件，可以加一个暂停/继续的快捷键，或者用滑块控制时间流速。这部分不难，但会让代码超出150行的限制。我的建议是：先保证核心版本跑通，再根据需求逐步加功能，不要一开始就追求大而全。

7.4 我在这个项目里踩过的三个坑

第一个坑是单位制混乱。一开始我混用了米和AU，结果地球的初始速度算出来是每秒几米，行星直接掉进太阳里。后来统一用归一化单位才解决。第二个坑是matplotlib的blit与文字渲染冲突。如果你在图上加了行星名称的文字标签，blit=True时文字不会更新，需要把文字对象也加入返回值。第三个坑是pygame和matplotlib的窗口冲突。在某些系统上，pygame初始化窗口后会抢占matplotlib的显示，解决办法是只初始化mixer不初始化display。

这个项目最让我满意的地方是它的"性价比"——150行代码换来的不仅是一个能跑的动画，更是对数值积分、坐标变换、动画渲染、音频集成四个知识点的串联实践。如果你能独立把这150行写出来并调通，说明你已经具备了科学计算可视化的基本工程能力。后续想深入的话，建议从"加入行星间引力"这个方向入手，那时候你就真正在解N体问题了，代码量会翻倍，但理解也会深一个层次。