DCTS/.agents/skills/tlusty-synspec-test/scripts/check_temperature.py

#!/usr/bin/env python3
"""
check_temperature.py — 对 fort.7/{stage}.7 做温度结构物理判定。

精确复刻 DCTS 权威判据 check_temperature_structure（conv_check.rs:463）。
退出码：0=物理（三项全过）/1=不物理 /2=文件缺失或无法解析（跳过，不判失败）。

判据（valid = 三项全过）：
  1. 表层 T < max_factor × Teff   （默认 max_factor=3.0）
  2. 每个深度 T ∈ [temp_floor, temp_ceiling]   （默认 [10, 1e8] K）
  3. 无 NaN/Inf

用法：
    python3 check_temperature.py outputs/nl.7 --teff 60000
    python3 check_temperature.py outputs/nc.7 --teff 55000 --max-factor 3.0
"""
import argparse
import math
import re
import sys

# 复刻 conv_check.rs:44 parse_fortran_float
# Fortran 指数≥100 时挤掉 E，如 -1.35E+118 → -1.35+118
_NO_E_EXP_RE = re.compile(r"^([+-]?[\d.]+)([+-]\d+)$")


def parse_fortran_float(s):
    """返回 float（可能为 inf/-inf）；无法解析返回 None。"""
    s = s.strip()
    try:
        return float(s)
    except ValueError:
        pass
    m = _NO_E_EXP_RE.match(s)
    if m:
        try:
            return float(f"{m.group(1)}E{m.group(2)}")  # 含 inf/-inf（溢出）
        except ValueError:
            pass
    return None


def check_temperature_structure(path, teff, max_factor=3.0, temp_floor=10.0, temp_ceiling=1e8):
    """复刻 conv_check.rs:463。返回 (temps, violations) 或 (None, None)=无法解析。"""
    try:
        with open(path, "r", errors="replace") as f:
            lines = f.readlines()
    except OSError:
        return None, None

    if not lines:
        return None, None

    # 1. 第1行 → nd, numpar
    tokens = lines[0].split()
    if len(tokens) < 2:
        return None, None
    try:
        nd = int(float(tokens[0]))
        numpar = int(float(tokens[1]))
    except ValueError:
        return None, None
    if nd <= 0 or numpar <= 0:
        return None, None

    # 2. 行数计算
    dm_lines = (nd + 5) // 6        # FORMAT 502: 6/行
    block_lines = (numpar + 4) // 5  # FORMAT 503: 5/行

    # 3. 跳过 DM 列阵
    idx = 1  # 第1行已读
    for _ in range(dm_lines):
        if idx >= len(lines):
            return None, None  # 截断 → 跳过
        idx += 1

    # 4. nd 个深度块，每块第一行第一token = TEMP
    temps = []
    for _ in range(nd):
        if idx >= len(lines):
            break  # 截断：已解析的深度保留
        first_line_tokens = lines[idx].split()
        if not first_line_tokens:
            idx += block_lines
            continue
        t = parse_fortran_float(first_line_tokens[0])
        temps.append(t if t is not None else float("nan"))
        idx += block_lines  # 跳到下一块（含本块剩余 block_lines-1 行）

    if not temps:
        return None, None

    # 5. 三项判据
    violations = []
    surface_ratio = temps[0] / teff if teff > 0 else float("inf")
    if surface_ratio > max_factor:
        violations.append(
            f"表层 T={temps[0]:.4g} > {max_factor}×Teff({teff})，ratio={surface_ratio:.4g}"
        )
    for i, t in enumerate(temps):
        if not math.isfinite(t):
            violations.append(f"深度 {i+1} T={t}（非有限，NaN/Inf）")
            break
        if t < temp_floor or t > temp_ceiling:
            violations.append(
                f"深度 {i+1} T={t:.4g} 越界 [{temp_floor}, {temp_ceiling}]"
            )
            break
    return temps, violations


def main():
    ap = argparse.ArgumentParser(description="fort.7 温度结构物理判定（DCTS 权威判据）")
    ap.add_argument("fort7", help="fort.7 或 {stage}.7 路径")
    ap.add_argument("--teff", type=float, required=True, help="Teff (K)")
    ap.add_argument("--max-factor", type=float, default=3.0, help="表层/Teff 上限（默认 3.0）")
    ap.add_argument("--floor", type=float, default=10.0, dest="temp_floor", help="温度下限 K（默认 10）")
    ap.add_argument("--ceiling", type=float, default=1e8, dest="temp_ceiling", help="温度上限 K（默认 1e8）")
    args = ap.parse_args()

    import os
    if not os.path.exists(args.fort7):
        print(f"⚠ 文件不存在: {args.fort7}（跳过，不判失败）")
        return 2

    temps, violations = check_temperature_structure(
        args.fort7, args.teff, args.max_factor, args.temp_floor, args.temp_ceiling
    )
    if temps is None:
        print(f"⚠ 无法解析 {args.fort7}（缺失/截断/格式错）—— 跳过，不判失败")
        return 2

    surface_ratio = temps[0] / args.teff if args.teff > 0 else float("inf")
    print(f"=== 温度结构判定: {args.fort7} (Teff={args.teff}) ===")
    print(f"  深度数:       {len(temps)}")
    print(f"  表层 T:       {temps[0]:.4g} K  ({surface_ratio:.4g}×Teff，Eddington 参考 0.84)")
    print(f"  底层 T:       {temps[-1]:.4g} K")
    nan_count = sum(1 for t in temps if not math.isfinite(t))
    print(f"  NaN/Inf 深度: {nan_count}")

    if not violations:
        print(f"✅ 物理（三项全过：表层<{args.max_factor}×Teff、各深度∈[{args.temp_floor},{args.temp_ceiling}]、无 NaN）")
        return 0
    print("❌ 不物理 —— 违规:")
    for v in violations:
        print(f"  - {v}")
    print("\n参考 references/physics-checks.md 的判据说明与修正历史。")
    return 1


if __name__ == "__main__":
    sys.exit(main())