计算迁移时间，再重建 Kepler-223 角度

先算一张局部迁移时间图

固定四地球质量的圆轨道行星，位于一天文单位，盘厚比零点零五五，面密度尺度万分之一。公开脚本计算二十五档温度斜率、三十一档对数黏性，共七百七十五格；每格保存两项力矩、有符号漂移率与局部迁移时间。

冻结网格有六百七十八格内迁、九十七格外迁。圆点标出二百二十二格：在采用的两项力矩各上下扰动百分之二十的包络内，它们的净力矩可以变号。这是确定性的敏感性检验，不是置信区间或已标定的误差分布。物理律与单位沿用前面三章。

两项接近抵消时，局部时间很大，恰好零力矩时发散。它是当前半径除以当前漂移率的绝对值，不是积分后的抵达时间。半径变化、盘演化或动力学力矩都会改变后续路径，静态零点提供了检验额外物理的位置。

\tau_{\rm local}=\frac{a}{|\dot a|},\qquad \nu=10^x\,\mathrm{AU}^2/t_0

横坐标是黏性指数，颜色记录方向与局部时间的对数。恰好零漂移时用空值明确表示无穷时间。

            python3 migration-timescale-map.py --output my-map.json --figure my-map.svg

对应阅读：Paardekooper · Cresswell & Nelson

先确认手里是什么观测

输入包含 Mills 等 arXiv:1612.07376v1 扩展数据表1的全部四十四个季度测时偏差：b、c、d、e 四颗行星各十一档。它们由真实 Kepler 测光拟合而来，不是合成凌日时刻、逐次凌日表或后验抽样。

固定作者稿共四百二十二万九千五百零九字节，完整 SHA-256 随下载元数据保存，表在 PDF 第三十页。这里只转录有引用的事实数值和自行编写的脚本，不复制论文文字、图或作者代码。表写的是 BJD 减去二四五四九零零，没有在此明确动力学时标；不擅自补成 UTC 或 TDB。

对应阅读：Mills et al. · Izidoro et al.

保留原表的线性星历

每个偏差都相对于特定参考星历定义。若换用另一张表里的动力学拟合周期，诊断便改变了。元数据保存扩展数据表1的四个精确参考周期与零凌日历元。

季度按时间顺序匹配，但各行星实际平均观测时刻不同。共同诊断时刻取四个实时时刻的算术平均；同一季度内把测时偏差暂时固定，再分别在最早和最晚平均时刻重算，得到对齐敏感性。

行星	参考周期（日）	零凌日（原表 BJD 偏移）
b	7.3840154	70.49489
c	9.8487130	71.37624
d	14.7883997	109.76775
e	19.7213435	68.10686

对应阅读：Mills et al. · REBOUND 5.2.2 and IAS15

晚到的凌日，对应落后的相位

在近圆、近侧视轨道下，依原文约定，凌日对应平经度的四分之一圈。凌日若晚到，在同一参考时刻的相位就落后于未偏移预测，这决定了测时项的负号。

用近圆平经度近似构造两个相邻三行星拉普拉斯角，再构造四体组合。第三个在代数上依赖前两个，不能当成第三项独立发现。

\lambda_i(t)\simeq\frac\pi2+\frac{2\pi}{P_i}(t-T_{0,i}-\delta T_i)

这是近圆近似，不是完整偏心光动力学变换。

\phi_1=-\lambda_b+2\lambda_c-\lambda_d,\qquad\phi_2=\lambda_c-3\lambda_d+2\lambda_e

整数系数让共同经度零点抵消。

\phi_3=-3\phi_1+2\phi_2=3\lambda_b-4\lambda_c-3\lambda_d+4\lambda_e

应直接从原始测时系数组合计算，避免把相关角当成独立量。

对应阅读：Mills et al. · Izidoro et al.

在空文件夹里运行公开材料

下载列出的脚本、CSV、来源元数据与两组教学轨道参考。基础复算只用 Python 标准库，不导入网站、私有仓库或 Lab Engine。它会打印力矩网格、实测角度范围，以及对一段合成回放的独立短程 RK4 检查。

清单校验科学输入和脚本的固定字节。哈希说明字节一致；输出则按数值容差比较，因为不同 Python 版本可能在浮点末位不同。基础脚本不声称重算整段漫长 IAS15 历史。

            python3 reproduce-migration.py --directory . --output my-results.json
python3 analyze-kepler223.py --directory . --output my-angles.json

对应阅读：Mills et al. · REBOUND 5.2.2 and IAS15

读出实际结果与敏感性

十一季度内重建角度落在下列范围。共同时间从最早改到最晚行星平均时刻，诊断变化都小于一度。这项对齐敏感性较小，却不会消除近圆近似与有限观测窗的局限。

原表的测时区间不对称，以相对最优值的带符号偏移给出，对应注明的卡方变化。脚本把四个个别区间组成一个矩形，再求线性组合的极值。这个包络既不是联合置信区间，也不是可信区间；它没有虚构独立性，也不产生后验概率。

组合	观测窗范围（度）	最大季度对齐敏感性（度）
内侧三行星	172.440–190.844	0.216
外侧三行星	47.075–75.434	0.786
四体组合	294.286–316.566	0.924

对应阅读：Mills et al. · Izidoro et al.

分清两种复算

测时计算复算真实观测的有限诊断。合成力矩网格单独检验 Lab 的局部处方，独立编写的四阶龙格—库塔积分则比较八圈回放与 IAS15。代码时间最大步长由零点零四减到零点零二、再减到零点零一，最大笛卡尔差从约十万分之一降到约五千万分之一量级。

完整合成重算另下载可选生成器、轨道与盘模块、打包脚本和 C 力源文件，在独立环境安装固定 REBOUND 5.2.2。生成器也能使用 Python 力回调；可选编译源文件只为加速长积分，网站不分发编译后的可执行文件。比较新运行时应记录积分器、参数、时长、相位采样与停止原因。

对应阅读：REBOUND 5.2.2 and IAS15 · Paardekooper

说清它还没有建立什么

短时间角度有界，本身不证明长期共振稳定。原文解释还使用相互作用轨道拟合与稳定性计算；偏心率、投影、先验和参数协方差都有作用，季度表不能重建全部信息。

同样，两颗合成行星的气体盘模型不能推断这个真实四行星系统的诞生位置与盘参数。检验形成历史，需要把前向模型与观测相连，并明确比较竞争历史。

对应阅读：Mills et al. · Izidoro et al.

离开网站，也能下载重算

合成局部迁移时间图：内迁、外迁与符号敏感格点 — 固定四地球质量、一天文单位的圆轨道模型。圆点表示对明确分量扰动的敏感性，不作置信解释。

由真实 Kepler-223 季度测时计算的角度，以及个别轮廓区间的矩形包络 — 每点使用四个实测季度偏差；误差棒是个别轮廓区间构成的矩形包络，不作联合置信解释。

把下列基础文件放进一个空目录。第一条命令只用 Python 标准库；下载文件按哈希检查字节，跨 Python 版本的输出按数值容差比较。

                python3 migration-timescale-map.py --output my-map.json --figure my-map.svg
python3 reproduce-migration.py --directory . --output my-results.json
python3 analyze-kepler223.py --directory . --output my-angles.json

可选：重算完整合成历史

这些额外科学源文件独立于浏览器 Engine。在独立环境安装固定 REBOUND 5.2.2；不加可选本地库时，生成器使用较慢的 Python 力回调，方程相同。

                  python3 -m venv migration-env
. migration-env/bin/activate
python3 -m pip install rebound==5.2.2
python3 generate-reference.py --orbits 12000 --samples 2400 --output my-full-history.json
# Same parameters, longer physical duration:
python3 generate-reference.py --orbits 50000 --samples 5000 --output my-later-history.json

可选打包脚本重建末段密回放时需要本地力库；可把随附 C 源码对照已安装的 REBOUND 头文件编译，生成器也可用 --native-library 加载它。上方材料包含数值方法许可与完整来源哈希。

# Optional native force and dense-clip regeneration:
REB_HEADER=$(python3 -c 'import pathlib,rebound; print(pathlib.Path(rebound.__file__).resolve().parent.parent/"src")')
cc -O3 -shared -fPIC -I "$REB_HEADER" migration-force.c -o migration-force.so
python3 generate-reference.py --orbits 12000 --samples 2400 --native-library ./migration-force.so --output my-full-history.json
python3 pack-reference.py --input my-full-history.json --id capture --native-library ./migration-force.so --output my-capture-trajectory.json
# Other ready-preset parameters:
# outward: --orbits 50000 --samples 5000 --slope 1.2
# saturated: --orbits 22000 --samples 2200 --slope 1.2 --viscosity 1e-9
# transition: --orbits 50000 --samples 5000 --transition
# gravity: --orbits 50000 --samples 5000 --density 0

亲手观察

把实测与合成问题分开

先运行公开测时计算，再开教学 Lab。
比较四行星拉普拉斯组合和 Lab 的两行星临界角。
各写一项有支持的推断，并注明观测或模型边界。

练习与答案

为何不能把前两个角的误差包络当独立误差，再得到第三个？

两角重复使用相同行星的测时偏差。四体组合必须从四个原始测量与系数传播；原表没有联合协方差或后验，所以结果仍只能称矩形包络。

档案后来给出更精确周期，是否应直接换入？

若仍保留旧测时偏差，就不应直接替换。偏差相对固定参考星历定义，换参考必须一致变换偏差，并标明这是一项新分析。

来源与计算范围

Paardekooper, Baruteau & Kley, arXiv:1007.4964v1 (2010; MNRAS 2011)
第3至5.7节：力矩尺度、扩散函数、有效绝热指数与修正拟合；式3至7、19、22至23、30至31、40、45至53，图2至3、6至8。二维低质量拟合，软化长度有明确取值。
Cresswell & Nelson, arXiv:0811.4322v1 (2008)
偏心率阻尼、迁移时间与小偏心率波时间处方；这里只用共面极限，并限制偏心率小于盘厚比的零点七倍。
Mills et al., arXiv:1612.07376v1, Nature 533, 509 (2016)
Methods 中的测时与拉普拉斯角重建，扩展数据表1（PDF第30页）及图3（第35页）。教程冻结全部44个季度测时值及原线性星历。
Izidoro et al., arXiv:1703.03634v2 (24 May 2017)
第2.2至2.3、5.2、6.1、8节，图7、14。研究把迁移、共振链与后续失稳联系起来；本页两行星演示不复现其总体分布。
REBOUND 5.2.2 and IAS15
IAS15 参考积分，精度参数为十的负十次方。可选生成脚本固定 REBOUND 5.2.2，公开教程不复制该包代码。