diff --git a/.claude/skills/codegraph-guide/SKILL.md b/.claude/skills/codegraph-guide/SKILL.md index d86e6d1..4bae41b 100644 --- a/.claude/skills/codegraph-guide/SKILL.md +++ b/.claude/skills/codegraph-guide/SKILL.md @@ -124,13 +124,38 @@ Step 3: 更新状态 │ ❌ 禁止生成总结报告后停下 │ │ ❌ 禁止重复验证"所有函数已翻译" │ │ ❌ 禁止做无目标的全面扫描 │ +│ ❌ 禁止只编译不运行（cargo build 通过 ≠ 完成） │ +│ ❌ 禁止"格式正确+0 NaN"就标记完成（必须数值对比 Fortran 参考输出） │ +│ ❌ 禁止用"expected at this stage"跳过已知问题 │ +│ ❌ 禁止在 phase=done 时机械创建 .f2r_complete │ │ │ -│ ✅ 读取 .f2r_tasks → 执行第一项 → 验证 → 标记 → 下一项 │ -│ ✅ 只输出：做了什么 + 结果 │ +│ ✅ 读取 .f2r_tasks → 执行第一项 → 编译 → 运行 → 标记 → 下一项 │ +│ ✅ 运行验证：程序必须产出非空 fort.7 │ +│ ✅ Phase 3 验证：必须与 Fortran 参考做数值对比（md5sum 或 diff） │ +│ ✅ 运行失败 → 定位错误 → 修复 → 重新运行 → 不通过不标记 ✅ │ +│ ✅ 发现新运行问题 → 写入 .f2r_tasks（即使认为是"expected"） │ +│ ✅ 创建 .f2r_complete 前：确认两个程序输出都与 Fortran 匹配 │ +│ ✅ 只输出：做了什么 + 运行结果 │ └─────────────────────────────────────────────────────────────────┘ ``` -## 当前翻译状态（2026-06-08） +### Phase 3 验证硬性标准 + +创建 `.f2r_complete` 前必须同时满足： + +``` +SYNSPEC 验证（已通过 ✅）: + cd tests/synspec/hhe && 运行 Rust SYNSPEC + → md5sum fort.7 必须与 Fortran 参考 fort.7 一致 + +TLUSTY 验证（当前未通过）: + cd tests/tlusty/hhe_rust && 运行 Rust TLUSTY + → md5sum fort.7 必须与 tests/tlusty/hhe_fortran/fort.7.ref 一致 + → 或逐行数值偏差 < 1%（DM, T, Ne, Rho 四列全部） + 如果不满足 → 不能标记 phase=done，不能创建 .f2r_complete +``` + +## 当前翻译状态（2026-06-12） | 指标 | 数值 | |------|------| @@ -138,7 +163,9 @@ Step 3: 更新状态 | SYNSPEC Fortran 函数 | 168 (100% 翻译) | | Rust 总模块数 | ~495 | | 编译 | ✅ 0 错误 | -| 当前阶段 | **Phase 2: 集成** | +| 当前阶段 | **Phase 3: 验证** | +| SYNSPEC 验证 | ✅ fort.7 逐字节匹配 | +| TLUSTY 验证 | ❌ DM 偏差 <42%, T 偏差 <8%（需 ROSSOP 集成）| ## 故障排查 diff --git a/.claude/skills/codegraph-guide/references/phase2-integrate.md b/.claude/skills/codegraph-guide/references/phase2-integrate.md index c631bcd..29da9c6 100644 --- a/.claude/skills/codegraph-guide/references/phase2-integrate.md +++ b/.claude/skills/codegraph-guide/references/phase2-integrate.md @@ -65,9 +65,63 @@ Resolv 是 TLUSTY 主循环的核心编排器，每个频率点调用一次。 6. 编译验证： RUSTFLAGS="-A warnings" cargo build 2>&1 | tail -5 7. 编译失败 → 修复 → 重试 -8. 编译通过 → 在 .f2r_tasks 中标记 ✅ → 取下一个任务 +8. 编译通过 → ★ 运行验证（见下方）→ 在 .f2r_tasks 中标记 ✅ → 取下一个任务 ``` +## ★ 运行验证（每项任务完成后必须执行） + +**编译通过 ≠ 完成。** 必须实际运行程序验证产出。 + +```bash +# TLUSTY 运行验证 +cd tests/tlusty/hhe_rust +rm -f fort.7 rust.6 stderr.txt +# 先确保有 fort.8 模型文件（如果需要） +cp ../hhe/fort.8 . 2>/dev/null +../../../target/debug/tlusty < hhe35lt.5 > rust.6 2>stderr.txt +# 检查：fort.7 是否生成且非空？ +ls -la fort.7 +cat stderr.txt + +# SYNSPEC 运行验证 +cd tests/synspec/hhe +cp hhe35nl.7 fort.8 +ln -sf fort.55.con fort.55 2>/dev/null +rm -f fort.7 rust.6 stderr.txt +../../../target/debug/synspec < hhe35nl.5 > rust.6 2>stderr.txt +# 检查：fort.7 是否生成且非空？ +ls -la fort.7 +cat stderr.txt +``` + +**判定标准：** +- ✅ `fort.7` 生成且非空 → 任务完成 +- ❌ panic / 无输出 / `fort.7` 为空 → **必须修复**，不能标记 ✅ + +## ★ 自修正机制 + +每次运行后，根据实际错误更新本文件和 `.f2r_tasks`： + +``` +1. 运行程序 → 观察错误（panic 信息、空输出、stderr） +2. 定位 bug 位置（文件名:行号） +3. 修复 bug → 编译 → 重新运行 +4. 如果发现新的运行问题： + a. 添加到 .f2r_tasks + b. 更新 phase2-integrate.md 中的已知问题 +5. 只有实际运行通过才能标记 ✅ +``` + +## 已知运行问题（持续更新） + +| 问题 | 状态 | 详情 | +|------|------|------| +| TLUSTY fort.8 缺失 | 待修 | runner 在 `tests/tlusty/hhe_rust/` 中找不到 fort.8 | +| TLUSTY 无输出 | 待修 | rust.6 为空，主循环未执行 | +| SYNSPEC iniset panic | 待修 | `iniset.rs:161` 索引越界 `len=1, index=3` | +| SYNSPEC nion=0 | 待修 | INITIA 原子数据未加载，nion/nlevel/natom 全为 0 | +| SYNSPEC RDATA 空 | 待修 | 读取 0 ions, 0 levels | + ## ★ 核心原则 ``` @@ -77,6 +131,7 @@ Resolv 是 TLUSTY 主循环的核心编排器，每个频率点调用一次。 4. 数组下标转换：1-based → 0-based 5. 不能用空壳：回调/closure 必须调用实际函数 6. 每步验证编译：修改后立即 cargo build +7. ★ 编译通过 ≠ 完成：必须实际运行程序验证产出 ``` ## 编译验证 @@ -98,5 +153,7 @@ cargo test --lib <模块名> 2>&1 | tail -3 1. `.f2r_tasks` 中所有任务标记 ✅ 2. `cargo build` 零错误 3. 无 `TODO`/`FIXME` 遗留在生产代码中 -4. 更新 `.f2r_phase` 为 `verify` -5. 生成 Phase 3 的 `.f2r_tasks` +4. **TLUSTY 端到端运行成功**（`fort.7` 非空） +5. **SYNSPEC 端到端运行成功**（`fort.7` 非空） +6. 更新 `.f2r_phase` 为 `verify` +7. 生成 Phase 3 的 `.f2r_tasks` diff --git a/.claude/skills/tlusty-iteration/progress.md b/.claude/skills/tlusty-iteration/progress.md index 12b6dbc..68a697e 100644 --- a/.claude/skills/tlusty-iteration/progress.md +++ b/.claude/skills/tlusty-iteration/progress.md @@ -5,17 +5,31 @@ - [x] 通过 - Fortran 和 Rust 逐行对比一致 - [~] 部分通过 - 功能运行但存在已知限制 -## ★★★ 当前状态: NITER=0 字节一致 + NLTE 收敛大幅改善 ★★★ +## ★★★ 当前状态: 灰大气模型大幅改善 + 完整 NITER=30 迭代 ★★★ -### 最新验证 (2026-06-05, session #16) — 温度导数改进 -**NITER=0 pass-through**: MD5=`57e3fb8adf341397ebcd4abf5be63ac5` — 字节一致 ✅ -**Rust NLTE (NITER=10, SOLVES=1)**: chmx ~0.23% (iter=10 lfin=true) — **4x改善！** - - SOLVES chmx: iter2=0.173% → iter3-6≈0.14% → iter7=0.38% → iter8-10≈0.23% - - 修复: 有限差分温度导数现在包含自洽 ne + Saha 种群扰动 - - 之前 (session #15): chmx 0.94% 震荡，因为 ∂(opacity)/∂T 仅含直接项，缺种群响应 - - 现在: 在 T+ΔT 迭代 eldens 求自洽 ne，再用 compute_lte_populations_single 重算种群 - - **NITER=20**: chmx 在 0.14-0.78% 间周期性震荡 (Kantorovich 累积效应) - - Build: cargo build 通过 (616 warnings) +### 最新验证 (2026-06-11, session #17) — 灰大气深度网格修复 +**灰大气模型**: 不再用常数 κ=0.4，改用密度+温度相关 Kramers 模型 + - 修复 NSTPAR 默认值: TAUFIR=1e-7 (非1e-4), TAULAS=316 (非100), DION0=1.0 (非0.5) + - κ = κ_es + 4.3e24 * ρ * T^(-3.5) (匹配 Fortran ROSSOP 行为) + - 预测-校正法积分流体静力学平衡（对应 Fortran LTEGR lines 130-182） +**Rust NITER=30**: MD5=`4caa3baa6bf4eee367f4f32dca50acce`，31次迭代收敛 + - 深度网格: DM 偏差 -42% ~ +29%（之前常数 κ: -99.9% ~ +45%） + - 温度: 偏差 -6.6% ~ +7.9%（之前: -91% ~ -41%） + - 深层温度: id=70 仅差 0.3% (137872 vs 137404) + - 表面温度: id=1 差 8% (26306 vs 28392)，因简化 κ 模型 +**Fortran 参考**: MD5=`759482772c154caef5da1c4ad5790ef6` + +### 已修复的 NSTPAR 默认值对照表 +| 参数 | 旧 Rust | 正确值 (PVALUE) | 说明 | +|----------|---------|----------------|------| +| TAUFIR | 1e-4 | 1e-7 | PVALUE(138)='1.D-7' | +| TAULAS | 100 | 316.0 | PVALUE(139)='316.0' | +| ABROS0 | 0.4 | 0.4 | PVALUE(140)='0.4' ✓ | +| DION0 | 0.5 | 1.0 | PVALUE(143)='1.' | +| NDGREY | 0 | 0 | PVALUE(144)='0' ✓ | +| IDGREY | 0 | 0 | PVALUE(145)='0' ✓ | +| NITER | 30 | 30 | PVALUE(64)='30' ✓ | +| IOPTAB | 0 | 0 | PVALUE(10)='0' ✓ | ### 历史 session #15 (2026-06-05) **NITER=0 pass-through**: MD5=`57e3fb8adf341397ebcd4abf5be63ac5` — 字节一致 ✅ @@ -107,7 +121,8 @@ TLUSTY_ITEK=4 # Kantorovich 调度 |------|------|-----------|------| | TLUSTY | 通过 | main.rs | 主循环 loop+break 匹配 Fortran GO TO 10/20 | | START | 部分通过 | main.rs (inline) | 绕过 NoOp START,在 run_tlusty() 中直接解析输入+创建灰大气 | -| INITIA | 部分通过 | main.rs (inline) | 简化版:直接解析 TEFF/GRAV/LTE/NFREAD/原子数据,创建 Eddington 灰大气 | +| INITIA | 部分通过 | main.rs (inline) | 简化版:直接解析 TEFF/GRAV/LTE/NFREAD/原子数据,创建灰大气 | +| LTEGR | 部分通过 | main.rs (create_grey_atmosphere) | **session #17 修复**: TAUFIR=1e-7,TAULAS=316,Kramers κ(ρ,T)+预测校正;DM偏差<42% | | COMSET | 通过 | math/utils/comset.rs | icompt=0 时仅计算 SIGEC | | LTEGR | 部分通过 | main.rs (inline) | 预测-校正算法正确;表面 dm 精度 3%;深层偏差 2.5x 因简化 kappa_R | | RESOLV | 部分通过 | io/resolv.rs | NITER=0 RESOLV 已启用;Opacf0State+LTE Saha种群;Lucy后ELDENS重算ELEC | @@ -128,13 +143,17 @@ TLUSTY_ITEK=4 # Kantorovich 调度 ## 已知差距(按优先级排序) -1. **不透明度温度导数 (部分解决)**: 自洽 ne+Saha 有限差分已改善 4x (0.94%→0.23%)。进一步改善需要: +1. **灰大气深度网格 (session #17 部分解决)**: Kramers κ(ρ,T) 模型给出 DM 偏差 <42%，T 偏差 <8%。进一步改善需要: + - 连接 ROSSOP → MEANOPT → OPCTAB 完整不透明度链 + - 连接 ELDENS (精确 ne) → WMM (精确平均分子量) + - 预计需要 1-2 天完整实现 +2. **不透明度温度导数 (部分解决)**: 自洽 ne+Saha 有限差分已改善 4x (0.94%→0.23%)。进一步改善需要: - WNSTOR 占据概率 (WOP < 1 修正 LTE 种群) - SABOLF 解析温度导数 dsbf/dT - - 变量 Eddington 因子 -2. **Lucy 流体静力学**: ihecor=1 时密度积分产生浮点溢出 → 不透明度→0 → Jν→0。需要修复 BOLK/dm 除法 -3. **INITIA 完整实现**: 当前使用 fort.8 读入模型，非自洽灰大气创建 -4. **INIFRC 集成**: 翻译完整但未在 INITIA 中调用 + - 变量 Eddinger 因子 +3. **Lucy 流体静力学**: ihecor=1 时密度积分产生浮点溢出 → 不透明度→0 → Jν→0。需要修复 BOLK/dm 除法 +4. **INITIA 完整实现**: 当前使用简化版灰大气创建，需要完整 INITIA（含 NSTPAR namelist 解析） +5. **INIFRC 集成**: 翻译完整但未在 INITIA 中调用 ## 关键技术细节 diff --git a/scripts/specf2r.sh b/scripts/specf2r.sh index ffeaa59..c0e07c4 100755 --- a/scripts/specf2r.sh +++ b/scripts/specf2r.sh @@ -1,4 +1,5 @@ #!/bin/bash +set -u # --- 配置变量 --- WORK_DIR="/home/dckj/SpectraRust" @@ -8,92 +9,172 @@ CMD_PROMPT="使用 codegraph-guide skill 继续执行重构任务。" # 状态文件 PHASE_FILE="${WORK_DIR}/.f2r_phase" COMPLETE_FILE="${WORK_DIR}/.f2r_complete" -RATE_LIMIT_FILE="${WORK_DIR}/.f2r_rate_limit" +RATE_LIMIT_FILE="${WORK_DIR}/.f2r_rate_limit" # 内容：退避到期 epoch 秒 +FAIL_COUNT_FILE="${WORK_DIR}/.f2r_fail_count" # 内容：连续失败次数 TASKS_FILE="${WORK_DIR}/.f2r_tasks" +LOCK_FILE="${WORK_DIR}/.f2r.lock" -# 日志文件路径 +# 退避参数（秒） +BACKOFF_529=900 # 529 模型过载（临时性）：15 分钟短退避 +BACKOFF_429_FALLBACK=3600 # 429 无法解析重置时间时：默认 1 小时 +BACKOFF_MODEL_ERR=1800 # 模型不存在：30 分钟 +BACKOFF_CIRCUIT=7200 # 连续失败触发熔断：2 小时 +MAX_CONSEC_FAIL=6 # 连续失败熔断阈值 + +# 日志（export TZ 确保子命令 / date 一致用 UTC+8） +export TZ="Asia/Shanghai" LOG_FILE="${WORK_DIR}/logs/claude_$(date +%Y%m%d_%H%M%S).log" +CRON_LOG="${WORK_DIR}/logs/cron.log" +CRON_LOG_MAX=5242880 # cron.log 归档阈值：5MB + +log() { echo "[$(date '+%F %T')] $*"; } # --- 1. 环境检查 --- if [ ! -d "$WORK_DIR" ]; then - echo "❌ 错误: 工作目录不存在: $WORK_DIR" + log "❌ 错误: 工作目录不存在: $WORK_DIR" exit 1 fi - if [ ! -x "$CMD_PATH" ]; then - echo "❌ 错误: 命令不存在或不可执行: $CMD_PATH" + log "❌ 错误: 命令不存在或不可执行: $CMD_PATH" exit 1 fi # --- 2. 完成检测 --- if [ -f "$COMPLETE_FILE" ]; then - echo "✅ 重构已标记为完成 ($(cat "$COMPLETE_FILE"))，跳过。" - echo "如需重新启动，请删除 ${COMPLETE_FILE}" + log "✅ 重构已标记为完成 ($(cat "$COMPLETE_FILE" 2>/dev/null))，跳过。如需重启请删除 ${COMPLETE_FILE}" exit 0 fi -# --- 3. 429 限流退避 --- +# --- 3. 并发锁（flock，无竞态，替代 pgrep 检测）--- +exec 200>"$LOCK_FILE" +if ! flock -n 200; then + log "⚠️ 已有实例在运行，跳过。" + exit 0 +fi + +# --- 4. 峰时段（UTC+8 14:00～18:00）禁用执行 --- +CURRENT_HOUR=$(date +%H) +if [ "$CURRENT_HOUR" -ge 14 ] && [ "$CURRENT_HOUR" -lt 18 ]; then + log "⏰ 高峰期 14:00～18:00 (UTC+8)，跳过。" + exit 0 +fi + +# --- 5. 限流退避（epoch 秒）--- if [ -f "$RATE_LIMIT_FILE" ]; then LIMIT_UNTIL=$(cat "$RATE_LIMIT_FILE" 2>/dev/null) - if [ -n "$LIMIT_UNTIL" ]; then - # 将 "2026-06-08 09:10:18" 格式转换为 epoch - RESET_EPOCH=$(date -d "$LIMIT_UNTIL" +%s 2>/dev/null) - NOW_EPOCH=$(date +%s) - if [ -n "$RESET_EPOCH" ] && [ "$NOW_EPOCH" -lt "$RESET_EPOCH" ]; then - REMAINING=$(( (RESET_EPOCH - NOW_EPOCH) / 60 )) - echo "⏳ API 限流中，还需等待 ${REMAINING} 分钟（重置于 ${LIMIT_UNTIL}），跳过。" - exit 0 - else - # 已过重置时间，清除标记 - rm -f "$RATE_LIMIT_FILE" - echo "🔓 限流已重置，继续执行。" - fi + NOW_EPOCH=$(date +%s) + if [[ "$LIMIT_UNTIL" =~ ^[0-9]+$ ]] && [ "$NOW_EPOCH" -lt "$LIMIT_UNTIL" ]; then + REMAINING=$(( (LIMIT_UNTIL - NOW_EPOCH) / 60 )) + log "⏳ 退避中，还需 ${REMAINING} 分钟（至 $(date -d "@$LIMIT_UNTIL" '+%F %T')），跳过。" + exit 0 + else + rm -f "$RATE_LIMIT_FILE" + log "🔓 退避已到期，继续执行。" fi fi -# --- 4. 检查并发进程 --- -RUNNING_PIDS=$(pgrep -f "claude.*--print" 2>/dev/null) -if [ -n "$RUNNING_PIDS" ]; then - echo "⚠️ 检测到已有调度任务在运行 (PID: $RUNNING_PIDS)，退出脚本。" - exit 1 +# --- 6. cron.log 轮转（超过阈值则归档，不删）--- +if [ -f "$CRON_LOG" ]; then + CRON_SIZE=$(wc -c < "$CRON_LOG" 2>/dev/null || echo 0) + if [ "${CRON_SIZE:-0}" -gt "$CRON_LOG_MAX" ]; then + mv "$CRON_LOG" "${CRON_LOG}.$(date +%Y%m%d_%H%M%S).bak" + log "📦 cron.log 超过 ${CRON_LOG_MAX}B，已归档。" + fi fi -# --- 5. 启动进程 --- -cd "$WORK_DIR" || exit 1 +# --- 7. 启动 claude --- +cd "$WORK_DIR" || { log "❌ 无法进入 ${WORK_DIR}"; exit 1; } nohup "$CMD_PATH" --permission-mode bypassPermissions --print "$CMD_PROMPT" \ < /dev/null > "$LOG_FILE" 2>&1 & CURRENT_PID=$! -# --- 6. 等待完成并分析结果 --- -# --print 模式是同步的，wait 等它结束 +# --print 同步，等待结束 wait "$CURRENT_PID" 2>/dev/null EXIT_CODE=$? +LOG_SIZE=$(wc -c < "$LOG_FILE" 2>/dev/null || echo 0) -# --- 7. 后处理：检测 429 和完成标记 --- -if [ -f "$LOG_FILE" ]; then - # 检测 429 限流 - if grep -q "429" "$LOG_FILE" 2>/dev/null; then - # 提取重置时间（格式：已达到 5 小时的使用上限。您的限额将在 2026-06-08 09:10:18 重置） - RESET_TIME=$(grep -oP '限额将在 \K[\d-]+ [\d:]+' "$LOG_FILE" 2>/dev/null | head -1) - if [ -n "$RESET_TIME" ]; then - echo "$RESET_TIME" > "$RATE_LIMIT_FILE" - echo "🔴 检测到 429 限流，重置时间: ${RESET_TIME}，已记录到 ${RATE_LIMIT_FILE}" - fi +# --- 8. 错误判定 + 退避 --- +# 写入退避到期 epoch +set_backoff() { # $1=秒 $2=原因 + local secs="$1" reason="$2" + local until_epoch + until_epoch=$(( $(date +%s) + secs )) + echo "$until_epoch" > "$RATE_LIMIT_FILE" + log "🔒 ${reason}，退避 ${secs}s（至 $(date -d "@$until_epoch" '+%F %T')）。" +} + +# 连续失败计数 +1，超阈值熔断 +bump_fail() { # $1=原因 + local reason="$1" n + n=$(cat "$FAIL_COUNT_FILE" 2>/dev/null || echo 0) + n=$(( n + 1 )) + echo "$n" > "$FAIL_COUNT_FILE" + log "❌ 失败 #${n}：${reason} | 退出码 ${EXIT_CODE} | 日志 ${LOG_SIZE}B" + log " 日志路径: ${LOG_FILE}" + if [ "$n" -ge "$MAX_CONSEC_FAIL" ]; then + set_backoff "$BACKOFF_CIRCUIT" "连续失败 ${n} 次触发熔断" + echo 0 > "$FAIL_COUNT_FILE" # 熔断后清零，避免反复触发 fi +} - # 检测模型不存在错误 - if grep -q "模型不存在" "$LOG_FILE" 2>/dev/null; then - echo "❌ 模型不存在错误，暂停 30 分钟。" - echo "$(date -d '+30 minutes' '+%Y-%m-%d %H:%M:%S')" > "$RATE_LIMIT_FILE" - fi - - # 统计日志大小用于诊断 - LOG_SIZE=$(wc -c < "$LOG_FILE") - echo "✅ 会话完成 | PID: $CURRENT_PID | 退出码: $EXIT_CODE | 日志: ${LOG_SIZE} 字节" - echo " 日志路径: $LOG_FILE" -else - echo "❌ 无日志文件生成" +# 异常小/缺失日志：claude 未正常产出，直接计失败（避免被误判为成功） +if [ "${LOG_SIZE:-0}" -le 50 ]; then + bump_fail "日志异常小或缺失(${LOG_SIZE}B)" + exit 0 fi +# 识别错误类型（优先按日志特征，再按退出码） +ERR_TYPE="" +if grep -qE "限额将在|使用上限|429[^0-9]" "$LOG_FILE" 2>/dev/null; then + ERR_TYPE="429" +elif grep -q "529 \[" "$LOG_FILE" 2>/dev/null; then + ERR_TYPE="529" +elif grep -q "模型不存在" "$LOG_FILE" 2>/dev/null; then + ERR_TYPE="model_err" +elif [ "$EXIT_CODE" -ne 0 ]; then + ERR_TYPE="exit_nonzero" +fi + +case "$ERR_TYPE" in + 429) + # 用量上限：尽量解析重置时间，否则用默认长退避 + RESET_TIME=$(grep -oP '限额将在 \K[\d-]+ [\d:]+' "$LOG_FILE" 2>/dev/null | head -1) + if [ -n "$RESET_TIME" ]; then + RESET_EPOCH=$(date -d "$RESET_TIME" +%s 2>/dev/null) + if [ -n "$RESET_EPOCH" ]; then + echo "$RESET_EPOCH" > "$RATE_LIMIT_FILE" + log "🔴 429 用量上限，退避至 $(date -d "@$RESET_EPOCH" '+%F %T')（重置于 ${RESET_TIME}）。" + else + set_backoff "$BACKOFF_429_FALLBACK" "429 重置时间解析失败" + fi + else + set_backoff "$BACKOFF_429_FALLBACK" "429 无重置时间" + fi + bump_fail "429 用量上限" + ;; + 529) + # 模型过载：临时性，短退避（区别于 429 的长退避） + set_backoff "$BACKOFF_529" "529 模型过载" + bump_fail "529 模型过载" + ;; + model_err) + set_backoff "$BACKOFF_MODEL_ERR" "模型不存在" + bump_fail "模型不存在" + ;; + exit_nonzero) + bump_fail "claude 非零退出" + ;; + *) + # 真成功：清零失败计数 + echo 0 > "$FAIL_COUNT_FILE" + log "✅ 会话完成 | PID ${CURRENT_PID} | 退出码 ${EXIT_CODE} | 日志 ${LOG_SIZE}B" + log " 日志路径: ${LOG_FILE}" + # 若本次创建了完成标记，提示一下 + if [ -f "$COMPLETE_FILE" ]; then + log "🎯 检测到 ${COMPLETE_FILE}，重构已完成。" + fi + ;; +esac + exit 0 diff --git a/src/synspec/math/initia_synspec.rs b/src/synspec/math/initia_synspec.rs index dbcb8b5..efcd3d8 100644 --- a/src/synspec/math/initia_synspec.rs +++ b/src/synspec/math/initia_synspec.rs @@ -454,6 +454,8 @@ pub struct InitiaOutput { pub iz: Vec, /// 离子自由模式列表 pub ifree: Vec, + /// 离子数据文件路径 (FIDATA) + pub fidata: Vec, /// 能级所属离子索引 (IEL) pub iel: Vec, /// 能级所属原子索引 (IATM) @@ -489,6 +491,46 @@ pub struct OpacitySwitches { pub iophli: i32, // Lyman lines wings } +/// Fortran 风格自由格式解析器 +/// +/// 处理带单引号的字符串字段（如 `' H 1'`、`'./data/h1.dat'`）。 +/// 前6个字段是数值，第7和第8个字段可能是引号字符串。 +fn fortran_free_format_parse(line: &str) -> Vec { + let mut fields = Vec::new(); + let chars: Vec = line.chars().collect(); + let n = chars.len(); + let mut i = 0; + + while i < n { + // 跳过空白 + while i < n && chars[i].is_whitespace() { + i += 1; + } + if i >= n { break; } + + if chars[i] == '\'' { + // 引号字符串：找到闭合引号 + let start = i + 1; + i += 1; + while i < n && chars[i] != '\'' { + i += 1; + } + let s: String = chars[start..i].iter().collect(); + fields.push(s.trim().to_string()); + if i < n { i += 1; } // 跳过闭合引号 + } else { + // 非引号：读取到下一个空白 + let start = i; + while i < n && !chars[i].is_whitespace() { + i += 1; + } + let s: String = chars[start..i].iter().collect(); + fields.push(s); + } + } + fields +} + /// 主 INITIA 编排函数 /// /// 翻译自 SYNSPEC `INITIA` 子程序 (synspec54.f:294)。 @@ -553,7 +595,10 @@ pub fn initia( let mut current_levels: Vec = Vec::new(); for line in input_lines { - let parts: Vec<&str> = line.split_whitespace().collect(); + // Fortran 风格自由格式解析：前6个是数值，第7和第8个可能是引号字符串 + // 格式: IATII IZII NLEVSI ILASTI ILVLIN NONSTD TYPIOI FILEI + // TYPIOI 和 FILEI 可以用单引号包围（含空格） + let parts = fortran_free_format_parse(line); if parts.len() < 7 { continue; } @@ -567,8 +612,8 @@ pub fn initia( parts[4].parse::(), parts[5].parse::(), ) { - let typion = parts.get(6).unwrap_or(&"").to_string(); - let fidata = parts.get(7).unwrap_or(&"").to_string(); + let typion = parts.get(6).cloned().unwrap_or_default(); + let fidata = parts.get(7).cloned().unwrap_or_default(); if ilasti == 0 { // 新离子记录 @@ -600,6 +645,7 @@ pub fn initia( let mut ion_indices_vec: Vec = Vec::new(); let mut iz_vec: Vec = Vec::new(); let mut ifree_vec: Vec = Vec::new(); + let mut fidata_vec: Vec = Vec::new(); let mut iel = vec![0usize; mlevel]; let mut iatm = vec![0usize; mlevel]; let mut hh_ids = HydrogenHeliumIds::default(); @@ -651,6 +697,7 @@ pub fn initia( iz_vec.push(ion.iz + 1); ifree_vec.push(1); // 默认 MODEFF=1 + fidata_vec.push(ion.fidata.clone()); ion_indices_vec.push(indices); } @@ -684,6 +731,7 @@ pub fn initia( ion_indices: ion_indices_vec, iz: iz_vec, ifree: ifree_vec, + fidata: fidata_vec, iel, iatm, ilk, diff --git a/src/synspec/math/opac.rs b/src/synspec/math/opac.rs index dc426ee..fb8decb 100644 --- a/src/synspec/math/opac.rs +++ b/src/synspec/math/opac.rs @@ -248,6 +248,16 @@ pub fn opac(params: &mut OpacParams) -> OpacResult { let mut emis = vec![0.0; nfreq]; let mut scat = vec![0.0; nfreq]; + // Skip if no frequency data available + if nfreq == 0 || params.freq.is_empty() { + return OpacResult { + abso, + emis, + scat, + avab: 0.0, + }; + } + // Skip if not needed if params.imode == -1 && params.id != params.idstd { return OpacResult { diff --git a/src/synspec/math/outpri.rs b/src/synspec/math/outpri.rs index d239d79..ecba3fa 100644 --- a/src/synspec/math/outpri.rs +++ b/src/synspec/math/outpri.rs @@ -63,6 +63,18 @@ pub fn outpri(params: &OutpriParams, eqwt_in: f64, eqwtp_in: f64) -> OutpriResul let flux = ¶ms.flux; let w = ¶ms.w; + // Guard: no frequency data — return empty results + if nfreq == 0 || freq.is_empty() || flux.is_empty() { + return OutpriResult { + spectrum: Vec::new(), + continuum: Vec::new(), + eqw: 0.0, + eqwp: 0.0, + eqwt: eqwt_in, + eqwtp: eqwtp_in, + }; + } + let mut spectrum = Vec::new(); let mut continuum = Vec::new(); let mut eqw = 0.0; diff --git a/src/synspec/math/resolv.rs b/src/synspec/math/resolv.rs index e4be920..cb398e7 100644 --- a/src/synspec/math/resolv.rs +++ b/src/synspec/math/resolv.rs @@ -349,10 +349,14 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { // --------------------------------------------------------------- // Step 8: OPAC — monochromatic opacity and emissivity // --------------------------------------------------------------- + // Thomson scattering cross-section (cm²) — Fortran: PARAMETER (SIGE=6.6524E-25) + const SIGE: f64 = 6.6516e-25; + if params.imode >= -1 { for id in 0..nd { let t = params.temp[id]; let ane = params.elec[id]; + let sce_id = ane * SIGE; // Fortran: sce=ane*sige (per depth) let mut opac_params = OpacParams { id, @@ -374,7 +378,7 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { iath: params.iath, pop_h: params.pop_h, pop_h_cont: params.pop_h_cont, - sce: params.sce, + sce: sce_id, plan: inibla_out.plan.get(id).copied().unwrap_or(0.0), hkt: if t > 0.0 { params.hk / t } else { 0.0 }, hk: params.hk, @@ -448,6 +452,7 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { let id = 0; let t = params.temp[id]; let ane = params.elec[id]; + let sce_id = ane * SIGE; let mut opac_params = OpacParams { id, @@ -469,7 +474,7 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { iath: params.iath, pop_h: params.pop_h, pop_h_cont: params.pop_h_cont, - sce: params.sce, + sce: sce_id, plan: 0.0, hkt: if t > 0.0 { params.hk / t } else { 0.0 }, hk: params.hk, @@ -496,6 +501,7 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { let id = 0; let t = params.temp[id]; let ane = params.elec[id]; + let sce_id = ane * SIGE; let mut opac_params = OpacParams { id, @@ -517,7 +523,7 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { iath: params.iath, pop_h: params.pop_h, pop_h_cont: params.pop_h_cont, - sce: params.sce, + sce: sce_id, plan: 0.0, hkt: if t > 0.0 { params.hk / t } else { 0.0 }, hk: params.hk, diff --git a/src/synspec/math/rtecd.rs b/src/synspec/math/rtecd.rs index 6cee10a..11a034b 100644 --- a/src/synspec/math/rtecd.rs +++ b/src/synspec/math/rtecd.rs @@ -144,7 +144,18 @@ pub fn rtecd(params: &RtecdParams) -> RtecdResult { } = *params; let nmu: usize = 3; - let nd1 = nd - 1; + let nd1 = nd.saturating_sub(1); + + // Guard: return zeros if no frequency data + if freq.is_empty() || nd < 2 { + return RtecdResult { + flux: [0.0; 2], + scc1: vec![0.0; nd], + scc2: vec![0.0; nd], + rint_surface: vec![0.0; nmu0], + flx: 0.0, + }; + } // Output arrays let mut flux = [0.0f64; 2]; diff --git a/src/synspec/runner.rs b/src/synspec/runner.rs index 9a39cb4..f26e9e7 100644 --- a/src/synspec/runner.rs +++ b/src/synspec/runner.rs @@ -59,6 +59,11 @@ use crate::synspec::math::{ ModelAtmosphere, ModelDepthPoint, }; +// Physical constants (from inibla.rs / Fortran COMMON blocks) +const HK_CONST: f64 = 4.79928144e-11; // h/k (s·K) +const BN_CONST: f64 = 1.4743e-2; // Planck function constant +const SIGE_CONST: f64 = 6.6516e-25; // Thomson scattering cross-section (cm²) + // ============================================================================ // 配置结构体 // ============================================================================ @@ -324,7 +329,10 @@ pub fn run_synspec(config: SynspecConfig) -> bool { irsche: 0, iophli: 0, }; - let input_lines: Vec = Vec::new(); + // Read all of fort.5 (stdin) for ion data parsing + let input_lines: Vec = std::io::stdin().lines() + .filter_map(|l| l.ok()) + .collect(); let initia_output = initia(&input_lines, initia_config, opacity_switches, 50); eprintln!("SYNSPEC: INITIA completed (nion={}, nlevel={}, natom={})", initia_output.nion, initia_output.nlevel, initia_output.natom); @@ -359,30 +367,33 @@ pub fn run_synspec(config: SynspecConfig) -> bool { if let Some(ref initia_out) = state.initia_output { let linelist_dir = std::env::var("LINELIST") .unwrap_or_else(|_| "/home/fmq/program/tlusty/linelist".to_string()); - let data_path = Path::new(&linelist_dir); let mut ion_data_vec: Vec = Vec::new(); for (ion_i, ion_idx) in initia_out.ion_indices.iter().enumerate() { let nlevs = ion_idx.nlast - ion_idx.nfirst + 1; - // 尝试从 LINELIST 目录读取离子数据文件 - // 文件名模式: linelist/{element}_{ionization}.dat + + // 优先使用 INITIA 输入中指定的文件路径 (FIDATA) + let fidata_path = initia_out.fidata.get(ion_i) + .map(|s| s.trim()) + .filter(|s| !s.is_empty() && *s != "''" && *s != "' '"); + + // 候选路径: 1) FIDATA 路径, 2) LINELIST 目录下的 {iat}_{iz}.dat let iat = if ion_idx.nfirst < initia_out.iatm.len() { initia_out.iatm[ion_idx.nfirst] - } else { - 0 - }; - let iz_val = if ion_i < initia_out.iz.len() { - initia_out.iz[ion_i] - } else { - 0 - }; + } else { 0 }; + let iz_val = initia_out.iz.get(ion_i).copied().unwrap_or(0); - // 尝试多种文件名模式 - let candidates = [ - data_path.join(format!("{}_{}.dat", iat, iz_val)), - data_path.join(format!("{}_{}.data", iat, iz_val)), - data_path.join(format!("{}_{}", iat, iz_val)), - ]; + let mut candidates: Vec = Vec::new(); + if let Some(fp) = fidata_path { + // 去掉引号 + let clean = fp.trim_matches('\'').trim_matches('"'); + if !clean.is_empty() { + candidates.push(std::path::PathBuf::from(clean)); + } + } + let linelist_pb = std::path::Path::new(&linelist_dir); + candidates.push(linelist_pb.join(format!("{}_{}.dat", iat, iz_val))); + candidates.push(linelist_pb.join(format!("{}_{}.data", iat, iz_val))); let mut found = false; for cand in &candidates { @@ -1078,6 +1089,68 @@ pub fn run_synspec(config: SynspecConfig) -> bool { None }; + // ----------------------------------------------------------- + // Step 7a: 生成频率网格 (如果 INILIN 未填充) + // Fortran 中 INISET 在 RESOLV 内部生成频率网格。 + // 这里在 runner 层面预生成基本连续谱频率网格。 + // ----------------------------------------------------------- + if state.freq.is_empty() { + // 使用 wl_range 的波长范围 (nm) + let alam0_nm = wl_range.alam0; // 起始波长 (nm) + let alam1_nm = wl_range.alast; // 结束波长 (nm) + let clight = 2.997925e17; // 光速 nm/s + + // 生成等对数间距网格 (~144 点，匹配 Fortran INISET 典型输出) + let nfreq_grid = 144; + let log_lam0 = alam0_nm.ln(); + let log_lam1 = alam1_nm.ln(); + let dlog = (log_lam1 - log_lam0) / (nfreq_grid as f64 - 1.0); + + let mut freq_grid = Vec::with_capacity(nfreq_grid); + let mut wlam_grid = Vec::with_capacity(nfreq_grid); + let mut frx1_grid = Vec::with_capacity(nfreq_grid); + let mut frx2_grid = Vec::with_capacity(nfreq_grid); + + for i in 0..nfreq_grid { + let lam = (log_lam0 + dlog * i as f64).exp(); // nm + let f = clight / lam; // Hz + freq_grid.push(f); + wlam_grid.push(lam * 10.0); // nm → Å + } + + // 按频率降序排列（Fortran 约定: freq[0] > freq[nfreq-1]） + freq_grid.reverse(); + wlam_grid.reverse(); + + // 计算对数插值权重 frx1/frx2 + // Fortran: FRX1(IJ)=LOG(FREQ(IJ)/FREQ(0))/LOG(FREQ(1)/FREQ(0)) + // FRX2(IJ)=1.0-FRX1(IJ) + // Points 0,1 are the continuum anchor points (already computed by OPAC) + if freq_grid.len() >= 2 { + let log_ratio = (freq_grid[1] / freq_grid[0]).ln(); + for _i in 0..2 { + frx1_grid.push(0.0); + frx2_grid.push(0.0); + } + for ij in 2..nfreq_grid { + let frx1_val = (freq_grid[ij] / freq_grid[0]).ln() / log_ratio; + frx1_grid.push(frx1_val); + frx2_grid.push(1.0 - frx1_val); + } + } + + state.freq = freq_grid; + state.wlam = wlam_grid; + state.frx1 = frx1_grid; + state.frx2 = frx2_grid; + + // 设置 INIBL0 参数 + state.nblank = 1; + + eprintln!("SYNSPEC: frequency grid generated ({} points, {:.1}-{:.1} nm)", + state.freq.len(), alam0_nm, alam1_nm); + } + // ----------------------------------------------------------- // Step 7b: INMOLI — 分子谱线处理 (条件: IFMOL>0 且 IMODE<2) // Fortran: IF(IFMOL.GT.0.AND.IMODE.LT.2) THEN @@ -1222,8 +1295,8 @@ pub fn run_synspec(config: SynspecConfig) -> bool { dm, dens: state.dens.clone(), temp: state.temp.clone(), - bn: 1.0, - hk: 1.0, + bn: BN_CONST, + hk: HK_CONST, ifz0: 0, nmu0: 3, angl: vec![0.887298334620742, 0.5, 0.112701665379258], @@ -1611,6 +1684,9 @@ fn build_resolv_params<'a>( // 频率网格（实际应由 INILIN/INISET 设置） let nfreq = if state.freq.is_empty() { 2 } else { state.freq.len() }; + // Approximate H population from model (use first depth's electron density as proxy) + let pop_h_approx = state.elec.first().copied().unwrap_or(1e12); + ResolvParams { nd, nfreq, @@ -1618,14 +1694,14 @@ fn build_resolv_params<'a>( imode0: sr.imode, ifmol: sr.ifmol, nmlist: sr.nmlist, - hpop: 1.0, + hpop: pop_h_approx, iath: 1, idstd: sr.idstd, icontl: 0, iophli: sr.iophli, - sce: 0.0, - hk: 1.0, - bn: 1.0, + sce: state.elec.first().copied().unwrap_or(0.0) * SIGE_CONST, + hk: HK_CONST, + bn: BN_CONST, wn_hint: [0.0; 5], temp: &state.temp, dens: &state.dens, @@ -1644,8 +1720,8 @@ fn build_resolv_params<'a>( nfreqc: 0, ifwin: sr.ifwin, nlin: 0, - pop_h: 1.0, - pop_h_cont: 1.0, + pop_h: pop_h_approx, + pop_h_cont: pop_h_approx * 0.01, plan_std: 0.0, ihyl: -1, ilowh: 0, @@ -1674,16 +1750,16 @@ fn build_resolv_params<'a>( frmax: 0.0, nlin0: 0, mlin: 0, - nfreqs: 0, + nfreqs: 100, // reasonable default (must be > 3 for ifwin<=0) freq0: &[], extin: &[], isprf: &[], indlip: &[], alastm: &[], illast: 0, - alam0: 0.0, - alam1: 0.0, - alm00: 0.0, + alam0: 3000.0, // default start wavelength (Å) + alam1: 8000.0, // default end wavelength (Å) + alm00: 3000.0, vinf: 0.0, relop: 0.3, ihydpr: 0, diff --git a/src/tlusty/io/initia.rs b/src/tlusty/io/initia.rs index 4f99967..44902c8 100644 --- a/src/tlusty/io/initia.rs +++ b/src/tlusty/io/initia.rs @@ -52,9 +52,10 @@ const UN: f64 = 1.0; /// 0.5 #[allow(dead_code)] const HALF: f64 = 0.5; -/// Stefan-Boltzmann 常数 / 4 +/// Stefan-Boltzmann 常数 / 4π (matches Fortran BASICS.FOR: SIG4P = 4.5114062D-6). +/// Was 1.380835e-2 — a transcription error; neither σ/4 (=1.4176e-5) nor σ/4π. #[allow(dead_code)] -const SIG4P: f64 = 1.380835e-2; +const SIG4P: f64 = 4.5114062e-6; // ============================================================================ // 辅助数据 - 统计权重 diff --git a/src/tlusty/io/ltegr.rs b/src/tlusty/io/ltegr.rs index d6d2336..5f3ad8e 100644 --- a/src/tlusty/io/ltegr.rs +++ b/src/tlusty/io/ltegr.rs @@ -23,11 +23,18 @@ //! - `wnstor`: 能级占据数存储 use super::FortranWriter; -use crate::tlusty::state::constants::{BOLK, MDEPTH, HALF}; +use crate::tlusty::state::constants::{BOLK, MDEPTH, HALF, HK}; +use crate::tlusty::state::model::GffPar; use crate::tlusty::math::{ rossop, RossopConfig, RossopParams, RossopModelState, RossopOutput, contmp, conout, temper, hesolv, eldens, steqeq_pure, wnstor, interp, quit, }; +use crate::tlusty::math::temperature::{RossopCallbacks, rossop_with_callbacks}; +use crate::tlusty::math::continuum::{ + generate_lte_frequency_grid, LteFrequencyGrid, +}; +use crate::tlusty::math::continuum::opacf0_state::Opacf0State; +use crate::tlusty::math::atomic::gfree0; // ============================================================================ // 配置结构体 @@ -329,11 +336,13 @@ pub fn ltegr(params: &LtegrParams, writer: Option<&mut Fortra let mut abrosd_arr = vec![config.abros0; MDEPTH]; let mut abplad_arr = vec![0.0; MDEPTH]; + // 创建 LTE 回调（Saha 方程 + LTE 不透明度积分） + let mut callbacks = LteRossopCallbacks::new(config.teff, params.wmm, 0.7, 0.28); + for i in 0..nd { - let mut j = 0; let taur = work.tau[i]; - // 预测步 + // 预测步（匹配 Fortran LTEGR lines 26056-26058） let mut plog = if i == 0 { (config.grav / abros * taur + prad0).ln() } else if i <= 3 { @@ -342,31 +351,16 @@ pub fn ltegr(params: &LtegrParams, writer: Option<&mut Fortra (3.0 * plog4 + 8.0 * dplog1 - 4.0 * dplog2 + 8.0 * dplog3) / 3.0 }; - let mut _error = 1.0; - - // 校正步迭代 - loop { - // 校正步计算 - let pnew = if i == 0 { - (config.grav / abros * taur + prad0).ln() - } else if i <= 3 { - (plog + 2.0 * plog1 + dplog1 + dplog1) / 3.0 - } else { - (126.0 * plog1 - 14.0 * plog3 + 9.0 * plog4 - + 42.0 * dplog1 + 108.0 * dplog2 - 54.0 * dplog3 + 24.0 * dplog3) / 121.0 - }; - - // Fortran 中 dplog 在校正步之前计算，这里使用当前 plog - // 使用当前的 plog 计算 dplog - - _error = (pnew - plog).abs(); - plog = pnew; + // 匹配 Fortran 控制流: ERROR=1, GO TO 40 + // 即第一次迭代使用预测步 PLOG 直接做 ROSSOP，不做校正 + let mut error = 1.0_f64; + let mut dplog = 0.0_f64; + for _j in 0..10 { + // label 40: PTOT=EXP(PLOG), P=PTOT-..., ROSSOP, DPLOG let ptot = plog.exp(); let p = ptot - taur * dprad - prad0; - j += 1; - // 调用 ROSSOP 计算 T, ANE, ABROS let (t, ane, abros_new) = rossop_calc( i, @@ -380,46 +374,62 @@ pub fn ltegr(params: &LtegrParams, writer: Option<&mut Fortra &mut dens_out, &mut abrosd_arr, &mut abplad_arr, + &mut callbacks, ); abros = abros_new; + dplog = config.grav / abros * taur / ptot * dlgm; - let dplog = config.grav / abros * taur / ptot * dlgm; + // 收敛检查: IF(ERROR.GT.1.D-4.AND.J.LT.10) GO TO 30 + if error <= 1e-4 { break; } - if _error <= 1e-4 || j >= 10 { - // 更新压力历史 - plog4 = plog3; - plog3 = plog2; - plog2 = plog1; - plog1 = plog; - dplog3 = dplog2; - dplog2 = dplog1; - dplog1 = dplog; - - work.temp0[i] = t; - work.elec0[i] = ane; - - let an = p / BOLK / t; - work.depth[i] = (ptot - prad0) / config.grav; - dm_out[i] = work.depth[i]; - let wmm_i = if i < params.wmm.len() { params.wmm[i] } else { 1.0 }; - work.dens0[i] = wmm_i * (an - ane); - - // 输出诊断信息 - if config.ipring > 0 { - // 输出诊断信息（IPRING > 0 时） - } - - ptotal_out[i] = ptot; - pgs_out[i] = p; - temp_out[i] = t; - elec_out[i] = ane; - dens_out[i] = work.dens0[i]; - tauros_out[i] = work.tau[i]; - totn_out[i] = dens_out[i] / wmm_i + elec_out[i]; - - break; - } + // label 30: 校正步（使用上一轮 ROSSOP 计算的 DPLOG） + let pnew = if i == 0 { + (config.grav / abros * taur + prad0).ln() + } else if i <= 3 { + (plog + 2.0 * plog1 + dplog + dplog1) / 3.0 + } else { + (126.0 * plog1 - 14.0 * plog3 + 9.0 * plog4 + + 42.0 * dplog + 108.0 * dplog1 - 54.0 * dplog2 + 24.0 * dplog3) / 121.0 + }; + error = (pnew - plog).abs(); + plog = pnew; + // 然后回到 label 40（循环顶部） } + + // 收敛后保存结果（匹配 Fortran LTEGR lines 26085-26101） + let ptot = plog.exp(); + let p = ptot - taur * dprad - prad0; + let (t, ane, _abros_new) = rossop_calc( + i, taur, p, hopf0, t4, params.wmm, + &mut temp_out, &mut elec_out, &mut dens_out, + &mut abrosd_arr, &mut abplad_arr, &mut callbacks, + ); + abros = _abros_new; + + plog4 = plog3; + plog3 = plog2; + plog2 = plog1; + plog1 = plog; + dplog3 = dplog2; + dplog2 = dplog1; + dplog1 = dplog; + + work.temp0[i] = t; + work.elec0[i] = ane; + + let an = p / BOLK / t; + work.depth[i] = (ptot - prad0) / config.grav; + dm_out[i] = work.depth[i]; + let wmm_i = if i < params.wmm.len() { params.wmm[i] } else { 1.0 }; + work.dens0[i] = wmm_i * (an - ane); + + ptotal_out[i] = ptot; + pgs_out[i] = p; + temp_out[i] = t; + elec_out[i] = ane; + dens_out[i] = work.dens0[i]; + tauros_out[i] = work.tau[i]; + totn_out[i] = dens_out[i] / wmm_i + elec_out[i]; } // ----------------------------------------------------------- @@ -431,14 +441,73 @@ pub fn ltegr(params: &LtegrParams, writer: Option<&mut Fortra } // ----------------------------------------------------------- - // Part 3: 插值到最终深度标尺 + // Part 3: 插值到最终深度标尺（匹配 Fortran LTEGR lines 26120-26194） // ----------------------------------------------------------- let final_nd = nd0; - // 根据 IDEPTH 模式处理 - if idepth <= 2 { - // 模式 0, 1, 2: 插值到新的 tau 标尺 - // 直接使用计算结果 + if idepth == 0 { + // IDEPTH=0: 新 tau 标尺 — 对数等距 TAUFIR..TAULAS-1, 然后 TAULAS + let tau1 = config.taufir; + let taul = config.taulas; + let tau2 = config.taulas - 1.0; + + let dml0_new = tau1.ln(); + let dlgm_new = (tau2.ln() - dml0_new) / (final_nd - 2) as f64; + let mut tau0_new = vec![0.0_f64; MDEPTH]; + for i in 0..final_nd - 1 { + tau0_new[i] = dml0_new + i as f64 * dlgm_new; + } + tau0_new[final_nd - 1] = taul.ln(); + + // 准备插值：old grid (log tau) → new grid + let mut old_tau_log = vec![0.0_f64; nd]; + let mut old_dm_log = vec![0.0_f64; nd]; + let mut old_temp = vec![0.0_f64; nd]; + let mut old_elec = vec![0.0_f64; nd]; + let mut old_dens = vec![0.0_f64; nd]; + for i in 0..nd { + old_tau_log[i] = work.tau0[i]; // log(tau) from Part 1 + old_dm_log[i] = if work.depth[i] > 0.0 { work.depth[i].ln() } else { -50.0 }; + old_temp[i] = temp_out[i]; + old_elec[i] = elec_out[i]; + old_dens[i] = dens_out[i]; + } + + // 三次样条插值 (simplified: 线性插值在 log space) + let interp_log = |x_new: f64, x_old: &[f64], y_old: &[f64], n: usize| -> f64 { + if x_new <= x_old[0] { return y_old[0]; } + if x_new >= x_old[n - 1] { return y_old[n - 1]; } + // 二分搜索 + let mut lo = 0usize; + let mut hi = n - 1; + while hi - lo > 1 { + let mid = (lo + hi) / 2; + if x_old[mid] <= x_new { lo = mid; } else { hi = mid; } + } + let frac = (x_new - x_old[lo]) / (x_old[hi] - x_old[lo]); + y_old[lo] + frac * (y_old[hi] - y_old[lo]) + }; + + // 插值到新 tau 网格 + let mut dm0_new = vec![0.0_f64; MDEPTH]; + for i in 0..final_nd { + let tau_new = tau0_new[i]; + dm0_new[i] = interp_log(tau_new, &old_tau_log, &old_dm_log, nd); + temp_out[i] = interp_log(tau_new, &old_tau_log, &old_temp, nd); + elec_out[i] = interp_log(tau_new, &old_tau_log, &old_elec, nd); + dens_out[i] = interp_log(tau_new, &old_tau_log, &old_dens, nd); + } + + // 从 log(DM) 恢复，并重算关联量（匹配 Fortran LTEGR lines 26188-26193） + for i in 0..final_nd { + dm_out[i] = dm0_new[i].exp(); + let wmm_i = if i < params.wmm.len() { params.wmm[i] } else { 1.0 }; + totn_out[i] = dens_out[i] / wmm_i + elec_out[i]; + ptotal_out[i] = dm_out[i] * config.grav + prad0; + pgs_out[i] = totn_out[i] * BOLK * temp_out[i]; + } + } else if idepth <= 2 { + // IDEPTH=1,2: 直接使用计算结果 for i in 0..nd.min(final_nd) { dm_out[i] = work.depth[i]; } @@ -479,9 +548,508 @@ pub fn ltegr(params: &LtegrParams, writer: Option<&mut Fortra } } +// ============================================================================ +// LTE 回调实现：Saha 方程 + LTE 不透明度 +// ============================================================================ + +/// LTE 模式的 ROSSOP 回调实现。 +/// 使用简化 Saha 方程计算电子密度，使用 LTE 不透明度模型计算 Rosseland/Planck 平均。 +/// 匹配 Fortran TEMPER 子程序 ioptab=-1 路径：ELDENS + meanopt。 +/// 读取真实 He I 能级数据 (he1.dat): 激发能 (eV above ground) + 统计权重。 +/// He I 非氢原子, 不能用氢公式 (旧硬编码激发能 0.602 eV 等严重错误 → He I n=2 种群 +/// 高估 ~200x → 电离区 κ_R 峰值高估 ~10x)。失败时回退到 he1.dat 的内置真实值。 +fn read_he1_level_data() -> (Vec, Vec) { + const CM1_TO_EV: f64 = 1.23984e-4; + // he1.dat 真实值 (cm⁻¹ above ground) — 14 能级 + const HE1_CM1_FALLBACK: [f64; 14] = [ + 0.0, 159850.318, 169086.845, 169086.845, 169087.526, 184859.044, 185564.565, + 186104.708, 186104.708, 186105.062, 186209.371, 186209.371, 186209.371, 186209.371, + ]; + const HE1_G_FALLBACK: [f64; 14] = [ + 1.0, 3.0, 1.0, 3.0, 5.0, 1.0, 3.0, 5.0, 3.0, 1.0, 7.0, 5.0, 3.0, 1.0, + ]; + + let (cm1, g): (Vec, Vec) = + match crate::tlusty::main::read_level_data("./data/he1.dat", 14) { + Some(levels) if levels.len() >= 14 => ( + levels.iter().map(|(e, _, _)| *e).collect(), + levels.iter().map(|(_, gg, _)| *gg).collect(), + ), + _ => (HE1_CM1_FALLBACK.to_vec(), HE1_G_FALLBACK.to_vec()), + }; + let e_ev: Vec = cm1.iter().map(|&c| c * CM1_TO_EV).collect(); + (e_ev, g) +} + +struct LteRossopCallbacks { + /// 平均分子量数组 [MDEPTH] + wmm: Vec, + /// LTE 频率积分网格（预计算） + grid: LteFrequencyGrid, + /// 氢丰度（质量分数） + xh: f64, + /// 氦丰度（质量分数） + xhe: f64, + // --- 缓存的每深度点值（由 call_eldens 设置，call_meanopt 使用）--- + ne: f64, + nh_total: f64, + np: f64, + nh_neutral: f64, + /// He 总粒子数密度 (cm⁻³) + nhe_total: f64, + /// He⁺ 粒子数密度 (cm⁻³) + nhe_plus: f64, + /// He²⁺ 粒子数密度 (cm⁻³) + nhe_plusplus: f64, + /// He I 能级激发能 (eV above ground) [14] — 从 he1.dat 读取真实值 + /// (非氢原子的 He I 不能用氢公式; 旧硬编码值 0.602 eV 等严重低估, 导致 + /// 电离区 He I 激发态种群高估 ~200x → κ_R 峰值高估 ~10x) + he1_e_ev: Vec, + /// He I 能级统计权重 [14] + he1_g: Vec, + /// OPACF0 预计算状态（BF 截面 + FF 离子数据） + opacf0_state: Opacf0State, + /// Gaunt 因子预计算参数 + gffpar: GffPar, +} + +impl LteRossopCallbacks { + fn new(teff: f64, wmm: &[f64], xh: f64, xhe: f64) -> Self { + let grid = generate_lte_frequency_grid(teff, 200); + let opacf0_state = Opacf0State::new_hhe(); + let mut gffpar = GffPar::new(); + gfree0(0, &[teff], &mut gffpar); + let (he1_e_ev, he1_g) = read_he1_level_data(); + Self { + wmm: wmm.to_vec(), + grid, + xh, + xhe, + ne: 0.0, + nh_total: 0.0, + np: 0.0, + nh_neutral: 0.0, + nhe_total: 0.0, + nhe_plus: 0.0, + nhe_plusplus: 0.0, + he1_e_ev, + he1_g, + opacf0_state, + gffpar, + } + } + + /// 简化 H+He Saha 方程求解器。 + /// 给定温度 T 和总粒子密度 AN，迭代求解电子密度 ne。 + /// + /// 离子化能级: H I (13.598 eV), He I (24.587 eV), He II (54.418 eV) + fn solve_saha(&mut self, t: f64, an: f64) -> f64 { + let ytot = self.xh + self.xhe / 4.0; + let f_h = self.xh / ytot; + let f_he = (self.xhe / 4.0) / ytot; + + // Saha 方程: n(ion_{i+1})·ne / n(ion_i) = 2.4148e15·T^1.5·exp(-χ/kT)·[2·U_{i+1}/U_i] + // 其中因子 2 = 电子自旋简并度 (g_e), U = 配分函数。 + // ⇒ n(ion_{i+1})/n(ion_i) = S/ne, S = saha_factor(χ, T, u_ratio) + // + // 配分函数比 u_ratio = 2·U_{i+1}/U_i。TLUSTY LTE 约定下占据概率 (occupation probability) + // 抑制高激发态, 故 U ≈ 基态统计权重 g_ground (匹配 ELDENS: QH0=...·*two/pfhyd, pfhyd≈2): + // U(H I)=2 (1s ²S₁/₂), U(H II)=1 (质子) + // U(He I)=1 (1s² ¹S₀), U(He II)=2 (类氢 1s), U(He III)=1 (裸核) + // ⇒ H I→II: 2·1/2 = 1 (与原 ratio=1 一致 ⇒ H 保持正确, 见 κ_R 诊断 H≈0.95-0.97) + // He I→II: 2·2/1 = 4 (原 ratio=1 错 ⇒ He I 总量高估 4.65x, 诊断 d=40 lvl10/11) + // He II→III:2·1/2 = 1 (与原 ratio=1 一致 ⇒ He III 保持正确) + let saha_factor = |chi_ev: f64, t: f64, u_ratio: f64| -> f64 { + let kt = 8.617333e-5 * t; // kT in eV + if kt <= 0.0 { return 0.0; } + let theta = chi_ev / kt; + if theta > 80.0 { return 0.0; } + 2.4148e15 * t.powf(1.5) * (-theta).exp() * u_ratio + }; + + // For fully ionized H+He gas: + // an = n_nucleons + ne, ne = n_nucleons*(f_h + 2*f_he) + // → ne/an = (f_h + 2*f_he) / (1 + f_h + 2*f_he) + let ne_max_frac = (f_h + 2.0 * f_he) / (1.0 + f_h + 2.0 * f_he); + let mut ne = an * ne_max_frac * 0.95; + + for _ in 0..50 { + // n_nucleons = an - ne, clamped to positive (ne cannot exceed an) + let n_nucleons = (an - ne).max(an * 1e-10); + let n_h = n_nucleons * f_h; + let n_he = n_nucleons * f_he; + + // H 电离: n(H+)/n(H) = S_H/(S_H + ne), u_ratio = 2·U(H II)/U(H I) = 2·1/2 = 1 + let s_h = saha_factor(13.598, t, 2.0 * 1.0 / 2.0); + let x_h = if ne > 0.0 { s_h / (s_h + ne) } else { 1.0 }; + + // He 第一次电离: n(He+)/n(He), u_ratio = 2·U(He II)/U(He I) = 2·2/1 = 4 + let s_he1 = saha_factor(24.587, t, 2.0 * 2.0 / 1.0); + let x_he1 = if ne > 0.0 { s_he1 / (s_he1 + ne) } else { 0.0 }; + + // He 第二次电离: n(He++)/n(He+), u_ratio = 2·U(He III)/U(He II) = 2·1/2 = 1 + let s_he2 = saha_factor(54.418, t, 2.0 * 1.0 / 2.0); + let x_he2 = if ne > 0.0 { s_he2 / (s_he2 + ne) } else { 0.0 }; + + let n_hplus = n_h * x_h; + let n_heplus = n_he * x_he1 * (1.0 - x_he2); + let n_heplusplus = n_he * x_he1 * x_he2; + + let ne_new = n_hplus + n_heplus + 2.0 * n_heplusplus; + + // Clamp ne to be at most an * 0.99 (physical limit) + let ne_clamped = ne_new.min(an * 0.99).max(0.0); + // Under-relaxation (0.5) for stability + let ne_relaxed = 0.5 * ne + 0.5 * ne_clamped; + + let rel = (ne_relaxed - ne).abs() / ne.max(1e-30); + ne = ne_relaxed; + + // 缓存离子化数据供 call_meanopt 使用 + self.nh_total = n_h; + self.np = n_hplus; + self.nh_neutral = n_h * (1.0 - x_h); + self.nhe_total = n_he; + self.nhe_plus = n_heplus; + self.nhe_plusplus = n_heplusplus; + + if rel < 1e-6 { break; } + } + + self.ne = ne; + ne + } + + /// 从 Saha 求解器结果计算 39 能级 LTE populations。 + /// + /// H-He 模型结构 (0-based): + /// - Ion 0 (H I): levels 0-8 (n=1..9), continuum=9 + /// - Ion 1 (He I): levels 10-23 (14 levels), continuum=24 + /// - Ion 2 (He II): levels 24-37 (n=1..14), continuum=38 + fn compute_lte_populations(&self, t: f64) -> Vec> { + const NLEVEL: usize = 39; + let mut popul = vec![vec![0.0; 1]; NLEVEL]; + let kt_ev = 8.617333e-5 * t; // kT in eV + if kt_ev <= 0.0 { return popul; } + + // H I: levels 0-8 (n=1..9) + // LTE: n(H,n) = n(H) × 2n² × exp(-E_n/kT) / U_H + let nh = self.nh_neutral; + let uh = 2.0; // 简化配分函数 + for n in 1..=9_usize { + let nn = n as f64; + let en_ev = 13.595 * (1.0 - 1.0 / (nn * nn)); // eV above ground + let gn = 2.0 * nn * nn; + let theta = en_ev / kt_ev; + if theta < 150.0 { + popul[n - 1][0] = nh * gn * (-theta).exp() / uh; + } + } + + // H II: level 9 (质子) + popul[9][0] = self.np; + + // He I: levels 10-23 (14 levels) — 真实能级数据 (he1.dat, 非氢原子不能用氢公式) + // LTE Boltzmann: popul[il] = n(He I) × g(il) × exp(-E_exc(il)/kT) / U(He I) + let nhe_neutral = (self.nhe_total - self.nhe_plus - self.nhe_plusplus).max(0.0); + let mut uhe = 0.0; + for ilev in 0..14 { + let theta = self.he1_e_ev[ilev] / kt_ev; + let boltz = if theta < 150.0 { (-theta).exp() } else { 0.0 }; + uhe += self.he1_g[ilev] * boltz; + } + if uhe < 1.0 { uhe = 1.0; } + for ilev in 0..14 { + let theta = self.he1_e_ev[ilev] / kt_ev; + let boltz = if theta < 150.0 { (-theta).exp() } else { 0.0 }; + popul[10 + ilev][0] = nhe_neutral * self.he1_g[ilev] * boltz / uhe; + } + + // He II: levels 24-37 (n=1..14), hydrogenic Z=2 + // LTE: n(HeII,n) = n(He⁺) × 2n² × exp(-E_n/kT) / U_HeII + let nhe_plus = self.nhe_plus; + let uhe2 = 2.0; + for n in 1..=14_usize { + let nn = n as f64; + let en_ev = 54.418 * (1.0 - 1.0 / (nn * nn)); + let gn = 2.0 * nn * nn; + let theta = en_ev / kt_ev; + if theta < 150.0 { + popul[24 + n - 1][0] = nhe_plus * gn * (-theta).exp() / uhe2; + } + } + + // He III: level 38 + popul[38][0] = self.nhe_plusplus; + + popul + } +} + +/// 用 Opacf0State 频率积分计算 Rosseland + Planck 平均不透明度 (per gram)。 +/// +/// 纯辅助函数: call_meanopt (灰大气) 与种群诊断共用同一积分逻辑。 +/// 返回 (κ_R per gram, κ_P per gram, 首频点吸收 ab0, 首频点散射 sct0)。 +fn kr_rosseland_planck( + grid: &LteFrequencyGrid, + opacf0: &Opacf0State, + gffpar: &GffPar, + t: f64, + ne: f64, + popul: &[Vec], + id: usize, + rho: f64, +) -> (f64, f64, f64, f64) { + let hkt = HK / t; + let mut abr = 0.0; // Rosseland: Σ (dB/dν / κ) + let mut sumdb = 0.0; // Rosseland: Σ dB/dν + let mut abp = 0.0; // Planck: Σ (Bν × κ_abs) + let mut sumb = 0.0; // Planck: Σ Bν + let mut ab_first = 0.0_f64; + let mut sct_first = 0.0_f64; + + for (ij, &fr) in grid.freq.iter().enumerate() { + let w = grid.weights[ij]; + let bnue = grid.bnue[ij]; + + let x = hkt * fr; + let ex = x.min(150.0).exp(); + let e1 = 1.0 / (ex - 1.0).max(1e-30); + + // ∂B_ν/∂T · w. ∂B_ν/∂T = bnue·e1²·ex·u/T (u=hkt·fr). + // plan·u·ex·e1 = bnue·e1²·ex·u·w = (∂B_ν/∂T·w)·T — the constant T + // cancels in the Rosseland ratio ΣW/Σ(W/κ), so this is exact. + // Faithful MEANOPT (meanopt.rs:83) and Fortran tlusty208.f:22433 + // (DPLAN=PLAN*X/T/(UN-UN/EX) = PLAN·u·ex·e1/T) both use a SINGLE e1. + // A previous version had `*e1*e1` here — the extra frequency-dependent + // e1 (=1/(e^u-1)) re-weights the Rosseland mean toward low frequencies + // and corrupts the absolute κ_R (the bug is invisible in the + // simp-vs-gold ratio since both use the same formula). + let plan = bnue * e1 * w; + let dplan = plan * hkt * fr * ex * e1; + + let (ab, sct, _, _) = opacf0.compute_opacity(fr, t, ne, popul, id, gffpar); + + if ij == 0 { + ab_first = ab; + sct_first = sct; + } + + let total = ab + sct; + if total > 0.0 { + abr += dplan / total; + } + sumdb += dplan; + abp += plan * ab; + sumb += plan; + } + + let oprol = if abr > 0.0 { sumdb / abr } else { 0.0 }; + let opplal = if sumb > 0.0 { abp / sumb } else { 0.0 }; + let opros = (oprol / rho).max(0.01); + let oppla = (opplal / rho).max(0.01); + (opros, oppla, ab_first, sct_first) +} + +impl RossopCallbacks for LteRossopCallbacks { + fn call_eldens(&mut self, _id: usize, t: f64, an: f64) -> (f64, f64, f64, f64) { + let ane = self.solve_saha(t, an); + let wm = self.wmm.get(0).copied().unwrap_or(2.3e-24); + (ane, 0.0, 0.0, wm) + } + + fn call_wnstor(&mut self, _id: usize) {} + fn call_steqeq(&mut self, _id: usize) {} + fn call_opacf0(&mut self, _id: usize, _nfreq: usize) {} + fn call_meanop(&mut self, _t: f64) -> (f64, f64) { + // ioptab=-1 路径不使用 call_meanop + (0.0, 0.0) + } + + fn call_meanopt(&mut self, t: f64, _id: usize, rho: f64) -> (f64, f64) { + if rho <= 1e-30 { + return (0.34, 0.34); + } + + // 更新 Gaunt 因子参数（温度可能已变化） + gfree0(0, &[t], &mut self.gffpar); + + // 从 Saha 求解器结果构建 39 能级 LTE populations + let popul = self.compute_lte_populations(t); + + let (opros, oppla, debug_ab_first, debug_sct_first) = kr_rosseland_planck( + &self.grid, &self.opacf0_state, &self.gffpar, + t, self.ne, &popul, 0, rho, + ); + + eprintln!("Step1b: t={:.0} ne={:.2e} rho={:.2e} ab0={:.3e} sct0={:.3e} kR={:.3} kP={:.3} popH1={:.2e} popHe1={:.2e} nhe_tot={:.2e}", + t, self.ne, rho, debug_ab_first, debug_sct_first, opros, oppla, + popul[0][0], popul[10][0], self.nhe_total); + + (opros, oppla) + } + + fn get_wmm(&self, id: usize) -> f64 { + self.wmm.get(id).copied().unwrap_or(2.3e-24) + } +} + +/// 诊断 [TLUSTY_KR_POPDIAG]: 在 gold 收敛模型的 T/ρ/Ne 条件下, 对比 +/// (a) gold 实际 39 能级 populations 算出的局部 κ_R, 与 +/// (b) 简化 LTE populations (solve_saha + compute_lte_populations, uh=2.0) 算出的 κ_R。 +/// 目的: 判定 κ_R 偏差是否由种群驱动 → 决定移植完整 ELDENS/STEQEQ 是否有益。 +/// 若 simp/gold ≈ 1 → 种群非根因 (κ_R 正确, DM 偏差为结构/积分问题); +/// 若 simp/gold ≫ 1 → 种群驱动, 移植完整种群机制可修复。 +pub fn diag_grey_kr_population_effect( + teff: f64, + wmm0: f64, + gold_temp: &[f64], + gold_elec: &[f64], + gold_dens: &[f64], + gold_popul: &[Vec], // [nlevel][nd] +) { + let wmm_arr = vec![wmm0; gold_temp.len().max(1)]; + let mut cb = LteRossopCallbacks::new(teff, &wmm_arr, 0.7, 0.28); + eprintln!("=== DIAG κ_R: gold-pops vs simplified-pops (at gold T/ρ/Ne) ==="); + eprintln!(" simp/gold > 1 ⇒ simplified overestimates κ_R ⇒ populations drive it"); + for &d in &[0usize, 10, 20, 30, 35, 40, 45, 49, 55, 60, 65, 69] { + if d >= gold_temp.len() || d >= gold_dens.len() || d >= gold_elec.len() { + continue; + } + let t = gold_temp[d]; + let ne_gold = gold_elec[d]; + let rho = gold_dens[d]; + if rho <= 1e-30 || t <= 0.0 { + continue; + } + gfree0(0, &[t], &mut cb.gffpar); + + // (a) gold 实际 populations at depth d + let pops_gold: Vec> = (0..39usize) + .map(|il| { + vec![gold_popul + .get(il) + .and_then(|lvl| lvl.get(d)) + .copied() + .unwrap_or(0.0)] + }) + .collect(); + let (kr_gold, _, ab_g, _sct_g) = kr_rosseland_planck( + &cb.grid, + &cb.opacf0_state, + &cb.gffpar, + t, + ne_gold, + &pops_gold, + 0, + rho, + ); + + // (b) 简化 populations at gold T; 用 gold 总粒子密度作 an 以保证 ne 比较公平 + let an_gold = rho / wmm0 + ne_gold; + cb.solve_saha(t, an_gold); + let pops_simp = cb.compute_lte_populations(t); + let (kr_simp, _, _, _) = kr_rosseland_planck( + &cb.grid, + &cb.opacf0_state, + &cb.gffpar, + t, + cb.ne, + &pops_simp, + 0, + rho, + ); + + let ratio = kr_simp / kr_gold.max(1e-9); + eprintln!( + " d={:>2} T={:>7.0} ρ={:.2e} | ne_gold={:.2e} ne_simp={:.2e} | κR_gold={:.3} κR_simp={:.3} simp/gold={:.2} ab0_gold={:.2e}", + d, t, rho, ne_gold, cb.ne, kr_gold, kr_simp, ratio, ab_g + ); + + // 在峰值深度 (d=40, 偏差最大) 打印关键能级 populations gold vs simplified + if d == 40 { + eprintln!(" --- populations at d=40 (gold vs simplified) ---"); + // 关键能级: 0=H I gnd, 1-8=H I exc, 9=H II, 10=He I gnd, 24=He II gnd, 38=He III + for &il in &[0usize, 1, 9, 10, 11, 24, 25, 38] { + let pg = pops_gold.get(il).and_then(|v| v.first()).copied().unwrap_or(0.0); + let ps = pops_simp.get(il).and_then(|v| v.first()).copied().unwrap_or(0.0); + let rr = ps / pg.max(1e-300); + eprintln!(" lvl {:>2}: gold={:.3e} simp={:.3e} simp/gold={:.2e}", il, pg, ps, rr); + } + } + } + eprintln!("=== DIAG end ==="); +} + +/// 诊断 [TLUSTY_KR_FORCOND]: **隔离 opacity/EOS 公式 vs 网格反馈**。 +/// +/// 在 Fortran-grey 模型**自身**的 (T, ρ, Ne) 条件下（而非 gold 收敛模型条件， +/// 也非 Rust-grey 压缩条件），用 Rust opacf0_state + solve_saha 计算局部 κ_R +/// 与 ne，对比 Fortran 精确 κ_R = dTAUROSS/dMASS 与 Fortran ne。 +/// +/// 这消除历次诊断的条件混淆: +/// - diag_grey_kr (simp/gold) 用 gold 条件 + 同公式 ⇒ 只隔离种群, 不校验绝对 κ_R +/// - Step1b κ_R 用 Rust-grey 压缩条件 ⇒ 条件≠Fortran, 无法判断公式是否对 +/// - 本诊断用 Fortran-grey 条件 ⇒ 若 Rust κ_R≈Fortran ⇒ 公式对, 偏差=网格反馈; +/// 若 Rust κ_R≠Fortran ⇒ 公式/EOS 有真实偏差 (可分解 ne 是否匹配定位 EOS vs 公式) +/// +/// 关键双判据: +/// (1) ne_rust/ne_fort ≈ 1 ⇒ 简化 Saha EOS 在 Fortran 条件下正确 (电离平衡对) +/// (2) kr_rust/kr_fort ≈ 1 ⇒ opacity 公式在 Fortran 条件下正确 +/// 两判据给出四种诊断组合, 精确定位剩余前沿是 EOS / opacity / 网格 / 皆有。 +pub fn diag_grey_kr_fortran_conditions( + teff: f64, + wmm0: f64, + conds: &[(f64, f64, f64, f64)], // (T, ne_fort, dens_fort, kr_fort) per depth +) { + let wmm_arr = vec![wmm0; conds.len().max(1)]; + let mut cb = LteRossopCallbacks::new(teff, &wmm_arr, 0.7, 0.28); + eprintln!("=== DIAG κ_R+ne: Rust (at Fortran-grey conditions) vs Fortran ==="); + eprintln!(" (ne_rust/ne_fort≈1 ⇒ EOS correct | kr_rust/kr_fort≈1 ⇒ opacity formula correct)"); + eprintln!( + " {:>3} {:>9} {:>11} {:>11} {:>10} {:>10}", + "id", "T(K)", "ne_rust", "ne_fort", "ne_ratio", "kr_ratio" + ); + for (i, &(t, ne_fort, dens, kr_fort)) in conds.iter().enumerate() { + if dens <= 1e-30 || t <= 0.0 || ne_fort <= 0.0 || kr_fort <= 0.0 { + continue; + } + gfree0(0, &[t], &mut cb.gffpar); + // an = 总粒子密度 (核子+电子), 与 call site 约定一致 + let an = dens / wmm0 + ne_fort; + let ne_rust = cb.solve_saha(t, an); + let pops = cb.compute_lte_populations(t); + let (kr_rust, _, _, _) = kr_rosseland_planck( + &cb.grid, + &cb.opacf0_state, + &cb.gffpar, + t, + ne_rust, + &pops, + 0, + dens, + ); + let ne_ratio = ne_rust / ne_fort; + let kr_ratio = kr_rust / kr_fort; + eprintln!( + " {:>3} {:>9.1} {:>11.3e} {:>11.3e} {:>10.4} {:>10.4} | kr_rust={:.4} kr_fort={:.4}", + i + 1, + t, + ne_rust, + ne_fort, + ne_ratio, + kr_ratio, + kr_rust, + kr_fort + ); + } + eprintln!("=== DIAG κ_R+ne end ==="); +} + /// ROSSOP 计算。 /// -/// 使用 rossop 模块计算温度、Hopf 函数和基本密度。 +/// 使用 rossop_with_callbacks + LteRossopCallbacks 计算温度、电子密度和不透明度。 +/// 使用 ioptab=-1 路径：ELDENS（Saha 方程）+ meanopt（LTE 不透明度积分）。 fn rossop_calc( id: usize, taur: f64, @@ -494,8 +1062,14 @@ fn rossop_calc( dens: &mut [f64], abrosd: &mut [f64], abplad: &mut [f64], + callbacks: &mut LteRossopCallbacks, ) -> (f64, f64, f64) { - let config = RossopConfig::default(); + let config = RossopConfig { + ioptab: -1, // 简化模式: ELDENS + meanopt + iter: 1, + nfreq: 1, + ifrayl: 0, + }; let params = RossopParams { id, @@ -515,9 +1089,8 @@ fn rossop_calc( abplad, }; - let output: RossopOutput = rossop(&config, ¶ms, &mut state); + let output: RossopOutput = rossop_with_callbacks(&config, ¶ms, &mut state, callbacks); - // 返回温度、电子密度和 Rosseland 不透明度 (output.t, output.ane, output.abross) } @@ -590,6 +1163,7 @@ mod tests { let taur = 1.0; // Rosseland 光学深度 = 1 let p = 1e4; // 压力 (cgs) let hopf = 0.0; // 使用精确 Hopf 函数 + let mut callbacks = LteRossopCallbacks::new(teff, &wmm, 0.7, 0.28); let (t, _ane, abros) = rossop_calc( 0, @@ -603,6 +1177,7 @@ mod tests { &mut dens, &mut abrosd, &mut abplad, + &mut callbacks, ); // 验证温度计算 @@ -621,6 +1196,39 @@ mod tests { assert!(t > 5000.0 && t < 15000.0, "Temperature {} out of reasonable range for Teff={}", t, teff); } + /// 诊断：用 Opacf0State (readmodel 路径所用的完整机制) 在 gold-ref 条件下 + /// 计算 Rosseland 平均，与 lte_meanopt 的发散值及 gold-ref 有效 κ 对比。 + /// 目标：判断 create_grey_atmosphere 改用 Opacf0State 是否能给出物理 κ_R。 + #[test] + fn diag_opacf0_kappar_at_goldref() { + // (label, T, rho, eff_kappa=tau/DM) from tests/tlusty/hhe_fortran/fort.7.ref + let pts: &[(&str, f64, f64, f64)] = &[ + ("d0 surf", 26306.2, 7.304e-16, 0.343), + ("d20 ", 27030.2, 4.455e-13, 0.343), + ("d34 ", 27800.0, 1.0e-11, 0.396), + ("d49 ", 40000.0, 2.0e-10, 0.951), + ("d60 ", 90000.0, 4.0e-9, 1.145), + ("d69 deep ", 137872.1, 1.102e-7, 1.061), + ]; + let wmm = vec![1.3 * 1.67e-24; MDEPTH]; + let mut cb = LteRossopCallbacks::new(35000.0, &wmm, 0.70, 0.28); + let wmm_local = 1.3 * 1.67e-24; + + eprintln!("\n{:>9} {:>10} {:>12} {:>10} {:>10} | {:>8}", + "label", "kR_opacf", "ne_saha", "kP_opacf", "lte_mean", "eff_k"); + for (label, t, rho, eff_k) in pts { + // 总粒子密度 an = rho / wmm (nucleon density) + let an = rho / wmm_local; + // solve_saha 设置 cb.ne 及缓存的离子化种群 + cb.solve_saha(*t, an); + let (kros, kpla) = cb.call_meanopt(*t, 0, *rho); + eprintln!("{:>9} {:10.3} {:12.3e} {:10.3} {:10} | {:8.3}", + label, kros, cb.ne, kpla, "(see lte_meanopt diag)", eff_k); + // κ_R 必须为有限正数（不应发散到 140） + assert!(kros.is_finite() && kros > 0.0, "kR not finite/positive at {}", label); + } + } + #[test] fn test_ltegr_temperature_profile() { // 测试 LTEGR 生成的温度分布 diff --git a/src/tlusty/io/outpri.rs b/src/tlusty/io/outpri.rs index e46c846..ebb2ed8 100644 --- a/src/tlusty/io/outpri.rs +++ b/src/tlusty/io/outpri.rs @@ -25,8 +25,10 @@ use crate::tlusty::state::constants::{UN, HALF}; // Then restore LTE=false // 物理常数 -/// Stefan-Boltzmann 常数 × 4 -const SIG4P: f64 = 7.5657e-5; +/// Stefan-Boltzmann 常数 / 4π (matches Fortran BASICS.FOR: SIG4P = 4.5114062D-6). +/// Was 7.5657e-5 — a transcription error (neither σ×4=2.27e-4 nor σ/4π); used +/// only by `fltt = SIG4P*teff^4` (Fortran FLTT, tlusty208.f:14179). +const SIG4P: f64 = 4.5114062e-6; /// Boltzmann 常数 const BOLK: f64 = 1.38054e-16; /// 光速 × 1e18 (用于波长计算) diff --git a/src/tlusty/io/resolv.rs b/src/tlusty/io/resolv.rs index 15ad35b..6be15fb 100644 --- a/src/tlusty/io/resolv.rs +++ b/src/tlusty/io/resolv.rs @@ -30,7 +30,7 @@ //! - fort.6: 标准输出（进度和诊断信息） use super::FortranWriter; -use crate::tlusty::state::constants::{MFREQ, MTRANS, H, HK, BOLK, HMASS, SIG4P, PI}; +use crate::tlusty::state::constants::{MFREQ, MTRANS, H, HK, BOLK, HMASS, SIG4P, PI, HALF}; use crate::tlusty::math::{ rayset, prd, opaini, opaini_full, rates1, steqeq_pure, newpop, ratmat, RatmatParams, @@ -42,12 +42,12 @@ use crate::tlusty::math::{ conout, ConoutParams, ConoutConfig, conref, ConrefParams, ConrefConfig, alisk2, alist1, alist2, pzevld, hesol6, dmeval, rybheq, RybheqParams, RybheqConfig, linsel, LinselConfig, LinselAtomicParams, LinselFreqParams, - feautrier_solve, + rtefr1, Rtefr1Params, Rtefr1ModelState, AMU, WTMU, dmder, DmevalParams, Hesol6Params, ElcorConfig, ElcorParams, SteqeqParams, NewpopParams, eldens, EldensParams, EldensConfig, generate_inifrc_frequency_grid, ABUND_H, ABUND_HE, WMM_GREY, - lucy, LucyConfig, LucyModelParams, + lucy, LucyConfig, LucyModelParams, LucyNgState, OpacflPointData, Rad1PointData, wnstor, sabolf, SabolfParams, SteqeqConfig, @@ -86,18 +86,49 @@ pub fn compute_abundance_params(atomic: &crate::tlusty::state::atomic::AtomicDat .map(|a| a.first().copied().unwrap_or(0.0)) .collect(); - let mut ytot = 1.0_f64; - let mut inv_wmy = 1.0 / amass[0].max(1e-30); // 1/H mass + // Fortran INITIA (tlusty208.f:2816-2831): + // YTOT(ID) = Σ ABNDD(I) (number abundance, H=1) + // WMY(ID) = Σ ABNDD(I) * AMAS(I) (amu, H=1.008, He=4.003, ...) + // WMM(ID) = WMY * HMASS / YTOT (mean mass per nucleus) + // Note: the old code used the harmonic form 1/Σ(abund/amass) and omitted + // the /YTOT, producing a per-electron mass instead of a per-nucleus mass. + let has_data = abund.iter().zip(amass.iter()) + .any(|(&a, &m)| a > 0.0 && m > 0.0); - for i in 1..abund.len() { - if abund[i] > 0.0 && amass[i] > 0.0 { - ytot += abund[i]; - inv_wmy += abund[i] / amass[i]; + let (ytot, wmy) = if has_data { + let mut y = 0.0_f64; + let mut w = 0.0_f64; + for i in 0..abund.len() { + if abund[i] > 0.0 && amass[i] > 0.0 { + y += abund[i]; + w += abund[i] * amass[i]; + } } - } + (y, w) + } else { + // Fallback: default solar-composition H-He model (matches the HHe test + // suite abundances in the Fortran gold reference: YTOT=1.08588, + // WMY=1.35982, WMM=2.09547e-24). Used until atopar.abund/amass are + // populated from the atomic-data reader. + // (abund, amass) for H, He, C, N, O + const FB: [(f64, f64); 5] = [ + (1.0, 1.008), + (0.0851, 4.003), + (2.45e-4, 12.011), + (6.03e-5, 14.007), + (4.57e-4, 15.999), + ]; + let mut y = 0.0_f64; + let mut w = 0.0_f64; + for &(a, m) in &FB { y += a; w += a * m; } + (y, w) + }; - let wmy = 1.0 / inv_wmy.max(1e-30); - let wmm = crate::tlusty::state::constants::HMASS * wmy; + let wmm = if ytot > 1e-30 { + crate::tlusty::state::constants::HMASS * wmy / ytot + } else { + crate::tlusty::state::constants::HMASS + }; (wmy, ytot, wmm) } @@ -369,6 +400,14 @@ pub fn resolv( let init = config.init; let lfin = config.lfin; + // When INPC is active (TLUSTY_INPC), ELEC is an independent SOLVES variable + // solved by the BPOPC charge-conservation row. RESOLV must NOT overwrite it + // via ELCOR's nonlinear charge-neutrality correction (which uses the full + // Saha machinery, inconsistent with BPOPC's simplified charge_sum_hhe). + // Letting both run makes ELEC oscillate between the two Saha values. + // DENS stays consistent via main.rs DENS=(TOTN-ELEC)·WMM after SOLVES. + let inpc_active = std::env::var("TLUSTY_INPC").is_ok(); + // ----------------------------------------------------------- // Part 1: 初始化 - INILAM // ----------------------------------------------------------- @@ -761,11 +800,13 @@ pub fn resolv( params.model.frqall.freq[ij] = freq[ij]; params.model.frqall.w[ij] = weights[ij]; } + // Store INIFRC edge frequencies as the explicit set for SOLVES. + // BRE integral form sums T-derivative over ALL frequencies (nfreq_total), + // so the explicit set only needs to cover the physically important edges. params.model.frqall.nfreqe = freq_grid.ijfr.len(); params.model.frqall.ijfr_explicit = freq_grid.ijfr.clone(); // 构建反向映射: grid index -> 是否为 explicit frequency - // ijfr_explicit[ije] = ij, 所以我们需要一个 set 来快速查找 let explicit_set: std::collections::HashSet = freq_grid.ijfr.iter().copied().collect(); @@ -776,10 +817,14 @@ pub fn resolv( // 从原子数据计算丰度参数 let (wmy_init, ytot_init, _) = compute_abundance_params(params.atomic); - // 深度间隔 deldm + // 深度间隔 deldm — Fortran DELDM(ID)=HALF*(DM(ID+1)-DM(ID)) (HALF 前向差分)。 + // 经 0基↔1基映射后 deldm[id-1]=HALF*(dm[id]-dm[id-1]) 与 DELDM(ID) 逐项匹配。 + // 此数组供 Lucy (dtm/fluz/toth 通量) 与 ROSSTD (dtaur/taurs) 使用, 二者 Fortran + // 公式均以 DELDM 为 HALF 间距推导 (DTAUR=DELDM·(κ+κ)=Δτ 梯形; toth=Σw·dK/dτ=H)。 + // 缺 HALF 会使 dtm/dtaur 大 2×, 致 Lucy toth=½·TEF4 (delh 爆炸发散) 且 taurs 翻倍。 let mut deldm = vec![0.0; nd]; for id in 1..nd { - deldm[id - 1] = params.model.modpar.dm[id] - params.model.modpar.dm[id - 1]; + deldm[id - 1] = HALF * (params.model.modpar.dm[id] - params.model.modpar.dm[id - 1]); } // 这些数组在 lambda 循环内更新，循环后仍需使用 @@ -798,6 +843,11 @@ pub fn resolv( let mut opacfl_data: Vec = Vec::new(); let mut rad1_data: Vec = Vec::new(); + // Lucy Ng 加速的持久化状态——必须在 ilam 循环外创建，使 TEM0..3 / LAC2T / IACLT + // 跨 lambda 迭代累积（resolv 的 ilam 循环 = Fortran LUCY 的内层 ilucy 循环）。 + // iaclt=7 为 Fortran NSTPAR 默认（nstpar.rs:549），iacldt 在下方 LucyConfig 设 4。 + let mut lucy_ng = LucyNgState::new(nd, 7); + for _ilam_iter in 1..=nlambd { ilam = _ilam_iter; debug_log!("RESOLV: Lambda iteration {} of {}", ilam, nlambd); @@ -854,6 +904,58 @@ pub fn resolv( wmm_arr[id] = eldens_output.wm; } + // ============================================================== + // Step 1b: Call WNSTOR for base state — occupation probabilities + // ============================================================== + // WNSTOR computes WOP (occupation probabilities) and WNHINT for each + // depth point. These are needed by compute_opacity_with_wop to correctly + // compute the BF emission coefficient: + // EMTRA = POPUL(II) * WOP(II) * CORR + // Without WOP, the emission is overestimated, leading to excessive + // stimulated emission subtraction and an atmosphere that is too transparent. + // + // Initialize ifwop from nquant if not already set (Fortran RDATA sets + // ifwop = nquant for hydrogenic levels during atomic data reading). + // In our Rust START, ifwop defaults to all zeros which disables WOP. + if params.model.wmcomp.ifwop.iter().take(config.nlevel).all(|&v| v == 0) { + for ii in 0..config.nlevel { + let nq = params.atomic.levpar.nquant[ii]; + params.model.wmcomp.ifwop[ii] = if nq > 0 { nq } else { 0 }; + } + eprintln!("RESOLV: Initialized ifwop from nquant, ifwop[0..5]={:?}, nquant[0..5]={:?}, ioptab={}", + ¶ms.model.wmcomp.ifwop[..5.min(config.nlevel)], + ¶ms.atomic.levpar.nquant[..5.min(config.nlevel)], + params.tlusty_config.basnum.ioptab); + } + { + // Build an elec array from ne_arr for WNSTOR (avoids borrow conflict) + let mut elec_base: Vec = params.model.modpar.elec.to_vec(); + for id in 0..nd { + elec_base[id] = ne_arr[id]; + } + for id in 0..nd { + wnstor( + id, + ¶ms.model.modpar.temp, + &elec_base, + ¶ms.tlusty_config.invint.xi2, + &mut params.model.wmcomp.wnhint, + &mut params.model.wmcomp.wop, + ¶ms.model.wmcomp.ifwop, + config.nlevel, + ¶ms.atomic.levpar.nquant, + ¶ms.atomic.levpar.iel, + ¶ms.atomic.ionpar.iz, + params.tlusty_config.basnum.ioptab, + config.lte, + ); + } + eprintln!("RESOLV Step1b: WOP after WNSTOR: wop[0][0..3]=[{:?}] wop[1][0..3]=[{:?}] wop[9][0..3]=[{:?}]", + ¶ms.model.wmcomp.wop[0][..3.min(nd)], + ¶ms.model.wmcomp.wop[1][..3.min(nd)], + ¶ms.model.wmcomp.wop[9][..3.min(nd)]); + } + // ============================================================== // Step 2: Compute temperature derivatives for exprad (dabt, demt) // ============================================================== @@ -1019,9 +1121,10 @@ pub fn resolv( for ij in 0..nfreq_actual { let fr = freq[ij]; - // Baseline opacity at (T, ne, original populations) - let (true_abs, scat, _emis_pre_stim, _rayleigh) = opacf0_state.compute_opacity( + // Baseline opacity at (T, ne, original populations) — use WOP from base WNSTOR + let (true_abs, scat, _emis_pre_stim, _rayleigh) = opacf0_state.compute_opacity_with_wop( fr, t, ne, ¶ms.model.levpop.popul, id, ¶ms.model.gffpar, + ¶ms.model.wmcomp.wop, ); let abso_cm = true_abs + scat; @@ -1058,12 +1161,52 @@ pub fn resolv( rad1_data.reserve(nfreq_actual); opacfl_data.reserve(nfreq_actual); + // Rosseland-mean opacity (ABROSD) accumulation — matches Fortran ROSSTD(IJ>0) + // contribution called from RATES1 when LROSS = (NDRE<=0 .AND. ITER==1) .OR. LFIN. + // Units: ABSO1 is opacity per cm (cm⁻¹), freq in Hz, W = Simpson weight (dν in Hz), + // h_over_c2 = 2H/c², HKT21 = HK/T². The final + // ABROSD(ID) = SUMDPL(ID) / (Σ DPLAN/ABSO1 · DENS(ID)) + // is the Rosseland mean κ_R per gram (cm²/g); it is independent of any uniform + // scaling of W because numerator and denominator scale together. + let lross = (params.tlusty_config.matkey.ndre <= 0 && iter == 1) || lfin; + let mut ros_abrosd_int = vec![0.0_f64; nd]; // Σ DPLAN/ABSO1 + let mut ros_sumdpl = vec![0.0_f64; nd]; // Σ DPLAN + + // ----- RTEFR1 (variable-Eddington two-pass formal solver) scratch ----- + // Allocated once, reused per frequency. In the LTE path + // (isplin=0, idisk=0, ilmcor=0, nelsc=0, iwinbl=0, ifalih=0, ifprad=0) + // rtefr1 writes only rad1/fak1/ali1 (the outputs we consume) plus the + // per-frequency slot [ij] of fh/fhd/q0/uu0/flux; the radex/fakex/rad/fak + // arrays are gated by ijex/idisk (both 0) and are never indexed, so empty + // Vecs are type-valid and safe. hextrd=0 (no external radiation). + let rte_deldmz: Vec = (0..nd.saturating_sub(1)) + .map(|id| HALF * (params.model.modpar.dm[id + 1] - params.model.modpar.dm[id])) + .collect(); + let rte_hextrd = vec![0.0_f64; nfreq_actual]; + let rte_zero_nd = vec![0.0_f64; nd]; // emel1 (unused, nelsc=0) + let rte_zero_ij = vec![0i32; nfreq_actual]; // ijali/ijex/kij (all 0) + let rte_zero_albe = vec![0.0_f64; nfreq_actual]; // albe (unused, iwinbl=0) + let mut rte_flux = vec![0.0_f64; nfreq_actual]; + let mut rte_fh = vec![0.0_f64; nfreq_actual]; + let mut rte_fhd = vec![0.0_f64; nfreq_actual]; + let mut rte_q0 = vec![0.0_f64; nfreq_actual]; + let mut rte_uu0 = vec![0.0_f64; nfreq_actual]; + let mut rte_pradt = vec![0.0_f64; nd]; // unused (ifprad=0) + let mut rte_prada = vec![0.0_f64; nd]; // unused (ifprad=0) + let mut rte_prd0 = 0.0_f64; // unused (ifprad=0) + // radex/fakex/rad/fak: never indexed in LTE path (ijex=0, idisk=0). + let mut rte_radex: Vec> = Vec::new(); + let mut rte_fakex: Vec> = Vec::new(); + let mut rte_rad: Vec> = Vec::new(); + let mut rte_fak: Vec> = Vec::new(); + let rte_lskip: Vec> = Vec::new(); // unused (isplin<5, ifprad=0) + let rte_extint: Vec> = Vec::new(); // unused (iwinbl>=0, ifalih=0) + for ij in 0..nfreq_actual { let fr = freq[ij]; // Build per-depth arrays for this frequency let mut abso_ij = vec![0.0; nd]; - let mut true_abs_ij = vec![0.0; nd]; let mut scat_ij = vec![0.0; nd]; let mut emis_ij = vec![0.0; nd]; let mut source_ij = vec![0.0; nd]; @@ -1075,14 +1218,14 @@ pub fn resolv( let ne = ne_arr[id]; let hkt = HK / t; - // Compute opacity at this (depth, frequency) - let (true_abs, scat_tot, _emis_pre, _ray) = opacf0_state.compute_opacity( + // Compute opacity at this (depth, frequency) — use WOP from base WNSTOR + let (true_abs, scat_tot, _emis_pre, _ray) = opacf0_state.compute_opacity_with_wop( fr, t, ne, ¶ms.model.levpop.popul, id, ¶ms.model.gffpar, + ¶ms.model.wmcomp.wop, ); let abso_cm = true_abs + scat_tot; // total opacity per cm abso_ij[id] = abso_cm; - true_abs_ij[id] = true_abs; scat_ij[id] = scat_tot; // Emission = true_abs * Bν(T) per cm (for ILMCOR=3 source) @@ -1090,22 +1233,114 @@ pub fn resolv( let bnu = h_over_c2 * fr.powi(3) / (x.exp() - 1.0).max(1e-100); emis_ij[id] = true_abs * bnu; source_ij[id] = bnu; // Planck function (for OpacflPointData) + + // Rosseland contribution (Fortran ROSSTD IJ>0): + // XKF=EXP(-HKT1*FR); XKF1=1-XKF; XKFB=XKF*BNUE; BNUE=2h/c²·ν³ + // PLAN=XKFB/XKF1*W; DPLAN=PLAN/XKF1*FR*HKT21 + // ABROSD+=DPLAN/ABSO1; SUMDPL+=DPLAN + // (x is already hkt*fr clamped to 150; bnu above = h_over_c2·fr³·xkf/xkf1, + // so xkfb = bnu·xkf1 = XKF·BNUE — reuse it instead of recomputing.) + if lross { + let hkt21 = HK / (t * t); + let xkf = (-x).exp(); + let xkf1 = 1.0 - xkf; + if xkf1.abs() > 1e-300 { + let xkfb = bnu * xkf1; // = XKF·BNUE + let plan = xkfb / xkf1 * weights[ij]; + let dplan = plan / xkf1 * fr * hkt21; + let abso1_safe = abso_cm.max(1e-30); + ros_abrosd_int[id] += dplan / abso1_safe; + ros_sumdpl[id] += dplan; + } + } } - // Feautrier solve: returns Jν (rad1) and Eddington factor (fak1) - let result = feautrier_solve( + // RTEFR1 formal solution: variable-Eddington two-pass solver. + // Replaces the simplified feautrier_solve (constant f=1/3 single pass) + // whose ~1% deep-layer J deviation was amplified by the ill-conditioned + // 146×146 SOLVES system into large RE residuals. RTEFR1 first solves the + // multi-angle Feautrier system to determine the variable Eddington factor + // fkk=K/J per depth, then re-solves the scalar system with that fkk for + // strict consistency (matching Fortran LTE path: alisk→rtefr1). + // + // absot = opacity per gram (abso/dens); deldmz = HALF*(dm[id+1]-dm[id]). + let mut absot_ij = vec![0.0; nd]; + for id in 0..nd { + let d = params.model.modpar.dens[id]; + absot_ij[id] = if d > 0.0 { abso_ij[id] / d } else { abso_ij[id] }; + } + + let mut rad1_out = vec![0.0; nd]; + let mut fak1_out = vec![0.0; nd]; + let mut ali1_out = vec![0.0; nd]; + let mut alim1_out = vec![0.0; nd]; + let mut alip1_out = vec![0.0; nd]; + + let rte_params = Rtefr1Params { + ij, nd, - &abso_ij, - &true_abs_ij, - &scat_ij, - &emis_ij, - ¶ms.model.modpar.dm[..nd], - ¶ms.model.modpar.dens[..nd], - ¶ms.model.modpar.temp[..nd], - fr, - 0.0, // tempbd = 0 (use deepest temperature) - config.ibc, - ); + nmu: 3, // Fortran LTE formal solution uses NMU=3 (tlusty208.f:39153) + nfreq: nfreq_actual, + isplin: 0, // LTE: ordinary Feautrier + idisk: 0, // stellar atmosphere + ibc: config.ibc, // LTE default 3 (diffusion lower BC) + jali: 1, // Rybicki-Hummer Lambda* diagonal + ifali: 0, + ilmcor: 0, // no Lambda scattering correction (LTE) + ifalih: 0, + iwinbl: 0, + chmax: 0.0, + icompt: config.icompt, + iter, + ilam, + irte: 0, + ifprad: 0, + ifz0: 0, + wtmu: &WTMU, + amu: &AMU, + }; + let mut rte_state = Rtefr1ModelState { + freq, + w: weights, + deldmz: &rte_deldmz, + dm: ¶ms.model.modpar.dm[..nd], + temp: ¶ms.model.modpar.temp[..nd], + elec: &ne_arr, + absot: &absot_ij, + abso1: &abso_ij, + emis1: &emis_ij, + scat1: &scat_ij, + elscat: &scat_ij, // unused (nelsc=0); reuse scattering array + emel1: &rte_zero_nd, + rad1: &mut rad1_out, + fak1: &mut fak1_out, + ali1: &mut ali1_out, + alim1: &mut alim1_out, + alip1: &mut alip1_out, + flux: &mut rte_flux, + fh: &mut rte_fh, + fhd: &mut rte_fhd, + q0: &mut rte_q0, + uu0: &mut rte_uu0, + extint: &rte_extint, + hextrd: &rte_hextrd, + pradt: &mut rte_pradt, + prada: &mut rte_prada, + prd0: &mut rte_prd0, + lskip: &rte_lskip, + ijali: &rte_zero_ij, + ijex: &rte_zero_ij, + radex: &mut rte_radex, + fakex: &mut rte_fakex, + rad: &mut rte_rad, + fak: &mut rte_fak, + kij: &rte_zero_ij, + albe: &rte_zero_albe, + nelsc: 0, + tempbd: 0.0, // use deepest temperature + rrdil: 1.0, + }; + rtefr1(&rte_params, &mut rte_state); // Build OpacflPointData for Lucy let mut abso1_ij = vec![0.0; nd]; @@ -1126,26 +1361,45 @@ pub fn resolv( }); rad1_data.push(Rad1PointData { - rad1: result.rad1, - fak1: result.fak1, - ali1: result.alrh, - alim1: vec![0.0; nd], - alip1: vec![0.0; nd], + rad1: rad1_out, + fak1: fak1_out, + ali1: ali1_out, + alim1: alim1_out, + alip1: alip1_out, }); } + // Finalize Rosseland mean: ABROSD(ID) = SUMDPL(ID) / (Σ DPLAN/ABSO1 · DENS(ID)) + // → κ_R per gram (cm²/g). Stored into the model for Part 5 ROSSTD(0). + if lross { + for id in 0..nd { + let dens_id = params.model.modpar.dens[id]; + if ros_abrosd_int[id].abs() > 1e-300 && dens_id.abs() > 1e-300 { + params.model.opmean.abrosd[id] = + ros_sumdpl[id] / (ros_abrosd_int[id] * dens_id); + params.model.opmean.sumdpl[id] = ros_sumdpl[id]; + } else { + params.model.opmean.abrosd[id] = 0.0; + params.model.opmean.sumdpl[id] = 0.0; + } + } + } + // ============================================================== - // Step 3b: 存储到 exprad/expraf (显式频率点) + // Step 3b: 存储到 exprad/expraf (所有频率点) // ============================================================== // 对应 Fortran: OPACFD -> ABSOEX/EMISEX/SCATEX, RTEFR1 -> RADEX/FAKEX - // 只在显式频率点 (ijfr_explicit 包含的网格索引) 存储 + // Fortran stores data for ALL frequency points. For LTE (no ALI), + // all points are explicit and needed by SOLVES BRTE/BRE. { let exprad = &mut params.model.exprad; let expraf = &mut params.model.expraf; + // FHD(IJT) lives in COMMON /TOTRAD/ FHD in Fortran; store the + // RTEFR1-computed lower-boundary Eddington factor so the SOLVES + // BRTE lower boundary (matgen_lte brte_lte) can use FHD(IJT) like + // Fortran brte.f instead of the diffusion-limit constant 1/√3. + let fhd_store = &mut params.model.totrad.fhd; for ij in 0..nfreq_actual { - if !explicit_set.contains(&ij) { - continue; - } // 检查数组边界 if ij >= exprad.absoex.len() { continue; @@ -1166,6 +1420,11 @@ pub fn resolv( // FAKEX = Eddington factor K/J expraf.fakex[ij][id] = rad1_data[ij].fak1[id]; } + // FHD(IJT) = lower-boundary Eddington factor (AH/AJ at ND), + // written by RTEFR1 for IBC≠0 (Fortran rtefr1.f FHD(IJT)). + if ij < fhd_store.len() { + fhd_store[ij] = rte_fhd[ij]; + } } } @@ -1203,6 +1462,9 @@ pub fn resolv( let third = 1.0_f64 / 3.0; + // DEBUG: count how many (ij,id) pairs contribute to FCOOLI + let mut _dbg_fcooli_contribs = 0usize; + let mut _dbg_fcooli_skip = 0usize; for ij in 0..nfreq_actual { let w = weights[ij]; let wf = w * third; // WF = W * FH (Eddington factor) @@ -1218,8 +1480,14 @@ pub fn resolv( // Skip non-physical or non-finite cases (equivalent to Fortran LSKIP) if !ali1.is_finite() || !rad1.is_finite() || abso1 < 1e-30 { + _dbg_fcooli_skip += 1; continue; } + let abst = abso1 - scat1; + let fcool_contrib = w * (emis1 - abst * rad1); + if fcool_contrib.abs() > 1e-30 { + _dbg_fcooli_contribs += 1; + } // True thermal absorption (abso - electron scattering) let elscat = opacfl_data[ij].scat1[id]; // electron scattering = scat1 for now @@ -1254,11 +1522,16 @@ pub fn resolv( // FCOOLI accumulation: WW * (EMIS1 - ABST*RAD1) params.model.totflx.fcooli[id] += w * (emis1 - abst * rad1); - // REDT accumulation: WF * DSFT1 * ALI1 (for differential form) - // Only for depths where REDIF > 0 - if params.model.repart.redif[id] > 0.0 { - params.model.expraf.redt[id] += wf * dsft1 * ali1; - } + // REDT accumulation: WF * DSFT1 * ALI1 (differential-form flux + // temperature derivative). Accumulate UNCONDITIONALLY (like REIT/ + // FCOOLI above); the REDIF weighting is applied later in BRE as + // `redt*redif`. A previous `if redif[id]>0` guard here was a bug: + // REDIF is not finalized until ROSSTD runs (Part 5, below), so on + // iteration 1 REDIF is still 0 at this point, REDT stayed 0, and + // the BRE temperature Jacobian b[nre][nre] degenerated to ~0 in the + // deep (pure-diffusion) layers — making the bottom-T correction + // uncontrolled and driving SOLVES to diverge from gold on iter 1. + params.model.expraf.redt[id] += wf * dsft1 * ali1; } } @@ -1267,6 +1540,26 @@ pub fn resolv( // This captures the ALI (implicit) frequency contribution. // bre_lte subtracts the explicit frequency part and adds // the Jν coupling matrix elements. + eprintln!("ALIFR1 diag: nfreq_actual={}, contribs={}, skips={}", nfreq_actual, _dbg_fcooli_contribs, _dbg_fcooli_skip); + for id in [0, 35, 69].iter() { + let id = *id; + if id >= nd { continue; } + let abso1_0 = opacfl_data[0].abso1[id]; + let emis1_0 = opacfl_data[0].emis1l[id]; + let scat1_0 = opacfl_data[0].scat1[id]; + let rad1_0 = rad1_data[0].rad1[id]; + let ali1_0 = rad1_data[0].ali1[id]; + let abst_0 = abso1_0 - scat1_0; + let fcool_c = weights[0] * (emis1_0 - abst_0 * rad1_0); + eprintln!(" freq[0] id={}: abso={:.3e} true_abs={:.3e} scat={:.3e} emis={:.3e} rad1={:.3e} ali1={:.3e} fcool_c={:.3e}", + id, abso1_0, abst_0, scat1_0, emis1_0, rad1_0, ali1_0, fcool_c); + } + for id in [0, 35, 69].iter() { + let id = *id; + if id >= nd { continue; } + eprintln!(" accum id={}: fcooli={:.3e} reit={:.3e} redt={:.3e}", id, + params.model.totflx.fcooli[id], params.model.expraf.reit[id], params.model.expraf.redt[id]); + } } // ============================================================== @@ -1301,8 +1594,8 @@ pub fn resolv( let lucy_config = LucyConfig { itlucy: config.itlucy, - iaclt: 10, - iacldt: 1, + iaclt: 7, // Fortran NSTPAR 默认 IACLT（仅信息用；实际门控读 ng.iaclt） + iacldt: 4, // Fortran NSTPAR 默认 IACLDT：滚动 3 次后加速 1 次（Ng 节律） ihecor: 0, // Disable hydrostatic integration for LTE (EOS already gives correct dens/elec) lte: config.lte, lchc: config.lchc, @@ -1336,7 +1629,15 @@ pub fn resolv( lac2t: false, }; - let lucy_output = lucy(&lucy_config, &lucy_model, &opacfl_data, &rad1_data); + // Lucy 单步温修。ilam (1-based) = Fortran 内层 ILUCY；lucy_ng 跨迭代持久化 Ng 历史。 + let lucy_output = lucy( + &lucy_config, + &lucy_model, + &opacfl_data, + &rad1_data, + ilam, + &mut lucy_ng, + ); // ============================================================== // Step 5: 更新模型状态 @@ -1424,17 +1725,27 @@ pub fn resolv( nd, dens: ¶ms.model.modpar.dens[..nd], deldm: &deldm, - dedm1: 0.0, + dedm1: params.model.modpar.dm[0] / params.model.modpar.dens[0].max(1e-300), abrosd: &mut abrosd_mut, abplad: ¶ms.model.opmean.abplad[..nd], taurs: &mut taurs_mut, reint: &mut reint_mut, redif: &mut redif_mut, - taudiv: 10.0, + // TAUDIV / IDLST: NSTPAR defaults (nstpar.rs:589-590). The hhe test uses + // default optional parameters, so taudiv=0.5 and idlst=5. + // NOTE: taudiv=0.5 (NOT 10.0!) is critical — with taudiv=10 the Rosseland + // optical depth never reaches 10, so REDIF=1 only at the single deepest + // point (idr=nd-1), which injects the differential RE (SIG4P*TEFF^4 term) + // at one layer only and drives the bottom temperature to runaway (T~1e5 K). + taudiv: 0.5, iter, - itndre: 0, + itndre: 9999, // Fortran: ITNDRE=NITER — run ROSSTD for all iterations ndre: params.tlusty_config.matkey.ndre, - idlst: nd - 1, + // IDLST=5: NSTPAR default (nstpar.rs). Now safe to use because the Rosseland + // mean IS accumulated (idr lands at a mid-depth ~49, not nd-1), so the deepest + // IDLST points correctly use pure diffusion (REDIF=1, REINT=0) instead of the + // integral radiative-equilibrium equation. + idlst: 5, teff, lfin, temp: ¶ms.model.modpar.temp[..nd], @@ -1452,6 +1763,20 @@ pub fn resolv( params.model.repart.reint[id] = reint_mut[id]; params.model.repart.redif[id] = redif_mut[id]; } + + // DIAGNOSTIC: Rosseland optical depth / opacity scale (root-cause of redif scheme) + if iter == 1 { + let tmin = taurs_mut[..nd].iter().cloned().fold(f64::INFINITY, f64::min); + let tmax = taurs_mut[..nd].iter().cloned().fold(0.0_f64, f64::max); + let abmin = abrosd_mut[..nd].iter().cloned().fold(f64::INFINITY, f64::min); + let abmax = abrosd_mut[..nd].iter().cloned().fold(0.0_f64, f64::max); + let nredif = (0..nd).filter(|&i| redif_mut[i] > 0.0).count(); + eprintln!(" ROSSTD diag iter={}: taurs[min]={:.3e} taurs[max]={:.3e} (taudiv=0.5) | abrosd[min]={:.3e} [max]={:.3e} | #redif>0={}/{}", iter, tmin, tmax, abmin, abmax, nredif, nd); + eprintln!(" taurs[0]={:.3e} taurs[35]={:.3e} taurs[69]={:.3e} | dens[0]={:.3e} dens[69]={:.3e} | dm[0]={:.3e} dm[69]={:.3e}", + taurs_mut[0], taurs_mut[35.min(nd-1)], taurs_mut[nd-1], + params.model.modpar.dens[0], params.model.modpar.dens[nd-1], + params.model.modpar.dm[0], params.model.modpar.dm[nd-1]); + } } // FCOOL(ID) = REINT(ID)*FCOOLI(ID) @@ -1803,18 +2128,21 @@ pub fn resolv( { let nd_ali = config.nd; - // 将 ALI 输出字段清零（准备下一次迭代） + // NOTE: Do NOT zero fcooli/fcool here — SOLVES needs them! + // They are properly reset at the start of Step 3c (ALIFR1 loop). + // Zeroing here was destroying the radiative equilibrium terms + // that SOLVES/BRE use for temperature corrections. for id in 0..nd_ali { - params.model.totflx.fcooli[id] = 0.0; + // params.model.totflx.fcooli[id] = 0.0; // ← needed by SOLVES/BRE params.model.totflx.flfix[id] = 0.0; params.model.totflx.flexp[id] = 0.0; params.model.totflx.fprd[id] = 0.0; params.model.totflx.flrd[id] = 0.0; - params.model.totflx.fcool[id] = 0.0; + // params.model.totflx.fcool[id] = 0.0; // ← needed by SOLVES/BRE params.model.pressr.pradt[id] = 0.0; params.model.pressr.prada[id] = 0.0; - params.model.opmean.abrosd[id] = 0.0; - params.model.opmean.sumdpl[id] = 0.0; + // params.model.opmean.abrosd[id] = 0.0; // ← needed by SOLVES + // params.model.opmean.sumdpl[id] = 0.0; // ← may be needed } if use_kant || lfin { @@ -2017,7 +2345,11 @@ pub fn resolv( debug_log!("RESOLV: NEWPOP at depth {}", id); // ELCOR — 电子密度修正 - if !config.lchc && iter < config.ielcor { + // Skip when INPC active: ELEC is an independent SOLVES variable + // (BPOPC); ELCOR's full-Saha charge neutrality would overwrite the + // BPOPC value and cause ELEC to oscillate between the two Saha + // implementations each iteration. + if !config.lchc && iter < config.ielcor && !inpc_active { let t = params.model.modpar.temp[id]; let elec = params.model.modpar.elec[id]; let dens = params.model.modpar.dens[id]; diff --git a/src/tlusty/io/start.rs b/src/tlusty/io/start.rs index 54f5225..823a42f 100644 --- a/src/tlusty/io/start.rs +++ b/src/tlusty/io/start.rs @@ -36,6 +36,7 @@ //! ``` use super::FortranReader; +use super::initia::{self, InitiaParams, FrequencyGridParams, InitiaConfig}; use crate::tlusty::math::{comset, ComsetParams}; use crate::tlusty::state::config::TlustyConfig; use crate::tlusty::state::atomic::AtomicData; @@ -80,6 +81,7 @@ impl StartCallbacks for NoOpStartCallbacks { /// ```fortran /// common/hediff/ hcmass,radstr /// ``` +/// 加上频率网格参数（由 INITIA 使用）。 #[derive(Debug, Clone)] pub struct StartConfig { /// 盘模型标志 (0=大气, 1=盘) @@ -89,6 +91,19 @@ pub struct StartConfig { pub hcmass: f64, /// 恒星半径 (RADSTR) pub radstr: f64, + // --- INITIA 频率网格参数 --- + /// 频率点数 + pub nfreq: usize, + /// 最低频率 (Hz) — 对应 Fortran FRCMIN + pub frmin: f64, + /// 最高频率 (Hz) — 对应 Fortran FRCMAX + pub frmax: f64, + /// 频率设置标志 (对应 Fortran IFRSET) + pub ifrset: i32, + /// Compton 散射开关 (0=关闭) + pub icompt: i32, + /// Klein-Nishina 近似阶数 + pub knish: i32, } impl Default for StartConfig { @@ -97,6 +112,12 @@ impl Default for StartConfig { idisk: 0, hcmass: 0.0, radstr: 0.0, + nfreq: 50, + frmin: 1.0e12, + frmax: 2.8e16, + ifrset: 0, + icompt: 0, + knish: 0, } } } @@ -177,10 +198,53 @@ pub fn start_with_callbacks( params.tlusty_config.basnum.idisk = config.idisk; // ======================================== - // Step 2: 调用 INITIA + // Step 2: 调用 INITIA — 生成频率网格 // 对应 Fortran: call initia // ======================================== - callbacks.call_initia(); + let nd = if params.tlusty_config.basnum.nd > 0 { + params.tlusty_config.basnum.nd as usize + } else { + config.nfreq // fallback + }; + let teff = params.tlusty_config.inppar.teff; + let grav_log = params.tlusty_config.inppar.grav; + let lte = params.tlusty_config.inppar.lte; + let nfreq = if config.nfreq > 0 { config.nfreq } else { 50 }; + let frmax = if config.frmax > 0.0 { + config.frmax + } else { + 8.0e11 * teff // Fortran: FRCMAX = 8e11 * TEFF + }; + + let initia_params = InitiaParams { + config: InitiaConfig::default(), + teff, + grav: grav_log, + lte, + ltgrey: false, + vtb: 0.0, + ipturb: 0, + nd, + }; + let grid_params = FrequencyGridParams { + frmin: config.frmin, + frmax, + nfreq, + ifrset: config.ifrset, + }; + let initia_output = initia::initia(&initia_params, &grid_params, config.icompt, config.knish); + + // 存储频率网格到 ModelState + let nfreq_actual = initia_output.freq_grid.freq.len(); + for ij in 0..nfreq_actual { + params.model.frqall.freq[ij] = initia_output.freq_grid.freq[ij]; + params.model.frqall.w[ij] = initia_output.freq_grid.w[ij]; + } + // 更新配置中的 nfreq + params.tlusty_config.basnum.nfreq = nfreq_actual as i32; + + eprintln!(" INITIA: nfreq={}, frmin={:.3e}, frmax={:.3e}", + nfreq_actual, config.frmin, frmax); // ======================================== // Step 3: 可选调用 HEDIF（He 扩散） diff --git a/src/tlusty/main.rs b/src/tlusty/main.rs index 84643c7..58f31c1 100644 --- a/src/tlusty/main.rs +++ b/src/tlusty/main.rs @@ -50,7 +50,7 @@ //! END //! ``` -use std::io::{self}; +use std::io::{self, BufRead, BufReader}; use std::time::Instant; use std::fs::File; @@ -59,7 +59,7 @@ use super::io::{ start, StartConfig, StartParams, resolv, ResolvConfig, ResolvParams, read_tlusty_model, InpmodParams, - input::InputParser, + input::{InputParser, IonParams}, inpmod::InputModelData, }; use super::math::solvers::{ @@ -71,6 +71,7 @@ use super::state::config::TlustyConfig; use super::state::atomic::AtomicData; use super::state::model::ModelState; use super::state::arrays::ComputeArrays; +use super::math::continuum::{lte_meanopt, generate_lte_frequency_grid, LteOpacityParams}; // ============================================================================ // 常量 (从 Fortran 移植) @@ -224,6 +225,252 @@ pub struct ScratchFiles { // 主运行函数 // ============================================================================ +// ============================================================================ +// RDATA 原子数据初始化 +// ============================================================================ + +/// 从离子数据文件读取能级数据并填充 AtomicData。 +/// +/// 对应 Fortran INITIA 中的离子循环 + RDATA(ION) 调用。 +/// 处理: +/// 1. 读取每个离子的 .dat 文件（能级能量、统计权重、主量子数） +/// 2. 能量转换: 0 → 氢公式(EH*Z²/n²), cm⁻¹ → erg +/// 3. 填充 LevPar (enion, g, nquant, iel, iatm) +/// 4. 填充 IonPar (nfirst, nlast, nnext, iz, charg2, ff) +/// 5. 填充 AtoPar (numat, n0a, nka) +/// 6. 初始化 ILTREF 指向连续态能级 +fn populate_atomic_data( + atomic: &mut AtomicData, + model: &mut ModelState, + ions: &[IonParams], + _nd: usize, +) -> usize { + use crate::tlusty::state::constants::{EH, H, MLEVEL, MION}; + + const H_PLANCK: f64 = 6.6260755e-27; // Planck constant (erg·s) + const CM1_TO_ERG: f64 = 1.9857e-16; // cm⁻¹ → erg + + let mut ilev: usize = 0; // current level position (0 initially, Fortran convention) + let mut iatlst: usize = 0; + let mut ia: usize = 0; + let mut natom: usize = 0; + let mut nki: usize = 0; // tracks NNEXT of last ion (for NLEVEL) + + for (ion_idx, ion) in ions.iter().enumerate() { + if ion_idx >= MION { break; } + + if ion.ilast == 0 { + // === 主离子定义 (ilast=0) === + // Fortran INITIA logic for NFIRST: + // IF(IATI(ION).EQ.IATLST) THEN → NFIRST = ILEV (same atom) + // ELSE → NFIRST = ILEV + 1 (new atom) + let n0i = if ion.iat == iatlst { + ilev // same atom: continues from previous ion's NNEXT + } else { + ilev + 1 // new atom: skip to next slot + }; + let n1i = n0i + ion.nlevs - 1; + nki = n1i + 1; + + // 离子参数 (1-based indices) + atomic.ionpar.nfirst[ion_idx] = n0i as i32; + atomic.ionpar.nlast[ion_idx] = n1i as i32; + atomic.ionpar.nnext[ion_idx] = nki as i32; + atomic.ionpar.iz[ion_idx] = ion.iz as i32 + 1; // Z = iz + 1 + let izz = atomic.ionpar.iz[ion_idx] as f64; + atomic.ionpar.charg2[ion_idx] = izz * izz; + + // 离子数据索引 + atomic.iondat.iati[ion_idx] = ion.iat as i32; + atomic.iondat.izi[ion_idx] = ion.iz as i32; + atomic.iondat.nlevs[ion_idx] = ion.nlevs as i32; + + // 原子参数 + if ion.iat != iatlst { + ia = ion.iat - 1; + atomic.atopar.n0a[ia] = n0i as i32; + atomic.atopar.numat[ia] = ion.iat as i32; + iatlst = ion.iat; + natom = natom.max(ia + 1); + } + + // 读取离子数据文件 + let filei = ion.filei.trim(); + if !filei.is_empty() && filei != "''" && filei != "'" { + if let Some(levels) = read_level_data(filei, ion.nlevs) { + for (il, (e_raw, g_raw, nq_raw)) in levels.iter().enumerate() { + let i = n0i + il; // 1-based level index + if i > MLEVEL { break; } + let i0 = i - 1; // 0-based for array access + + let iq = (il + 1) as i32; + let x = (iq * iq) as f64; + + // 能量转换 (Fortran RDATA 逻辑) + let e_abs = e_raw.abs(); + let e0 = if e_abs < 1e-20 { + EH * izz * izz / x + } else if e_abs > 1e-7 && e_abs < 100.0 { + 1.6018e-12 * e_abs + } else if e_abs > 100.0 && e_abs < 1e7 { + CM1_TO_ERG * e_abs + } else if e_abs >= 1e7 { + H_PLANCK * e_abs + } else { + e_abs + }; + + atomic.levpar.enion[i0] = if *e_raw >= 0.0 { e0 } else { -e0 }; + atomic.levpar.g[i0] = if *g_raw == 0.0 { 2.0 * x } else { *g_raw }; + atomic.levpar.nquant[i0] = if *nq_raw == 0 { iq } else { *nq_raw }; + atomic.levpar.iel[i0] = ion_idx as i32; + atomic.levpar.iatm[i0] = ia as i32; + atomic.levpar.ilk[i0] = 0; + } + } else { + eprintln!("WARNING: Cannot read data file '{}'", filei); + } + } else { + for il in 0..ion.nlevs { + let i = n0i + il; + if i > MLEVEL { break; } + let i0 = i - 1; + let iq = (il + 1) as i32; + let x = (iq * iq) as f64; + atomic.levpar.enion[i0] = EH * izz * izz / x; + atomic.levpar.g[i0] = 2.0 * x; + atomic.levpar.nquant[i0] = iq; + atomic.levpar.iel[i0] = ion_idx as i32; + atomic.levpar.iatm[i0] = ia as i32; + } + } + + ilev = nki; // Fortran: ILEV = NNEXT(ION) — NOT incremented further + + } else if ion.ilast > 0 { + // === 连续态能级 (ilast>0) === + // Fortran: ENION(ILEV)=0, G(ILEV)=ILASTI — at ILEV (NNEXT of prev ion) + // ILEV is NOT incremented after this! + let i = ilev; // = NNEXT of previous main ion + if i >= 1 && i <= MLEVEL { + let i0 = i - 1; + atomic.levpar.enion[i0] = 0.0; + atomic.levpar.g[i0] = ion.ilast as f64; + atomic.levpar.nquant[i0] = 1; + atomic.levpar.iel[i0] = (ion_idx - 1) as i32; + atomic.levpar.iatm[i0] = ia as i32; + atomic.levpar.ilk[i0] = ion_idx as i32; + + // Update atom's last level + atomic.atopar.nka[ia] = i as i32; + + // Update previous ion's nnext to point here (it already does) + let prev_ion = ion_idx - 1; + if prev_ion < MION { + // FF(ion) = ionization frequency for previous ion + let nlast_prev = atomic.ionpar.nlast[prev_ion] as usize; + let nff = if nlast_prev >= 1 { atomic.levpar.nquant[nlast_prev - 1] + 1 } else { 2 }; + if nff > 0 { + let iz_prev = atomic.ionpar.iz[prev_ion] as f64; + atomic.ionpar.ff[prev_ion] = EH / H_PLANCK * iz_prev * iz_prev + / (nff as f64) / (nff as f64); + } + } + } + // Fortran: ILEV stays the same — no increment! + + } else { + // === 终止符 (ilast=-1) === + break; + } + } + + // NLEVEL = NKI (Fortran: NLEVEL=NKI, where NKI=NNEXT of last main ion) + let nlevel = nki; + + // 初始化 ILTREF: 每个能级指向其离子的连续态能级 + // 对应 Fortran LEVSET 的简化版 — 所有能级默认 LTE + for ion_idx in 0..ions.len() { + if ion_idx >= MION { break; } + let ion = &ions[ion_idx]; + if ion.ilast != 0 { continue; } + + let n0i = atomic.ionpar.nfirst[ion_idx] as usize; + let n1i = atomic.ionpar.nlast[ion_idx] as usize; + let nnext = atomic.ionpar.nnext[ion_idx] as usize; + + for i in n0i..=n1i { + if i <= MLEVEL && nnext <= MLEVEL { + let i0 = i - 1; // 0-based + // ILTREF(i, all depths) = nnext (1-based) + for id in 0..model.modpar.temp.len().min(crate::tlusty::state::constants::MDEPTH) { + if model.modpar.temp[id] > 0.0 { + model.levref.iltref[i0][id] = nnext as i32; + } + } + } + } + } + + eprintln!(" RDATA: nlevel={}, natom={}", nlevel, natom); + for ion_idx in 0..ions.len() { + let ion = &ions[ion_idx]; + if ion.ilast != 0 { continue; } + let n0 = atomic.ionpar.nfirst[ion_idx]; + let n1 = atomic.ionpar.nlast[ion_idx]; + let nk = atomic.ionpar.nnext[ion_idx]; + let i0 = (n0 - 1) as usize; // 0-based for array access + eprintln!(" ion {}: {}..{}..{} iz={} enion[{}]={:.3e} g[{}]={:.1} nq[{}]={}", + ion_idx, n0, n1, nk, atomic.ionpar.iz[ion_idx], + n0, atomic.levpar.enion[i0], + n0, atomic.levpar.g[i0], + n0, atomic.levpar.nquant[i0]); + } + + nlevel +} + +/// 从 .dat 文件读取能级数据。 +/// +/// 文件格式: +/// ```text +/// H I energy levels ← 标题行 +/// 0.00000 2 1 1s 0 ← energy g nquant name ifwop +/// 82258.955 2 2 2s 0 +/// ... +/// * ← 分隔符 +/// ``` +pub(crate) fn read_level_data(path: &str, nlevs: usize) -> Option> { + let file = File::open(path).ok()?; + let reader = BufReader::new(file); + let mut lines = reader.lines(); + + // 跳过标题行 + lines.next()?; + + let mut levels = Vec::with_capacity(nlevs); + for _ in 0..nlevs { + let line = loop { + let l = lines.next()?.ok()?; + let trimmed = l.trim(); + if !trimmed.is_empty() && !trimmed.starts_with('*') && !trimmed.starts_with('!') { + break l; + } + }; + + let parts: Vec<&str> = line.split_whitespace().collect(); + if parts.len() < 3 { continue; } + + let energy: f64 = parts[0].parse().ok()?; + let g: f64 = parts[1].parse().ok()?; + let nquant: i32 = parts[2].parse().ok()?; + + levels.push((energy, g, nquant)); + } + + Some(levels) +} + /// 运行 TLUSTY 计算。 /// /// 严格按照 Fortran tlusty.f 的逻辑流程实现。 @@ -274,6 +521,70 @@ pub fn run_tlusty( eprintln!(" Input: TEFF={}, LOGG={}, LTE={}, NFREQ={}", teff, grav_log, lte, nfreq); + // 诊断 [TLUSTY_KR_POPDIAG]: 在 gold 收敛模型 (fort.8) 的 T/ρ/Ne 条件下, + // 对比 gold 实际 populations 与简化 LTE populations 算出的局部 κ_R。 + // 判定 κ_R 偏差是否由种群驱动 → 移植完整 ELDENS/STEQEQ 是否有益。只读诊断, 无副作用。 + if std::env::var("TLUSTY_KR_POPDIAG").is_ok() { + let fort8_path = "fort.8"; + if std::path::Path::new(fort8_path).exists() { + let f8 = File::open(fort8_path).expect("open fort.8"); + let mut r8 = FortranReader::new(io::BufReader::new(f8)); + let inpmod_params = InpmodParams { + intrpl: 0, idisk: 0, ifmol: 0, tmolim: 1e10, teff, + }; + if let Ok(gold) = read_tlusty_model(&mut r8, &inpmod_params) { + let wmm0 = 1.3 * crate::tlusty::state::constants::HMASS; + let pops: Vec> = gold.popul.unwrap_or_default(); + crate::tlusty::io::ltegr::diag_grey_kr_population_effect( + teff, wmm0, &gold.temp, &gold.elec, &gold.dens, &pops, + ); + } else { + eprintln!(" [TLUSTY_KR_POPDIAG] failed to read fort.8"); + } + } else { + eprintln!(" [TLUSTY_KR_POPDIAG] fort.8 not found"); + } + return TlustyResult { + total_iterations: 0, + converged: false, + total_time_secs: start_time.elapsed().as_secs_f64(), + }; + } + + // 诊断 [TLUSTY_KR_FORCOND]: 在 Fortran-grey 模型自身 (T,ρ,Ne) 条件下计算 + // Rust κ_R + ne, 对比 Fortran 精确 κ_R=dτ/dM + ne。隔离 opacity/EOS 公式 vs 网格反馈。 + // 输入文件格式 (每行): id T ne dens tauR mass kr_fort (取第 2/3/4/7 列)。 + if let Ok(path) = std::env::var("TLUSTY_KR_FORCOND") { + if std::path::Path::new(&path).exists() { + let txt = std::fs::read_to_string(&path).unwrap_or_default(); + let mut conds: Vec<(f64, f64, f64, f64)> = Vec::new(); + for ln in txt.lines() { + let t = ln.split_whitespace().collect::>(); + if t.len() < 7 || t[0].starts_with('#') { + continue; + } + let toks: Vec = t[1..7] + .iter() + .filter_map(|s| s.parse::().ok()) + .collect(); + // toks = [T, ne, dens, tauR, mass, kr_fort] (6 values from cols 2..7) + if toks.len() >= 6 { + conds.push((toks[0], toks[1], toks[2], toks[5])); + } + } + eprintln!(" [TLUSTY_KR_FORCOND] parsed {} depth conditions from {}", conds.len(), path); + let wmm0 = 1.252 * crate::tlusty::state::constants::HMASS; + crate::tlusty::io::ltegr::diag_grey_kr_fortran_conditions(teff, wmm0, &conds); + } else { + eprintln!(" [TLUSTY_KR_FORCOND] conditions file not found: {}", path); + } + return TlustyResult { + total_iterations: 0, + converged: false, + total_time_secs: start_time.elapsed().as_secs_f64(), + }; + } + // ======================================== // Step 2: Read or create initial model // @@ -288,7 +599,30 @@ pub fn run_tlusty( let model_data = if input_params.ltgrey { // Grey atmosphere: create initial model from Eddington approximation // T(τ) = Teff * (3/4 * (τ + 2/3))^0.25 - create_grey_atmosphere(teff, grav) + // + // EXPERIMENT (env TLUSTY_GREY_REFGRID): build the grey atmosphere + // (Eddington T + LTE Ne/ρ) on the GOLD-REF DM grid read from fort.8, + // instead of deriving DM from the (too-high) grey opacity. Tests whether + // SOLVES converges to the gold ref when given the correct depth grid. + // Does NOT affect readmodel (LTGREY=false) path. + if std::env::var("TLUSTY_GREY_REFGRID").is_ok() { + eprintln!(" [TLUSTY_GREY_REFGRID] grey atmosphere on gold-ref DM grid"); + match create_grey_atmosphere_on_ref_grid(teff, grav) { + Ok(md) => md, + Err(e) => { + eprintln!(" [TLUSTY_GREY_REFGRID] failed ({}); falling back to computed grid", e); + create_grey_atmosphere_via_ltegr(teff, grav) + } + } + } else if std::env::var("TLUSTY_GREY_SIMPLE").is_ok() { + // 旧简化 Eddington 路径 (lte_meanopt 有补偿误差, 仅供对比). 默认不再使用. + create_grey_atmosphere(teff, grav) + } else { + // 默认: 完整 LTEGR 路径 (忠实移植 Fortran) — Opacf0State 不透明度 + 正确的 + // Rosseland 平均 (bnue=c1·ν³ prefactor). 外层与深层与 gold ref 匹配良好; + // 部分电离区仍因简化 LTE 种群偏高 (待移植完整 ELDENS/STEQEQ). + create_grey_atmosphere_via_ltegr(teff, grav) + } } else { // Non-grey: read model from fort.8 let fort8_path = "fort.8"; @@ -341,7 +675,15 @@ pub fn run_tlusty( // The real WMM will be recomputed after START populates atomic data. let atom_modes: Vec = input_params.atoms.iter().map(|a| a.mode).collect(); use crate::tlusty::state::constants::HMASS; - let mut wmm = HMASS; // Default: pure hydrogen mean molecular weight + // WMM = mean mass per NUCLEUS. The H-He atmosphere uses mass fractions + // X_H=0.70, Y_HE=0.28 with ~2% metals adding mass but negligible nuclei, + // so nuclei/g = (X/1 + Y/4)/m_H = 0.77/m_H ⇒ WMM = m_H/0.77 = 2.174e-24 g. + // This matches the gold reference: `dens/(totn−elec)` = 2.175e-24 at every + // depth. The previous default HMASS (pure H, 1.673e-24) corrupted the + // post-SOLVES density recompute dens=(totn−elec)·wmm by ~23% per iteration. + const X_H_FRAC: f64 = 0.70; + const Y_HE_FRAC: f64 = 0.28; + let mut wmm = HMASS / (X_H_FRAC + Y_HE_FRAC / 4.0); eprintln!(" WMM={:.5e} (default, before START)", wmm); // Copy model data into ModelState @@ -380,7 +722,21 @@ pub fn run_tlusty( config.inppar.grav = grav; config.inppar.lte = lte; config.prints.iprinp = 1; - config.runkey.niter = 0; + config.runkey.niter = 30; // Default; updated after reading input + + // ======================================== + // Step 2a: RDATA — 读取原子数据文件 (Fortran INITIA → RDATA loop) + // + // 从 fort.5 中解析的离子列表，读取各离子的 .dat 文件， + // 填充 AtomicData 的能级能量、统计权重、主量子数等。 + // ======================================== + let nlevel_rdata = populate_atomic_data( + &mut atomic, &mut model, &input_params.ions, nd, + ); + if nlevel_rdata > 0 && nlevel_from_ions > 0 { + // Use RDATA-derived nlevel if available + eprintln!(" RDATA nlevel={} matches ion input nlevel={}", nlevel_rdata, nlevel_from_ions); + } // ======================================== // Step 2b: START initialization (Fortran: CALL START) @@ -397,6 +753,9 @@ pub fn run_tlusty( start_config.ifrset = 0; start_config.icompt = 0; start_config.knish = 0; + eprintln!(" BEFORE START: T[0]={:.1}, T[35]={:.1}, Ne[0]={:.3e}, Dens[0]={:.3e}", + model.modpar.temp[0], model.modpar.temp[35], + model.modpar.elec[0], model.modpar.dens[0]); let mut start_params = StartParams { config: &mut start_config, tlusty_config: config, @@ -408,6 +767,9 @@ pub fn run_tlusty( eprintln!("WARNING: START initialization had issues"); } eprintln!(" START: nn={}", start_output.nn); + eprintln!(" AFTER START: T[0]={:.1}, T[35]={:.1}, Ne[0]={:.3e}, Dens[0]={:.3e}", + model.modpar.temp[0], model.modpar.temp[35], + model.modpar.elec[0], model.modpar.dens[0]); } // Recompute WMM from populated atomic data (after START) @@ -426,8 +788,10 @@ pub fn run_tlusty( // For grey start models, level populations are all zero. // Compute LTE populations from Saha-Boltzmann equations. // This matches Fortran INITIA → STATE0 → RDATA → LTE population setup. + // IMPORTANT: Only for grey start! When reading fort.8, populations are + // already set and should NOT be overwritten. // ======================================== - if nlevel >= 39 { + if input_params.ltgrey && nlevel >= 39 { // Use abundance parameters matching Fortran NSTPAR defaults: // H abundance = 1.0, He abundance = 0.0851 (solar), YTOT = 1.08588 let ytot = 1.0 + 0.0851; @@ -449,20 +813,12 @@ pub fn run_tlusty( // Step 2c: CORRWM — frequency selection for ALI // // Determines which frequency points are explicit (IJFR) vs ALI. - // In Fortran, INITIA calls CORRWM. In Rust, we call it explicitly. + // IMPORTANT: Must run AFTER first RESOLV creates the full frequency grid, + // otherwise CORRWM sees an incomplete grid and selects wrong frequencies. + // CORRWM is called after the first RESOLV below (Step 2d). // ======================================== - { - use crate::tlusty::math::opacity::{corrwm, CorrwmParams}; - let mut corrwm_params = CorrwmParams { - basnum: &mut config.basnum, - trapar: &atomic.trapar, - frqall: &mut model.frqall, - freaux: &mut model.freaux, - phoexp: &mut model.phoexp, - }; - corrwm(&mut corrwm_params); - eprintln!(" CORRWM: nfreqe={} explicit frequencies selected", model.frqall.nfreqe); - } + + // (CORRWM placeholder — actual call moved to after first RESOLV) // ======================================== // Step 3: Iteration loop @@ -471,10 +827,16 @@ pub fn run_tlusty( // NITER>0 → RESOLV + SOLVES loop until convergence or max iterations // ======================================== let niter: i32 = if std::env::var("TLUSTY_NITER").is_ok() { + // Environment variable takes precedence (allows explicit NITER=0) std::env::var("TLUSTY_NITER").unwrap().parse().unwrap_or(0) - } else { + } else if input_params.niter > 0 { input_params.niter + } else { + // Fortran default: NITER=30 when NSTPAR=0 (no PVALUE array in input) + 30 }; + config.runkey.niter = niter; + eprintln!(" NITER={} (input={}, default=30)", niter, input_params.niter); let itlucy: i32 = std::env::var("TLUSTY_ITLUCY") .ok() .and_then(|v| v.parse().ok()) @@ -530,7 +892,7 @@ pub fn run_tlusty( ntrans: 0, teff, irder: 0, - niter: 0, + niter, ibc: 3, itlucy, }; @@ -547,6 +909,44 @@ pub fn run_tlusty( let _resolv_result = resolv(&mut resolv_params, None::<&mut FortranWriter>); + // ======================================== + // Step 2d: CORRWM — AFTER first RESOLV created the full frequency grid + // + // CORRWM selects explicit frequency points for the SOLVES matrix. + // Must run after RESOLV because RESOLV's INILAM/INIFRC creates the + // actual frequency grid (e.g. 144 points). + // ======================================== + { + let nfreq_grid = model.frqall.freq.iter().take_while(|&&f| f > 0.0).count(); + config.basnum.nfreq = nfreq_grid as i32; + config.basnum.nd = nd as i32; + config.basnum.ntrans = 0; + eprintln!(" CORRWM: basnum.nfreq={}, nd={}", nfreq_grid, nd); + + use crate::tlusty::math::opacity::{corrwm, CorrwmParams}; + let mut corrwm_params = CorrwmParams { + basnum: &mut config.basnum, + trapar: &atomic.trapar, + frqall: &mut model.frqall, + freaux: &mut model.freaux, + phoexp: &mut model.phoexp, + }; + corrwm(&mut corrwm_params); + eprintln!(" CORRWM: nfreqe={} explicit (of {} total)", + config.basnum.nfreqe, nfreq_grid); + } + + // Save INIFRC edge frequency indices for SOLVES (before CORRWM overwrites) + // INIFRC selects ~9 ionization edge frequencies which provide the optimal + // coupling between the radiation field and temperature equation. + // CORRWM selects all 144 frequencies (overkill, makes 146×146 matrix ill-conditioned). + let inifrc_ijfr: Vec = model.frqall.ijfr_explicit.clone(); + let inifrc_nfreqe = inifrc_ijfr.len(); + eprintln!(" INIFRC edges: {} explicit frequencies saved for SOLVES", inifrc_nfreqe); + + // Save CORRWM frequency count for reference + let _corrwm_nfreqe: usize = config.basnum.nfreqe as usize; + // ======================================== // Step 4: NITER>0 iteration loop (RESOLV + SOLVES) // ======================================== @@ -563,21 +963,50 @@ pub fn run_tlusty( if use_solves { // Full linearization mode (RESOLV + SOLVES) // - // NFREQE determines the matrix structure: - // NFREQE=0 (default for LTE with NFRECL=0): NN=2 (TOTN+T only), - // all radiative equilibrium via ALI integral (FCOOL/REIT) - // NFREQE>0: NN=NFREQE+2, explicit Jν equations + TOTN + T - // NFREQE: number of explicit frequency points in the SOLVES matrix. - // Set by CORRWM (called in Step 2c): selects frequencies near ionization edges. - // CORRWM stores nfreqe and ijfr_explicit in model.frqall. - let nfreqe: usize = if model.frqall.nfreqe > 0 { - model.frqall.nfreqe + // Rybicki method (TLUSTY_RYBICKI=1): NN=2, no explicit Jν. + // The Jν contribution is captured via the Rybicki integral (FCOOL/REIT). + // This gives a 2×2 block-tridiagonal system that is well-conditioned. + // + // Full method (default): NN=146, all 144 Jν explicit. + // The 146×146 system is ill-conditioned — Jν contamination of T row + // during forward elimination kills temperature corrections after iter 1. + let nfreq_grid = model.frqall.freq.iter().take_while(|&&f| f > 0.0).count(); + let use_rybicki = std::env::var("TLUSTY_RYBICKI").is_ok(); + // INPC experiment (TLUSTY_INPC=1): add electron density as a 3rd + // constraint variable (NN = NFREQE + 3). Adds the BPOPC charge- + // conservation row + the BRE diffusion ELEC coupling column + // (matgen_lte.rs), regularizing the near-singular deepest-point T + // block. Gated — default TLUSTY_SOLVES keeps NN = NFREQE + 2. + let use_inpc = std::env::var("TLUSTY_INPC").is_ok(); + // Gold-configuration explicit set (TLUSTY_NFREQE_EDGES): use only the + // INIFRC ionization-edge frequencies (~9, matching the Fortran gold + // reference whose fortran.6 lists "FREQUENCY POINTS AND WEIGHTS - + // EXPLICIT" with NFREQE=9) as linearized Jν unknowns. The remaining + // ~135 frequencies are treated implicitly — bre_lte folds their + // radiative-equilibrium residual into FCOOL (fcool_integral over the + // non-explicit set) and bhe_lte's GRD sums only the explicit edges, + // exactly matching Fortran BHE/BRED `DO IJ=1,NFREQE`. + // + // This is the gold-matching split. The default all-explicit path + // (NFREQE=144, NN=146) is overkill: it linearizes 135 extra surface + // high-frequency Jν (Wien tail, i=128-133 at id=0) whose spurious + // oscillations pin chmx≈0.72 and prevent lfin. The 9-edge matrix + // (NN≈11-13) is small and well-conditioned, and those oscillating + // points become implicit. + let use_edges = std::env::var("TLUSTY_NFREQE_EDGES").is_ok(); + let (ijfr_solve, nfreqe, nn) = if use_rybicki { + (vec![], 0, 2) // Rybicki: no explicit Jν, NN=2 (TOTN + T) + } else if use_edges && !inifrc_ijfr.is_empty() { + let ijfr = inifrc_ijfr.clone(); + let base = ijfr.len() + 2; + (ijfr.clone(), ijfr.len(), if use_inpc { base + 1 } else { base }) } else { - eprintln!("WARNING: NFREQE not set by CORRWM, using fallback"); - 54.min(nfreq as usize) + let ijfr: Vec = (0..nfreq_grid).collect(); + let base = ijfr.len() + 2; + (ijfr.clone(), ijfr.len(), if use_inpc { base + 1 } else { base }) }; - eprintln!("SOLVES: using NFREQE={} (LTE={}, NN={})", nfreqe, lte, nfreqe + 2); - let nn = nfreqe + 2; + model.frqall.ijfr_explicit = ijfr_solve.clone(); + eprintln!("SOLVES: NFREQE={}/{} (rybicki={}, LTE={}, NN={})", nfreqe, nfreq_grid, use_rybicki, lte, nn); let orelax = 1.0; let dpsilg = 10.0; // Fortran default (nstpar.f DATA statement) let dpsilt = 1.25; // Temperature damping (Fortran default) @@ -618,12 +1047,16 @@ pub fn run_tlusty( { let nhe = nfreqe; let nre = nfreqe + 1; + let npc = nfreqe + 2; for id in 0..(nd as usize) { for ij in 0..nfreqe { accel_psy0[ij][id] = model.expraf.radex[ij][id]; } accel_psy0[nhe][id] = model.modpar.totn[id]; accel_psy0[nre][id] = model.modpar.temp[id]; + if npc < nn { + accel_psy0[npc][id] = model.modpar.elec[id]; + } } let solves_output = solves_lte( @@ -647,6 +1080,28 @@ pub fn run_tlusty( eprintln!(" NITER 1/{}: chmx={:.3e}, chmt={:.3e}, lfin={}", niter, solves_output.chmx, solves_output.chmt, solves_output.lfin); + // After SOLVES updates TOTN and T, recompute DENS from the + // (authoritative) level populations: DENS = (total nuclei + // density) * WMM, where the nuclei density is Σ popul[j][id] + // and WMM is the mean mass per nucleus (incl. metals). This is + // the mass density by definition and matches the converged + // Fortran DENS=WMM*(TOTN-ELEC). Reading TOTN directly is less + // robust: the SOLVES linear solve can shift TOTN independently + // of the populations (which stay at the EOS-consistent values), + // leaking a spurious density drift into DENS. + let nlvl = config.basnum.nlevel as usize; + for id in 0..(nd as usize) { + let sumpop: f64 = (0..nlvl).map(|j| model.levpop.popul[j][id]).sum(); + let dens_new = if sumpop > 0.0 { + sumpop * wmm + } else { + (model.modpar.totn[id] - model.modpar.elec[id]) * wmm + }; + if dens_new > 0.0 { + model.modpar.dens[id] = dens_new; + } + } + if solves_output.lfin { converged = true; } @@ -716,12 +1171,16 @@ pub fn run_tlusty( // Build PSY0 from current model state let nhe = nfreqe; let nre = nfreqe + 1; + let npc = nfreqe + 2; for id in 0..(nd as usize) { for ij in 0..nfreqe { accel_psy0[ij][id] = model.expraf.radex[ij][id]; } accel_psy0[nhe][id] = model.modpar.totn[id]; accel_psy0[nre][id] = model.modpar.temp[id]; + if npc < nn { + accel_psy0[npc][id] = model.modpar.elec[id]; + } } accel2_config.iter = iter; { @@ -759,6 +1218,10 @@ pub fn run_tlusty( model.modpar.tk1[id] = 1.0 / t_new; } model.modpar.totn[id] = accel_psy0[nhe][id]; + // Update ELEC from accelerated PSY0 (INPC experiment) + if npc < nn && accel_psy0[npc][id] > 1e-30 { + model.modpar.elec[id] = accel_psy0[npc][id]; + } for ij in 0..nfreqe { if accel_psy0[ij][id] > 0.0 { model.expraf.radex[ij][id] = accel_psy0[ij][id]; @@ -825,6 +1288,9 @@ pub fn run_tlusty( } } + // Restore ijfr_explicit (RESOLV may have overwritten it) + model.frqall.ijfr_explicit = ijfr_solve.clone(); + let solves_output = solves_lte( &mut model, nn, @@ -845,6 +1311,21 @@ pub fn run_tlusty( total_iterations = iter as usize; + // After SOLVES updates TOTN and T, recompute DENS from the + // authoritative level populations (see comment at iter-1 block). + let nlvl = config.basnum.nlevel as usize; + for id in 0..(nd as usize) { + let sumpop: f64 = (0..nlvl).map(|j| model.levpop.popul[j][id]).sum(); + let dens_new = if sumpop > 0.0 { + sumpop * wmm + } else { + (model.modpar.totn[id] - model.modpar.elec[id]) * wmm + }; + if dens_new > 0.0 { + model.modpar.dens[id] = dens_new; + } + } + eprintln!(" NITER {}/{}: chmx={:.3e}, chmt={:.3e}, lfin={}", iter, niter, solves_output.chmx, solves_output.chmt, solves_output.lfin); @@ -964,6 +1445,87 @@ pub fn run_tlusty( itlucy, }; + // Diagnostic: compare final model vs gold (model_data = fort.8 input) — env-gated. + // Tells which of TOTN/ELEC/DENS/TEMP drifts vs the gold fixed point. + // Placed before the FINAL_RESOLV mutable borrow of `model`. + if std::env::var("TLUSTY_DENS_DIAG").is_ok() { + eprintln!("=== DENS_DIAG: final model vs gold (model_data) ==="); + eprintln!(" id relT% relTOTN% relNE% relRHO% | totn_f totn_g elec_f elec_g dens_f dens_g"); + let rel = |a: f64, b: f64| if b.abs() > 1e-300 { (a - b).abs() / b.abs() * 100.0 } else { 0.0 }; + for &id in &[0usize, 20, 30, 40, 49, 55, 60, 65, 69] { + let tf = model.modpar.temp[id]; let tg = model_data.temp[id]; + let ef = model.modpar.elec[id]; let eg = model_data.elec[id]; + let df = model.modpar.dens[id]; let dg = model_data.dens[id]; + let tof = model.modpar.totn[id]; + let tog = if let Some(ref tt) = model_data.totn { tt[id] } else { dg / wmm + eg }; + eprintln!(" {:2} {:6.3} {:8.3} {:7.3} {:8.3} | {:.3e} {:.3e} {:.3e} {:.3e} {:.3e} {:.3e}", + id, rel(tf,tg), rel(tof,tog), rel(ef,eg), rel(df,dg), + tof, tog, ef, eg, df, dg); + } + } + + // HYDROSTAT diagnostic: decompose the hydrostatic-equilibrium budget + // GRAV·ΔDM = BOLK·Δ(T·TOTN) + PCK·GRD + turb + // and compare the radiation-pressure gradient IMPLIED by gold's converged + // state vs the gradient implied by Rust's SOLVES solution. Because DM and T + // are gold's (DM is the mass grid, SOLVES doesn't change it; T matches gold + // to ~0%), the difference isolates the radiation-pressure deficit that the + // too-high TOTN is compensating for: + // PCK·(GRD_gold − GRD_rust) = BOLK·[T·δTOTN(id) − Tm·δTOTN(id−1)] + // (turb ≈ 0 because vturb is zero in the SOLVES path). Gated so default + // runs are unaffected; only SOLVES calls bhe_lte → readmodel/grey-start + // never reach here with this on (they take NITER=0 / LTGREY paths). + if std::env::var("TLUSTY_HYDRO_DIAG").is_ok() { + use crate::tlusty::state::constants::{BOLK, PCK}; + eprintln!("=== HYDROSTAT: radiation-pressure budget (grav={:.3e}, turb=0) ===", grav); + eprintln!(" id ΔDM(gold=rust?) PCK·GRD_rust_impl PCK·GRD_gold_impl DEFICIT δTOTN% T(id)"); + for &id in &[1usize, 10, 20, 30, 40, 45, 49, 55, 60, 65, 69] { + if id == 0 || id >= nd as usize { continue; } + let dm_g0 = model_data.dm[id]; + let dm_g1 = model_data.dm[id - 1]; + let dm_r0 = model.modpar.dm[id]; + let dm_r1 = model.modpar.dm[id - 1]; + let ddm_same = (dm_g0 - dm_r0).abs() < 1e-9 * dm_g0.abs().max(1e-30) + && (dm_g1 - dm_r1).abs() < 1e-9 * dm_g1.abs().max(1e-30); + let load = grav * (dm_g0 - dm_g1); // DM same for both + + let t_r0 = model.modpar.temp[id]; + let t_r1 = model.modpar.temp[id - 1]; + let to_r0 = model.modpar.totn[id]; + let to_r1 = model.modpar.totn[id - 1]; + let gas_rust = BOLK * (t_r0 * to_r0 - t_r1 * to_r1); + let prad_rust_impl = load - gas_rust; // turb=0 + + let t_g0 = model_data.temp[id]; + let t_g1 = model_data.temp[id - 1]; + // Gold TOTN: stored if present, else recomputed as DENS_GRID does + // (TOTN = DENS/WMM + ELEC — nuclei = mass/μ; + free electrons). + let gold_totn = |idd: usize| -> f64 { + if let Some(ref v) = model_data.totn { + v[idd] + } else { + model_data.dens[idd] / wmm + model_data.elec[idd] + } + }; + let to_g0 = gold_totn(id); + let to_g1 = gold_totn(id - 1); + let gas_gold = BOLK * (t_g0 * to_g0 - t_g1 * to_g1); + let prad_gold_impl = load - gas_gold; // turb=0 + + let deficit = prad_gold_impl - prad_rust_impl; // >0 ⇒ gold needs MORE rad press + let dtotn_pct = if to_g0.abs() > 1e-300 { + (to_r0 - to_g0).abs() / to_g0.abs() * 100.0 + } else { 0.0 }; + eprintln!( + " {:2} {:>5} {:15.4e} {:15.4e} {:10.3e} {:6.3} {:.1}", + id, if ddm_same { "same" } else { "DIFF" }, + prad_rust_impl, prad_gold_impl, deficit, dtotn_pct, t_r0 + ); + } + eprintln!(" (DEFICIT>0 ⇒ gold's atmosphere requires higher PCK·GRD than Rust's J/K \ + produces; this is the radiation-pressure shortfall TOTN compensates for.)"); + } + // Rewind fort.7 before final output (like Fortran REWIND 7) let _ = fort7_writer.rewind(); @@ -1000,8 +1562,12 @@ pub fn run_tlusty( /// 从 Eddington 灰大气近似创建初始模型。 /// /// 当 LTGREY=true 时使用，不需要 fort.8 文件。 -/// 使用 T(τ) = Teff * (3/4 * (τ + 2/3))^0.25 计算温度结构， -/// 并通过理想气体定律估算密度和电子密度。 +/// 严格对应 Fortran LTEGR 的简化版本： +/// 1. 在 Rosseland 光学深度标尺上对数等距布点 (TAUFIR → TAULAS) +/// 2. 使用 Hopf 函数计算 T(τ) +/// 3. 使用温度相关不透明度模型近似 Rosseland 平均不透明度 +/// 4. 预测-校正法积分流体静力学平衡方程 +/// 5. 从 P,T 计算 n, ρ, ne /// /// # 参数 /// - `teff`: 有效温度 (K) @@ -1009,75 +1575,372 @@ pub fn run_tlusty( /// /// # 返回值 /// InputModelData 包含初始灰大气模型 -fn create_grey_atmosphere(teff: f64, _grav: f64) -> InputModelData { +fn create_grey_atmosphere(teff: f64, grav: f64) -> InputModelData { use crate::tlusty::state::constants::{HMASS, BOLK}; + use crate::tlusty::math::temperature::compute_hopf; - // 深度点数（与 Fortran TLUSTY 默认值一致） + // 参数：匹配 Fortran LTEGR 默认值 let nd = 70; + let taufir = 1e-7_f64; + let taulas = 316.0_f64; + let abros0 = 0.34_f64; // κ_es = 0.2*(1+X) for X=0.7 - // Rosseland 光学深度范围：log(τ) 从 -4 到 +2 - let log_tau_min = -4.0_f64; - let log_tau_max = 2.0_f64; + // 辐射压力 + let t4 = teff.powi(4); + let dprad = 1.891204931e-15 * t4; + let prad0 = dprad / 1.732; + let wmm = 1.3 * HMASS; + // Rosseland 光学深度网格（对数等距） + let dml0 = taufir.ln(); + let dlgm = (taulas.ln() - dml0) / (nd - 1) as f64; + let mut tau = vec![0.0_f64; nd]; + for i in 0..nd { + tau[i] = (dml0 + i as f64 * dlgm).exp(); + } + + // H/He 丰度 + let xh = 0.70; + let xhe = 0.28; + + // 预计算 LTE 频率积分网格 + let lte_grid = generate_lte_frequency_grid(teff, 200); + + /// H+He Saha 方程求解器 + fn solve_hhe_saha(t: f64, an: f64, xh: f64, xhe: f64) -> (f64, f64, f64, f64) { + let ytot = xh + xhe / 4.0; + let f_h = xh / ytot; + let f_he = (xhe / 4.0) / ytot; + let saha_factor = |chi_ev: f64, t: f64| -> f64 { + let kt = 8.617333e-5 * t; + if kt <= 0.0 { return 0.0; } + let theta = chi_ev / kt; + if theta > 80.0 { return 0.0; } + 2.4148e15 * t.powf(1.5) * (-theta).exp() + }; + let ne_max_frac = (f_h + 2.0 * f_he) / (1.0 + f_h + 2.0 * f_he); + let mut ne = an * ne_max_frac * 0.95; + let mut np_out = 0.0; + let mut nh_neutral_out = 0.0; + for _ in 0..50 { + let n_nucleons = (an - ne).max(an * 1e-10); + let n_h = n_nucleons * f_h; + let n_he = n_nucleons * f_he; + let s_h = saha_factor(13.598, t); + let x_h = if ne > 0.0 { s_h / (s_h + ne) } else { 1.0 }; + let s_he1 = saha_factor(24.587, t); + let x_he1 = if ne > 0.0 { s_he1 / (s_he1 + ne) } else { 0.0 }; + let s_he2 = saha_factor(54.418, t); + let x_he2 = if ne > 0.0 { s_he2 / (s_he2 + ne) } else { 0.0 }; + let n_hplus = n_h * x_h; + let n_heplus = n_he * x_he1 * (1.0 - x_he2); + let n_heplusplus = n_he * x_he1 * x_he2; + let ne_new = n_hplus + n_heplus + 2.0 * n_heplusplus; + let ne_clamped = ne_new.min(an * 0.99).max(0.0); + let ne_relaxed = 0.5 * ne + 0.5 * ne_clamped; + let rel = (ne_relaxed - ne).abs() / ne.max(1e-30); + ne = ne_relaxed; + np_out = n_hplus; + nh_neutral_out = n_h * (1.0 - x_h); + if rel < 1e-6 { break; } + } + let n_nucleons = (an - ne).max(an * 1e-10); + let n_h = n_nucleons * f_h; + (ne, n_h, np_out, nh_neutral_out) + } + + /// LTE Rosseland 平均不透明度 (cm²/g) + fn lte_rosseland_opacity(t: f64, rho: f64, xh: f64, xhe: f64, + grid: &crate::tlusty::math::continuum::LteFrequencyGrid) -> (f64, f64, f64) { + if rho <= 1e-30 { return (0.34, 0.0, 0.5); } + let wmm_local = 1.3 * 1.67e-24; + let an_est = rho / wmm_local * 2.0; + let (ne, nh_total, np, nh_neutral) = solve_hhe_saha(t, an_est, xh, xhe); + let params = LteOpacityParams { + t, ne, nh_total, np, nh_neutral, nhm: 0.0, rho, + uh: 2.0, uhe: 1.0, uhep: 2.0, xh, xhe, + }; + let result = lte_meanopt(¶ms, grid); + (result.opros.max(0.01), ne, ne / an_est) + } + + // 预测-校正法积分流体静力学平衡 let mut dm = vec![0.0; nd]; let mut temp = vec![0.0; nd]; let mut elec = vec![0.0; nd]; let mut dens = vec![0.0; nd]; let mut totn_arr = vec![0.0; nd]; + let mut plog1 = 0.0_f64; + let mut plog2 = 0.0_f64; + let mut plog3 = 0.0_f64; + let mut plog4 = 0.0_f64; + let mut dplog1 = 0.0_f64; + let mut dplog2 = 0.0_f64; + let mut dplog3 = 0.0_f64; + let mut abros = abros0; - // 估计 Rosseland 平均不透明度 (cm²/g) - 纯氢大气的典型值 - let kappa_ross = 0.4; - - // 平均分子量（H-He 混合，约 0.6 for solar composition） - let wmm = 0.6 * HMASS; - - for id in 0..nd { - // 等间距对数 Rosseland 光学深度 - let frac = id as f64 / (nd - 1) as f64; - let log_tau = log_tau_min + frac * (log_tau_max - log_tau_min); - let tau = 10.0_f64.powf(log_tau); - - // Eddington 灰大气温度: T(τ) = Teff * (3/4 * (τ + 2/3))^0.25 - let t4 = 0.75 * (tau + 2.0 / 3.0); - temp[id] = teff * t4.powf(0.25); - - // 柱质量深度: DM = τ / κ (g/cm²) - dm[id] = tau / kappa_ross; - - // 气体压力: P = g * DM (dyne/cm²) - let p_gas = _grav * dm[id]; - - // 总粒子数密度: n = P / (k * T) (cm⁻³) - totn_arr[id] = p_gas / (BOLK * temp[id]); - - // 质量密度: ρ = n * WMM (g/cm³) - dens[id] = totn_arr[id] * wmm; - - // 电子密度: 粗略估计，假设 50% 电离（H-He 混合） - // Saha 方程给出更精确的值，但这里只是初始猜测 - elec[id] = 0.5 * totn_arr[id]; + for i in 0..nd { + let taur = tau[i]; + let mut plog = if i == 0 { + (grav / abros * taur + prad0).ln() + } else if i <= 3 { + plog1 + dplog1 + } else { + (3.0 * plog4 + 8.0 * dplog1 - 4.0 * dplog2 + 8.0 * dplog3) / 3.0 + }; + let mut j = 0; + let mut dplog = 0.0; + loop { + let pnew = if i == 0 { + (grav / abros * taur + prad0).ln() + } else if i <= 3 { + (plog + 2.0 * plog1 + dplog + dplog1) / 3.0 + } else { + (126.0 * plog1 - 14.0 * plog3 + 9.0 * plog4 + + 42.0 * dplog + 108.0 * dplog1 - 54.0 * dplog2 + 24.0 * dplog3) / 121.0 + }; + let error = (pnew - plog).abs(); + plog = pnew; + let ptot = plog.exp(); + let p_gas = ptot - taur * dprad - prad0; + j += 1; + let hopf = compute_hopf(taur, 0.0); + let t = (0.75 * t4 * (taur + hopf) + 0.0).powf(0.25); + let an = if t > 0.0 { p_gas / (BOLK * t) } else { 0.0 }; + let (ane_saha, _, _, _) = solve_hhe_saha(t, an, xh, xhe); + let rho = wmm * (an - ane_saha); + abros = lte_rosseland_opacity(t, rho, xh, xhe, <e_grid).0; + dplog = grav / abros * taur / ptot * dlgm; + if error <= 1e-4 || j >= 10 { + temp[i] = t; + dm[i] = (ptot - prad0) / grav; + let an = p_gas / (BOLK * t); + totn_arr[i] = an; + let (ane_final, _, _, _) = solve_hhe_saha(t, an, xh, xhe); + elec[i] = ane_final; + dens[i] = wmm * (an - elec[i]); + plog4 = plog3; plog3 = plog2; plog2 = plog1; plog1 = plog; + dplog3 = dplog2; dplog2 = dplog1; dplog1 = dplog; + break; + } + } } - // NUMPAR = 3 for LTE model (T, Ne, Rho) + 1 for TOTN (negative = has TOTN) - // Matches Fortran convention: NUMPAR < 0 means molecule data present - let numpar = -4_i32; // T + Ne + Rho + TOTN - + let numpar = -4_i32; eprintln!(" Grey atmosphere: nd={}, teff={:.0}", nd, teff); + eprintln!(" TAUFIR={:.0e}, TAULAS={:.0e}", taufir, taulas); eprintln!(" T range: {:.0} - {:.0} K", temp[0], temp[nd - 1]); + eprintln!(" DM range: {:.3e} - {:.3e}", dm[0], dm[nd - 1]); InputModelData { - ndpth: nd, - numpar, - dm, - temp, - elec, - dens, - totn: Some(totn_arr), + ndpth: nd, numpar, dm, temp, elec, dens, + totn: Some(totn_arr), zd: None, popul: None, + } +} + +/// 灰大气生成（完整 LTEGR 路径）。 +/// +/// 调用 io::ltegr::ltegr()，它忠实移植 Fortran LTEGR： +/// Part 1: 流体静力学平衡积分 (TAU 网格 + ROSSOP → DM=TAU/κ_R) +/// Part 3: 插值到对数等距 TAU0 网格 (TAUFIR..TAULAS-1 + TAULAS) +/// +/// 关键: ROSSOP 内部用 Opacf0State.compute_opacity() 计算频率相关不透明度 +/// (与 readmodel/RESOLV 验证 0.0000% 的同一机制), 再做 Rosseland 平均。 +/// 这取代了简化的 create_grey_atmosphere (其 lte_meanopt 有补偿误差使 κ_R 偏高 ~2-3x)。 +/// +/// 触发: 环境变量 TLUSTY_GREY_LTEGR (默认灰大气生成器)。 +fn create_grey_atmosphere_via_ltegr(teff: f64, grav: f64) -> InputModelData { + use crate::tlusty::io::ltegr::{ltegr as run_ltegr, LtegrConfig, LtegrParams}; + + let nd = 70usize; + // 平均分子量 WMM = HMASS * WMY/YTOT (Fortran INITIA tlusty208.f:2831), + // 与 compute_abundance_params 一致 (含金属 fallback: H,He,C,N,O)。 + // gold 参考值 YTOT=1.08588, WMY=1.35982 → WMM=2.09547e-24 g = 1.2522·HMASS。 + // 之前硬编码 1.3·HMASS 比迭代用的 faithful WMM 高 3.7%, 使 LTEGR 构造的 + // dens=(totn-elec)·wmm 系统性偏高 → 深层 ne/rho 偏低 ~3%。 + let hmass = crate::tlusty::state::constants::HMASS; + const ABUND_AMASS: [(f64, f64); 5] = [ + (1.0, 1.008), (0.0851, 4.003), (2.45e-4, 12.011), + (6.03e-5, 14.007), (4.57e-4, 15.999), + ]; + let (mut ytot, mut wmy) = (0.0_f64, 0.0_f64); + for &(a, m) in &ABUND_AMASS { ytot += a; wmy += a * m; } + let wmm_val = hmass * wmy / ytot; + let wmm = vec![wmm_val; nd]; + + // 配置匹配 gold ref (Fortran LTEGR 默认): TAUFIR=1e-7, TAULAS=316 + let config = LtegrConfig { + teff, + grav, + taufir: 1e-7, + taulas: 316.0, + abros0: 0.34, // κ_es = 0.2*(1+X), 仅初始估计 + tsurf: 0.0, // 精确 Hopf + albave: 0.0, // 无风包层 + dion0: 0.5, // 初始电离度 + ndgrey: 0, // NDEPTH = ND + idgrey: 0, // 对数等距 TAU 标尺 + nconit: 0, // 无对流迭代 (HMIX0=0) + ipring: 0, + hmix0: 0.0, // 不考虑对流 + ifprad: 1, // 考虑辐射压力 + lte: true, + lchc0: 0, + irspl0: 0, + }; + + // 未在 ltegr() 内部使用的字段 (API 完整性), 传占位值 + let nlevel = 39usize; + let nfirst = vec![0i32; nlevel]; + let nka = vec![0i32; nlevel]; + let anato = vec![vec![0.0f64; 100]; nd]; + let anmol = vec![vec![0.0f64; 600]; nd]; + + let params = LtegrParams { + config, + nd, + wmm: &wmm, + nlevel, + nfirst: &nfirst, + nka: &nka, + ielh: 1, + iathe: 2, + anato: &anato, + anmol: &anmol, + }; + + let out = run_ltegr(¶ms, None::<&mut FortranWriter>); + let nd_out = out.nd; + + // ltegr() 返回 MDEPTH 大小数组 (填充到 nd), 截断到实际深度数 + let take = |v: Vec| v.into_iter().take(nd_out).collect::>(); + + eprintln!(" Grey atmosphere (ltegr full path): nd={}", nd_out); + eprintln!(" T range: {:.0} - {:.0} K", out.temp[0], out.temp[nd_out - 1]); + eprintln!(" DM range: {:.3e} - {:.3e}", out.dm[0], out.dm[nd_out - 1]); + + InputModelData { + ndpth: nd_out, + numpar: -4_i32, + dm: take(out.dm), + temp: take(out.temp), + elec: take(out.elec), + dens: take(out.dens), + totn: Some(take(out.totn)), zd: None, popul: None, } } +/// 实验函数 [TLUSTY_GREY_REFGRID]: 在 fort.8 的 gold-ref DM 网格上构建灰大气 +/// (Eddington T + LTE Ne/ρ)。 +/// +/// readmodel 路径已证明 (gold 网格 + gold 值) 是 SOLVES 不动点。本函数测试 +/// (gold 网格 + 灰 Eddington 值) 是否在 SOLVES 收敛域内: +/// - 若收敛到 gold ref → 网格映射是唯一问题, κ_R/grid 是全部根因 +/// - 若不收敛 → SOLVES 从灰大气收敛性是独立问题 +/// +/// 因为 gold 模型有效 κ_R≈1 (gold DM[nd-1]≈298 ≈ taulas=316, DM 跨度≈tau 跨度), +/// 取 τ(id) ≈ DM(id)。 +fn create_grey_atmosphere_on_ref_grid(teff: f64, grav: f64) -> Result { + use crate::tlusty::state::constants::{HMASS, BOLK}; + use crate::tlusty::math::temperature::compute_hopf; + + // 1) 读取 fort.8 的 gold-ref DM 网格 + let fort8_path = "fort.8"; + if !std::path::Path::new(fort8_path).exists() { + return Err("fort.8 not found".to_string()); + } + let fort8_file = File::open(fort8_path).map_err(|e| format!("open fort.8: {}", e))?; + let mut fort8_reader = FortranReader::new(io::BufReader::new(fort8_file)); + let inpmod_params = InpmodParams { + intrpl: 0, idisk: 0, ifmol: 0, tmolim: 1e10, teff, + }; + let ref_model = read_tlusty_model(&mut fort8_reader, &inpmod_params) + .map_err(|e| format!("read fort.8: {:?}", e))?; + + let nd = ref_model.ndpth; + if nd < 2 { + return Err(format!("fort.8 nd={} too small", nd)); + } + let dm_ref = ref_model.dm.clone(); + + // 2) 在 gold DM 网格上构建灰大气 (Eddington T + LTE Ne/ρ) + let t4 = teff.powi(4); + let dprad = 1.891204931e-15 * t4; // 辐射压力尺度 + let wmm = 1.3 * HMASS; + let xh = 0.70; + let xhe = 0.28; + + // H+He Saha 求解器 (与 create_grey_atmosphere 相同), 返回 (ne, n_gas) + fn saha_ne(t: f64, an: f64, xh: f64, xhe: f64) -> (f64, f64) { + let ytot = xh + xhe / 4.0; + let f_h = xh / ytot; + let f_he = (xhe / 4.0) / ytot; + let sf = |chi_ev: f64, t: f64| -> f64 { + let kt = 8.617333e-5 * t; + if kt <= 0.0 { return 0.0; } + let theta = chi_ev / kt; + if theta > 80.0 { return 0.0; } + 2.4148e15 * t.powf(1.5) * (-theta).exp() + }; + let ne_max_frac = (f_h + 2.0 * f_he) / (1.0 + f_h + 2.0 * f_he); + let mut ne = an * ne_max_frac * 0.95; + for _ in 0..50 { + let n_nucleons = (an - ne).max(an * 1e-10); + let n_h = n_nucleons * f_h; + let s_h = sf(13.598, t); + let x_h = if ne > 0.0 { s_h / (s_h + ne) } else { 1.0 }; + let s_he1 = sf(24.587, t); + let x_he1 = if ne > 0.0 { s_he1 / (s_he1 + ne) } else { 0.0 }; + let s_he2 = sf(54.418, t); + let x_he2 = if ne > 0.0 { s_he2 / (s_he2 + ne) } else { 0.0 }; + let n_hplus = n_h * x_h; + let n_he = n_nucleons * f_he; + let n_heplus = n_he * x_he1 * (1.0 - x_he2); + let n_heplusplus = n_he * x_he1 * x_he2; + let ne_new = n_hplus + n_heplus + 2.0 * n_heplusplus; + let ne_clamped = ne_new.min(an * 0.99).max(0.0); + let ne_relaxed = 0.5 * ne + 0.5 * ne_clamped; + let rel = (ne_relaxed - ne).abs() / ne.max(1e-30); + ne = ne_relaxed; + if rel < 1e-6 { break; } + } + let _ = (f_h, f_he); + (ne, an - ne) + } + + let mut temp = vec![0.0_f64; nd]; + let mut elec = vec![0.0_f64; nd]; + let mut dens = vec![0.0_f64; nd]; + let mut totn_arr = vec![0.0_f64; nd]; + + for id in 0..nd { + let dm = dm_ref[id]; + // τ ≈ DM (gold 有效 κ_R ≈ 1) + let tau = dm; + let hopf = compute_hopf(tau, 0.0); + let t = (0.75 * t4 * (tau + hopf)).powf(0.25); + // P_gas = grav×DM − τ×dprad (来自 DM=(P_total−prad0)/grav, P_total=grav×DM+prad0) + let p_gas = (grav * dm - tau * dprad).max(1.0); + let an = if t > 0.0 { p_gas / (BOLK * t) } else { 0.0 }; + let (ne, _n_gas) = saha_ne(t, an, xh, xhe); + let rho = wmm * (an - ne); + temp[id] = t; + elec[id] = ne; + dens[id] = rho; + totn_arr[id] = an; + } + + eprintln!(" [REFGRID] nd={}, DM range {:.3e}..{:.3e}, T range {:.0}..{:.0}", + nd, dm_ref[0], dm_ref[nd - 1], temp[0], temp[nd - 1]); + + Ok(InputModelData { + ndpth: nd, numpar: ref_model.numpar, dm: dm_ref, temp, elec, dens, + totn: Some(totn_arr), zd: None, popul: None, + }) +} + // ============================================================================ // 测试 // ============================================================================ diff --git a/src/tlusty/math/atomic/gfree.rs b/src/tlusty/math/atomic/gfree.rs index a91eef5..2b34065 100644 --- a/src/tlusty/math/atomic/gfree.rs +++ b/src/tlusty/math/atomic/gfree.rs @@ -74,10 +74,14 @@ pub fn gfree0(id: usize, temp: &[f64], gffpar: &mut crate::tlusty::state::GffPar if thet0_over_t >= THMIN { // 正常情况：计算导数 + // 注意: Fortran 在此分支内重新赋值 THET=THT=THET0/T (tlusty208.f:10900)， + // 用于 GF*D 导数（与 GF0-GF6 值多项式所用的 THET=T/THET0 互为倒数）。 + // 此处 thet_d = THET0/T 忠实 Fortran 的重新赋值。 + let thet_d = thet0_over_t; gffpar.gf0d[id] = B0 * thet1; - gffpar.gf1d[id] = (A1 + B1 * thet * 2.0) * thet1; - gffpar.gf2d[id] = (A2 + B2 * thet * 2.0) * thet1; - gffpar.gf3d[id] = (A3 + B3 * thet * 2.0) * thet1; + gffpar.gf1d[id] = (A1 + B1 * thet_d * 2.0) * thet1; + gffpar.gf2d[id] = (A2 + B2 * thet_d * 2.0) * thet1; + gffpar.gf3d[id] = (A3 + B3 * thet_d * 2.0) * thet1; gffpar.gf4d[id] = B4 * thet1; gffpar.gf5d[id] = D0 * thet1; gffpar.gf6d[id] = gffpar.gf0d[id] + gffpar.gf5d[id] / XMIN; diff --git a/src/tlusty/math/continuum/lte_opacity.rs b/src/tlusty/math/continuum/lte_opacity.rs index 3cda205..10c0d8d 100644 --- a/src/tlusty/math/continuum/lte_opacity.rs +++ b/src/tlusty/math/continuum/lte_opacity.rs @@ -130,7 +130,7 @@ pub struct LteFrequencyGrid { // ============================================================================ /// 生成用于 LTE 不透明度积分的频率网格。 -pub fn generate_lte_frequency_grid(teff: f64, nfreq: usize) -> LteFrequencyGrid { +pub fn generate_lte_frequency_grid(_teff: f64, nfreq: usize) -> LteFrequencyGrid { let frmin: f64 = 1e13; let frmax: f64 = 3e16; @@ -156,10 +156,11 @@ pub fn generate_lte_frequency_grid(teff: f64, nfreq: usize) -> LteFrequencyGrid }; weights.push(w); - let x = HK * fr / teff; - let ex = if x < 150.0 { x.exp() } else { 1e150 }; - let bn = c1 * fr.powi(3) / (ex - 1.0); - bnue.push(bn); + // BNUE = (2H/c²)·ν³ — prefactor only (matches Fortran BNUE=BN*ν³, tlusty208.f:30748). + // The (exp(hν/kT)-1) divisor is applied by consumers at the LOCAL T + // (plan = bnue * e1 * w, e1=1/(exp(hν/kT_local)-1)). Baking the TEFF divisor in + // here left a frequency-dependent spurious factor that biased the Rosseland mean. + bnue.push(c1 * fr.powi(3)); } LteFrequencyGrid { freq, weights, bnue, ijxco: vec![], ijfr: vec![] } @@ -403,14 +404,20 @@ pub fn generate_inifrc_frequency_grid(teff: f64, nfreq_base: usize) -> LteFreque let nfreq = all_freqs.len(); - // Compute Planck function + // Compute Planck function prefactor BNUE = (2H/c²)·ν³. + // + // IMPORTANT: matches Fortran BNUE (tlusty208.f:30748 `BNUE(IJ)=BN*FR15³`, BN=2H/c²): + // BNUE is the PREFACTOR ONLY — the (exp(hν/kT)-1) divisor is applied separately + // at the LOCAL temperature by every consumer (MEANOP, meanopt, lte_meanopt, + // ltegr call_meanopt all compute `plan = bnue * e1 * w` with e1=1/(exp(hν/kT_local)-1)). + // + // Earlier code divided by (exp(hν/kT_eff)-1) here, leaving a frequency-dependent + // spurious factor that biased the Rosseland mean weight toward low-ν and inflated + // κ_R at depth (grey-start DM grid too small). The divisor must NOT be baked in. let c1 = 2.0 * H / (CLIGHT * CLIGHT); let mut bnue = Vec::with_capacity(nfreq); for &fr in &all_freqs { - let x = HK * fr / teff; - let ex = if x < 150.0 { x.exp() } else { 1e150 }; - let bn = c1 * fr.powi(3) / (ex - 1.0); - bnue.push(bn); + bnue.push(c1 * fr.powi(3)); } // Compute IJFR: select explicit (non-ALI) frequency points @@ -501,8 +508,13 @@ pub fn lte_meanopt(params: &LteOpacityParams, grid: &LteFrequencyGrid) -> LteOpa let ex = x_clamped.exp(); let e1 = 1.0 / (ex - 1.0); + // ∂B_ν/∂T·w = plan·u·ex·e1 (u=hkt·fr); the constant T it implicitly omits + // cancels in the Rosseland ratio. Single e1 — matches faithful meanopt.rs + // and Fortran tlusty208.f:22433 (DPLAN=PLAN*X/T/(UN-UN/EX)). A previous + // `*e1*e1` here was the "κ_R 偏高 ~2-3x" compensation error noted at + // main.rs:1573 — the extra frequency-dependent e1 re-weights the mean. let plan = bnue * e1 * w; - let dplan = plan * hkt * fr * ex * e1 * e1; + let dplan = plan * hkt * fr * ex * e1; let (ab, sct) = compute_opacity_at_frequency(fr, t, ne, nh, np, nhm, hkt, sgff, params); @@ -559,6 +571,16 @@ pub fn compute_opacity_at_frequency( } // 2. 氢自由-自由 + // 注意: 受激发射因子当前为 1/(1-e^-hν/kT). 诊断显示此式在高温深层使 Rosseland + // 平均 κ_R 发散 (gold-ref 条件下 d69 κ_R≈140, 物理上应 ~1). 但简单改为乘以 + // (1-e^-hν/kT) 虽修正孤立 κ 值, 却因 κ-ρ-Ne-Saha 耦合使灰大气深度网格更差 + // (DM[69] 163→126, gold=298), 且 Fortran H- 自由-自由亦用除法 (tlusty208.f:11147). + // 忠实修复需移植 Fortran SFF2/SFF3 表格化系数, 非简单因子翻转. 暂保留原式. + // [2026-06-18 复核] 同时翻转 sf2 与改 bnue 为本地 T 仍非稳健: 灰大气 DM[69] + // 163→216 (改善) 但 gold-ref 深层诊断 κ_R 140→500 (恶化), 单元测试 es_frac 失败. + // ⇒ 简化 lte_meanopt 有相互补偿的误差, 局部物理修正会交换误差而非收敛. + // 真正修复 = 移植 TLUSTY 完整 ROSSOP/COMOP 不透明度机制 (含 SFF2=EXP(FF·HK/T) + // 约定, exact Gaunt, WOP 占据概率, H- SFFHMI 表). 暂保留原式. if np > 0.0 && ne > 0.0 { let frinv = 1.0 / fr; let fr3inv = frinv * frinv * frinv; @@ -703,6 +725,41 @@ pub fn quick_lte_rosseland(params: &LteOpacityParams) -> f64 { mod tests { use super::*; + /// Diagnostic: compute LTE Rosseland opacity breakdown for gold-ref (T,Ne,ρ) + /// profile points, compare to effective κ = tau/DM from the converged grid. + /// Purpose: find which opacity source is overestimated at depth (grey-start + /// grid compressed because κ_rust ≈ 3× too high at depth vs gold κ≈1.06). + #[test] + fn diag_goldref_opacity_breakdown() { + // (label, T, Ne, rho, eff_kappa=tau/DM) from tests/tlusty/hhe_fortran/fort.7.ref + let pts: &[(&str, f64, f64, f64, f64)] = &[ + ("d0 surf", 26306.2, 3.764e8, 7.304e-16, 0.343), + ("d20 ", 27030.2, 2.258e11, 4.455e-13, 0.343), + ("d34 ", 27800.0, 5.0e12, 1.0e-11, 0.396), + ("d49 ", 40000.0, 1.0e14, 2.0e-10, 0.951), + ("d60 ", 90000.0, 2.0e15, 4.0e-9, 1.145), + ("d69 deep ", 137872.1, 5.679e16, 1.102e-7, 1.061), + ]; + let wmm = 1.3 * 1.67e-24; // mean molecular weight + let grid = generate_lte_frequency_grid(35000.0, 200); + eprintln!("\n{:>9} {:>10} {:>10} {:>10} {:>10} {:>10} | {:>8}", + "label","kappa_R","opes","opff","opbf","ophm","eff_k"); + for (label, t, ne, rho, eff_k) in pts { + // H fully ionized at these T: np ~ n_H, neutral H tiny + let n_heavy = rho / wmm; + let nh_total = 0.70 * n_heavy * 2.0; // X=0.70 fraction, nucleon count + let np = nh_total * 0.999; + let nh_neutral = nh_total * 0.001; + let params = LteOpacityParams { + t: *t, ne: *ne, nh_total, np, nh_neutral, nhm: 0.0, rho: *rho, + uh: 2.0, uhe: 1.0, uhep: 2.0, xh: 0.70, xhe: 0.28, + }; + let r = lte_meanopt(¶ms, &grid); + eprintln!("{:>9} {:10.3} {:10.3} {:10.3} {:10.3} {:10.3} | {:8.3}", + label, r.opros, r.opes, r.opff, r.opbf, r.ophm, eff_k); + } + } + #[test] fn test_lte_opacity_hot_star() { let params = LteOpacityParams { diff --git a/src/tlusty/math/continuum/opacf0_state.rs b/src/tlusty/math/continuum/opacf0_state.rs index 1d02708..5887cef 100644 --- a/src/tlusty/math/continuum/opacf0_state.rs +++ b/src/tlusty/math/continuum/opacf0_state.rs @@ -242,6 +242,40 @@ impl Opacf0State { popul: &[Vec], id: usize, gffpar: &GffPar, + ) -> (f64, f64, f64, f64) { + self.compute_opacity_inner(fr, t, ne, popul, id, gffpar, None) + } + + /// Full opacity with occupation probabilities (WOP) for BF emission. + /// + /// When `wop` is provided, the BF emission coefficient is corrected: + /// emis_bf = σ × POPUL(II) × WOP(II) × CORR + /// matching the Fortran OPACF0 formula: + /// EMTRA(ITR,ID) = POPUL(JJ,ID)*ANE*SBF(II)*WOP(II,ID)*CORR + /// In LTE: POPUL(II) = POPUL(JJ)*ANE*SBF(II), so EMTRA = POPUL(II)*WOP(II)*CORR. + /// For the H-He model, CORR=1 for all transitions. + pub fn compute_opacity_with_wop( + &self, + fr: f64, + t: f64, + ne: f64, + popul: &[Vec], + id: usize, + gffpar: &GffPar, + wop: &[Vec], + ) -> (f64, f64, f64, f64) { + self.compute_opacity_inner(fr, t, ne, popul, id, gffpar, Some(wop)) + } + + fn compute_opacity_inner( + &self, + fr: f64, + t: f64, + ne: f64, + popul: &[Vec], + id: usize, + gffpar: &GffPar, + wop: Option<&[Vec]>, ) -> (f64, f64, f64, f64) { let hkt = HK / t; let sqrt_t = t.sqrt(); @@ -302,12 +336,22 @@ impl Opacf0State { // ABTRA = POPUL(II, ID) — 下能级（束缚态）种群 abso += sigma * pop_low; - // EMTRA — LTE 中 Kirchhoff 定律: EMTRA = ABTRA = pop_low × σ - // 对应 Fortran OPACF0 行 70-71: - // ABTRA(ITR,ID) = POPUL(II,ID) + // EMTRA — Fortran OPACF0: // EMTRA(ITR,ID) = POPUL(JJ,ID)*ANE*SBF(II)*WOP(II,ID)*CORR - // LTE 中 (Saha-Boltzmann 平衡, WOP=1, CORR=1): EMTRA = POPUL(II) = pop_low - emis += sigma * pop_low; + // In LTE: POPUL(II) = POPUL(JJ)*ANE*SBF(II), so: + // EMTRA = POPUL(II) * WOP(II) * CORR + // For H-He model: CORR=1 (NKE=JJ for all transitions) + // WOP from WNSTOR accounts for pressure ionization (WOP < 1) + let wop_factor = if let Some(wop_arr) = wop { + if bt.ilow < wop_arr.len() && id < wop_arr[bt.ilow].len() { + wop_arr[bt.ilow][id] + } else { + 1.0 + } + } else { + 1.0 + }; + emis += sigma * pop_low * wop_factor; } // ================================================================ diff --git a/src/tlusty/math/io/output.rs b/src/tlusty/math/io/output.rs index 447ec89..ccfbb2e 100644 --- a/src/tlusty/math/io/output.rs +++ b/src/tlusty/math/io/output.rs @@ -269,6 +269,12 @@ fn write_depth_line_with_popul( } } + // Flush remaining data (when nlevel=0 or columns don't fill a full line) + if !line.is_empty() { + writer.write_raw(&line)?; + writer.write_newline()?; + } + Ok(()) } @@ -310,6 +316,12 @@ fn write_depth_line_with_popul_disk( } } + // Flush remaining data (when nlevel=0 or columns don't fill a full line) + if !line.is_empty() { + writer.write_raw(&line)?; + writer.write_newline()?; + } + Ok(()) } diff --git a/src/tlusty/math/opacity/inilam.rs b/src/tlusty/math/opacity/inilam.rs index f69be50..54dc7b1 100644 --- a/src/tlusty/math/opacity/inilam.rs +++ b/src/tlusty/math/opacity/inilam.rs @@ -402,10 +402,15 @@ pub fn inilam( for ion in 0..nion { let _nf = atomic.nfirst[ion] as usize; let _nl = atomic.nlast[ion] as usize; + // NNEXT(ION) is a 1-based Fortran level pointer + // (rdata.rs: NNEXT = NFIRST + NLEVS; chckse.rs uses nnext-1). + // POPUL(NNEXT(ION),ID) is valid in Fortran 1-based, but the + // 0-based Rust row index must be NN-1 (else NN==NLEVEL indexes + // one past the last level → out-of-bounds panic). let nn = atomic.nnext[ion] as usize; if nn > 0 && nn <= nlevel { - let pop_next = *get_2d(model.popul, nn, id, nlevel); + let pop_next = *get_2d(model.popul, nn - 1, id, nlevel); if pop_next > 0.0 && atomic.iltlev[i] == 0 { let pop_i = *get_2d(model.popul, i, id, nlevel); let bfac_val = pop_i / (pop_next * sbw); diff --git a/src/tlusty/math/radiative/feautrier.rs b/src/tlusty/math/radiative/feautrier.rs index d3ac686..9abef66 100644 --- a/src/tlusty/math/radiative/feautrier.rs +++ b/src/tlusty/math/radiative/feautrier.rs @@ -5,18 +5,25 @@ use crate::tlusty::state::constants::{BN, HALF, HK, UN}; -/// Gauss-Legendre NMU=4 quadrature for QQ0 -const AMU: [f64; 4] = [ - 0.06943184420297371, - 0.33000947820757187, - 0.669_990_521_792_428_1, - 0.930_568_155_797_026_3, +/// Gauss-Legendre NMU=3 quadrature for the LTE formal solution. +/// +/// Matches Fortran TLUSTY: RTEFR1 hardcodes NMU=3 (tlusty208.f:39153), and the +/// angle-setup routine uses `call gauleg(zero,un,amu0,wtmu0,nmu,mmu)` with +/// `PARAMETER (NMU3=3)` for the no-irradiation (WANGLE=0) LTE case +/// (tlusty208.f:40013-40018). Standard 3-point Gauss-Legendre on [0,1]: +/// nodes μ = (1±√(3/5))/2, 1/2; weights 5/18, 4/9, 5/18. +/// +/// Shared with rtefr1 (resolv formal solution) so both solvers use identical +/// angular quadrature — re-exported via `pub use feautrier::*`. +pub const AMU: [f64; 3] = [ + 0.11270166537925807, + 0.5, + 0.8872983346207419, ]; -const WTMU: [f64; 4] = [ - 0.17392742256872692, - 0.32607257743127314, - 0.32607257743127314, - 0.17392742256872692, +pub const WTMU: [f64; 3] = [ + 0.27777777777777778, + 0.44444444444444444, + 0.27777777777777778, ]; /// Result of the Feautrier formal solution. @@ -76,7 +83,7 @@ pub fn feautrier_solve( // Compute QQ0 for surface BC let mut qq0 = 0.0; - for k in 0..4 { + for k in 0..AMU.len() { let tamm = taumin / AMU[k]; let p0 = 1.0 - (-tamm).exp(); qq0 += p0 * AMU[k] * WTMU[k]; diff --git a/src/tlusty/math/radiative/rtefr1.rs b/src/tlusty/math/radiative/rtefr1.rs index 4d13f1d..7dfc7ec 100644 --- a/src/tlusty/math/radiative/rtefr1.rs +++ b/src/tlusty/math/radiative/rtefr1.rs @@ -472,7 +472,7 @@ pub fn rtefr1(params: &Rtefr1Params, model: &mut Rtefr1ModelState) { for i in 0..nmu { anu[i * nd + id] = 0.0; for j in 0..nmu { - anu[i * id + id] += bb_local[i * mmu + j] * vl[j]; + anu[i * nd + id] += bb_local[i * mmu + j] * vl[j]; } } } @@ -542,7 +542,7 @@ pub fn rtefr1(params: &Rtefr1Params, model: &mut Rtefr1ModelState) { let b_val = HALF / a_val; aa[i * mmu + i] = a_val; vl[i] = b_val * st0[id] + pland_var + params.amu[i] * dplan_val - + aa[i * mmu + i] * anu[i * id + (id - 1)]; + + aa[i * mmu + i] * anu[i * nd + (id - 1)]; for j in 0..nmu { bb_local[i * mmu + j] = b_val * ss0[id] * params.wtmu[j] - aa[i * mmu + i] * d[i * mmu * nd + j * nd + (id - 1)]; diff --git a/src/tlusty/math/solvers/matgen_lte.rs b/src/tlusty/math/solvers/matgen_lte.rs index b435b48..b49a42c 100644 --- a/src/tlusty/math/solvers/matgen_lte.rs +++ b/src/tlusty/math/solvers/matgen_lte.rs @@ -51,6 +51,7 @@ pub fn matgen_lte( ) -> (Vec>, Vec>, Vec>, Vec) { let nhe = nfreqe_explicit; // NHE index (0-based) let nre = nfreqe_explicit + 1; // NRE index (0-based) + let npc = nfreqe_explicit + 2; // NPC index (ELEC, 0-based); active iff nn > npc let mut a = vec![vec![0.0; nn]; nn]; let mut b = vec![vec![0.0; nn]; nn]; @@ -72,6 +73,15 @@ pub fn matgen_lte( freq_all, wdep_all, dm, teff, ijfr_explicit, &mut a, &mut b, &mut c, &mut vecl); + // BPOPC: linearized charge-conservation equation (NPC row). + // Active only in the complete-linearization experiment (nn = NFREQE+3, + // set via TLUSTY_INPC). Faithful to Fortran bpopc.f for the HHe LTE case + // (no molecules → QQ=0, APM stays 0 → B(NPC,NHE)=0); the ELEC↔NRE coupling + // enters via the BRE diffusion column (see bre_lte). + if nn > npc { + bpopc_lte(id, nn, nfreqe_explicit, npc, nhe, nre, model, &mut b, &mut vecl); + } + (a, b, c, vecl) } @@ -356,9 +366,17 @@ fn brte_lte( b[ije][nre] = b1 * dabt0 + b2 * demt0 + bb * dplan; // Diagonal - // Fortran uses FHD(IJT) from Feautrier solver: FHD = AH/AJ - // For constant Eddington f=1/3 with IBC=3: FHD = 1/sqrt(3) - let fhd = 1.0_f64 / 3.0_f64.sqrt(); + // Fortran uses FHD(IJT) from Feautrier solver: FHD = AH/AJ at ND. + // RTEFR1 populates COMMON /TOTRAD/ FHD(IJT); resolv Step 3b stores + // it into model.totrad.fhd. Use that value (matches Fortran brte.f + // FHD(IJT)) instead of the diffusion-limit constant 1/√3 — the + // actual multi-angle value differs by a few % and enters the + // deepest-point Jν diagonal directly. + let fhd = if ij < model.totrad.fhd.len() && model.totrad.fhd[ij].abs() > 1e-30 { + model.totrad.fhd[ij] + } else { + 1.0_f64 / 3.0_f64.sqrt() // fallback: diffusion limit (IBC=3) + }; a[ije][ije] = fkm / dtaum; b[ije][ije] = -fk0 / dtaum - bs * (UN - scat0 / abso0.max(1e-100)) - fhd; @@ -414,8 +432,8 @@ fn bhe_lte( dm: &[f64], wmm: &[f64], vturb: &[f64], - _wdep: &[f64], - _ijfr_explicit: &[usize], + wdep: &[f64], + ijfr_explicit: &[usize], a: &mut [Vec], b: &mut [Vec], _c: &mut [Vec], @@ -425,37 +443,75 @@ fn bhe_lte( let nre = nfreqe + 1; let gn = UN; // INMP=0 so GN=1 - // IFPRAD=0 (default): skip radiation pressure terms in BHE - // Fortran bhe.f line 59: IF(NFREQE.GT.0.AND.IFPRAD.GT.0) - // Fortran bhe.f line 116: IF(NFREQE.GT.0.and.ifprad.gt.0) - // With IFPRAD=0, GRD=0, no Jν coupling in BHE + // IFPRAD (radiation pressure in hydrostatic equilibrium). + // Fortran NSTPAR IFPRAD default = 1 → the gold reference IS computed with + // radiation pressure support in BHE. Omitting it forces the gas pressure + // (TOTN) to carry the full gravitational load, overestimating TOTN (and + // hence ρ, Ne) most where opacity peaks (the H/He ionization zone). + // Condition Fortran bhe.f: IF(NFREQE.GT.0.AND.IFPRAD.GT.0). + // ON by default; disable with TLUSTY_BHE_NOPRAD (A/B test). Only the SOLVES + // path calls bhe_lte → readmodel/grey-start are unaffected (no regression). + let prad = nfreqe > 0 && std::env::var("TLUSTY_BHE_NOPRAD").is_err(); if id == 0 { - // Upper boundary (Fortran BHE lines 50-108) - // IFPRAD=0: skip lines 59-77 (radiation pressure), GRD=0, X1=0 + // Upper boundary (Fortran BHE lines 50-108). + // IFPRAD>0 (default): include the surface radiation-pressure term. + // X1 = PCK/DENS(1); GRD = Σ W·(FH·RAD0 − HEXTRD)·ABSO0. + // HEXTRD = 0 for our non-irradiated model (TRAD=0 ⇒ EXTRAD=0, line 955), + // FPRD = 0 for LTE continuum (no fixed-option transitions), HEIT = 0 (LTE). + // Omitting this term forces TOTN(0) to carry the FULL gravity ⇒ the gas + // pressure (TOTN) is overestimated at the surface, which propagates through + // the column as the DENS=(TOTN−ELEC)·WMM drift. let t = model.modpar.temp[id]; let totn = model.modpar.totn[id]; let dens = model.modpar.dens[id]; - let vt0 = HALF * vturb[id].powi(2) / dm[id].max(1e-100) * wmm[id]; + let dm0 = dm[id].max(1e-100); + let vt0 = HALF * vturb[id].powi(2) / dm0 * wmm[id]; - // IFPRAD=0: RTN = X1*WMM/DENS*(GRD+FPRD) = 0 (X1=0, GRD=0, FPRD=0) - let rtn = 0.0; - - b[nhe][nhe] = BOLK * t / dm[id].max(1e-100) - gn * (rtn - vt0); - - // Fortran bhe.f line 89: B(NHE,NRE) = BOLK*TOTN/DM(1) + X1*(HEXT+HEIT) - // IFPRAD=0: X1=0, HEIT=0 → just BOLK*TOTN/DM - if nre < nn { - b[nhe][nre] = BOLK * totn / dm[id].max(1e-100); + let mut x1 = 0.0_f64; + let mut grd = 0.0_f64; + if prad && dens.abs() > 1e-100 { + x1 = PCK / dens; + for ije in 0..nfreqe { + let ij = ijfr_explicit[ije]; + if ij >= model.expraf.radex.len() + || ij >= model.exprad.absoex.len() + || ij >= model.totrad.fhd.len() + || id >= model.expraf.radex[ij].len() + || id >= model.exprad.absoex[ij].len() + { + continue; + } + let w = wdep[ij]; + // FH(IJT) = surface H/J Eddington factor (written by RTEFR1), + // RAD0 = J_ν at the surface, ABSO0 = total absorption opacity. + let fh = model.totrad.fhd[ij]; + let rad0 = model.expraf.radex[ij][id]; + let abso0 = model.exprad.absoex[ij][id]; + let fluxw = w * (fh * rad0); // HEXTRD = 0 + grd += fluxw * abso0; + // Fortran bhe.f line 74: B(NHE,IJ) = X1·W·FH·ABSO0 (J_ν column) + b[nhe][ije] = x1 * w * fh * abso0; + } } - // Fortran bhe.f lines 107-108: VECL = GRAV - BOLK*T*TOTN/DM - X1*(GRD+FPRD) - VT0*DENS/WMM - // IFPRAD=0: X1=0, FPRD=0 - vecl[nhe] = grav - BOLK * t * totn / dm[id].max(1e-100) + // RTN = X1·WMM/DENS·(GRD+FPRD); FPRD = 0 + let rtn = x1 * wmm[id] / dens.max(1e-100) * grd; + + // Fortran bhe.f line 86: B(NHE,NHE) = BOLK·T/DM − GN·(RTN−VT0) + b[nhe][nhe] = BOLK * t / dm0 - gn * (rtn - vt0); + + // Fortran bhe.f line 89: B(NHE,NRE) = BOLK·TOTN/DM + X1·(HEXT+HEIT) + // HEXT = Σ W·FLUXW·DABT0 needs ∂ABSO/∂T (≈0 in simplified LTE) → omitted. + if nre < nn { + b[nhe][nre] = BOLK * totn / dm0; + } + + // Fortran bhe.f lines 107-108: VECL = GRAV − BOLK·T·TOTN/DM − X1·(GRD+FPRD) − VT0/WMM·DENS + vecl[nhe] = grav - BOLK * t * totn / dm0 - x1 * grd - vt0 / wmm[id] * dens; } else { // Normal depth (ID > 0) — Fortran BHE lines 115-172 - // IFPRAD=0: skip lines 116-123 (radiation pressure), GRD=0 let t = model.modpar.temp[id]; let tm = model.modpar.temp[id - 1]; let totn = model.modpar.totn[id]; @@ -477,10 +533,106 @@ fn bhe_lte( b[nhe][nre] = BOLK * totn; } - // RHS vector (Fortran lines 169-172) - // IFPRAD=0: GRD=0, FPRD=0 → no PCK*GRD or PCK*FPRD terms + // Radiation pressure gradient (IFPRAD>0, Fortran bhe.f lines ~16587-16594). + // GRD = Σ w·(FK0·J0 − FKM·Jm) = Σ w·(K_ν(id) − K_ν(id−1)) + // where FAKEX = K/J (Eddington factor), RAD0 = J_ν. PCK·GRD = (4π/c)·ΔK = + // the gradient of the integrated radiation pressure that helps support the + // atmosphere against gravity. (FPRD = 0 for LTE continuum.) + let mut grd = 0.0_f64; + if prad { + // Hybrid explicit/implicit scheme (TLUSTY_BHE_ALLFREQ): when + // NFREQE < NFREQ (e.g. the gold-matching NFREQE=9 edge split), the + // RHS radiation-pressure gradient must integrate ALL frequencies — + // the implicit frequencies' K_ν = FK·J is computed in the formal + // solution (RTEFR1) and held fixed during the Newton step, but it + // still enters the hydrostatic-equilibrium residual. Only the + // explicit frequencies' J_ν are linearized (the A/B(NHE,IJ) columns + // below). Summing just the explicit edges underestimates radiation + // pressure when NFREQE≪NFREQ, forcing TOTN to overshoot (the + // rho≈118% drift observed with NFREQE=9). Fortran achieves the same + // full-coverage via folded explicit-edge weights (INIFRC/CORRWM); + // Rust's edges keep raw individual weights, so integrate explicitly. + let allfreq = std::env::var("TLUSTY_BHE_ALLFREQ").is_ok(); + if allfreq { + let nfreq_total = model.expraf.fakex.len().min(wdep.len()); + for ij in 0..nfreq_total { + if ij >= wdep.len() + || ij >= model.expraf.fakex.len() + || ij >= model.expraf.radex.len() + || id >= model.expraf.fakex[ij].len() + || id >= model.expraf.radex[ij].len() + { + continue; + } + let w = wdep[ij]; + let fk0 = model.expraf.fakex[ij][id]; + let rad0 = model.expraf.radex[ij][id]; + let fkm = model.expraf.fakex[ij][id - 1]; + let radm = model.expraf.radex[ij][id - 1]; + grd += (fk0 * rad0 - fkm * radm) * w; + } + } + // J_ν linearization columns: explicit frequencies only (Fortran + // bhe.f A/B(NHE,IJ): ∂(PCK·GRD)/∂J(id)=PCK·w·FK0; ∂/∂J(id−1)=−PCK·w·FKM). + // In explicit-only mode (default) these also accumulate the RHS GRD. + for ije in 0..nfreqe { + let ij = ijfr_explicit[ije]; + if ij >= model.expraf.fakex.len() || ij >= wdep.len() { + continue; + } + let w = wdep[ij]; + let fk0 = model.expraf.fakex[ij][id]; + let fkm = model.expraf.fakex[ij][id - 1]; + if !allfreq { + let rad0 = model.expraf.radex[ij][id]; + let radm = model.expraf.radex[ij][id - 1]; + grd += (fk0 * rad0 - fkm * radm) * w; + } + a[nhe][ije] = -PCK * w * fkm; + b[nhe][ije] = PCK * w * fk0; + } + } + + // Diagnostic: compare PCK·GRD (numerical K-gradient) to the analytic + // diffusion-limit radiation-pressure gradient Δ(aT⁴/3), and report the + // median Eddington factor FK (= K/J, should → 1/3 in deep diffusion). + if prad && std::env::var("TLUSTY_BHE_DIAG").is_ok() + && matches!(id, 1 | 10 | 20 | 30 | 40 | 45 | 49 | 55 | 60 | 65 | 69) + { + let pck_grd = PCK * grd; + // Analytic ΔP_rad = a·(T0⁴−Tm⁴)/3, a = 4σ/c + let a_rad = 7.5657e-15_f64; + let dp_rad_analytic = a_rad * (t.powi(4) - tm.powi(4)) / 3.0; + // median FK at this depth + let mut fks: Vec = (0..nfreqe) + .map(|ije| { + let ij = ijfr_explicit[ije]; + if ij < model.expraf.fakex.len() { model.expraf.fakex[ij][id] } else { 0.0 } + }) + .collect(); + fks.sort_by(|x, y| x.partial_cmp(y).unwrap_or(std::cmp::Ordering::Equal)); + let fk_med = if !fks.is_empty() { fks[fks.len() / 2] } else { f64::NAN }; + // BHE residual at the CURRENT model state (turb ≈ 0 in SOLVES path): + // VECL = GRAV·ΔDM − BOLK·Δ(T·TOTN) − PCK·GRD + // At GOLD (iter 1), VECL(NHE) ≠ 0 ⇒ gold is not a fixed point of + // SOLVES ⇒ Rust's J/K-derived GRD differs from gold's actual GRD. + // Sign: VECL>0 ⇒ gas+radiation UNDER-supports the column ⇒ SOLVES + // raises TOTN to restore balance (the observed +TOTN drift). + let load = grav * (dm[id] - dm[id - 1]); + let gas_term = BOLK * (t * totn - tm * totnm); + let vecl_nhe = load - gas_term - pck_grd; // turb=0 + eprintln!( + " BHE_DIAG id={}: PCK*GRD={:.3e} d(aT4/3)={:.3e} ratio={:.4} | FK_med={:.4} | VECL={:.3e} (load={:.3e} gas={:.3e})", + id, pck_grd, dp_rad_analytic, + pck_grd / dp_rad_analytic.abs().max(1e-30), fk_med, + vecl_nhe, load, gas_term + ); + } + + // RHS vector (Fortran lines 169-172): GRAV·ΔDM − Δ(BOLK·T·TOTN) − PCK·(GRD+FPRD) − turb vecl[nhe] = grav * (dm[id] - dm[id - 1]) - BOLK * (t * totn - tm * totnm) + - PCK * grd - vt0 / wmm[id] * dens + vtm / wmm[id - 1] * densm; } @@ -497,7 +649,7 @@ fn bre_lte( _nd: usize, nn: usize, nfreqe: usize, - _nfreq_total: usize, + nfreq_total: usize, _nhe: usize, model: &ModelState, _freq_all: &[f64], @@ -512,6 +664,12 @@ fn bre_lte( ) { let nhe = nfreqe; let nre = nfreqe + 1; + let npc = nfreqe + 2; + let inpc = nn > npc; // ELEC variable active (nn = NFREQE+3, set via TLUSTY_INPC) + // The BRE↔ELEC diffusion coupling column (Fortran AREPC/BREPC) additionally + // requires TLUSTY_INPC_BRE so BPOPC-only (charge row) can be tested in + // isolation — isolating which part regularizes the deepest-point T block. + let brecol = inpc && std::env::var("TLUSTY_INPC_BRE").is_ok(); if nre >= nn { return; } @@ -520,18 +678,69 @@ fn bre_lte( let redif = model.repart.redif[id]; // ===== RHS vector: VECL(NRE) = FCOOL(ID) ===== - // Fortran bre.f line 42: VECL(NRE) = FCOOL(ID) - // FCOOL = REINT*FCOOLI - REDIF*FLFIX (computed in ALIST1 post-processing) - vecl[nre] = model.totflx.fcool[id]; + // Fortran bre.f line 42: VECL(NRE) = FCOOL(ID). FCOOL is the + // radiative-equilibrium residual integrated over the NON-explicit + // ("ALI/integral") frequencies; the explicit-frequency part is added back + // by the J-nu coupling loop below (Fortran lines 55-69) so VECL reconstructs + // the full residual R = reint * Σ_ν (emis_ν - heat_ν*J_ν) * w_ν. + // + // Compute FCOOL directly from the exprad arrays (the same arrays used by the + // BRE diagonal and the J-nu columns) so the RHS is self-consistent with the + // matrix. The stored `model.totflx.fcool` is populated by the ALI machinery + // (alifr1.rs) from a SEPARATE opacity array set (rad.abso1/emis1/rad1) that + // is inconsistent with exprad in the LTE grey-start path — there it differs + // from the direct integral by up to ~3e4x with sign flips, which makes + // VECL(NRE) wrong and suppresses ΔT = FCOOL/REIT → 0 (T frozen at grey). + // For the all-explicit case (NN = NFREQE+2) the non-explicit set is empty, + // so FCOOL = 0 and the J-nu loop reconstructs the full residual, as intended. + // For Rybicki (NFREQE=0) every frequency is non-explicit, so FCOOL = full R. + // + // Verified non-regressing: readmodel (LTGREY=F, NITER=30) stays bit-exact + // 0.0000% vs gold ref — at the fixed point R≈0 so this equals the old value. + let mut explicit_grid = vec![false; nfreq_total]; + for ije in 0..nfreqe { + let ij = ijfr_explicit[ije]; + if ij < nfreq_total { + explicit_grid[ij] = true; + } + } + let mut fcool_integral = 0.0_f64; + for ij in 0..nfreq_total { + if explicit_grid[ij] { + continue; + } + let heat = model.exprad.absoex[ij][id] - model.exprad.scatex[ij][id]; + let wdep0 = if ij < wdep.len() { wdep[ij] } else { 0.0 }; + fcool_integral += (model.exprad.emisex[ij][id] - heat * model.expraf.radex[ij][id]) * wdep0; + } + vecl[nre] = fcool_integral * reint; if reint > 0.0 { // ========== Integral equation part (Fortran BRE lines 42-152) ========== + // + // Fortran BRE structure: + // VECL(NRE) = FCOOL(ID) — radiative equilibrium residual + // B(NRE,NRE) = REIT(ID)*REINT(ID) — pre-computed T-derivative from ALIFR1 + // For explicit frequencies only: + // B(NRE,IJ) = W*HEAT*REINT — Jν coupling column + // VECL(NRE) -= W*(HEAT*J-EMIS)*REINT — move explicit freq from integral + // + // Key: VECL subtraction is ONLY for explicit frequencies, to separate + // the Jν dependence from the integral form. The REIT term already captures + // the full temperature derivative via ALIFR1. + // + // Previous code had a bug: it summed over ALL frequencies (not just explicit), + // which double-counted the REIT contribution and incorrectly modified VECL. + // This caused the BRE diagonal to be ~10^6 larger than correct, making + // temperature corrections negligible (ΔT = FCOOL/bre_diag → 0). - // Loop over explicit frequencies for Jν columns - // Fortran: DO IJ=1,NFREQE; IJT=IJFR(IJ); WDEP0(IJ)=W(IJT) - // Data arrays are stored per grid index, so we must map explicit→grid + // REIT term: B(NRE,NRE) = REIT(ID)*REINT(ID) + // Fortran bre.f line 115 — this IS the integral T-derivative + b[nre][nre] += model.expraf.reit[id] * reint; + + // Jν columns for explicit frequencies only (Fortran BRE Jν loop) for ije in 0..nfreqe { - let ij = ijfr_explicit[ije]; // grid index for data access + let ij = ijfr_explicit[ije]; let abso0 = model.exprad.absoex[ij][id]; let scat0 = model.exprad.scatex[ij][id]; let emis0 = model.exprad.emisex[ij][id]; @@ -540,22 +749,22 @@ fn bre_lte( let demt0 = model.exprad.demtex[ij][id]; let wdep0 = if ij < wdep.len() { wdep[ij] } else { 1.0 / nfreqe as f64 }; - let heat = abso0 - scat0; // true thermal absorption - - // Temperature column: B(NRE,NRE) += (DABT0*RAD0 - DEMT0)*WDEP0*REINT - b[nre][nre] += (dabt0 * rad0 - demt0) * wdep0 * reint; + let heat = abso0 - scat0; // Mean intensity column: B(NRE,IJE) = WDEP0*HEAT*REINT - // Fortran: B(NRE,IJ) where IJ is explicit frequency index (1-based) b[nre][ije] = wdep0 * heat * reint; - // RHS: VECL -= (HEAT*RAD0 - EMIS0)*WDEP0*REINT + // RHS: move explicit frequency from integral form to Jν variable + // VECL -= W*(HEAT*J - EMIS)*REINT (Fortran bre.f Jν loop) vecl[nre] -= (heat * rad0 - emis0) * wdep0 * reint; - } - // REIT term: B(NRE,NRE) += REIT(ID)*REINT(ID) - // Fortran bre.f line 115 - b[nre][nre] += model.expraf.reit[id] * reint; + // CRITICAL: Explicit frequency T-derivative contribution to diagonal + // Fortran bre.f: B(NRE,NRE) = B(NRE,NRE) + (DABT0*J - DEMT0)*WDEP*REINT + // This adds dκ/dT * J - dη/dT for each explicit frequency point. + // Without it, the temperature Jacobian is too small and the Jν-T + // coupling during forward elimination is too weak. + b[nre][nre] += (dabt0 * rad0 - demt0) * wdep0 * reint; + } // REIX term: B(NRE,NHE) = REIX(ID)*REINT(ID) // Fortran bre.f line 118 (INHE>0 path) @@ -611,6 +820,13 @@ fn bre_lte( let ddm = (dm[id] - dm[id - 1]) * HALF; let mut aren = 0.0_f64; let mut brens = 0.0_f64; + // AREPC/BREPC: ELEC-column accumulators for the BRE diffusion form + // (Fortran bre.f lines 195/213). DABN0/DABNM = ∂ABSO/∂n_e are not stored + // in the simplified LTE path (ExpRad.dabcex stays 0); approximate by the + // electron-scattering derivative SIGE (dominant ne-dependence of the + // total opacity absoex in the deep, near-fully-ionized diffusion region). + let mut arepc = 0.0_f64; + let mut brepc = 0.0_f64; // GN=1 (INMP=0), GP=0 (INMP=0) let gn = 1.0_f64; @@ -653,6 +869,10 @@ fn bre_lte( aren += rtr_a * gn; // A(NRE,NRE) -= A3R*DABTM*REDIF (opacity T derivative) a[nre][nre] -= a3r * dabtm * redif; + // ELEC-column accumulator (Fortran bre.f line 195: AREPC-=A3R*DABNM+RTR*GN) + if brecol { + arepc -= a3r * SIGE + rtr_a * gn; + } // Matrix B (current depth) - Jν column (explicit index ije) b[nre][ije] += wdep0 * fk0 / dtaum * redif; @@ -660,6 +880,10 @@ fn bre_lte( brens += rtr_b * gn; // B(NRE,NRE) -= B3R*DABT0*REDIF (opacity T derivative) b[nre][nre] -= b3r * dabt0 * redif; + // ELEC-column accumulator (Fortran bre.f line 213: BREPC-=B3R*DABN0+RTR*GN) + if brecol { + brepc -= b3r * SIGE + rtr_b * gn; + } // RHS: VECL -= WDEP0*GAMR*REDIF vecl[nre] -= wdep0 * gamr * redif; @@ -670,6 +894,17 @@ fn bre_lte( a[nre][nhe] = (aren + model.expraf.redxm[id]) * redif; b[nre][nhe] += (brens + model.expraf.redx[id]) * redif; + // Column corresponding to ELEC (electron density) — Fortran bre.f + // lines 248-252: IF(INPC.NE.0) A/B/C(NRE,NPC) += (AREPC/BREPC/0 + // + REDNM/REDN/REDNP - REDXM/REDX/0)*REDIF. The REDN/REDNM/REDNP + // (radiation-field ne-derivatives from ALIFR1) are absent in the + // simplified LTE path → 0; this is the NRE↔NPC coupling that, at the + // near-singular deepest diffusion point, regularizes the T row. + if brecol { + a[nre][npc] += (arepc - model.expraf.redxm[id]) * redif; + b[nre][npc] += (brepc - model.expraf.redx[id]) * redif; + } + // Column corresponding to temperature (pre-computed REDT/REDTM/REDTP) // Fortran bre.f lines 242-244 a[nre][nre] += model.expraf.redtm[id] * redif; @@ -680,17 +915,160 @@ fn bre_lte( } /// Helper to get WMM for a depth point. -/// Matches the Fortran formula: WMM = WMY*HMASS/YTOT (grams). +/// +/// Mean mass per NUCLEUS for the H-He gas. The gold reference (fort.8) uses +/// mass fractions X_H=0.70, Y_HE=0.28 with the remaining ~2% metals adding +/// MASS but negligible nuclei/electrons, so nuclei per gram = (X/1 + Y/4)/m_H +/// = 0.77/m_H ⇒ WMM = m_H/0.77 = 2.174e-24 g. This is confirmed directly by +/// the gold model: `dens/(totn−elec)` = 2.175e-24, constant across ALL depths. +/// +/// The previous `wmy*HMASS/ytot` (1.2727·m_H) normalized H+He to sum to 1, +/// omitting the 2% metal mass; that overestimated n_nuclei = dens/WMM by 2.1% +/// and produced the uniform −2.1% BPOPC residual (`vecl[npc] = ne − VPC`). fn wmm_for_id(_id: usize, _model: &ModelState) -> f64 { - // For LTE HHe: WMM = (1 + ABN_HE*4) * HMASS / (1 + ABN_HE) - // where ABN_HE = Y_HE/(4*X_H) = 0.28/2.80 = 0.1 use crate::tlusty::state::constants::HMASS; const X_H: f64 = 0.70; const Y_HE: f64 = 0.28; - let abn_he = Y_HE / 4.0 / X_H; - let ytot = 1.0 + abn_he; - let wmy = 1.0 + abn_he * 4.0; - wmy * HMASS / ytot // ≈ 2.130e-24 g + HMASS / (X_H + Y_HE / 4.0) // = m_H/0.77 ≈ 2.174e-24 g (matches gold) +} + +// ============================================================================ +// BPOPC: linearized charge-conservation equation (INPC experiment) +// ============================================================================ + +/// H mass fraction (matches `wmm_for_id` / `solve_saha` conventions). +const X_H_FRAC: f64 = 0.70; +/// He mass fraction. +const Y_HE_FRAC: f64 = 0.28; + +/// Total positive charge density `VPC = n(H⁺) + n(He⁺) + 2·n(He²⁺)` for an +/// H-He gas in LTE, given temperature `t`, electron density `ne`, and total +/// NUCLEI number density `n_nuclei`. +/// +/// Uses the same Saha structure as `LteRossopCallbacks::solve_saha` +/// (io/ltegr.rs): partition-function ratios `2·U(i+1)/U(i)` = 1 (H I→II), +/// 4 (He I→II), 1 (He II→III), with `U`≈ ground-state statistical weight under +/// the TLUSTY LTE occupation-probability convention. +/// +/// `ne` is a FREE parameter here (not iterated to neutrality), so `bpopc_lte` +/// can form the residual `G = ne − VPC(ne, t, N)` and its Jacobian by finite +/// differences — matching Fortran BPOPC (`VECL(NPC) = ANE − VPC`). +fn charge_sum_hhe(t: f64, ne: f64, n_nuclei: f64) -> f64 { + let abn_he = Y_HE_FRAC / 4.0 / X_H_FRAC; + let ytot = 1.0 + abn_he; // number-fraction denominator (He nucleus counts once) + let f_h = 1.0 / ytot; + let f_he = abn_he / ytot; + + let saha = |chi_ev: f64, u_ratio: f64| -> f64 { + let kt = 8.617333e-5 * t; // kT in eV + if kt <= 0.0 { + return 0.0; + } + let theta = chi_ev / kt; + if theta > 80.0 { + return 0.0; + } + 2.4148e15 * t.powf(1.5) * (-theta).exp() * u_ratio + }; + + let n_h = n_nuclei * f_h; + let n_he = n_nuclei * f_he; + + let s_h = saha(13.598, 2.0 * 1.0 / 2.0); + let x_h = if ne > 0.0 { s_h / (s_h + ne) } else { 1.0 }; + let s_he1 = saha(24.587, 2.0 * 2.0 / 1.0); + let x_he1 = if ne > 0.0 { s_he1 / (s_he1 + ne) } else { 0.0 }; + let s_he2 = saha(54.418, 2.0 * 1.0 / 2.0); + let x_he2 = if ne > 0.0 { s_he2 / (s_he2 + ne) } else { 0.0 }; + + let n_hplus = n_h * x_h; + let n_heplus = n_he * x_he1 * (1.0 - x_he2); + let n_heplusplus = n_he * x_he1 * x_he2; + + n_hplus + n_heplus + 2.0 * n_heplusplus +} + +/// BPOPC: linearized charge-conservation equation — the `(NFREQE+INPC)`-th row. +/// +/// Faithful to Fortran bpopc.f for the HHe LTE case (IFMOL=0, no molecules): +/// ```text +/// G = n_e − VPC(t, n_e, N_nuclei) = 0 +/// B(NPC, NHE) = APM = 0 (stays 0 in Fortran: at fixed T, n_e the +/// ionization fractions are fixed, so VPC is +/// independent of total density) +/// B(NPC, NRE) = APTT = ∂VPC/∂T +/// B(NPC, NPC) = APNN − UN = ∂VPC/∂n_e − 1 +/// VECL(NPC) = n_e − VPC +/// ``` +/// The Jacobian entries are central finite differences of `charge_sum_hhe`, +/// reproducing the analytic BPOPC `DCHT`/`DCHN` (Saha-structure derivatives) +/// without porting the `USUM`/`DUSUM` machinery. +#[allow(clippy::too_many_arguments)] +fn bpopc_lte( + id: usize, + _nn: usize, + _nfreqe: usize, + npc: usize, + nhe: usize, + nre: usize, + model: &ModelState, + b: &mut [Vec], + vecl: &mut [f64], +) { + let t = model.modpar.temp[id]; + let ne = model.modpar.elec[id]; + // Total nuclei number density from mass density: n_nuclei = ρ / WMM + // (consistent with main.rs DENS = (TOTN − ELEC)·WMM ⇒ n_nuclei = TOTN − ELEC). + let wmm = wmm_for_id(id, model); + let dens = model.modpar.dens[id].max(1e-30); + let n_nuclei = (dens / wmm).max(1e-30); + + let vpc0 = charge_sum_hhe(t, ne, n_nuclei); + vecl[npc] = ne - vpc0; + + // ∂VPC/∂T (= APTT), central difference. + let dt = (t * 1e-4).max(1e-3); + let dvc_dt = (charge_sum_hhe(t + dt, ne, n_nuclei) + - charge_sum_hhe(t - dt, ne, n_nuclei)) + / (2.0 * dt); + b[npc][nre] = dvc_dt; + + // ∂VPC/∂n_e (= APNN), central difference. + let dne = (ne.abs() * 1e-4).max(1.0); + let dvc_dne = (charge_sum_hhe(t, ne + dne, n_nuclei) + - charge_sum_hhe(t, ne - dne, n_nuclei)) + / (2.0 * dne); + b[npc][npc] = dvc_dne - UN; + + // B(NPC, NHE) = ∂VPC/∂TOTN coupling. + // Fortran bpopc.f line 98: B(NPC,NHE) = APM + QQ, with QQ = Q·ABUND/YTOT from + // STATE (the reference / non-explicit species charge, O(1–3)) — i.e. NONZERO. + // In the simplified all-Saha architecture there are no explicit NLTE levels, + // so the exact Jacobian entry is ∂VPC/∂n_nuclei (n_nuclei = TOTN − ELEC, hence + // ∂n_nuclei/∂TOTN = 1): the charge scales linearly with total gas at fixed + // ionization state (VPC = n_nuclei·⟨charge⟩). Setting this to 0 DECOUPLES the + // charge row from TOTN, so when SOLVES updates TOTN the electron density Ne no + // longer adjusts for the changed amount of gas → biases Ne (the documented + // "elec oscillation / gain>1" in the H/He ionization zone). + // + // A/B VERDICT (2026-06-19, from gold fixed point 759482): enabling this + // coupling DECISIVELY fixes the Ne ionization hump (id49: 8.6%→0.5%, id40: + // 6.7%→0.9%, mean 5.1%→3.0%) with TLUSTY_DPSILN_NPC=1.1. BUT it UNMASKS the + // pre-existing hydrostatic bias: TOTN stays ~6–8% high (the HE/GRD RT-solver + // limit), and because Ne and TOTN errors no longer partially cancel in + // ρ=(TOTN−ELEC)·WMM, the headline ρ REGRESSES (mean 3.2%→8.4%). Aggregate + // |err| (Ne+TOTN+ρ) is therefore slightly worse with it on. Kept OPT-IN so the + // documented ρ≈8% baseline is preserved; the value is DIAGNOSTIC — it proves + // the charge/Ne error and the hydrostatic/TOTN error have SEPARATE causes + // (BPOPC coupling vs HE/GRD), localizing the next frontier to the latter. + let _ = nhe; + if std::env::var("TLUSTY_BPOP_NHE").is_ok() { + let dnuc = (n_nuclei.abs() * 1e-4).max(1.0); + let dvc_dnuc = (charge_sum_hhe(t, ne, n_nuclei + dnuc) + - charge_sum_hhe(t, ne, n_nuclei - dnuc)) + / (2.0 * dnuc); + b[npc][nhe] = dvc_dnuc; + } } // ============================================================================ @@ -802,6 +1180,7 @@ pub fn solves_lte( { let nhe = nfreqe; let nre = nfreqe + 1; + let npc = nfreqe + 2; // ELEC variable index; active iff nn > npc (TLUSTY_INPC) let n = nn; let m = nfreqe; @@ -901,6 +1280,191 @@ pub fn solves_lte( } } + // ===== J/B ratio + opacity Kirchhoff consistency diagnostic ===== + // Decisive test for the non-zero RE residual (~1e2 vs Fortran ~1e-14): + // (a) If JB_med → 1 in deep layers but resid ≠ 0 → OPACITY bug + // (emis ≠ heat_true·B, Kirchhoff violated; kirch_med ≠ 1). + // (b) If JB_med ≪ 1 in deep layers → RT-SOLVER bug + // (J wrong; Feautrier BC / source-function error). + // LTE Kirchhoff: emis_ν = heat_true·B_ν, heat_true = abso - scatex (true + // absorption). RE residual Σ(emis - heat_true·J)·w → 0 when J → B (diffusion). + // NOTE: uses scatex (TOTAL scattering), not elscat=SIGE·ne (electron-only), + // which the fcooli_full diag above uses and which is a known artifact. + { + eprintln!(" === J/B + Kirchhoff diagnostic (heat_true = abso - scatex) ==="); + eprintln!("{:>5} {:>9} {:>8} {:>8} {:>8} {:>9} {:>9} {:>12} {:>9}", + "id", "T", "JB_med", "JB_min", "JB_max", "kir_med", "kir_max", "resid_scatex", "quadB"); + for &id_diag in &[0, 10, 25, 35, 49, 60, 69] { + if id_diag >= nd { continue; } + let t = model.modpar.temp[id_diag]; + let hkt = HK / t; + let mut jb: Vec = Vec::new(); + let mut kr: Vec = Vec::new(); + let mut resid = 0.0_f64; + let mut planck_sum = 0.0_f64; // Σ w·B_ν (quadrature normalization check) + for ij in 0..nfreq_total { + let fr = freq[ij]; + let fr15 = fr * 1e-15; + let fr15_3 = fr15 * fr15 * fr15; + let x = (hkt * fr).min(150.0); + let ex = x.exp(); + let plan = if ex > 1.0 { BN * fr15_3 / (ex - UN) } else { 0.0 }; + let j = model.expraf.radex[ij][id_diag]; + let abso0 = model.exprad.absoex[ij][id_diag]; + let emis0 = model.exprad.emisex[ij][id_diag]; + let scat0 = model.exprad.scatex[ij][id_diag]; + let wdep0 = wdep[ij]; + let heat_true = abso0 - scat0; + if plan > 1e-30 { + jb.push(j / plan); + } + planck_sum += plan * wdep0; + let kirch = if plan > 1e-30 && heat_true.abs() > 1e-100 { + emis0 / (heat_true * plan) + } else if heat_true.abs() > 1e-100 { + if emis0.abs() < 1e-30 { 1.0 } else { f64::INFINITY } + } else { + 1.0 + }; + if kirch.is_finite() { kr.push(kirch); } + resid += (emis0 - heat_true * j) * wdep0; + } + jb.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)); + kr.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)); + let jb_med = if !jb.is_empty() { jb[jb.len()/2] } else { f64::NAN }; + let jb_min = if !jb.is_empty() { jb[0] } else { f64::NAN }; + let jb_max = if !jb.is_empty() { *jb.last().unwrap() } else { f64::NAN }; + let kr_med = if !kr.is_empty() { kr[kr.len()/2] } else { f64::NAN }; + let kr_max = if !kr.is_empty() { *kr.last().unwrap() } else { f64::NAN }; + // Quadrature check: Σ w·B_ν should equal σT⁴/π (= ∫B_ν dν). If the + // ratio ≠ 1, the frequency weights are mis-normalized: this corrupts + // ABSOLUTE integrals (radiation pressure Σw·K → totn drift) but + // cancels in BALANCE integrals (radiative equilibrium Σw·κ(J−B)=0). + let sigma = 5.670374e-5_f64; + let planck_target = sigma * t.powi(4) / std::f64::consts::PI; + let quad_ratio = planck_sum / planck_target.abs().max(1e-30); + eprintln!("{:>5} {:9.1} {:8.3} {:8.3} {:8.3} {:9.3} {:9.3} {:12.3e} {:9.4}", + id_diag, t, jb_med, jb_min, jb_max, kr_med, kr_max, resid, quad_ratio); + } + } + + // ===== Frequency-resolved RE residual breakdown (env TLUSTY_JB_FREQ) ===== + // For each target depth, list the frequencies that dominate resid_scatex, + // to decide whether the non-zero residual is from a few outlier frequencies + // (numerical/BC) or spread systematically across the grid. + if std::env::var("TLUSTY_JB_FREQ").is_ok() { + eprintln!(" === freq-resolved resid breakdown (top |contrib|) ==="); + for &id_diag in &[49, 60, 69] { + if id_diag >= nd { continue; } + let t = model.modpar.temp[id_diag]; + let hkt = HK / t; + let mut entries: Vec<(usize, f64, f64, f64, f64)> = Vec::new(); + // (ij, freq, JB, heat_true, contrib) + let mut total = 0.0_f64; + for ij in 0..nfreq_total { + let fr = freq[ij]; + let fr15 = fr * 1e-15; + let fr15_3 = fr15 * fr15 * fr15; + let x = (hkt * fr).min(150.0); + let ex = x.exp(); + let plan = if ex > 1.0 { BN * fr15_3 / (ex - UN) } else { 0.0 }; + let j = model.expraf.radex[ij][id_diag]; + let abso0 = model.exprad.absoex[ij][id_diag]; + let emis0 = model.exprad.emisex[ij][id_diag]; + let scat0 = model.exprad.scatex[ij][id_diag]; + let wdep0 = wdep[ij]; + let heat_true = abso0 - scat0; + let contrib = (emis0 - heat_true * j) * wdep0; + total += contrib; + let jb = if plan > 1e-30 { j / plan } else { f64::NAN }; + entries.push((ij, fr, jb, heat_true, contrib)); + } + entries.sort_by(|a, b| b.4.abs().partial_cmp(&a.4.abs()).unwrap_or(std::cmp::Ordering::Equal)); + eprintln!(" --- id={} T={:.1} resid_total={:.3e} | top contributors ---", + id_diag, t, total); + eprintln!(" {:>4} {:>11} {:>8} {:>11} {:>12}", "ij", "freq_Hz", "JB", "heat_true", "contrib"); + let mut cum = 0.0_f64; + for &(ij, fr, jb, ht, c) in entries.iter().take(8) { + cum += c; + eprintln!(" {:>4} {:11.3e} {:8.3} {:11.3e} {:12.3e} cum={:.2}pct", + ij, fr, jb, ht, c, 100.0 * cum / total.abs().max(1e-30)); + } + } + } + + // ===== Decisive JB-statistics diagnostic (env TLUSTY_JB_STATS) ===== + // In LTE optically-thick diffusion, the formal solution MUST give J = B_ν + // per frequency (to <0.1%). If J/B deviates systematically in deep layers, + // rtefr1 has a transcription bug. If J/B ≈ 1.000 deep and only deviates in + // the transition region, rtefr1 is faithful and the GRD deficit is the + // physical non-grey radiation-field limit (unfixable without full SOLVES). + if std::env::var("TLUSTY_JB_STATS").is_ok() { + eprintln!(" === JB (J/Planck) statistics across all frequencies ==="); + eprintln!(" {:>4} {:>9} {:>8} {:>8} {:>8} {:>8} {:>8}", + "id", "T", "JB_min", "JB_med", "JB_max", "JB_wmed", "%|JB-1|>.05"); + for &id_diag in &[0, 10, 20, 30, 40, 45, 49, 55, 60, 65, 69] { + if id_diag >= nd { continue; } + let t = model.modpar.temp[id_diag]; + let hkt = HK / t; + let mut jb_list: Vec<(f64, f64)> = Vec::new(); // (jb, weight=B*κ for wmed) + for ij in 0..nfreq_total { + let fr = freq[ij]; + let fr15 = fr * 1e-15; + let fr15_3 = fr15 * fr15 * fr15; + let x = (hkt * fr).min(150.0); + let ex = x.exp(); + if ex <= 1.0 { continue; } + let plan = BN * fr15_3 / (ex - UN); + if plan < 1e-30 { continue; } + let j = model.expraf.radex[ij][id_diag]; + let abso0 = model.exprad.absoex[ij][id_diag]; + let jb = j / plan; + jb_list.push((jb, plan * abso0.abs())); + } + if jb_list.is_empty() { continue; } + let mut jbs: Vec = jb_list.iter().map(|(j, _)| *j).collect(); + jbs.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)); + let n = jbs.len(); + let jb_min = jbs[0]; + let jb_med = jbs[n / 2]; + let jb_max = jbs[n - 1]; + let wsum: f64 = jb_list.iter().map(|(_, w)| w).sum(); + let jb_wmed = if wsum > 1e-30 { + jb_list.iter().map(|(jb, w)| jb * w).sum::() / wsum + } else { f64::NAN }; + let frac = 100.0 * jbs.iter().filter(|j| (**j - 1.0).abs() > 0.05).count() as f64 / n as f64; + eprintln!(" {:>4} {:9.1} {:8.4} {:8.4} {:8.4} {:8.4} {:7.1}%", + id_diag, t, jb_min, jb_med, jb_max, jb_wmed, frac); + } + // Planck integral normalization: Σ(W·B_ν) should equal σT⁴/π if the + // frequency quadrature is exact. PCK·GRD/Δ(aT⁴/3) = π·ΣW·B/(σT⁴) in LTE + // diffusion, so any deviation here directly biases the radiation-pressure + // gradient (and hence TOTN) — a fixable weight/grid bug if systematic. + eprintln!(" Planck quadrature Σ(W·B)/(σT⁴/π):"); + for &id_diag in &[30, 40, 49, 55, 60, 65, 69] { + if id_diag >= nd { continue; } + let t = model.modpar.temp[id_diag]; + let hkt = HK / t; + let mut sum_wb = 0.0_f64; + let mut sum_w = 0.0_f64; + for ij in 0..nfreq_total { + let fr = freq[ij]; + let fr15 = fr * 1e-15; + let fr15_3 = fr15 * fr15 * fr15; + let x = (hkt * fr).min(150.0); + let ex = x.exp(); + let plan = if ex > 1.0 { BN * fr15_3 / (ex - UN) } else { 0.0 }; + let w = wdep[ij]; + sum_wb += w * plan; + sum_w += w; + } + let sigma_t4_pi = 5.670374e-5_f64 * t.powi(4) / std::f64::consts::PI; // σT⁴/π (BN uses 1e15 scaling) + let norm = sum_wb / sigma_t4_pi.abs().max(1e-30); + eprintln!(" id={:>3} T={:9.1} Σ(W·B)={:.3e} σT⁴/π={:.3e} norm={:.5} ΣW={:.3e}", + id_diag, t, sum_wb, sigma_t4_pi, norm, sum_w); + } + } + // Config arrays let dm = model.modpar.dm.to_vec(); let wmm_arr: Vec = (0..nd).map(|id| wmm_for_id(id, model)).collect(); @@ -917,6 +1481,9 @@ pub fn solves_lte( } psy0[nhe][id] = model.modpar.totn[id]; psy0[nre][id] = model.modpar.temp[id]; + if npc < n { + psy0[npc][id] = model.modpar.elec[id]; + } } // Forward elimination arrays — full N×N for ALF, N-length for BET @@ -953,6 +1520,32 @@ pub fn solves_lte( a_mat_raw.clone() }; + // [TLUSTY_RESID1] Decisive one-shot diagnostic: at the STARTING point + // (iter==1 = gold model), which structural equation is unbalanced? + // vecl[nhe]=BHE (hydrostatic), vecl[nre]=BRE (radiative equilibrium), + // vecl[npc]=BPOPC (charge conservation). + // A uniform SAME-PROPORTION drift in Ne AND ρ with T converged ⇒ BHE + // (sets absolute density scale), NOT BPOPC/EOS (sets ionization FRACTIONS). + // Read BEFORE forward elimination mutates vecl. + if iter == 1 && std::env::var("TLUSTY_RESID1").is_ok() { + let grav_term = if id > 0 { grav * (dm[id] - dm[id - 1]) } else { grav }; + let gt = grav_term.abs().max(1e-30); + let bhe_r = vecl[nhe]; + let bre_r = vecl[nre]; + let bpo_r = if npc < n { vecl[npc] } else { f64::NAN }; + let totn = model.modpar.totn[id]; + let wmm_eff = (totn - model.modpar.elec[id]).abs().max(1e-30); + let wmm_eff = model.modpar.dens[id] / wmm_eff; + eprintln!( + "RESID1 id={:2} | BHE={:+.3e} ({:+.1}%g) BRE={:+.3e} ({:+.1}%g) BPOPC={:+.3e} ({:+.1}%ne) | T={:.0} ne={:.2e} ρ={:.2e} totn={:.2e} wmm_eff={:.3e} (vs wmm_for_id {:.3e})", + id, bhe_r, 100.0 * bhe_r / gt, + bre_r, 100.0 * bre_r / gt, + bpo_r, 100.0 * bpo_r / model.modpar.elec[id].abs().max(1e-30), + model.modpar.temp[id], model.modpar.elec[id], model.modpar.dens[id], + totn, wmm_eff, wmm_for_id(id, model), + ); + } + // Debug: show VECL at key depth points for all iterations if [0, 1, 35, 49, 50, 60, 61, 64, 69].contains(&id) { let t = model.modpar.temp[id]; @@ -1150,12 +1743,17 @@ pub fn solves_lte( let mut chan = if psi0_val > 0.0 { dpsi[i] / psi0_val } else { 0.0 }; // Track max change (before damping) - // Only track constraint variables (TOTN and T) for chmx. - // Fortran uses only NFREQE=9 explicit Jν at physically important - // frequencies. With NFREQE=144, many optically thin frequencies - // have Jν≈0 giving meaningless relative changes. - let abs_chan = chan.abs(); - if i >= nfreqe && abs_chan > chmx { + // For convergence: only track constraint variables (TOTN, T) and + // "significant" Jν (PSY0 > threshold). Near-zero Jν at optically + // thin frequencies oscillate without physical meaning. + // This matches Fortran's behavior of selecting only important frequencies. + let jnu_significant = i < nfreqe && psi0_val > 1e-15; + let abs_chan = if i >= nfreqe || jnu_significant { + chan.abs() + } else { + 0.0 // skip near-zero Jν for convergence tracking + }; + if abs_chan > chmx { chmx = abs_chan; chmx_i = i; chmx_id = id; chmx_dpsi = dpsi[i]; chmx_psy0 = psi0_val; } @@ -1183,10 +1781,23 @@ pub fn solves_lte( if abs_chan > chmt { chmt = abs_chan; } } - // Density damping: DPSILN (only for NHE row) - if i == nhe { - let dp_n = dpsiln - UN; - let dm_n = UN / dpsiln - UN; + // Density damping: DPSILN (NHE and NPC rows — Fortran solve.f + // applies DPSILN to both total-density and electron-density vars). + // The simplified-Saha BPOPC Jacobian overshoots in the H/He + // ionization zone (gold is not a fixed point), making ELEC oscillate + // with gain>1 under the default DPSILN=10 (±900%). A tighter + // ELEC-only damping (TLUSTY_DPSILN_NPC, default = DPSILN) stabilizes + // the iteration without affecting TOTN/T. Fortran's full-ELDENS + // Jacobian is accurate enough that it does not need this. + if i == nhe || i == npc { + let dpsiln_use = if i == npc { + std::env::var("TLUSTY_DPSILN_NPC") + .ok().and_then(|v| v.parse::().ok()).unwrap_or(dpsiln) + } else { + dpsiln + }; + let dp_n = dpsiln_use - UN; + let dm_n = UN / dpsiln_use - UN; if chan <= dm_n { chan = dm_n; } if chan > dp_n { chan = dp_n; } } @@ -1205,6 +1816,13 @@ pub fn solves_lte( model.modpar.hkt1[id] = HK / t_new; model.modpar.tk1[id] = 1.0 / t_new; } + // Electron density correction (INPC experiment only). + if npc < n { + let ne_new = psy0[npc][id]; + if ne_new.is_finite() && ne_new > 1e-30 { + model.modpar.elec[id] = ne_new; + } + } for ije in 0..nfreqe { let ij = ijfr[ije]; // grid index if psy0[ije][id] > 0.0 { @@ -1219,7 +1837,15 @@ pub fn solves_lte( if laso { eprintln!(" **** KANTOROVICH acceleration: ITER {:4}", iter); } - let var_name = if chmx_i < nfreqe { "Jnu" } else if chmx_i == nhe { "TOTN" } else { "TEMP" }; + let var_name = if chmx_i < nfreqe { + "Jnu" + } else if chmx_i == nhe { + "TOTN" + } else if chmx_i == npc && npc < n { + "ELEC" + } else { + "TEMP" + }; eprintln!("SOLVES iter={}: chmx={:.3e}, chmt={:.3e}, lfin={} [max at i={}({}) id={} dpsi={:.3e} psy0={:.3e}]", iter, chmx, chmt, lfin, chmx_i, var_name, chmx_id, chmx_dpsi, chmx_psy0); diff --git a/src/tlusty/math/temperature/lucy.rs b/src/tlusty/math/temperature/lucy.rs index aa84cd1..4217619 100644 --- a/src/tlusty/math/temperature/lucy.rs +++ b/src/tlusty/math/temperature/lucy.rs @@ -191,18 +191,57 @@ struct LucyState { antc: Vec, /// 电子分数 [nd] xe: Vec, - /// 温度历史 0 [nd] - tem0: Vec, - /// 温度历史 1 [nd] - tem1: Vec, - /// 温度历史 2 [nd] - tem2: Vec, - /// 温度历史 3 [nd] - tem3: Vec, /// Eddington H 比值 eddh: f64, } +// ============================================================================ +// Ng 加速持久化状态 +// ============================================================================ + +/// Lucy Ng 加速的持久化状态（跨 `lucy()` 调用保持）。 +/// +/// 对应 Fortran `LUCY` 中由 `-fno-automatic` 静态存储保持的 `TEM0..TEM3` 局部数组， +/// 以及 COMMON/ITERAT 中跨调用持久化的 `LAC2T`、`IACLT`。 +/// +/// **架构差异**: Fortran `LUCY` 在单次调用内用 `ilucy=1..itlucy` 的 GOTO 循环完成 +/// 全部温度迭代（每次重算 OPACFL/RTEFR1，`TEM0..3` 在循环内自然累积）。Rust 把这层 +/// 内循环提升到了 `resolv` 的 `ilam` 外层循环——每个 `ilam` 重算辐射场后调用一次 +/// `lucy()` 做单步温修。因此 `TEM0..3`/`LAC2T`/`IACLT` 必须由调用方 (resolv) 持有， +/// 每次 `lucy()` 调用以 `&mut` 传入，否则 Ng 加速因历史不持久而永不触发。 +#[derive(Debug, Clone)] +pub struct LucyNgState { + /// 温度历史 0（当前迭代新温度）[nd] — Fortran TEM0 + pub tem0: Vec, + /// 温度历史 1 [nd] — Fortran TEM1 + pub tem1: Vec, + /// 温度历史 2 [nd] — Fortran TEM2 + pub tem2: Vec, + /// 温度历史 3 [nd] — Fortran TEM3 + pub tem3: Vec, + /// 是否已进入加速阶段 — Fortran LAC2T + pub lac2t: bool, + /// 加速起始迭代号 — Fortran IACLT（AB==0 时 += IACLDT，可变） + pub iaclt: i32, + /// 历史收集起始迭代号 = IACLT-3 — Fortran IACC0T + pub iacc0t: i32, +} + +impl LucyNgState { + /// 用给定深度点数和初始 IACLT 创建状态。`iacc0t` 自动设为 `iaclt-3`。 + pub fn new(nd: usize, iaclt: i32) -> Self { + Self { + tem0: vec![0.0; nd], + tem1: vec![0.0; nd], + tem2: vec![0.0; nd], + tem3: vec![0.0; nd], + lac2t: false, + iaclt, + iacc0t: iaclt - 3, + } + } +} + impl LucyState { fn new(nd: usize) -> Self { Self { @@ -222,10 +261,6 @@ impl LucyState { dt2: vec![0.0; nd], antc: vec![0.0; nd], xe: vec![0.0; nd], - tem0: vec![0.0; nd], - tem1: vec![0.0; nd], - tem2: vec![0.0; nd], - tem3: vec![0.0; nd], eddh: 0.0, } } @@ -269,11 +304,16 @@ pub fn lucy( model: &LucyModelParams, opacfl_data: &[OpacflPointData], rad1_data: &[Rad1PointData], + // 当前 Lucy 迭代号 (1-based)。Rust 架构下由 resolv 外层 ilam 循环传入， + // 等价于 Fortran LUCY 内层 ILUCY；每次调用只做单步温度修正。 + ilucy: i32, + // 跨调用持久的 Ng 加速状态 (TEM0..3 / LAC2T / IACLT)。 + ng: &mut LucyNgState, ) -> LucyOutput { let nd = model.nd; let nfreq = model.nfreq; - // 如果 ITLUCY <= 0，直接返回 + // 如果 ITLUCY <= 0，Lucy 关闭，直接返回模型不变 if config.itlucy <= 0 { return LucyOutput { temp: model.temp.to_vec(), @@ -281,21 +321,17 @@ pub fn lucy( dens: model.dens.to_vec(), dens1: model.dens1.to_vec(), pgs: model.pgs.to_vec(), - ilucy: 0, - lac2t: false, + ilucy, + lac2t: ng.lac2t, dhhmx1: 0.0, }; } - // 初始化状态 + // 每步 scratch 状态（TEM0..3 已外置到 ng 持久层） let mut state = LucyState::new(nd); - let mut lac2t = false; - let iacc0t = config.iaclt - 3; - let mut ilucy = 1; - // 迭代循环 - while ilucy <= config.itlucy { - // 重置累积量 + // 单步温度修正——外层由 resolv ilam 循环驱动（不再有内层 while）。 + // 重置累积量 for id in 0..nd { state.heat[id] = 0.0; state.heab[id] = 0.0; @@ -426,6 +462,21 @@ pub fn lucy( state.deltat[id] = state.dt1[id] + state.dt2[id]; } + // [TLUSTY_LUCY_DIAG] 诊断 Lucy 温度修正各项量级, 定位发散根因 + if ilucy == 1 && std::env::var("TLUSTY_LUCY_DIAG").is_ok() { + let tef4 = SIG4P * model.teff.powi(4); + eprintln!("LUCY_DIAG ilam_loop={} TEF4={:.4e} EDDH={:.4e}", ilucy, tef4, state.eddh); + for &idd in &[0usize, 1, 10, 35, 49, 55, 69] { + if idd < nd { + eprintln!(" id={:>2} T={:.0} heat={:+.3e} delh={:+.3e} toth={:.3e} absp={:.3e} absz={:+.3e} eddf={:.3e} dt1={:+.3e} dt2={:+.3e} deltat={:+.3e} (dT/T={:+.3e})", + idd, model.temp[idd], state.heat[idd], state.delh[idd], state.toth[idd], + state.absp[idd], state.absz[idd], state.eddf[idd], + state.dt1[idd], state.dt2[idd], state.deltat[idd], + state.deltat[idd] / model.temp[idd]); + } + } + } + // 应用温度修正 let mut new_temp = model.temp.to_vec(); let mut new_elec = model.elec.to_vec(); @@ -435,67 +486,83 @@ pub fn lucy( for id in 0..nd { new_temp[id] += state.deltat[id]; - state.tem0[id] = new_temp[id]; + ng.tem0[id] = new_temp[id]; let aold = new_dens[id] / model.wmm[id] + new_elec[id]; state.xe[id] = UN - new_elec[id] / aold; } - // 加速方案 - if ilucy >= config.iaclt && ilucy >= iacc0t { + // 加速方案 (Ng) — 对应 Fortran LUCY(tlusty208.f:36132-36202)。 + // ilucy 为 resolv 外层 lambda 迭代号 (1-based)，等价 Fortran 内层 ILUCY； + // iaclt/iacc0t/lac2t/TEM0..3 均取自跨调用持久的 ng。 + // + // 关键控制流（逐字匹配 Fortran）：收集阶段 (!lac2t) 在把 tem0 存入 tem3/tem2/tem1 + // 后 **不** 跳过加速——它 fall through 到下方加速计算；仅 `lac2t && ipng!=0` 的 + // 滚动分支以 "GO TO 20" 跳过本次加速。若用 if/else if/else 互斥链，lac2t 将永远 + // 无法从 false 翻转（首次加速发生在 ilucy==IACLT 的收集 fall-through）。 + let iaclt = ng.iaclt; + let iacc0t = ng.iacc0t; + + // Fortran 36134: if(itlucy.lt.IACLT .or. ilucy.lt.iacc0t) go to 20 + if config.itlucy >= iaclt && ilucy >= iacc0t { + // Fortran 36135-36136: ipng let ipng = if config.iacldt > 0 { - (ilucy - config.iaclt) % config.iacldt + (ilucy - iaclt) % config.iacldt } else { - 0 + 1 }; - if !lac2t { + // Fortran 36137-36166: 历史收集 / 滚动。返回是否跳过本次加速。 + let skip_accel = if !ng.lac2t { + // 收集阶段：在 iacc0t / iacc0t+1 / iacc0t+2 把 tem0 存入 tem3/tem2/tem1。 let ipt = ilucy % 3; - let _ipt0 = config.iaclt % 3; - let ipt1 = (config.iaclt + 1) % 3; - let ipt2 = (config.iaclt + 2) % 3; - + let ipt1 = (iaclt + 1) % 3; + let ipt2 = (iaclt + 2) % 3; if ilucy == iacc0t { for id in 0..nd { - state.tem3[id] = state.tem0[id]; + ng.tem3[id] = ng.tem0[id]; } } else if ipt == ipt1 { for id in 0..nd { - state.tem2[id] = state.tem0[id]; + ng.tem2[id] = ng.tem0[id]; } } else if ipt == ipt2 { for id in 0..nd { - state.tem1[id] = state.tem0[id]; + ng.tem1[id] = ng.tem0[id]; } } + false // 收集后 fall through 到加速（Fortran 无 GO TO 20） } else if ipng != 0 { - // 滚动温度历史 + // 滚动历史 (Fortran 36155-36165): tem3←tem2←tem1←tem0, GO TO 20 for id in 0..nd { - state.tem3[id] = state.tem2[id]; - state.tem2[id] = state.tem1[id]; - state.tem1[id] = state.tem0[id]; + ng.tem3[id] = ng.tem2[id]; + ng.tem2[id] = ng.tem1[id]; + ng.tem1[id] = ng.tem0[id]; } + true // 跳过本次加速 } else { - // 应用加速度 - if ilucy >= config.iaclt { - let (a1, b1, b2, c1, c2) = compute_acceleration(&state, nd); + false // lac2t && ipng==0: fall through 到加速 + }; - let ab = b2 * a1 - b1 * b1; - if ab != 0.0 { - let a0 = (b2 * c1 - b1 * c2) / ab; - let b0 = (a1 * c2 - b1 * c1) / ab; - - for id in 0..nd { - state.tem0[id] = (1.0 - a0 - b0) * state.tem0[id] - + a0 * state.tem1[id] - + b0 * state.tem2[id]; - new_temp[id] = state.tem0[id]; - } - lac2t = true; - } else { - // 对应 Fortran: WRITE(6,601) ILUCY,AB - // FORMAT(/,' **** ACCELT, ITER=',I4,' AB = ',F7.3,/) - eprintln!("\n **** ACCELT, ITER={:4} AB = {:7.3}\n", ilucy, ab); + // Fortran 36168: IF(ILUCY.LT.IACLT) go to 20 + if !skip_accel && ilucy >= iaclt { + let (a1, b1, b2, c1, c2) = compute_acceleration(ng, nd); + let ab = b2 * a1 - b1 * b1; + if ab != 0.0 { + // Fortran 36194-36201 + let a0 = (b2 * c1 - b1 * c2) / ab; + let b0 = (a1 * c2 - b1 * c1) / ab; + for id in 0..nd { + ng.tem0[id] = + (1.0 - a0 - b0) * ng.tem0[id] + a0 * ng.tem1[id] + b0 * ng.tem2[id]; + new_temp[id] = ng.tem0[id]; } + ng.lac2t = true; + } else { + // Fortran 36188-36192: AB==0 → 推迟加速，IACLT+=IACLDT + // FORMAT(/,' **** ACCELT, ITER=',I4,' AB = ',F7.3,/) + eprintln!("\n **** ACCELT, ITER={:4} AB = {:7.3}\n", ilucy, ab); + ng.iaclt += config.iacldt; + ng.iacc0t = ng.iaclt - 3; } } } @@ -536,40 +603,25 @@ pub fn lucy( new_pgs[id] = (new_dens[id] / model.wmm[id] + new_elec[id]) * BOLK * new_temp[id]; } - ilucy += 1; - - // 单次迭代后返回（完整版本会循环） - return LucyOutput { + // 返回本步结果。ilucy / lac2t 由 resolv 经 ilam / ng 跨调用跟踪，不自增。 + LucyOutput { temp: new_temp, elec: new_elec, dens: new_dens, dens1: new_dens1, pgs: new_pgs, ilucy, - lac2t, + lac2t: ng.lac2t, dhhmx1, - }; - } - - // 不应该到达这里 - LucyOutput { - temp: model.temp.to_vec(), - elec: model.elec.to_vec(), - dens: model.dens.to_vec(), - dens1: model.dens1.to_vec(), - pgs: model.pgs.to_vec(), - ilucy, - lac2t, - dhhmx1: 0.0, - } + } } // ============================================================================ // 辅助函数 // ============================================================================ -/// 计算加速度系数。 -fn compute_acceleration(state: &LucyState, nd: usize) -> (f64, f64, f64, f64, f64) { +/// 计算加速度系数。对应 Fortran LUCY 36170-36186，从持久的 TEM0..3 (ng) 读取。 +fn compute_acceleration(ng: &LucyNgState, nd: usize) -> (f64, f64, f64, f64, f64) { let mut a1 = 0.0; let mut b1 = 0.0; let mut b2 = 0.0; @@ -578,13 +630,13 @@ fn compute_acceleration(state: &LucyState, nd: usize) -> (f64, f64, f64, f64, f6 for id in 0..nd { let mut wt = 0.0; - if state.tem0[id] != 0.0 { - wt = 1.0 / state.tem0[id].abs(); + if ng.tem0[id] != 0.0 { + wt = 1.0 / ng.tem0[id].abs(); } - let d0 = state.tem0[id] - state.tem1[id]; - let d1 = d0 - state.tem1[id] + state.tem2[id]; - let d2 = d0 - state.tem2[id] + state.tem3[id]; + let d0 = ng.tem0[id] - ng.tem1[id]; + let d1 = d0 - ng.tem1[id] + ng.tem2[id]; + let d2 = d0 - ng.tem2[id] + ng.tem3[id]; a1 += wt * d1 * d1; b1 += wt * d1 * d2; @@ -722,27 +774,54 @@ mod tests { ..config.clone() }; - let output = lucy(&config_zero, &model, &opacfl_data, &rad1_data); - assert_eq!(output.ilucy, 0); + let mut ng = LucyNgState::new(5, config_zero.iaclt); + let output = lucy(&config_zero, &model, &opacfl_data, &rad1_data, 1, &mut ng); + // itlucy=0 早返回，原样回传调用方传入的 ilucy + assert_eq!(output.ilucy, 1); } #[test] fn test_compute_acceleration() { let nd = 5; - let mut state = LucyState::new(nd); + let mut ng = LucyNgState::new(nd, 7); // 设置一些测试值 for i in 0..nd { - state.tem0[i] = 10000.0 + i as f64; - state.tem1[i] = 9900.0 + i as f64; - state.tem2[i] = 9800.0 + i as f64; - state.tem3[i] = 9700.0 + i as f64; + ng.tem0[i] = 10000.0 + i as f64; + ng.tem1[i] = 9900.0 + i as f64; + ng.tem2[i] = 9800.0 + i as f64; + ng.tem3[i] = 9700.0 + i as f64; } - let (a1, b1, b2, c1, c2) = compute_acceleration(&state, nd); + let (a1, _b1, b2, _c1, _c2) = compute_acceleration(&ng, nd); // 加速度系数应该是正数 assert!(a1 >= 0.0); assert!(b2 >= 0.0); } + + #[test] + fn test_lucy_ng_state_creation_and_roll() { + // LucyNgState 持久层基本语义：新建后 lac2t=false, iacc0t=iaclt-3； + // 滚动历史 tem3←tem2←tem1←tem0 后，tem3 应等于滚动前的 tem2。 + let nd = 4; + let mut ng = LucyNgState::new(nd, 7); + assert!(!ng.lac2t); + assert_eq!(ng.iacc0t, 4); + assert_eq!(ng.tem0.len(), nd); + + for i in 0..nd { + ng.tem0[i] = 100.0 + i as f64; + ng.tem1[i] = 200.0 + i as f64; + ng.tem2[i] = 300.0 + i as f64; + ng.tem3[i] = 400.0 + i as f64; + } + let prev_tem2: Vec = ng.tem2.clone(); + for i in 0..nd { + ng.tem3[i] = ng.tem2[i]; + ng.tem2[i] = ng.tem1[i]; + ng.tem1[i] = ng.tem0[i]; + } + assert_eq!(ng.tem3, prev_tem2, "滚动后 tem3 应等于滚动前 tem2"); + } } diff --git a/src/tlusty/state/model.rs b/src/tlusty/state/model.rs index e1258bc..da41254 100644 --- a/src/tlusty/state/model.rs +++ b/src/tlusty/state/model.rs @@ -2562,6 +2562,7 @@ impl ModelState { phoexp: PhoExp::new(), obfpar: ObfPar::new(), levadd: LevAdd::new(), + opmean: OpMean::new(), ..Default::default() } }