数值求出并检查临界风解

先把输入存好

把下方的脚本、输入、出处说明、运行说明与固定输出下载到同一普通文件夹。输入定义了五组单组分合成大气：默认风、强束缚尾部、稀薄大气、XUV 功率不足，以及临界点低于储气边界。这里没有实测行星剖面，也没有 MESA 演化输出。

JSON 固定 SI 常数、边界半径、密度、给定温度、XUV 辐照度与可用功率比例。脚本把输入的精确字节哈希写入输出。第一次先原样运行，需要改参数时另存副本。出处文件列明字节数和 SHA-256，并区分自行给定的教学输入与方程来源。

python3 -W error replay-wind.py --input escape-input.json --out wind-replay-rerun.json
# 可选解析参考，在虚拟环境中安装固定依赖：
python3 -m venv .venv
source .venv/bin/activate
python3 -m pip install -r requirements-reference.txt
python3 reference-lambert.py --out lambert-reference-rerun.json

对应阅读：Cranmer 2004, New views of the solar wind with the Lambert W function, v1 · Volkov et al. 2011, Thermally driven atmospheric escape, arXiv v1

把分子变零的位置括起来

对 $2c^2/r-GM/r^2$ 在 $0.01R_{\rm b}$ 到 $100R_{\rm b}$ 间做二分：每次保留两端异号的半区间。如果没有异号，脚本直接报错。固定五个例子均在范围内；换新输入时要重新确认。

把根与 $GM/(2c^2)$ 比较，默认得到 $2.41574013\times10^7\,\mathrm m$ 。找到这个位置，只是让分子变零，速度还必须沿加速分支经过 $u=c$ 。若直接在该点做除法，分母也为零。

对应阅读：Cranmer 2004, New views of the solar wind with the Lambert W function, v1

从正则斜率出发

积分变量取 $q=\ln\mathcal M$ ，这样极小的亚声速仍能保持有限的对数值。在 $x=1\pm10^{-5}$ 用局部展开给初值，再以四阶 Runge–Kutta 向内或向外积分。离声速点不足 $10^{-4}$ 的目标，直接使用展开。

最大径向步长为 $0.0015\min(x,1)$ 。每个目标半径都从声速附近独立积分。输出还在 $0.03r_{\rm s}$ 至 $12r_{\rm s}$ 的七个位置，把最大步长减半，保存对数 Mach 数的变化。这个检查衡量积分的步长敏感性；解析函数提供另一条求解路径，可以直接比较同一半径的 Mach 数。

\frac{dq}{dx}=\frac{2(1/x-1/x^2)}{e^{2q}-1},\quad q=\delta-\delta^2+\frac{13}{12}\delta^3+O(\delta^4)

远离精确临界点时，专门的指数减一函数减小分母相消误差；临界邻域本身由展开处理。

对应阅读：Cranmer 2004, New views of the solar wind with the Lambert W function, v1

与两个解析分支比较

设 $D=x^{-4}e^{3-4/x}$ 。Lambert-W 函数由 $W(z)e^{W(z)}=z$ 定义，在这里有两个实分支。SciPy 将它们记为零与负一。

可选脚本在声速内取主分支，外侧取另一实分支；精确 $x=1$ 直接指定正则极限。它在十一个半径与 ODE 比较。在两倍声速半径，如果错取主分支，会得到亚声速的小根。两根都可能让隐式方程残差很小，因此还必须检查分支与加速方向。

\mathcal M^2=-W_0(-D)\quad(x<1),\qquad\mathcal M^2=-W_{-1}(-D)\quad(x>1)

整条加速跨声速解，在两个实分支共同的声速端点换支。

对应阅读：Cranmer 2004, New views of the solar wind with the Lambert W function, v1

分别检查数值与物理条件

从底部边界恢复密度，逐点计算 $4\pi r^2\rho u$ 。另算等温动量积分 $u^2/2+c^2\ln(n/n_{\rm b})-GM/r$ 。前者检查密度归一，后者对独立积分的速度形状误差敏感。

固定默认解的质量通量相对极差约为十的负十五次方，动量积分极差以声速平方为单位小于十的负十次方。这些是数值残差，没有包含温度、碰撞截面或忽略冷却造成的不确定度。

下方表格与曲线直接读取公开重跑输出。强束缚例子的 Jeans 示例约为每秒十毫克，形式 Parker 声速区的 Knudsen 数却极大。功率不足例子在数值上平滑，所需机械功率约为可用 XUV 功率的752倍。临界点埋入底部的例子，从储气边界起就已超声速。输出保留这些解释失败的案例。

独立正则 ODE 的默认合成结果。横轴为半径与储气半径比，采用对数刻度；虚线为数值找到的临界半径。

replay-wind.py 的实际合成结果：固定类中性氢质量与有效碰撞截面，给定温度，储气半径为地球的两倍，外边界为储气半径的十二倍。输入列明每例全部设置；即使物理条件失效，表中仍列形式风解流量。
情形	声速与底部半径比	形式 Parker 流量 · 千克每秒	声速区 Knudsen 数	所需与可用功率比
默认风	$1.9$	$7.70\times10^{6}$	$4.02\times10^{-3}$	$2.28\times10^{-1}$
罕见尾部	$1.52\times10^{1}$	$6.77\times10^{-4}$	$1.83\times10^{8}$	$2.05\times10^{-9}$
稀薄气体	$5.06$	$7.90\times10^{2}$	$9.05\times10^{1}$	$2.80\times10^{-4}$
热量不足	$1.26$	$1.84\times10^{7}$	$1.37\times10^{-3}$	$7.52\times10^{2}$
临界点低于底部	$2.11\times10^{-1}$	$7.11\times10^{7}$	$7.26\times10^{-5}$	$3.78$

对应阅读：Volkov et al. 2011, Thermally driven atmospheric escape, arXiv v1 · Murray-Clay, Chiang & Murray2009, Atmospheric escape from hot Jupiters, v2

只改一个条件试试

复制输入，只把默认例子的底部数密度加倍。运行前先预测风流量、声速区 Knudsen 数及所需与可用功率比。三者应分别乘二、除二、乘二。速度方程只含引力与声速，因此速度形状保持不变。

再把给定外半径从十二倍底部半径改为二十四倍，重跑机械能比较。气体走得更远，还要继续爬升并加速，需要供给更多能量。引用等温风的热量需求时，要带上这个外边界。

公开文件可按 README 中的科学教学材料许可复用。若要研究某颗具体行星，下一步应加入能量与化学方程，并为它建立合适的内外边界。

对应阅读：Murray-Clay, Chiang & Murray2009, Atmospheric escape from hot Jupiters, v2 · Cranmer 2004, New views of the solar wind with the Lambert W function, v1

数值求出并检查临界风解

先把输入存好

把分子变零的位置括起来

从正则斜率出发

与两个解析分支比较

分别检查数值与物理条件

只改一个条件试试

练习与解释性答案

动手检验预测

把保存的风解与交互计算比较

来源与计算范围