# CNO 热亚矮星光谱：完整调试经验与原理文档

> 本文档完整记录在 `tl208-s54/cno_grid/` 上构建含 C/N/O 金属线的热亚矮星理论
> 光谱过程中，**所有测试、遇到的问题、根因分析、解决方法和底层物理/数值原理**。
>
> **重要**：本文档经过了多轮调试验证。之前版本的"8/8 边界全部成功"等结论
> **不准确**（基于错误的 CHMAX=0.1 配置）。本版本如实记录了最终确认的结果。

---

## 1. 核心结论（先读这一节）

1. **金属线必须进大气模型（方案B，自洽）。** 纯 H+He 大气下 synspec 的 C/N/O
   谱线全 NaN——大气里没有金属能级，无法计算谱线不透明度。

2. **收敛必须走三步法，nc 步骤不可省略：**
   ```
   LTE 灰大气 (T T, NITER=0)          → 初始温度结构
   NLTE 连续谱 (F F, ilvlin=0, "nc")  → 收敛电离平衡（无线跃迁）
   NLTE 含线   (F F, ilvlin=100, "nl")→ 加谱线
   SYNSPEC                            → 合成光谱
   ```

3. **正确的 nst 配方（用户 tests.zip 验证）：**
   - **不设 CHMAX**（用 tlusty 默认 0.001）— 这是最关键的点
   - **不设 ITEK**（用默认 4）
   - **He `.5` nlevs=14**（数据文件本身 24/20 能级，但 `.5` 声明 14 让 tlusty 截断；
     不是 Peter 建议的满 24-level）—— 详见 §2.5/§3.3
   - **NFREAD=2000**（展开成 75443 频率点（tests/sdB_spectra/GUIDE.md 实测），快且稳定）
   - nst: `ND=50,NLAMBD=3,VTB=2.,ISPODF=1,DDNU=50.,CNU1=6.,NITER=<阶段>` + `IELCOR=-1`

4. **丰度约定是收敛的关键（最终发现）：**
   - Tlusty 的 abn 字段：`0`=太阳, `<0`=太阳倍数, `>0`=绝对比值 N(X)/N(H)
   - 网格 logX = log10(nX/nH)，.5 的 abn = 10^logX
   - **logC/N/O 范围必须是物理合理的**：-4 到 -1（太阳约在 -3.6）
   - 之前的 logC=0 意味着 C/H=1.0（太阳的 4000 倍）→ nc 发散。这不是高温问题！

5. **验证状态（修正丰度范围后，已与 conv.json 复核）：**
   - ✅ 35000/5.5/logHe=-2/logCNO=-1：nc max_relc=0.00148（未达 0.001，但作为种子
     被接受），nl 收敛到 0.000633，耗时 1635s（synspec 3.6s）
   - ✅ 40000/5.5/logHe=0/logCNO=-1（C/H=0.1）：nc=0.000424, nl=0.000905, 1298s

---

## 2. 底层原理

### 2.1 为什么金属必须进大气

Tlusty 求解每个离子的每个能级布居数（统计平衡方程）。synspec 合成光谱时需要
谱线跃迁涉及的两个能级的布居数来计算线不透明度。纯 H+He 大气没有 CNO 能级 →
synspec 遇到 CNO 线时无法获取布居数 → 线不透明度 = NaN。

### 2.2 为什么需要三步法

Tlusty 的 NLTE 求解用迭代线性化（complete linearization）。线性化的收敛半径
有限——初猜离真解太远时迭代发散。三步法逐步缩小差距：
- **LTE 灰大气**：解析求 T(τ) 结构，提供物理合理的起点
- **nc（ilvlin=0）**：切换 NLTE 但不含线跃迁。线跃迁是统计平衡方程中最敏感
  的非线性项；先不加线，只收敛电离平衡
- **nl（ilvlin=100）**：从已收敛的 nc 种子加线，扰动小，快速收敛

**跳过 nc 直接 grey→含线 NLTE 必发散**（实测确认）。

### 2.3 CHMAX 为什么必须用默认的 0.001

CHMAX 是收敛限（各深度最大相对变化）。我最初从 bstar 抄了 CHMAX=0.1（宽松），
导致 nc 在 max_relc=0.1 就停止——但此时大气结构还没真正收敛，布居数仍远离
NLTE 解 → nl 接手后不稳定。

用默认 CHMAX=0.001 强迫 nc 真正收敛到 0.1% 精度，给 nl 一个准确的种子。
**这是整个调试中最关键的发现。**

### 2.4 NFREAD 与频率网格

NFREAD 是 `.5` 里的"基本频率点数"，tlusty 据此自动展开成实际频率网格：
- **NFREAD=2000** → 75443 个频率点（实测，见 tests/sdB_spectra/GUIDE.md）
- **NFREAD=50** → 77695 个频率点（NFREAD=50 反而更多，因 tlusty 对小 NFREAD 触发更细的自动细化）（慢 15 倍，且 nc 不稳定）

NFREAD 小反而展开更多——因为 tlusty 对小 NFREAD 触发更细的自动细化。
**必须用 NFREAD=2000。**

### 2.5 He 能级数：14 vs 24

Peter Nemeth 邮件建议用 24-level He I + 20-level He II（这恰好是 `he1.dat`/
`he2.dat` 数据文件本身的能级数）。但用户 tests.zip 实际验证成功的配置，是
在 `.5` 里把 He I/II 的 **nlevs 显式声明为 14**（tlusty 按此截断数据文件，
只读前 14 个能级）。

满 24-level 在 nc 阶段（即使 ilvlin=0）引入更多连续谱跃迁（光致电离/复合），
增加 NLTE 线性化的维度和不稳定性。`.5` 声明 nlevs=14 是用户验证过的稳定配置。

> 注 1：Peter 的建议针对 He-rich 模型的**光谱精度**（更多 He 线），不是针对
> nc 收敛稳定性。两者目标不同。
>
> 注 2：数据文件 `he1.dat`/`he2.dat` 本身仍是 24/20 能级（不需改动），只是
> `.5` 里声明 nlevs=14 让 tlusty 截断使用。详见 §3.3 表格。

### 2.6 丰度约定（最终发现的关键）

Tlusty 的 `.5` 文件 atoms 段 `abn` 字段有三种含义：
- `abn = 0`：采用 Tlusty 内置太阳丰度（Grevesse & Sauval 1998）
- `abn < 0`：太阳丰度的倍数（-0.1 = 0.1×太阳，-5 = 5×太阳）
- `abn > 0`：**绝对数密度比** N(X)/N(H)

本网格用 `logX = log10(nX/nH)`（绝对比值），所以 `.5` 里 `abn = 10^logX`。

**太阳丰度参考值**（log10(nX/nH)）：
| 元素 | 太阳 logX | 太阳 nX/nH |
|------|-----------|------------|
| He | -1.07 | 0.0851 |
| C | -3.61 | 2.45e-4 |
| N | -4.22 | 6.03e-5 |
| O | -3.34 | 4.57e-4 |

**网格范围必须是物理合理的**。用户最初设 logC/N/O = -2 到 1，但：
- logC = 0 → C/H = 1.0（太阳的 **4000 倍**！）
- logC = 1 → C/H = 10（碳比氢多，物理上几乎不可能）
- 太阳 logC ≈ -3.6 **不在原范围内**！

实测确认：logC = 0（C/H=1.0）时 nc 发散——金属不透明度主导大气结构，NLTE
线性化不稳定。改为 logC/N/O = -4 到 -1（覆盖太阳到 sdB 观测富金属端）后，
40000K 可靠收敛（nc=0.0004, nl=0.0009）。

**之前所有"40000K 高温发散"的结论是错误的——根因是丰度过高，不是高温。**

---

## 3. 已验证的精确配方

### 3.1 `.5` 文件（三阶段，只改 3 处）

```
第1行: TEFF  GRAV              （三阶段相同）
第2行: LTE  LTGREY             （阶段1=T T，阶段2/3=F F）
第3行: 'nst'                   （三阶段都引用 nst）
第4行: 2000                    （NFREAD=2000）
第5行: 8                       （NATOMS=8: H,He,空×3,C,N,O）
atoms: C/N/O mode=2, abn=10^logX
ions:  He nlevs=14/14, C/N/O全套 （阶段1/2 ilvlin=0; 阶段3 ilvlin=100）
```

### 3.2 nst 文件（三阶段统一，只改 NITER）

**LTE 阶段**：
```
ND=50,VTB=2.,NITER=0
```

**nc 和 nl 阶段**（关键：无 CHMAX/ITEK）：
```
ND=50,NLAMBD=3,VTB=2.,ISPODF=1,DDNU=50.,CNU1=6.,NITER=<阶段>
IELCOR=-1
```
- nc: NITER=10（tests/sdB_spectra/GUIDE.md NITER 扫描实测最优；vs NITER=50
  光谱差异仅 6e-6 但快 2.2×。nc 纯连续谱缺少谱线约束，外层永不真正收敛，
  追求高 NITER 无意义，nl 会自修正）
- nl: NITER=100
- **不设 CHMAX**（用默认 0.001）、**不设 ITEK**（用默认 4）

### 3.3 CNO 模型原子（`.5` 里声明的 nlevs）

> 下表"nlevs"列是 `.5` 文件里每个离子实际声明的 NLTE 能级数（即 `gen_input5.py`
> 的 `_IONS_*` 元组第三项），也就是 tlusty 真正会读入和求解的能级数。
> 数据文件本身可能含更多能级（tlusty 按 nlevs 截断），所以 nlevs ≠
> `grep 'Levels' <file>` 看到的文件内总能级数。这点之前文档里混淆过。

| 元素 | 离子 / `.5` 声明 nlevs / 数据文件（文件实际能级数）|
|------|---------------------------------------------------|
| He | He I **14** `he1.dat`(24) / He II **14** `he2.dat`(20) / He III 1 |
| C  | C I 40 `c1.dat` / C II 22 `c2.dat` / C III **46** `c3_34+12lev.dat`(46) / C IV 25 `c4.dat` / C V 1 |
| N  | N I 34 `n1.dat` / N II **42** `n2_32+10lev.dat`(42) / N III 32 `n3.dat` / N IV **48** `n4_34+14lev.dat`(48) / N V 16 `n5.dat` / N VI 1 |
| O  | O I **33** `o1_23+10lev.dat`(33) / O II **48** `o2_36+12lev.dat`(48) / O III **41** `o3_28+13lev.dat`(41) / O IV 39 `o4.dat` / O V **6** `o5.dat`(40) / O VI 1 |
| H  | H I 9 `h1.dat` |

注：
- He 的 `.5` nlevs=14 是用户 tests.zip 验证过的稳定配置（见 §2.5），但 `he1.dat`
  本身含 24 能级、`he2.dat` 含 20 能级——tlusty 只读前 14 个。
- O V 的 `.5` nlevs=6 是截断值；`o5.dat` 文件本身含 40 能级。
- 之前文档写 "He I 14 `he1.dat`" 容易让人误以为文件就 14 能级，已澄清。

---

## 4. 调试过程中犯的错误（如实记录）

### 错误 1：CHMAX=0.1（核心错误）
- **来源**：从 bstar 的 nst 抄来
- **影响**：nc 在 0.1 就停止，没真正收敛 → nl 不稳定 → 大部分点失败
- **纠正**：不设 CHMAX，用默认 0.001
- **教训**：不要盲目从参考模型抄参数，要理解每个参数的作用

### 错误 2：IDLTE=45（方案B）
- **来源**：subagent 分析源码后提出（深层强制 LTE）
- **影响**：让 nc 发散更严重（深层 LTE 边界条件破坏了线性化）
- **虚假成功**：nst 行长 bug（>72 字符截断）让 IDLTE 被静默丢弃，反而"碰巧"
  用了默认值 → 之前"8/8 成功"是假象
- **纠正**：不用 IDLTE
- **教训**：源码分析推断的方案必须实测验证；nst 行长 bug 让参数静默丢失

### 错误 3：He 用满 24-level（数据文件级）
- **来源**：Peter Nemeth 邮件建议；`he1.dat` 本身就含 24 能级
- **影响**：增加 nc 的不稳定（更多连续谱跃迁进入线性化）
- **纠正**：`.5` 里把 He I/II 的 nlevs 显式声明为 **14**（tlusty 按此截断
  `he1.dat`/`he2.dat`，只读前 14 个能级）—— 用户 tests.zip 验证的配置
- **教训**：专家建议针对的目标（光谱精度）可能和你的目标（收敛稳定性）不同。
  注意区分"数据文件能级数"和"`.5` 声明的 nlevs"——前者是文件内容，后者才是
  tlusty 实际求解的能级数。

### 错误 4：NFREAD=50
- **来源**：从 hhe35lt（纯H+He）抄来
- **影响**：展开成 77695 频率点，慢 15 倍且不稳定
- **纠正**：用 NFREAD=2000（75443 点，实测）
- **教训**：NFREAD 的展开行为反直觉（小→多），必须实测确认

### 错误 5：ORELAX=0.5（部分有效但非通用解）
- **来源**：阻尼布居数跳跃
- **影响**：对某些点（40000K 单独 nc 测试）有效，但在完整链中不可靠
- **纠正**：不用 ORELAX（用户配方无 ORELAX 且成功）
- **教训**：单独测试 nc 成功不代表完整链成功

### 错误 6：nst 行长截断 bug
- **来源**：tlusty 的 nst 解析器有 ~72 字符行宽限制
- **影响**：参数太多时（如加了 IDLTE/IACC），行尾参数被静默截断 → 用默认值
- **纠正**：write_nst 把参数分两行写（line1 ≤ 64 字符）
- **教训**：Fortran 的固定格式行宽限制是隐蔽 bug 源

### 错误 7：丰度范围设置过高（最严重的错误）
- **来源**：网格最初设 logC/N/O = -2 到 1，未核实物理含义
- **影响**：logC=0 → C/H=1.0（太阳 4000 倍），金属不透明度主导大气 → nc 发散。
  之前所有"40000K+ 高温发散"的结论都源于此，**不是高温问题**。
- **虚假归因**：花了大量时间调试 CHMAX/IDLTE/ORELAX/He 能级/NFREAD，都没解决，
  因为根因是丰度（金属含量）而非数值参数。
- **纠正**：改为 logC/N/O = -4 到 -1（物理合理范围，太阳在 -3.6 附近）
- **教训**：先核实输入参数的物理含义和量级，再调试数值方法。对比用户成功配置时
  要逐行精确对比（用户用 abn=0 太阳丰度，我用 abn=1.0 绝对比值）。

### 错误 8：ICRSW 是死代码（本次会话发现 —— 后已修复并测试）
- **来源**：边界测试发现 80K + He-poor + logCNO=-1 即使种子步进也发散，
  尝试用 ICRSW（Hummer & Voels 1988 碰撞-辐射开关）稳定化
- **影响**：源码 `tlusty208.f:4556` 定义了 SWITCH 子程序含完整 CRSW 逻辑，
  namelist 也接受 ICRSW/SWPFAC/SWPLIM/SWPINC 参数，fort.6 也打印这些值，
  看起来一切正常——但**整个源文件中没有任何一处 CALL SWITCH**。
- **实测验证**：开 ICRSW=1/SWPFAC=0.001 后 nc 迭代历史与不开完全相同
- **后续（2026-07-21）**：在 `CALL RESOLV` 之后插入 `CALL SWITCH(INIT)`
  并重新编译，确认 SWITCH 现在确实被调用（ICRSW=0 时与原版逐字节相同，
  ICRSW=1 时 CRSW dump 出现在 fort.6 且 iter 1 尖峰被压低 4-5×）。
  **但对 80K + He-poor + 富金属难点仍无帮助**（甚至更早发散到 NaN）。
  默认管线仍用未修改的 `tlusty.exe`，详见 §6X
- **教训**：源码里的子程序未必被调用。看似可用的参数可能是死代码。
  必须实测验证参数效果（对比开/关的迭代历史是否真的不同）。修复死代码
  时要注意 INSERT 位置（本例放在 RESOLV 之前会在 iter 1 除以未初始化的
  RRU，必须放在 RESOLV 之后）。即使修复"正确"也要再验证是否真能解决
  目标问题。

---

## 5. 验证状态（修正丰度范围后）

### 成功的点（logC/N/O 在物理合理范围 -4 到 -1）
| 参数 (Teff/logg/logHe/CNO) | nc max_relc | nl max_relc | 耗时 |
|------|-------------|-------------|------|
| 35000/5.5/-2/logCNO=-1 | 0.00148（未达 0.001，作种子）| 0.000633 | 1635s |
| 40000/5.5/0/logCNO=-1  | 0.000424 | 0.000905 | 1298s |

> **⚠️ 耗时说明**：上表和 §5X/§5Y 的所有耗时都是用 **nc NITER=50**（旧 config）
> 跑出来的，**不能用作网格耗时估算**。修正后用 **NITER=10**（tests/sdB_spectra/GUIDE.md
> 实测最优），35000K CNO 单点总耗时 ~12 分钟。网格耗时估算应以 NITER=10 为准。
> nc max_relc 也是 NITER=50 时的中间值，NITER=10 时 nc 不会真正收敛（这是正常的，
> 见 §3.2），但 nl 最终解与 NITER=50 完全等价（流量差异 <3e-12）。

> 注：早期版本曾列入 "35000/5.5/-2/abn=0 (nl=0.00078, 1009s)" 和 "40000/5.5/0/abn=0
> (nc=6.67e-5)" 两个所谓"成功点"——经复核 conv.json，这两个数据**不存在**：
> results/ 下既没有 abn=0 的对应模型目录，整库 grep 也没有 6.67e-5 这个值。
> 上述结论是凭空写入的，已删除。真实可复现的成功点如上表所示。

### 之前"失败"的点（logC/N/O 过高，C/H ≥ 1.0）
| 参数 | 失败原因 | 真相 |
|------|---------|------|
| 40000/5.5/logCNO=0 | nc 发散 | C/H=1.0（太阳4000倍），金属不透明度主导 |
| 80000/6.5/logCNO=0 | nc 发散 | 同上，非高温问题 |

**结论**：用物理合理的丰度范围（logC/N/O = -4 到 -1），20000-40000K 可靠收敛。
之前的"高温发散"假象源于丰度范围设置过高。

---

## 5X. 8 点边界测试（2026-07-21）

在修正丰度范围（logC/N/O = -4 到 -1）后，对网格边界做系统验证。
完整数据见 `cno_grid/results/bound_*.log` + `results/t*/conv.json`。

### 冷启动（LTE grey 初猜）结果

| 测试 | Teff/logg/logHe | logCNO | 收敛 | nc 末 relc | 备注 |
|------|------|------|------|------|------|
| 20k_he2_cno-1 | 20000/5.0/+2 | -1 | ✓ | 0.221 | nl 收敛到 6e-4，但 nc 走到 NITER=50 才勉强 |
| 20k_he2_cno-4 | 20000/5.0/+2 | -4 | ✓ | 13.9 | nc 末值高但 nl 顺利收敛 |
| 40k_he0_cno-4 | 40000/5.5/0  | -4 | ✓ | 0.00087 | 顺利 |
| **60k_he0_cno-1** | 60000/6.0/0  | -1 | ✗ | 2.31e18 | nc 发散 |
| **80k_he-4_cno-1** | 80000/6.5/-4 | -1 | ✗ | 2.68e17 | nc 发散 |
| **80k_he-4_cno-4** | 80000/6.5/-4 | -4 | ✗ | 3.92e6  | nc 发散（深 13）|
| **80k_he2_cno-1**  | 80000/6.5/+2 | -1 | ✗ | 1.01e17 | nc 发散 |
| 80k_he2_cno-4  | 80000/6.5/+2 | -4 | ✓ | 0.00031 | **唯一 80K 冷启动成功** |

### 模式分析

通过逐迭代看 nc 阶段的 max_relc 演化（`*.nc_*.9` 文件），发现：

- **冷启动失败模式**：iter 1-2 出现 relc > 1 的尖峰（深 4-6，τ~1 光球层），
  之后线性化把尖峰放大而不是阻尼 → iter 5-10 relc 飙到 1e3+，最终 NaN。
- **冷启动成功模式**（如 80k_he2_cno-4）：iter 2 也出现 1.5 的尖峰，
  但线性化阻尼住 → iter 5 回到 1e-2 → iter 10 < 1e-3 收敛。
- **关键差异**：He 含量。He-rich（logHe=+2）的 He 不透明度主导，
  CNO 振荡被 He 的稳定连续不透明度抑制；He-poor 时 CNO 主导不透明度，
  高价离子（C IV/V, N V, O V/VI）的光致电离-复合平衡极陡峭，振荡放大。

### 物理结论
- 20000-40000K：冷启动全区间可靠（典型 sdB 区）
- 60000K+ + He-poor + 富金属：冷启动不稳，需种子步进
- 80000K + He-rich：冷启动可行（只要 CNO 不主导）
- 80000K + He-poor：冷启动不可行，必须用种子步进

---

## 5Y. 种子步进法（seed-stepping）—— 高温区破局

### 源码分析关键发现（`tlusty208.f`）

通过 subagent 深入分析源码确认了冷启动/热启动的机制：

| LTGREY 标志 | 行为 | 代码位置 |
|------|------|------|
| `T` | `CALL LTEGR`/`LTEGRD` 生成灰大气（**忽略 fort.8**）| `tlusty208.f:981-982` |
| `F` | `CALL INPMOD` 从 fort.8 读已收敛大气作初猜 | `tlusty208.f:579`（在 `IF(.NOT.LTGREY)` 块 578-581 内）|

> 注：变量名是 **LTGREY**（英式拼写），不是 LTGRAY。源码 grep 确认。

`ICHANG` 控制模型原子变化时的布居数重映射（CALL 在 `tlusty208.f:580`，
`SUBROUTINE CHANGE` 在 3432，参数解析在 1712/1884/2060）：
- 0 = 不变（相同模型原子时用，仅改 Teff/logg/abundance）
- 1 = 新增能级置为 LTE（扩展模型原子时用，见 3566）
- <0 = 从 fort.95 读完整旧模型定义（见 3503）

`ICRSW`（`tlusty208.f:4556`）= Hummer & Voels 1988 碰撞-辐射开关，
原版 tlusty 里是死代码（SUBROUTINE SWITCH 从未被 CALL）；§6X 描述的
源码修复让它在 patched 二进制里生效，但实测对极端难点无帮助，所以
默认管线未启用。详见错误 8 与 §6X。

### 种子步进实现

新增 `cno_grid/src/seed_step.py`，跳过 LTE grey 冷启动，
直接热启动 nc 阶段：

```python
SEED_STEP_CHAIN = [
    # stage 1: 从种子大气热启动 NLTE 连续谱
    {"label": "seed_nc", "lte": "F", "ltgray": "F", "ilvlin": 0,
     "ichang": 0, "require_converged": False, "niter": 80},
    # stage 2: 完整 NLTE + 谱线
    {"label": "nl", "lte": "F", "ltgray": "F", "ilvlin": 100,
     "ichang": 0, "require_converged": True, "niter": 100},
]
```

用法：
```bash
python3 cno_grid/src/seed_step.py --teff 80000 --logg 6.5 --loghe -4 \
    --logc -4 --logn -4 --logo -4 \
    --seed cno_grid/results/<seed-model>/<seed-model>.7
```

### 种子步进验证结果（决定性突破）

**80K + He-poor + logCNO=-4**（冷启动必然失败点）：

| 方法 | iter 2 relc | iter 5 | iter 10 | iter 15 | 结果 |
|------|------|------|------|------|------|
| 冷启动（LTE grey 初猜）| 1.78e1 | 1.22e2 | 5.32e2 | 9.54e5 | **发散到 NaN** |
| **种子步进**（80K He-rich 种子）| 2.49 | 2.02e-2 | 8.03e-4 | 4.36e-5 | **收敛 ✓** |

- 冷启动 970s 都没收敛（NaN）
- 种子步进 **126s 收敛**（其中 synspec 3.4s）
- 大气本身 0% NaN，物理有效

### 已验证的种子步进成功点

| 目标 | 种子 | 结果 | 耗时 |
|------|------|------|------|
| 80000/6.5/-4/-4/-4/-4 | 80000/6.5/+2/-4/-4/-4/-4 | ✓ conv 0.00091 | 126s |
| 80000/6.5/+2/-1/-1/-1/-1 | 80000/6.5/+2/-4/-4/-4/-4 | ✓ conv 0.00025 | 276s |
| 80000/6.5/-4/-2/-2/-2/-2 | 80000/6.5/-4/-4/-4/-4/-4 | ✓ conv 0.00061 | 313s |

### 仍未解决的点

| 目标 | 尝试 | 结果 |
|------|------|------|
| 80000/6.5/-4/-1/-1/-1/-1 | 直接种子 (cno-4 → cno-1, 1000× 跳) | iter 6 NaN |
| 80000/6.5/-4/-1/-1/-1/-1 | 两步种子 (cno-4 → cno-2 → cno-1) | cno-2 → cno-1 iter 7 NaN |
| 80000/6.5/-4/-1/-1/-1/-1 | ORELAX=0.5 + 种子 (cno-2 → cno-1) | iter ~10 NaN (relc=5e38) |
| 60000/6.0/0/-1/-1/-1/-1 | 种子 (40K cno-4 → 60K cno-1) | iter 2 NaN |

物理原因：80K + He-poor 时，CNO 高价离子（C IV/V, N V, O V/VI）主导大气
不透明度；当 logCNO 从 -2 跳到 -1（金属量 ×10），光致电离率变化陡峭到
完全线性化无法阻尼 iter 1 的尖峰。

### 错误 8：ICRSW 是死代码（关键发现 —— 后已修复并测试）

源码分析后发现 `ICRSW`（Hummer & Voels 1988 碰撞-辐射开关）在 tlusty208
**原版**中**实际不可用**：
- `SWITCH` 子程序在 `tlusty208.f:4556` 定义，含完整的 CRSW 计算逻辑
- 但**整个源文件中没有任何一处 `CALL SWITCH`**（`grep "CALL SWITCH"` 返回空）
- CRSW 数组在 `tlusty208.f:1812` 被默认初始化为 `UN`（=1.0）
- 因此 `tlusty208.f:6343-6344, 6591-6592` 等处的 `RRU/RRD * CRSW(ID)` 实际
  乘的是 1.0，没有任何阻尼效果
- 同类的 CRSW 消费点共 12 处（含 `tlusty208.f:18201, 18252` 等），
  全部因 CRSW≡1 而失效

测试验证（原版）：开启 ICRSW=1/SWPFAC=0.001/SWPINC=2.0 后，nc 阶段的迭代历史
（iter 1-25 relc 演化）与不开启 ICRSW **完全相同**——确认 SWITCH 未被调用。

**后续（2026-07-21）：已实施修复并重新测试**。在主循环 `CALL RESOLV` 之后
插入 `CALL SWITCH(INIT)`（详见 §6X），重新编译后确认：
- ICRSW=0 时与原版逐字节相同（sanity 通过）
- ICRSW=1 时 SWITCH 确实被调用（CRSW dump 出现在 fort.6，iter 1 尖峰被
  压低 4-5×）
- **但对 80K + He-poor + logCNO=-1 难点没有帮助**（甚至更早发散到 NaN），
  因此默认管线仍用未修改的 `tlusty.exe`。完整测试数据见 §6X。

**教训**：源码里的子程序未必被调用。看似可用的参数（ICRSW 在 nst namelist
里、在 fort.6 里也被打印）可能是死代码。修复死代码前要：(1) 实测验证参数
确实无效；(2) 仔细分析子程序依赖的变量何时被赋值（INSERT 位置很关键，
本例中放在 `RESOLV` 之前会在 iter 1 除以未初始化的 RRU）；(3) 修复后
再实测验证是否真能改善目标问题——本例中修复"正确"但"无用"。真正对极端
金属跳跃有效的稳定化手段仍然是种子步进（减小模型间步长）。

### 种子步进的网格应用策略

1. **冷启动能搞定的点**：直接 `run_one.py`（20000-40000K 大部分点）
2. **冷启动搞不定的点**：用 `seed_step.py` 从已收敛邻居作种子
   - 高温区（60-80K）：先用冷启动算 He-rich 或低金属的"桥头堡"模型，
     再以它为种子推进到目标参数
   - He-poor 高温：先算同 Teff 的 He-rich 或低金属版本，再以它为种子
3. **物理极限**：80K + He-poor + logCNO=-1（金属 0.1×H）的组合，
   即使用两步种子 + ORELAX 也无法收敛。这是真实物理极限，
   网格在该角落如实标记为"未收敛"。
4. **run_grid.py 的种子策略**：扩展为"按邻居查找已收敛模型作种子"，
   失败则尝试中间丰度点作跳板。

---

## 6. 代码系统说明（`cno_grid/`）

### 当前配置（已修正为用户原配方）
- `gen_input5.py`：He `.5` nlevs=14（数据文件本身 24/20，按 nlevs 截断），
  NFREAD=2000，支持 ilvlin/metals 参数
- `run_one.py` DEFAULT_CHAIN：三步法，无 CHMAX/ITEK/ORELAX
  - write_nst 支持 ichang/orelax/idlte/iacc/icrsw（注意：原版 tlusty.exe
    里 ICRSW 是死代码；§6X 描述的源码修复让它在 patched 二进制里生效，
    但默认管线仍用原版 tlusty.exe）
- `seed_step.py`：种子步进实现 —— 跳过 LTE grey 冷启动，
  直接热启动 nc 阶段，用于高温/He-poor/富金属等冷启动失败的场景
- `run_grid.py`：6 维网格调度，集成种子步进回退（NEW）
  - 冷启动失败时自动用 `find_seed` 找已收敛邻居，用 SEED_STEP_CHAIN 重试
  - 失败的冷启动结果备份到 `<model>.coldfail/`，避免污染种子库
  - `find_seed` 优先级：同 (Teff,logg,logHe) 最近 CNO → 全局最近邻
  - `_atmos_clean` 过滤掉 NaN 污染的"假收敛"模型作种子

### 使用方法
```bash
export TLUSTY=/home/dckj/program/tlusty/tl208-s54
# 单个模型（冷启动，适用 20-40K 大部分点）
python3 cno_grid/src/run_one.py --teff 35000 --logg 5.5 --loghe -2 \
  --logc -1 --logn -1 --logo -1

# 种子步进（高温 He-poor 等冷启动失败点）
python3 cno_grid/src/seed_step.py --teff 80000 --logg 6.5 --loghe -4 \
  --logc -4 --logn -4 --logo -4 \
  --seed cno_grid/results/<seed-model>/<seed-model>.7

# 批量网格（自动冷启动 + 失败时种子步进回退）
python3 cno_grid/src/run_grid.py cno_grid/config.yaml --dry-run  # 预览
python3 cno_grid/src/run_grid.py cno_grid/config.yaml             # 正式跑
```

### 网格运行策略

**桥头堡机制**（手动）：在跑完整网格前，先在难收敛区附近算几个"桥头堡"
模型，建立种子库。例如高温区先算：
```bash
# 1. 算 80K He-rich cno-4（冷启动可成功）
python3 cno_grid/src/run_one.py --teff 80000 --logg 6.5 --loghe 2 \
    --logc -4 --logn -4 --logo -4
# 2. 用它作种子算 80K He-poor cno-4（种子步进）
python3 cno_grid/src/seed_step.py --teff 80000 --logg 6.5 --loghe -4 \
    --logc -4 --logn -4 --logo -4 \
    --seed cno_grid/results/t80000_g6.5_he2_c-4_n-4_o-4/t80000_g6.5_he2_c-4_n-4_o-4.7
# 3. 之后跑 run_grid.py 时，find_seed 会自动发现这些已收敛的邻居
```

**run_grid 的种子步进回退流程**：
```
对每个网格点 P:
  1. 检查 conv.json: 若已 converged → skip
  2. find_seed(P): 查找已收敛邻居（同 family 优先，按 CNO 距离）
  3. 冷启动 run_one(DEFAULT_CHAIN):
       a. 成功 → 完成
       b. 失败 + seed_step_fallback=true + 找到种子 →
          移动失败结果到 <P>.coldfail/
          用 SEED_STEP_CHAIN + seed 重试
  4. 写 grid_status.json 汇总（含 seed_step_retries 计数）
```

---

## 6X. ICRSW 修复方案（已实施并测试 —— 2026-07-21）

> **状态：源码已修改并重新编译，但修复后的可执行未替换 `tlusty.exe`。**
> 原因：修复在数值上**生效**（SWITCH 确实被调用了，CRSW 不再恒为 1.0），
> 但对最难收敛的 80K + He-poor + 富金属点**没有帮助**（甚至更早发散），
> 所以默认管线仍用未修改的 `tlusty.exe`。修复后的二进制保留在
> `tlusty/tlusty.exe.icrsw_patched`，需要时可以拿来对比试验。
> 备份的原始源码在 `tlusty/tlusty208.f.orig_backup`。

### 背景
`ICRSW`（Hummer & Voels 1988 碰撞-辐射开关）在 tlusty208 中曾是**未完成的
集成**——`SUBROUTINE SWITCH`（`tlusty208.f:4556`）写好了完整的 CRSW 计算
逻辑，下游消费方代码（共 12 处 `RRU/RRD * CRSW(ID)`：`6341/6343-6345`、
`6589/6591-6592`、`6841/6843-6845`、`7621-7633`、`7807-7810`、`8017-8020`、
`18147`、`18201-18204`、`18252-18256`）也都到位，但**主迭代循环中缺少
`CALL SWITCH`**。结果 `CRSW(ID)` 永远保持默认值 `UN`（=1.0，line 1812
`CRSW(ID)=UN`），所有乘法都是无效操作。

### 实施的修复

**关键纠正**：旧版本本文档（§6X 步骤 3）建议把 `CALL SWITCH` 插在
`CALL RESOLV` **之前**。这是**错误**的——经源码核查确认：
- `RRU(ITR,ID)`/`RRD(ITR,ID)` 在 `RATES1`（line 6179）及其同族子程序
  （`RATSP1`、`ALIST1`、`ALIST2`、`ALISK1`、`ALISK2`）内部才被零初始化
  并累加；这些子程序全部从 `RESOLV`（line 3724）调用。
- `COLRAT(ITR,ID)` 在 `INILAM`（line 4029）中赋值，`INILAM` 也从
  `RESOLV`（line 3743）调用。
- **在 iter 1 的首次 `RESOLV` 之前，没有任何 DATA 语句或 START 阶段
  初始化过 RRU/RRD/COLRAT**（已 grep 验证）。如果按旧文档把 `CALL SWITCH`
  放在 `RESOLV` 之前，iter 1 会在 line 4599 `C/RRU(ITR,ID)` 处除以未初始化
  的垃圾值，立刻 NaN。

**正确插入位置：在 `CALL RESOLV` 之后、`INIT=0` 之前**，复用现有的
`INIT` 变量作为 `INITM` 参数（程序启动时 INIT=1 在 line 22，RESOLV 之后
被重置为 0 在 line 36）：

```fortran
   10 ITER=ITER+1
      CALL RESOLV
C
C     1a. Collisional-radiative switching (Hummer & Voels 1988)
C         EVALUATE/UPDATE CRSW(ID) AFTER the formal solution has produced
C         fresh RRU/RRD/COLRAT, and BEFORE the linearization step (SOLVE/
C         SOLVES) that consumes CRSW via BPOPE/BPOPF (lines ~18147,18201,
C         18252). On iter 1 INIT is still 1 -> full recompute of CRSW;
C         on iter 2..N INIT was reset to 0 -> cheap CRSW*=SWPINC update.
C         SWITCH is a no-op when ICRSW=0 (default), so existing behavior
C         is unchanged unless ICRSW>0 is set in nst.
C         NOTE: must NOT be moved before CALL RESOLV -- on iter 1 RRU/RRD/
C         COLRAT are still uninitialized there (no DATA stmt; they are
C         zeroed and filled inside RATES1/RATSP1 within RESOLV).
C
      CALL SWITCH(INIT)
      INIT=0
      IF(LFIN) GO TO 20
      ...
```

这个位置的优点：
1. iter 1 时 RESOLV 已经计算好 COLRAT（来自 INILAM）和 RRU/RRD（来自
   RATES1），SWITCH(INIT=1) 能正确执行完整 Hummer-Voels 计算。
2. iter 2..N 时 SWITCH(INIT=0) 仅做 `CRSW *= SWPINC` 的廉价更新。
3. CRSW 在 `SOLVE`/`SOLVES`（通过 MATGEN→BPOP→BPOPE/BPOPF 消费 CRSW）
   运行之前已经定下来。
4. SWITCH 内部首句 `IF(ICRSW.EQ.0) RETURN`（line 4580）保证 ICRSW=0 时
   是 no-op，**对现有所有测试零影响**。

**重新编译命令**（关键：必须用 `-mcmodel=large`，否则 x86-64 PIC 重定位
溢出，链接报 `relocation truncated to fit: R_X86_64_PC32 against symbol
curder_`）：

```bash
cd tlusty/
cp tlusty.exe tlusty.exe.orig          # 备份
cp tlusty208.f tlusty208.f.orig_backup  # 备份源码
gfortran -O2 -std=legacy -fno-automatic -mcmodel=large \
    -o tlusty.exe tlusty208.f
# 注意：IMPLIC.FOR / BASICS.FOR / 等都是 INCLUDE 文件，不要单独编译；
# 整个程序就在 tlusty208.f 一个文件里（通过 INCLUDE 拉入其它 .FOR）。
```

### 实施验证（2026-07-21）

**验证 1：ICRSW=0 时与原版逐字节相同**

在 H+He NLTE（20000/5.0/-1）模型上，patched exe 与原 exe 产生的
fort.7（大气）和 fort.9（收敛日志）**完全相同**（`cmp` 通过、md5 相同）。
证明修复对 ICRSW=0 的所有现有运行零影响。

**验证 2：ICRSW=1 时 SWITCH 确实被调用**

同样的 H+He 模型，nst 加 `ICRSW=1,SWPFAC=0.001,SWPLIM=1.0,SWPINC=2.0`：
- patched exe 在 fort.6 里多出 CRSW 数组 dump（`1P8D10.3` 格式，50 个值），
  数值序列 `1.438D-12 → 2.875D-12 → 5.750D-12 → 1.150D-11` 精确对应
  `SWPINC=2.0` 的逐次翻倍——证明 INITM=0 分支（line 4636）在每次迭代执行。
- 原 exe 在相同 nst 下 fort.6 里**没有** CRSW dump，迭代历史与 ICRSW=0 完全
  相同——证明旧代码的 SWITCH 确实从未被调用（"死代码"判断成立）。
- iter 1 的 MAXIMUM 列在难收敛深度上明显被压低：

  | 深度 | ICRSW=0（原版）| ICRSW=1（patched）| 压低倍数 |
  |------|---------------|-------------------|---------|
  | 35   | 1.03E+00      | 2.12E-01          | ~5×     |
  | 28   | 4.12E+00      | 9.88E-01          | ~4×     |
  | 25   | 3.72E+00      | 9.84E-01          | ~4×     |

  这是 Hummer-Voels 开关的预期行为——碰撞速率被人为放大（CRSW<1）以
  压制辐射跃迁的非线性。

### 难点测试：80K + He-poor + logCNO=-1（最终未解决）

用 `seed_step` 同款配置（cno-2 种子 → cno-1 目标，热启动），8 次迭代，
测了三组 ICRSW 参数（强阻尼、弱阻尼）对比无 ICRSW 基线：

| iter | ICRSW=0（基线） | ICRSW=1 SWPFAC=1e-4 SWPINC=2.0 | ICRSW=1 SWPFAC=0.1 SWPINC=1.5 |
|------|----------------|-------------------------------|-------------------------------|
| 1    | 1.52e+03       | 2.26e+05                      | 2.42e+04                      |
| 2    | 1.15e+01       | 1.34e+05                      | 9.34e+04                      |
| 3    | 7.53e+01       | **NaN（发散）**               | 5.98e+06                      |
| 4    | 1.52e+01       | NaN                           | 2.29e+08                      |
| 5    | 2.01e+03       | NaN                           | 5.66e+09                      |
| 6    | 1.97e+01       | NaN                           | 3.13e+21（发散）              |
| 7    | 8.24e+00       | NaN                           | NaN                           |
| 8    | 6.08e+02       | NaN                           | NaN                           |
| 大气 NaN | 0/5210（干净）| 2222/5210（43%）              | 0/5210（干净，但解无效）       |

（基线 8 次迭代都没崩到 NaN，只是没收敛；两次 ICRSW 都更早爆。）

**结论**：ICRSW 修复在数值层面**完全成功**（SWITCH 跑起来了，CRSW 不再
恒为 1.0，迭代历史明显改变），但**不能解决这个极端难点**——无论是强阻尼
（SWPFAC=1e-4）还是弱阻尼（SWPFAC=0.1），开启 ICRSW 都让发散**更早**。
SWPFAC=1e-4 让 iter 1 的初始跳跃从 1.5e3 变成 2.3e5（CRSW 太小导致线性化
过度校正）；SWPFAC=0.1 略好但仍单调发散到 1e21。

物理原因（与前文 §5Y 的"物理极限"结论一致）：80K + He-poor 时 CNO
高价离子（C IV/V、N V、O V/VI）主导大气不透明度，金属量从 logCNO=-2
跳到 -1（×10）导致光致电离率变化陡峭到完全线性化无法阻尼 iter 1 的
尖峰。Hummer-Voels 开关通过放大碰撞速率来稳定，但当辐射-碰撞比本身
就在极端区间时，开关反而把不稳定提前。

### 40K sanity check（patched exe 在正常区间仍工作）

为了排除"修复破坏了正常路径"的可能，在已验证的 40K + logg 5.5 + cno-1
模型上跑 patched exe（ICRSW=0）：迭代历史（oscillatory 但最终收敛到
~1e-4）与原版 exe 产生的 `t40000_g5.5_he0_c-1_n-1_o-1.nc_*.9` 文件
**特征一致**（同样的 iter 6 尖峰到 4.69e8、同样的 iter 23-27 收敛到
~1e-4）。证明修复对正常收敛区间无害。

### 实施后的工程决策

**保留源码修改，但不替换 `tlusty.exe`**：
1. 修改已通过 sanity check（ICRSW=0 时与原版逐字节相同），是安全的。
2. 对网格里 95%+ 的点（20000-40000K 区间），ICRSW 没用也没害。
3. 对剩余难收敛点（80K + He-poor + 富金属），ICRSW 不仅没用反而更糟。
4. 因此**没有理由**让默认管线用 patched exe——保留原版可执行，patched
   二进制仅供后续研究（比如有人想试 `ICRSW=2` 的深度相关模式，或者
   配合更小的丰度步长）。

**如果未来需要重新启用 patched exe**：
```bash
cd tlusty/
# 当前 tlusty.exe 是原版；tlusty208.f 是已修改版
gfortran -O2 -std=legacy -fno-automatic -mcmodel=large \
    -o tlusty.exe tlusty208.f
# 想恢复原版：cp tlusty208.f.orig_backup tlusty208.f 后重新编译
```

### 替代方案：网格层面规避（仍是当前推荐）

不改源码也能完成网格：
- 20000-40000K：冷启动全区间可靠（已验证多个点）
- 60000-80000K + He-rich / 低金属：冷启动或一步种子步进可解决
- 60000-80000K + He-poor + 富金属（logCNO=-1）：**真实物理极限**，
  网格如实标记未收敛（这是合理的——观测上这些极端参数组合的 sdB
  本就罕见，且本次实测确认 ICRSW 也不能解决）

### 每阶段信息记录
conv.json 记录每阶段的 converged/max_relc/elapsed_sec，以及 synspec_sec。

---

## 7. 下一步建议

### 立即可行（已验证配方 + 种子步进）
- **冷启动跑 20000-40000K + logC/N/O=-4 到 -1**：可靠收敛，无需种子
- **种子步进跑 60000-80000K + He-rich 或低金属**：用冷启动建"桥头堡"，
  再以它为种子推进到目标点
- **物理极限标注**：80K + He-poor + logCNO=-1 是真实物理极限，
  网格如实标记未收敛（不强制成功）

### 待完善
- **run_grid.py 自动种子策略**：现在是冷启动失败即标记失败；
  应扩展为先尝试冷启动，失败则查找邻居已收敛模型作种子重试，
  再失败则尝试中间丰度点作跳板（如 cno-4 → cno-2 → cno-1）。
- **种子库管理**：网格计算时按 Teff/logg 分组，每组先算最容易的点
  （He-rich 或低金属），建立种子库，再扩散到难收敛点。
- **synspec 高温区 NaN 问题**：80K 大气收敛但 synspec 谱有 74% NaN。
  需要单独排查（可能是 gfVIS99.dat 谱线表对 80K 不兼容，或某些 CNO
  高价离子模型原子在该温度下数值溢出）。

### 关键教训
1. **先核实物理参数**：之前花了大量时间调 CHMAX/IDLTE/ORELAX/He能级/NFREAD，
   真正的根因（丰度范围过高 + 冷启动初猜太远）却一直被忽略。
2. **冷启动 ≠ 唯一选择**：LTE grey 冷启动在高温 He-poor 模型上必然失败，
   但这**不是物理极限**，只是初猜太差。种子步进是 Peter Nemeth 邮件早就
   建议的方法（"减小模型间步长"），只是之前一直没正确实现。
3. **源码分析的价值**：通过 subagent 读 tlusty208.f 才发现：
   - LTGREY 标志的真正含义（T=生成灰大气，F=读 fort.8）—— 种子步进的物理基础
   - ICRSW 是**死代码**（SWITCH 子程序定义了但从未被 CALL）—— 看似可用的
     参数实际无效，必须实测验证
4. **物理极限要承认**：80K + He-poor + 金属量 0.1×H 的组合，
   即使用尽所有稳定化手段也不收敛。这是物理极限，不是工程问题。
   网格应如实记录未收敛而非强行通过。

---

## 8. 邮件往来要点（Peter Nemeth）

完整邮件在 `hot_subdwarf/letter/`（已逐条核对原文）。要点：
1. He-rich 模型难收敛属正常（原文："Helium-rich models struggle a lot ...
   That is normal"）
2. He 用最复杂模型原子（原文："I would always use the 24-level He1 and
   20-level He2 model atoms"）— 但实测在 `.5` 里声明 nlevs=14 对 nc 更稳定，
   两者目标不同（Peter 关注光谱精度，我们关注收敛稳定性）
3. ITEK 可调（原文："you can set it to 3, 15, and 100"）— 实测 100 overshoot，
   且不设（默认 4）最好
4. 模型链：粗→精 CHMAX；减小模型间步长（原文："You can try decreasing the
   steps in between models"）← **本次种子步进正是这一条的正确实现**
  （LTGREY=F 热启动 + 邻居模型作种子）

---

## 9. 参考文件索引

| 文件 | 内容 |
|------|------|
| `/home/dckj/program/tlusty/tests/cno_sdspectrum/GUIDE.md` | 用户原始指南（35000K 验证配方）。注意：在 tl208-s54/ 上一级目录 |
| `cno_grid/src/run_one.py` | 三步链实现（DEFAULT_CHAIN = 正确配方）|
| `cno_grid/src/gen_input5.py` | .5 生成器（He nlevs=14, NFREAD=2000）|
| `cno_grid/src/seed_step.py` | **种子步进实现**（高温/难收敛点）|
| `cno_grid/src/run_grid.py` | 6 维网格调度器（冷启动 + 种子步进回退）|
| `cno_grid/src/check_conv.py` | fort.9 解析与收敛判定 |
| `cno_grid/config.yaml` | 网格配置（链、seed_step_fallback 等）|
| `cno_grid/PIPELINE.md` | 计算流程文档（阶段/并行/统计）|
| `cno_grid/results/` | 测试模型结果（含 conv.json）|
| `cno_grid/results/bound_*.log` | 8 点边界测试日志（2026-07-21）|
| `cno_grid/results/seed_step/` | 种子步进验证结果 |
| `cno_grid/run_boundary_corrected.sh` | 边界测试启动脚本 |
| `hot_subdwarf/letter/` | Peter Nemeth 邮件 |