从保存状态重算一条棒势轨道

运行前先看清小包里有什么

科学ZIP只包含自有Python程序、数值输入、来源说明和解释文字。清单逐项列出自身以外的成员，重放入口在计算前检查它们；改坏一个字节，或者指定已经存在的输出目录，都会拒绝继续。

引力势为 $\Phi=\ln R+Af(R)\cos[2(\phi-\Omega_pt)]$ ，取 $A=0.01$ 、 $\Omega_p=1$ 、 $R_b=0.5$ 。示踪星群状态用惯性坐标，周期母轨道则用随棒坐标。能量、角动量都按单位质量计，没有拟合银河系的长度或质量尺度。

与交互画面比较时，请选择任意一个母轨道预设，它们使用这份相同引力场。主默认状态的强度是0.04，属于另外一次释放。

python3 -W error replay.py --output /absolute/new/empty/run

第一步对二百五十六颗示踪星分别走五千一百二十步和一万零二百四十步，保存细网格结果。第二步用极坐标引力的独立DOP853方法对照十个原始测试初态，再计算解析圆轨道、小本轮振动和规定的振幅增长。第三、四步分别求内外两组母轨道根，实际走完各自完整周期。

第五步读取已完成的文件计算诊断；第六步从这些保存数组画图。包里没有放完成轨迹来代替积分。六条命令严格串行，每条最多运行一百二十秒。

作者第一次从空目录运行时，六步连同首次绘图和字体处理共用 $23.37\,\mathrm{s}$ ，输出约 $1.90\,\mathrm{MB}$ 。这是本地实测时间，不是所有电脑都一样的承诺；收条中分开记录了每条命令。

打开 population/population-record.npz，其中有时间、四分量惯性状态、等权采样权重，以及细步长力矩积分。每帧用 $L_z=xv_y-yv_x$ 算角动量，减去初值，再与保存的力矩时间积分比较。

图中选择的示踪星，末态 $\Delta L_z\simeq-0.0596448$ 。因 $v_0=1$ ，轴对称参考导引半径减少相同数值；这既不是当前半径，也还没有指定永久迁移后的终点。

保存的力矩差约为 $4.7\times10^{-14}$ 。KDK速度步和配套力矩求和共享一个离散角动量恒等式，因此接近舍入误差的吻合，首先检查的是内部一致性。独立十星轨迹对照的误差更大，按预定缩放范数约为 $2.44\times10^{-5}$ ，它才提供另一种轨迹精度证据。

整个示踪星群的中细差除以三，最大约 $2.53\times10^{-5}$ ，小于原定 $2\times10^{-4}$ 门槛；每次细网格漂移检查得到的最小半径约 $0.7814$ 。这些量都不是只从二百五十七个保存时刻推出来的。

六条新求母轨道的参考返回误差最大约 $1.3\times10^{-14}$ 。另外三种KDK步长还会与完整参考轨迹比较；求根返回极小，不意味着另一种离散积分也有同样小的轨迹误差。

示踪星群统计量 $Q_2$ 对双倍方位角的复指数做等权平均，幅度表示这份样本的角向不均匀。它没有按面积、质量、恒星年龄或光度加权，不能拿来代替星系亮度图。

Q_2(t)=\sum_{i=1}^{256}\frac{1}{256}e^{2\mathrm{i}[\phi_i(t)-\Omega_pt]}

复平均把双倍相对方位角的余弦和正弦分别平均，幅度表示这份等权样本的角向不均匀。

在初始半径1.45处，原宽括区两端的返回速度都为负。重放会记录实际残差，保留这次括区失败；成功的六条母轨道使用的是另行声明的窄区间。

快而暖的释放取 $A=0.04,\Omega_p=1.4,\sigma_R=0.08$ ，会穿过允许分支的内边界 $R=0.5$ 。读者程序在首次非法漂移就停止，并保存这个状态，没有忽略越界来拼出一条有效的完整轨迹。

增长例子先保留“Jacobi应当不变”这一假设的失败，再检查包含振幅功的正确关系。外层程序成功重现反例，不会让反例中的物理假设变成成功。

读新曲线前，先写参考系、单位、初始采样和完整时间变化。诊断使用的权重要与样本一致；除了守恒量，还要用独立积分或解析极限比较完整轨迹。

改变一个物理条件前，先规定误差与适用范围怎样检查。新轨道若进入排除的内侧分支，就留下失败输入；扰动若增长或消退，就先加显式功再判断守恒。末态属于这段历史，不能替代这段历史。