# fort.14（`.emflux`）全 NaN 现象分析

> 调研日期：2026-08-08
> 范围：TLUSTY 208 输出文件 fort.14（出射谱 λ-Fλ，runner 快照为 `.emflux`）中出现全 NaN 的物理与数值含义。
> 核心问题：**`.emflux` 全 NaN 意味着什么？是否代表 TLUSTY 大气模型收敛失败？**

---

## 0. TL;DR

**`.emflux` 全文件 NaN = TLUSTY 大气模型结构（TEMP/POPUL/ELEC）数值发散到不可用程度的可靠标志。** 源码机制与真实数据双重佐证：

- **源码层面**：fort.14 的流量 `FLUX` 由辐射转移求解器（RTEFR1）基于大气结构（TEMP/POPUL）计算。大气发散 → 不透明度计算产生 Inf/0 → 源函数 Inf/Inf=NaN → 整列 FLUX 全 NaN → OUTPRI 在耗尽 NITER 上限后写进 fort.14。全 NaN 只能由大气发散引起，**不存在"大气收敛良好但辐射转移单独失败"的路径**。
- **数据层面**：162 个网格点中，**62/62（100%）的 emflux-NaN 点 `converged=false`**；`converged=true` 且 emflux 含 NaN 的点 = **0**。没有任何局部 NaN 样本。

**对校验设计的结论**：emflux bolometric 校验把"全 NaN 判失败"是**正确的**，语义应表述为"**大气收敛失败（辐射转移求解发散）**"，而非仅"辐射转移失败"。它与 `atmosphere_has_nan` 同级、互为冗余，且能捕获后者的检测盲区。

---

## 1. fort.14 的写出机制

### 1.1 写出位置与格式

fort.14 由 `SUBROUTINE OUTPRI`（`tlusty208.f:14156`）写出，唯一活跃写出语句在 **`tlusty208.f:14188`**：

```fortran
14181      DO IJ=1,NFREQ
14182         IJP=IJ
14183         IF(ispodf.eq.0) IJP=JIK(IJ)
14184         IF(IJX(IJP).NE.-1) THEN
14185            WRITE(13,602) FREQ(IJP),FLUX(IJP),FH(IJP)
14186            TOTF=TOTF+FLUX(IJP)*W(IJP)
14187            FLAM=FLUX(IJP)*FREQ(IJP)*FREQ(IJP)/2.997925E18
14188            write(14,614) 2.997925e18/freq(ijp),slam
14189         END IF
14190      END DO
```

- **FORMAT** `614 FORMAT(F15.3,1pe15.3)`（`tlusty208.f:14406`）—— 第一列波长 λ(Å)，第二列 Fλ。
- **写出变量**：`FLUX(IJP)`（COMMON/SURFAC/，`MODELQ.FOR:135`），不是局部量。`FLAM` 是从 Fν 到 Fλ 的换算（`Fλ = Fν · ν² / c`）。
- **关键**：`FLUX` 是直接从 COMMON 块读已存的值，OUTPRI 本身**不重新求解辐射转移**。

### 1.2 写出时机

OUTPRI 仅在 `IF(LFIN)` 块内被调用（`tlusty208.f:3887-3889`）：

```fortran
3887:      IF(LFIN) THEN
3888:         CALL OUTPRI
3889:      END IF
```

`LFIN` 的判定（`tlusty208.f:47596`）：

```fortran
47596:      LFIN=ABS(CHMX).LE.CHMAX.OR.ITER.GE.NITER
```

即**最后一次迭代**（收敛 `CHMX≤CHMAX` 或耗尽迭代上限 `ITER≥NITER`）才写 fort.14。这意味着：fort.14 反映的是**最终大气状态**下的辐射转移解。

---

## 2. NaN 的产生链（源码追踪）

### 2.1 FLUX 的计算路径

`FLUX(IJ)` 的赋值有三处（都在辐射转移求解子程序中）：

| 位置 | 子程序 | 语句 |
|------|--------|------|
| `tlusty208.f:31800` | `RTEFR1`（Feautrier 形式解） | `FLUX(IJ)=AHout*half` |
| `tlusty208.f:32225` | `RTEDF2` | `FLUX(IJ)=AH` |
| `tlusty208.f:32848` | `RTEDF1` | `FLUX(IJ)=FH(IJ)*RAD1(1)-half*hextrd(ij)` |

这些在每次迭代的频率循环中被调用（如 `tlusty208.f:4269-4271`）：

```fortran
4269:      DO IJ=1,NFREQ
4270:         CALL OPACF1(IJ)    ← 不透明度计算
4271:         CALL RTEFR1(IJ)    ← 辐射转移求解
```

### 2.2 NaN 的具体产生点

**第一步：不透明度计算 OPACF1（`tlusty208.f:4800`）依赖大气结构**

| 源码行 | 语句 | 发散时 |
|--------|------|--------|
| `tlusty208.f:4857` / `:45375` | `HKT1(ID)=HK/T`（T=TEMP(ID)） | TEMP→0 时 HKT1→Inf |
| `tlusty208.f:4858` | `XKF(ID)=EXP(-HKT1(ID)*FR)` | HKT1·FR 极大时下溢到 0 |
| `tlusty208.f:31444` / `:33455` | `ABTRA(ITR,ID)=POPUL(II,ID)` | POPUL 含 NaN → ABTRA 含 NaN |
| `tlusty208.f:31446` | `EMTRA(ITR,ID)=POPUL(JJ,ID)*ANE*SBF(II)*...` | POPUL 极大值 × 其他 → Inf |

**第二步：辐射转移 RTEFR1（`tlusty208.f:32265`）源函数计算**

| 源码行 | 语句 | 发散时 |
|--------|------|--------|
| `tlusty208.f:32352` | `ST0(ID)=EMIS1(ID)/AB0(ID)` | **两个无穷大相除 = NaN**（典型源头） |
| `tlusty208.f:32361` | `SS0(ID)=-SCAT1(ID)/AB0(ID)` | 同上 |

NaN 一旦进入 `RAD1`（辐射转移矩阵解），就通过 `AH`/`AHout` 传播到 `FLUX(IJ)`，然后波及整列所有频率点。

### 2.3 为什么发散后不会提前 STOP

TLUSTY 有一个发散护栏（`tlusty208.f:47587-47591`）：

```fortran
47587:      IF(ITER.NE.1 .AND. ABS(CHMX).GT.1.D16) THEN
47588:         WRITE(6,610) ITER,CHMX
47590:         STOP
47591:      END IF
```

但 **CHMX 是温度相对变化的最大值**（`tlusty208.f:47468` `CHT=CHANGT(ID)/TEMP(ID)`），当 TEMP 含 NaN 时 CHMX 也变 NaN。而 **NaN 与任何数比较都返回 `.FALSE.`**：
- `ABS(NaN).GT.1.D16` 为假 → **STOP 不触发**（护栏对 NaN 失效）。
- `ABS(NaN).LE.CHMAX` 为假 → `LFIN` 只能靠 `ITER.GE.NITER`（耗尽上限）成立。

因此程序**一路跑到 NITER，然后带着含 NaN 的 TEMP/POPUL 走完最后一次频率循环，把整列 FLUX 全写成 NaN，最后 OUTPRI 把全 NaN 的 FLUX 倒进 fort.14**。

---

## 3. 两种 NaN 场景的区分

fort.14 的 NaN 理论上有两种可能来源：

| 场景 | 描述 | 实际是否存在 |
|------|------|-------------|
| **A. 大气模型发散** | TEMP/POPUL 含 NaN → FLUX 全 NaN | ✅ 存在（62 例） |
| **B. 大气收敛良好，局部辐射转移溢出** | 某些波长点不透明度极值 → 局部 NaN | ❌ 数据中不存在 |

**源码判定**：场景 B 在机制上极不可能。FLUX 全 NaN 要求所有频率点的源函数都变为 NaN，这只能由大气结构整体发散引起。局部 NaN（少数频率点）需要大气结构有效但个别频率点溢出——但 TLUSTY 的频率循环对每个频率点用同一套 TEMP/POPUL，大气有效时不会让某个频率点单独产生 NaN（除非该频率点恰好触发了吸收系数的精确除零，但这在实测数据中从未出现）。

**数据判定**：在 162 个有 `.emflux` 的网格点中，**所有 62 个 emflux-NaN 点都是 100% 全 NaN**（如 t55000_g5.0：4514/4514 行全 NaN；t60000_g6.5：4420/4420 行全 NaN）。**没有任何一个点是局部/特定波长 NaN**。这直接排除了场景 B。

---

## 4. 真实数据验证

### 4.1 全量统计（162 个网格点）

| 情形 | 数量 | conv.json converged |
|------|------|---------------------|
| `.emflux` 全 NaN + `.7` 含 NaN | **53** | 全部 false |
| `.emflux` 全 NaN + `.7` 检测无 NaN | **9** | 全部 false |
| `.emflux` 正常 + `.7` 含 NaN | **0** | — |
| 两者都正常（0 NaN） | **100** | 混合 |

**决定性事实**：
1. **62/62（100%）的 emflux-NaN 点 `converged=false`**。
2. `converged=true` 且 emflux 含 NaN 的点 = **0**。
3. 全 NaN 点都是 100% 整文件 NaN，无局部 NaN。

### 4.2 典型全 NaN 点

`t55000_g5.0_he-2_c-1_n-4_o-4`：
- `.7` 从第 11 行（depth 1 的 TEMP/ELEC/DENS/POPUL）开始全是 NaN（5200/5210 行 NaN）。
- `conv.json`：`converged:false`、`atmosphere_has_nan:true`、`best_max_relc:39100`（seed_nc 第 3 次迭代发散）。
- `.nl.err`：`IEEE_INVALID_FLAG IEEE_UNDERFLOW_FLAG IEEE_DENORMAL` 浮点异常。

### 4.3 检测器盲区：9 个"`.7` 检测无 NaN 但 emflux 全 NaN"的点

这 9 个点最关键——它们决定了 emflux 校验是否有独立价值：

- `.7` 不含字面 NaN/Inf（`atmosphere_has_nan` 检测器通过），但实际已发散成**巨大但有限的非物理值**：
  - TEMP ≈ 1.015e-15（趋近 0）
  - POPUL 达 5e12（极大）
  - `final_max_relc` 在 1e16 ~ 1e27（灾难性发散）
- 这些值不是字面 NaN，但在 OPACF1 里 `HK/TEMP`→巨大、POPUL 乘积→Inf、`EMIS1/AB0`→Inf/Inf=NaN。
- `synspec_error` 为 `"spec 含 NaN/Inf/溢出行 (共 303554 行)"` —— 发散大气喂给 SYNSPEC 后整谱也是 NaN/溢出。
- **原因**：`atmosphere_has_nan`（`conv_check.rs:18`）的正则 `(?i)(\bnan\b|\binf(?:inity)?\b|\*{3,}|[eE]\+(?:3\d{2}|[4-9]\d{2,}))` 只匹配 `nan/inf/***/E+300+`，**抓不到"巨大但有限的发散值"**（如 5e12 POPUL、1e-15 TEMP）。

**这 9 个点靠 `converged=false`（max_relc 极大）+ emflux 全 NaN + synspec NaN 三重否决才被正确标记。emflux 全 NaN 校验在这类点上恰好补位了 `atmosphere_has_nan` 的盲区。**

---

## 5. 对校验设计的影响

### 5.1 emflux bolometric 校验的语义修正

现有 `check_emflux_bolometric`（`conv_check.rs`）把"全 NaN（占比>50%）"判失败。根据本调研：

- **判定正确**：全 NaN 是大气发散的充分（且几乎必要的）标志。
- **语义需精确化**：错误描述从 ~~"辐射转移可能失败"~~ 改为 **"大气收敛失败（辐射转移求解发散）"**。

```rust
// 修正后（conv_check.rs）
error: Some(format!(
    "emflux NaN/不可解析占比 {:.0}% 过高（>{:.0}%），大气收敛失败（辐射转移求解发散）",
    nan_ratio * 100.0, 50.0
)),
```

### 5.2 与 `atmosphere_has_nan` 的关系

emflux 全 NaN 校验与 `atmosphere_has_nan` **同级、互为冗余**：

| 校验 | 覆盖的发散形态 | 盲区 |
|------|---------------|------|
| `atmosphere_has_nan` | `.7` 含字面 NaN/Inf/`***`/E+300+ | 巨大但有限的发散值（如 5e12 POPUL、1e-15 TEMP） |
| `check_emflux_bolometric`（NaN 分支） | 大气发散导致的辐射转移全 NaN | 需 TLUSTY 启用且 fort.14 存在 |
| `check_fort9`（max_relc） | max_relc 极大（收敛失败） | 无法区分真收敛/假收敛 |

三者交叉覆盖，任一独立即可拦截大部分发散点。emflux 全 NaN 校验的独特价值在于捕获那 9 个 `atmosphere_has_nan` 盲区点。

### 5.3 不影响的部分

- **局部 NaN（如 `.spec` 中 Balmer 线区域的 NaN）**：这是 SYNSPEC 处理 Balmer 线 Stark 加宽时的正常行为（与 `tests/synspec/hhe/results/hhe35nc.spec` 一致），不是收敛失败。emflux 校验不涉及 `.spec`，且阈值设在 50%（容忍局部 NaN，只拦全 NaN）。

---

## 6. 结论

1. **`.emflux` 全 NaN 的物理含义**：TLUSTY 大气模型结构（TEMP/POPUL/ELEC）已数值发散到不可用程度，导致最后一次频率循环里 OPACF1→RTEFR1 把整列 FLUX 全部污染为 NaN，OUTPRI 在耗尽 NITER 上限后写入 fort.14。

2. **是否等于 TLUSTY 收敛失败？—— 是**。源码机制（CHMX 变 NaN 后 STOP 护栏与收敛判定双双失效，全 NaN 只能在大气发散时产生）+ 数据（62/62 全 NaN 点 converged=false，converged=true 点 0 个有 NaN）双重佐证。**全 NaN emflux 是大气收敛失败的充分（且几乎必要的）标志。**

3. **"`.7` 大气无 NaN（按检测器）且 converged=true，但 `.emflux` 全 NaN"是否可能？—— 不可能**。162 点中此组合为 0 例。出现"`.7` 检测无 NaN + emflux 全 NaN"的 9 例，都是 `.7` 实际已发散成巨大有限值/近零温度（检测器盲区），converged 全为 false。

4. **emflux 全 NaN 校验的独立价值**：能捕获 `atmosphere_has_nan` 正则盲区（巨大有限发散值），起到交叉验证补位作用。

---

## 附录：关键源码行号

均在 `/home/fmq/program/tlusty/tl208-s54/tlusty/tlusty208.f`：

| 功能 | 行号 |
|------|------|
| fort.14 写出（OUTPRI） | 14188 |
| FORMAT 614 | 14406 |
| FLUX 来源 COMMON | `MODELQ.FOR:135` |
| FLUX 赋值（辐射转移解） | 31800, 32225, 32848 |
| NaN 产生（不透明度 HK/TEMP） | 4857, 45375 |
| NaN 产生（POPUL 布居） | 31444, 31446, 33455 |
| NaN 传播（源函数 EMIS1/AB0） | 32352, 32361 |
| 发散 STOP 护栏（对 NaN 失效） | 47587-47591 |
| LFIN 收敛判定 | 47596 |
| OUTPRI 调用时机（LFIN） | 3887-3889 |