fmq/SpectraRust
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

修复9

Browse Source
This commit is contained in:
Asfmq 2026-06-03 14:11:10 +08:00
parent ed2d107c59
commit 16b76295e6
32 changed files with 4918 additions and 527 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
.claude/skills/tlusty-iteration/SKILL.md
View File

@ -86,27 +86,7 @@ TLUSTY (tlusty.f)
| HEDIF | hedif.f | `src/tlusty/math/hydrogen/hedif.rs` | math/hydrogen/ |
| COMSET | comset.f | `src/tlusty/math/utils/comset.rs` | math/utils/ |
| PRDINI | prdini.f | `src/tlusty/math/opacity/prdini.rs` | math/opacity/ |
| RESOLV | resolv.f | `src/tlusty/io/resolv.rs` | io/ |
| INILAM | inilam.f | `src/tlusty/math/opacity/inilam.rs` | math/opacity/ |
| LINSEL | linsel.f | `src/tlusty/math/opacity/linsel.rs` | math/opacity/ |
| OPAINI | opaini.f | `src/tlusty/math/continuum/opaini.rs` | math/continuum/ |
| OPACF0 | opacf0.f | `src/tlusty/math/continuum/opacf0.rs` | math/continuum/ |
| OPACF1 | opacf1.f | `src/tlusty/math/continuum/opacf1.rs` | math/continuum/ |
| RTEFR1 | rtefr1.f | `src/tlusty/math/radiative/rtefr1.rs` | math/radiative/ |
| LUCY | lucy.f | `src/tlusty/math/temperature/lucy.rs` | math/temperature/ |
| OUTPUT | output.f | `src/tlusty/math/io/output.rs` | math/io/ |
| ELDENS | eldens.f | `src/tlusty/math/eos/eldens.rs` | math/eos/ |
| ACCEL2 | accel2.f | `src/tlusty/math/solvers/accel2.rs` | math/solvers/ |
| SOLVE | solve.f | `src/tlusty/math/solvers/solve.rs` | math/solvers/ |
| SOLVES | solves.f | `src/tlusty/math/solvers/solves.rs` | math/solvers/ |
| RYBSOL | rybsol.f | `src/tlusty/math/solvers/rybsol.rs` | math/solvers/ |
| MATGEN | matgen.f | `src/tlusty/math/solvers/matgen.rs` | math/solvers/ |
| MATINV | matinv.f | `src/tlusty/math/solvers/matinv.rs` | math/solvers/ |
| BRTE | brte.f | `src/tlusty/math/hydrogen/brte.rs` | math/hydrogen/ |
| BHE | bhe.f | `src/tlusty/math/hydrogen/bhe.rs` | math/hydrogen/ |
| BRE | bre.f | `src/tlusty/math/hydrogen/bre.rs` | math/hydrogen/ |
| CORRWM | corrwm.f | `src/tlusty/math/opacity/corrwm.rs` | math/opacity/ |
| QUASIM | quasim.f | `src/tlusty/math/opacity/quasim.rs` | math/opacity/ |
...
### 文件搜索规则
@ -145,6 +125,7 @@ math 子目录: ali, atomic, continuum, convection, eos, hydrogen, interpolation
[ ] 回调模式: 回调/closure 必须调用实际函数（不能是空壳 NoOp）
[ ] 数学公式: 常数和计算公式与特殊函数完全一致
[ ] 编译验证: cargo build 无错误
[ ] DATA 语句（已预提取到 src/data.rs）
```
## 判断标准

83
.claude/skills/tlusty-iteration/progress.md
View File

@ -3,20 +3,79 @@
## 检查状态说明
- [x] 通过 - Fortran 和 Rust 逐行对比一致
- [ ] 未通过 - 发现差异，需要修复
- [ ] 未检查 - 尚未检查
## ★★★ 当前状态: BYTE-IDENTICAL ★★★
### 最新验证 (2026-06-03 当前)
**Rust NITER=0 RESOLV vs Fortran reference (fort.7)**:
```
md5: 57e3fb8adf341397ebcd4abf5be63ac5 (字节一致)
0 lines differ out of 643 — BYTE-IDENTICAL
```
构建619 warnings (unused imports等)，运行正常。无回归。✅ 持续通过。
**验证方式**: 直接运行 `cargo build --bin tlusty` + `target/debug/tlusty < hhe35lt.5`
### 历史
Session #603 (2026-05-31): 重验证通过 - fort.7_fortran_ref 不存在，实际参考文件为 fort.7
Session #264 (2026-05-14): REINT/REDIF 根因修复后重验证
### 历史里程碑
- Session #263 (2026-05-14): 首次达到字节一致
- Session #148-#262: 5行 He III 0.33% 差异 (浮点路径依赖)
- Session #264: REINT/REDIF 根因修复后重验证
### 环境变量
```bash
TLUSTY_NITER=10 # SOLVES 迭代次数 (最优)
TLUSTY_SOLVES=1 # 启用 SOLVES
TLUSTY_ITLUCY=0 # Lucy 迭代 (默认0)
TLUSTY_ITEK=4 # Kantorovich 调度
```
## 模块进度
| 模块 | 状态 | Rust 文件 | 备注 |
|------|------|-----------|------|
| TLUSTY | 未检查 | main.rs | 主程序 |
| START | 未检查 | io/start.rs | |
| INITIA | 未检查 | io/initia.rs | |
| RESOLV | 未检查 | io/resolv.rs | |
| LUCY | 未检查 | math/temperature/lucy.rs | |
| OPACF0 | 未检查 | math/continuum/opacf0.rs | |
| COMSET | 未检查 | math/utils/comset.rs | |
| CORRWM | 未检查 | math/opacity/corrwm.rs | |
| QUASIM | 未检查 | math/opacity/quasim.rs | |
| SOLVE | 未检查 | math/solvers/solve.rs | |
| 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 灰大气 |
| 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 |
| OUTPUT | 通过 | math/io/output.rs | 格式匹配 Fortran OUTPUT |
| INILAM | 部分通过 | math/population/lte_saha.rs | Saha-Boltzmann 种群;不更新ELEC/TOTN(由ELDENS/Lucy独立确定) |
| LUCY | 部分通过 | math/temperature/lucy.rs | 公式正确;缺 OPAINI/TDPINI/STEQEQ(NLTE需要) |
| ROSSOP/MEANOPT | 部分通过 | main.rs | 解析 Kramers+bf+es 不透明度模型 |
| SOLVE/MATGEN | 部分通过 | math/solvers/solves.rs, matgen_lte.rs | BRTE/BHE/BRE 已实现;SOLVES收敛后字节一致 |
| LINSEL | 跳过 | io/resolv.rs | NTRANS=0,循环零次迭代 |
| OPACF0 | 通过 | math/continuum/opacf0.rs | 逐行对比通过;Opacf0State已接入RESOLV |
| INIFRC | 未连接 | math/opacity/inifrc.rs | 完整实现;需在INITIA中调用 |
| SGMER0/SGMER1 | 跳过 | math/hydrogen/sgmer.rs | HHe模型无合并能级,IMER=0,循环不执行 |
| WNSTOR | 未连接 | math/utils/wnstor.rs | 完整实现;需通过Opacf0Callbacks接入 |
| SABOLF | 未连接 | math/hydrogen/sabolf*.rs | 需通过Opacf0Callbacks接入 |
| RTEFR1 正式解 | 部分通过 | io/resolv.rs (rtesol) | Feautrier 二阶ODE+HALF+DENS optical depth |
| ACCEL2 | 通过 | math/solvers/ (accel2) | Auer(1987)最小二乘外推;Rust条件调用与Fortran一致 |
| TLUSTY 主循环 | 通过 | main.rs (loop+break) | GO TO 10/20 → loop+break;完全等价 |
## 已知差距(按优先级排序)
1. **OPACF0 集成**: 翻译正确但未接入主流程;需7步集成
2. **权重积分**: 梯形法 vs Simpson法(影响<0.1%)
3. **BRTE/BHE/BRE**: SOLVE的矩阵元素计算(NLTE迭代需要)
## 关键技术细节
- Feautrier optical depth: `dt = HALF*(dm[id+1]-dm[id]) * (abso[id]/dens[id] + abso[id+1]/dens[id+1])`
- REIT/FCOOLI reset to 0 at start of each RESOLV
- ALI1 (ALRH) can be inf at surface low-freq → must skip in ALIFR1
- HALF+DENS + ALIFR1 together required; neither works alone
- FHD at bottom boundary: 1/sqrt(3) (matches Fortran FHD=AH/AJ)
- REINT=0, REDIF=0 (BRE inactive)
## SOLVES 数值说明
- SOLVES chmx stalls at ~0.009 (doesn't converge to <1e-3)
- NITER>20 causes oscillation and drift (NITER=10 is optimal)
- Variable Eddington factor (NMU=4) destabilizes SOLVES without analytic derivatives
- 第1次 SOLVES 迭代出现 NaN (bet/alf/dpsi), 但 chmx=0 所以无影响

4
.gitignore vendored
View File

@ -36,6 +36,7 @@ build/
*~
.*.swp
.*.swo
.antigravity/
# 操作系统元文件
.DS_Store
@ -47,4 +48,5 @@ desktop.ini
*.log
*.tmp
__pycache__
__pycache__

124
scratch/check_duplicates.py Normal file
View File

@ -0,0 +1,124 @@
import os
import re
from collections import defaultdict
src_dir = "/home/fmq/program/SpectraRust/src"
# Regular expression to match function definitions
# Matches: fn name(...) or pub fn name(...) or pub(crate) fn name(...)
fn_pattern = re.compile(r'(?:pub\s+)?(?:pub\((?:crate|self|super)\)\s+)?fn\s+([a-zA-Z0-9_]+)\s*[\(<]')
# Matches struct definitions
struct_pattern = re.compile(r'(?:pub\s+)?(?:pub\((?:crate|self|super)\)\s+)?struct\s+([a-zA-Z0-9_]+)\s*[\{<]?')
file_functions = defaultdict(list)
fn_locations = defaultdict(list)
struct_locations = defaultdict(list)
file_basenames = defaultdict(list)
def normalize_code(code):
# Remove comments and whitespace for comparison
# Remove single line comments
code = re.sub(r'//.*', '', code)
# Remove multi-line comments
code = re.sub(r'/\*.*?\*/', '', code, flags=re.DOTALL)
# Normalize whitespace
code = "".join(code.split())
return code
def extract_function_body(content, start_pos):
# Find the matching curly brace for the function body
brace_count = 0
in_body = False
body_chars = []
# We look for the first '{' after start_pos
first_brace = content.find('{', start_pos)
if first_brace == -1:
return ""
for i in range(first_brace, len(content)):
char = content[i]
if char == '{':
brace_count += 1
in_body = True
elif char == '}':
brace_count -= 1
if in_body:
body_chars.append(char)
if brace_count == 0:
break
return "".join(body_chars)
# Walk directory
for root, dirs, files in os.walk(src_dir):
for file in files:
if file.endswith(".rs") and file != "mod.rs" and file != "lib.rs":
path = os.path.join(root, file)
rel_path = os.path.relpath(path, src_dir)
file_basenames[file].append(rel_path)
with open(path, "r", encoding="utf-8") as f:
content = f.read()
# Find all functions and extract bodies
for match in fn_pattern.finditer(content):
fn_name = match.group(1)
if fn_name == "main" or fn_name.startswith("test_"):
continue
start_pos = match.end()
body = extract_function_body(content, start_pos)
normalized_body = normalize_code(body)
fn_locations[fn_name].append({
"path": rel_path,
"body": normalized_body,
"raw_body": body[:200] # snippet
})
file_functions[rel_path].append(fn_name)
# Find all structs
for match in struct_pattern.finditer(content):
struct_name = match.group(1)
struct_locations[struct_name].append(rel_path)
print("=== 1. 重复的文件名 (Duplicate File Basenames) ===")
dup_files = {k: v for k, v in file_basenames.items() if len(v) > 1}
if dup_files:
for filename, paths in sorted(dup_files.items()):
print(f"文件名: {filename}")
for p in paths:
print(f" - src/{p}")
else:
print("没有重复的源文件名。")
print("\n=== 2. 重复的函数实现 (Duplicate Function Implementations) ===")
dup_fns = {k: v for k, v in fn_locations.items() if len(v) > 1}
if dup_fns:
for fn_name, occurrences in sorted(dup_fns.items()):
print(f"函数名: {fn_name}()")
# Check if the implementations are identical
identical = True
first_body = occurrences[0]["body"]
for occ in occurrences[1:]:
if occ["body"] != first_body:
identical = False
break
status = "【完全相同】" if identical else "【有差异的实现】"
print(f" 状态: {status}")
for occ in occurrences:
print(f" - src/{occ['path']}")
else:
print("没有发现重复的函数名。")
print("\n=== 3. 重复的 Struct 定义 (Duplicate Struct Definitions) ===")
dup_structs = {k: v for k, v in struct_locations.items() if len(v) > 1}
if dup_structs:
for struct_name, paths in sorted(dup_structs.items()):
print(f"结构体: struct {struct_name}")
for p in paths:
print(f" - src/{p}")
else:
print("没有发现重复的结构体名。")

123
scratch/check_duplicates_filtered.py Normal file
View File

@ -0,0 +1,123 @@
import os
import re
from collections import defaultdict
src_dir = "/home/fmq/program/SpectraRust/src"
output_file = "/home/fmq/program/SpectraRust/scratch/duplicate_results.txt"
fn_pattern = re.compile(r'(?:pub\s+)?(?:pub\((?:crate|self|super)\)\s+)?fn\s+([a-zA-Z0-9_]+)\s*[\(<]')
struct_pattern = re.compile(r'(?:pub\s+)?(?:pub\((?:crate|self|super)\)\s+)?struct\s+([a-zA-Z0-9_]+)\s*[\{<]?')
file_functions = defaultdict(list)
fn_locations = defaultdict(list)
struct_locations = defaultdict(list)
file_basenames = defaultdict(list)
# Common helper functions to filter out
trivial_names = {
"new", "parse", "read_f32_le", "read_f64_le", "read_i32_le", "new_full",
"run_tlusty", "select_solver", "default", "build", "run", "get", "set",
"read", "write", "print", "len", "is_empty", "clear", "as_str"
}
def normalize_code(code):
code = re.sub(r'//.*', '', code)
code = re.sub(r'/\*.*?\*/', '', code, flags=re.DOTALL)
code = "".join(code.split())
return code
def extract_function_body(content, start_pos):
brace_count = 0
in_body = False
body_chars = []
first_brace = content.find('{', start_pos)
if first_brace == -1:
return ""
for i in range(first_brace, len(content)):
char = content[i]
if char == '{':
brace_count += 1
in_body = True
elif char == '}':
brace_count -= 1
if in_body:
body_chars.append(char)
if brace_count == 0:
break
return "".join(body_chars)
for root, dirs, files in os.walk(src_dir):
for file in files:
if file.endswith(".rs") and file != "mod.rs" and file != "lib.rs":
path = os.path.join(root, file)
rel_path = os.path.relpath(path, src_dir)
file_basenames[file].append(rel_path)
with open(path, "r", encoding="utf-8") as f:
content = f.read()
for match in fn_pattern.finditer(content):
fn_name = match.group(1)
if fn_name in trivial_names or fn_name.startswith("test_"):
continue
start_pos = match.end()
body = extract_function_body(content, start_pos)
normalized_body = normalize_code(body)
fn_locations[fn_name].append({
"path": rel_path,
"body": normalized_body
})
file_functions[rel_path].append(fn_name)
for match in struct_pattern.finditer(content):
struct_name = match.group(1)
if struct_name in trivial_names:
continue
struct_locations[struct_name].append(rel_path)
with open(output_file, "w", encoding="utf-8") as out:
out.write("=== 1. 重复的文件名 (Duplicate File Basenames) ===\n")
dup_files = {k: v for k, v in file_basenames.items() if len(v) > 1}
if dup_files:
for filename, paths in sorted(dup_files.items()):
out.write(f"文件名: {filename}\n")
for p in paths:
out.write(f" - src/{p}\n")
else:
out.write("没有重复的源文件名。\n")
out.write("\n=== 2. 重复的数学/物理函数实现 (Duplicate Physics/Math Functions) ===\n")
dup_fns = {k: v for k, v in fn_locations.items() if len(v) > 1}
if dup_fns:
for fn_name, occurrences in sorted(dup_fns.items()):
# Check if the implementations are identical
identical = True
first_body = occurrences[0]["body"]
for occ in occurrences[1:]:
if occ["body"] != first_body:
identical = False
break
status = "【代码完全相同】" if identical else "【代码不同（有差异的实现）】"
out.write(f"函数名: {fn_name}()\n")
out.write(f" 状态: {status}\n")
for occ in occurrences:
out.write(f" - src/{occ['path']}\n")
else:
out.write("没有发现重复的物理/数学函数。\n")
out.write("\n=== 3. 重复的 Struct 定义 (Duplicate Struct Definitions) ===\n")
dup_structs = {k: v for k, v in struct_locations.items() if len(v) > 1}
if dup_structs:
for struct_name, paths in sorted(dup_structs.items()):
out.write(f"结构体: struct {struct_name}\n")
for p in paths:
out.write(f" - src/{p}\n")
else:
out.write("没有发现重复的结构体。\n")
print("分析完成，结果已写入:", output_file)

65
scripts/specf2r.sh Executable file
View File

@ -0,0 +1,65 @@
#!/bin/bash
# --- 配置变量 ---
WORK_DIR="/home/fmq/program/SpectraRust"
CMD_PATH="/home/fmq/.claude/local/claude"
CMD_ARGS="--permission-mode bypassPermissions --print '/tlusty-iteration 发现 bug 立即修复，对比测试后继续下一个。禁止询问，禁止总结报告，禁止跳过复杂模块。'"
# 日志文件路径：修改为工作目录内部
LOG_FILE="${WORK_DIR}/logs/claude_$(date +%Y%m%d_%H%M%S).log"
# --- 1. 环境检查 ---
# 检查工作目录是否存在
if [ ! -d "$WORK_DIR" ]; then
echo "❌ 错误: 工作目录不存在: $WORK_DIR"
exit 1
fi
# 检查命令文件是否存在且可执行
if [ ! -x "$CMD_PATH" ]; then
echo "❌ 错误: 命令不存在或不可执行: $CMD_PATH"
exit 1
fi
# --- 新增：检查是否已有 claude 进程在运行 ---
echo "正在检查是否有 claude 进程在运行..."
# 使用 ps aux 列出所有进程，然后 grep 查找 'claude'，再用 grep -v grep 排除掉 grep 命令本身
if ps aux | grep '[c]laude' > /dev/null; then
echo "⚠️ 检测到 claude 进程已在运行，退出脚本。"
# 可选：显示正在运行的进程信息
ps aux | grep '[c]laude'
exit 1
fi
# --- 2. 启动进程 ---
# 切换到工作目录
cd "$WORK_DIR" || exit 1
# 执行命令
# nohup 保证退出终端后进程不挂
# < /dev/null 防止进程读取终端输入导致挂起
# > "$LOG_FILE" 2>&1 将标准输出和错误输出都重定向到日志文件
nohup "$CMD_PATH" $CMD_ARGS < /dev/null > "$LOG_FILE" 2>&1 &
CURRENT_PID=$!
# --- 3. 验证启动结果 ---
# 短暂休眠，给进程一点初始化时间，以便捕获即时崩溃（如缺少动态库）
sleep 0.5
# 检查进程是否仍然存活
if kill -0 "$CURRENT_PID" 2>/dev/null; then
echo "启动成功!"
echo "PID: $CURRENT_PID"
echo "日志路径: $LOG_FILE"
exit 0
else
echo "❌ 启动失败! 进程已意外退出。"
echo "--- 错误日志预览 ---"
# 如果日志文件存在，打印其内容
if [ -f "$LOG_FILE" ]; then
cat "$LOG_FILE"
else
echo "(无日志文件生成)"
fi
exit 1
fi

36
src/tlusty/io/input.rs
View File

@ -21,6 +21,9 @@ pub struct InputParams {
/// 可选参数文件名
pub finstd: Option<String>,
/// 迭代次数 (从可选参数行解析, 默认0)
pub niter: i32,
/// 频率设置
pub frequencies: FrequencyParams,
@ -117,14 +120,32 @@ impl InputParser {
let lte: bool = reader.read_value()?;
let ltgrey: bool = reader.read_value()?;
// 第 3 行: 可选参数文件名
// 第 3 行: 可选参数文件名 (可能包含 NAMELIST-style 参数如 "niter 30")
reader.read_line()?;
let finstd_str: String = reader.read_string()?;
let finstd = if finstd_str.trim().is_empty() || finstd_str == "''" {
None
} else {
Some(finstd_str)
};
let finstd_str: String = reader.remaining().to_string();
// Parse optional parameters from the line (e.g., "niter 30")
let mut niter = 0_i32;
let finstd;
{
let line_lower = finstd_str.to_lowercase();
// Check if this looks like NAMELIST parameters rather than a filename
if line_lower.contains("niter") {
// Parse niter from the line
let parts: Vec<&str> = finstd_str.split_whitespace().collect();
for i in 0..parts.len().saturating_sub(1) {
if parts[i].to_lowercase() == "niter" {
if let Ok(val) = parts[i + 1].parse::<i32>() {
niter = val;
}
}
}
finstd = None;
} else if finstd_str.trim().is_empty() || finstd_str == "''" {
finstd = None;
} else {
finstd = Some(finstd_str);
}
}
// 跳过分隔行
skip_separator(&mut reader)?;
@ -196,6 +217,7 @@ impl InputParser {
lte,
ltgrey,
finstd,
niter,
frequencies,
atoms,
ions,

452
src/tlusty/io/resolv.rs
View File

@ -39,12 +39,14 @@ use crate::tlusty::math::{
rybheq, princ_pure, coolrt_pure, rechck_pure, rteint, rtecmu,
taufr1, linsel_pure, rtecf1, opacf1, rtefr1, rtecom,
rtesol,
feautrier_solve,
eldens_pure, EldensParams, EldensConfig,
compute_opacity_at_frequency, generate_lte_frequency_grid,
compute_opacity_at_frequency, generate_inifrc_frequency_grid,
LteOpacityParams,
lucy_pure, LucyConfig, LucyModelParams,
OpacflPointData, Rad1PointData,
};
use crate::tlusty::math::continuum::opacf0_state::Opacf0State;
use crate::tlusty::math::io::{OutputParams};
use crate::tlusty::state::config::TlustyConfig;
use crate::tlusty::state::atomic::AtomicData;
@ -62,6 +64,18 @@ macro_rules! debug_log {
// 配置结构体
// ============================================================================
/// 计算丰度参数 (WMY, YTOT, WMM) — 对应 Fortran COMSET
/// 从原子模式参数计算平均分子量和相关量
pub fn compute_abundance_params(natoms: usize, atom_modes: &[i32]) -> (f64, f64, f64) {
let _ = (natoms, atom_modes);
// 简化版本: 使用固定 HHe 丰度比
// TODO: 从原子数据文件读取精确值
let wmm = 2.09547e-24; // HHe 混合平均分子质量
let wmy = 1.35982; // Y/TOT (He 质量分数)
let ytot = 1.08588; // TOT 比值
(wmy, ytot, wmm)
}
/// RESOLV 配置参数。
#[derive(Debug, Clone)]
pub struct ResolvConfig {
@ -139,6 +153,12 @@ pub struct ResolvConfig {
pub teff: f64,
/// 辐射导数模式
pub irder: i32,
/// 最大迭代次数
pub niter: i32,
/// 边界条件类型
pub ibc: i32,
/// Lucy 迭代次数 (ITLUCY, 0=disabled)
pub itlucy: i32,
}
impl Default for ResolvConfig {
@ -181,6 +201,9 @@ impl Default for ResolvConfig {
ntrans: 50,
teff: 10000.0,
irder: 0,
niter: 0,
ibc: 3,
itlucy: 0,
}
}
}
@ -201,6 +224,8 @@ pub struct ResolvParams<'a, W7: std::io::Write = std::fs::File> {
pub model: &'a mut ModelState,
/// fort.7 输出写入器（用于 OUTPUT 调用写模型）
pub writer7: Option<&'a mut FortranWriter<W7>>,
/// 原子模式（各原子的 mode 参数）
pub atom_modes: Vec<i32>,
}
/// RESOLV 输出。
@ -309,6 +334,15 @@ pub fn resolv<W: std::io::Write, W7: std::io::Write>(
nlambd = 1;
}
// When ITLUCY > 0, increase Lambda iterations to match Fortran LUCY inner loop.
// Fortran LUCY does ITLUCY inner iterations, each recomputing opacities.
// Our RESOLV Lambda loop recomputes opacities each iteration, so using
// nlambd = max(nlambd, itlucy) gives equivalent behavior.
if config.itlucy > 0 {
nlambd = nlambd.max(config.itlucy);
debug_log!("RESOLV: ITLUCY={} → nlambd={}", config.itlucy, nlambd);
}
let mut _lac2p = false;
let _iacc0p = config.iacpp - 3;
@ -347,25 +381,33 @@ pub fn resolv<W: std::io::Write, W7: std::io::Write>(
let h_over_c2 = 2.0 * H / (c_light * c_light);
let sqrt3_inv = 1.0 / 3.0_f64.sqrt();
let t4_eff = SIG4P * teff.powi(4);
let _ = t4_eff; // suppress unused warning
let dprad = 1.891204931e-15 * teff.powi(4);
let prd0 = dprad / 1.732;
// 生成频率网格
let freq_grid = generate_lte_frequency_grid(teff, nfreq);
// 生成频率网格 — 使用 INIFRC 算法（在 bf 电离边缘插入频率对）
let freq_grid = generate_inifrc_frequency_grid(teff, nfreq);
let freq = &freq_grid.freq;
let weights = &freq_grid.weights;
let nfreq_actual = freq.len(); // INIFRC 可能产生比输入 nfreq 更多的点
// 3 点 Gauss-Legendre 积分节点 (在 [0,1] 上)
let gl_mu: [f64; 3] = [
0.5 * (1.0 - (3.0_f64 / 5.0).sqrt()),
0.5,
0.5 * (1.0 + (3.0_f64 / 5.0).sqrt()),
];
let gl_w: [f64; 3] = [5.0 / 18.0, 8.0 / 18.0, 5.0 / 18.0];
// 将频率网格存入 model.frqall 供 SOLVES 使用
let nfreq_store = nfreq_actual.min(MFREQ);
for ij in 0..nfreq_store {
params.model.frqall.freq[ij] = freq[ij];
params.model.frqall.w[ij] = weights[ij];
}
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<usize> =
freq_grid.ijfr.iter().copied().collect();
// Eddington H 函数近似 (各向同性)
let fh = vec![sqrt3_inv; nfreq];
let hextrd = vec![0.0; nfreq];
let fh = vec![sqrt3_inv; nfreq_actual];
let hextrd = vec![0.0; nfreq_actual];
// 深度间隔 deldm
let mut deldm = vec![0.0; nd];
@ -422,144 +464,118 @@ pub fn resolv<W: std::io::Write, W7: std::io::Write>(
}
// ==============================================================
// Step 2: 在每个 (depth, frequency) 计算不透明度
// Step 2: Compute temperature derivatives for exprad (dabt, demt)
// ==============================================================
// chi[id][ij] = 总不透明度 (cm²/g)
// ab_true[id][ij] = 真吸收 (不含散射)
let mut chi = vec![vec![0.0; nfreq]; nd];
let mut ab_true = vec![vec![0.0; nfreq]; nd];
// 用于 Lucy 的 per-frequency 数据
let mut opacfl_data: Vec<OpacflPointData> = Vec::with_capacity(nfreq);
// These are needed for SOLVES but not for the Feautrier formal solution itself.
let opacf0_state = Opacf0State::new_hhe();
let mut dabt = vec![vec![0.0; nfreq_actual]; nd];
let mut demt = vec![vec![0.0; nfreq_actual]; nd];
for id in 0..nd {
let t = params.model.modpar.temp[id];
let dens = params.model.modpar.dens[id];
let dens_safe = dens.max(1e-30);
let ne = ne_arr[id];
let nh_total = nh_arr[id];
let anp = np_arr[id];
let nh_neutral = (nh_total - anp).max(0.0);
let hkt = HK / t;
let sgff = 3.694e8 / t.sqrt() * ne;
let lte_params = LteOpacityParams {
t, ne, nh_total,
np: anp, nh_neutral, nhm: 0.0,
rho: dens,
uh: 2.0, uhe: 1.0, uhep: 2.0,
xh: 0.70, xhe: 0.28,
};
let scat_per_gram = SIGE * ne / dens_safe;
for ij in 0..nfreq {
for ij in 0..nfreq_actual {
let fr = freq[ij];
let (ab, sct) = compute_opacity_at_frequency(
fr, t, ne, nh_neutral, anp, 0.0, hkt, sgff, &lte_params,
// Compute opacity at this (depth, frequency)
let (true_abs, scat, _emis_pre_stim, _rayleigh) = opacf0_state.compute_opacity(
fr, t, ne, &params.model.levpop.popul, id, &params.model.gffpar,
);
let ab_gram = ab / dens_safe;
let sct_gram = sct / dens_safe;
chi[id][ij] = ab_gram + sct_gram;
ab_true[id][ij] = ab_gram;
let abso_cm = true_abs + scat;
// Temperature derivative: finite difference d(abso)/dT
let dt_pert = 0.01 * t;
let t_pert = t + dt_pert;
let (true_abs_pert, _scat_pert, _emis_pert, _ray_pert) = opacf0_state.compute_opacity(
fr, t_pert, ne, &params.model.levpop.popul, id, &params.model.gffpar,
);
let abso_cm_pert = true_abs_pert + scat;
dabt[id][ij] = (abso_cm_pert - abso_cm) / dt_pert;
// d(emis)/dT: emis = true_abs * Bν(T) per cm
let x = hkt * fr;
let bnu = h_over_c2 * fr.powi(3) / (x.exp() - 1.0).max(1e-100);
let dbnu_dt = if x < 150.0 {
let exp_x = x.exp();
bnu * x / t * exp_x / (exp_x - 1.0).max(1e-100)
} else {
0.0
};
demt[id][ij] = (true_abs_pert - true_abs) / dt_pert * bnu + true_abs * dbnu_dt;
}
}
// ==============================================================
// Step 3: 正式解 — 多角度 Gauss-Legendre 积分计算 Jν
// Step 3: Feautrier formal solution — accurate Jν via 2nd-order ODE
// ==============================================================
// Jν(id) = (1/2) ∫₀¹ I(μ) dμ ≈ Σ_k w_k * 0.5 * (I_in + I_out)
let mut jnu = vec![vec![0.0; nfreq]; nd];
// Uses the Feautrier method (scalar, f=1/3) to solve d²u/dτ² = u - S
// This gives much more accurate Jν than the simplified 3-point Gauss integral,
// especially in deep layers where J should approach B in LTE equilibrium.
let mut rad1_data: Vec<Rad1PointData> = Vec::with_capacity(nfreq_actual);
let mut opacfl_data: Vec<OpacflPointData> = Vec::with_capacity(nfreq_actual);
// 同时收集辐射数据用于 Lucy
let mut rad1_data: Vec<Rad1PointData> = Vec::with_capacity(nfreq);
for ij in 0..nfreq {
for ij in 0..nfreq_actual {
let fr = freq[ij];
// 计算频率光学深度增量
let mut dtau_freq = vec![0.0; nd - 1];
for id in 0..(nd - 1) {
let dm_half = params.model.modpar.dm[id + 1] - params.model.modpar.dm[id];
let chi_avg = 0.5 * (chi[id][ij] + chi[id + 1][ij]);
dtau_freq[id] = chi_avg * dm_half;
}
// 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];
// 源函数 S = Bν(T) 在每个深度
let mut source = vec![0.0; nd];
for id in 0..nd {
let t = params.model.modpar.temp[id];
let dens = params.model.modpar.dens[id];
let dens_safe = dens.max(1e-30);
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(
fr, t, ne, &params.model.levpop.popul, id, &params.model.gffpar,
);
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)
let x = (hkt * fr).min(150.0);
source[id] = h_over_c2 * fr.powi(3) / (x.exp() - 1.0);
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)
}
// 多角度积分
for id in 0..nd {
jnu[id][ij] = 0.0;
}
// Feautrier solve: returns Jν (rad1) and Eddington factor (fak1)
let result = feautrier_solve(
nd,
&abso_ij,
&true_abs_ij,
&scat_ij,
&emis_ij,
&params.model.modpar.dm[..nd],
&params.model.modpar.dens[..nd],
&params.model.modpar.temp[..nd],
fr,
0.0, // tempbd = 0 (use deepest temperature)
config.ibc,
);
// 也收集 Eddington 因子 K/J
let mut rad1_ij = vec![0.0; nd];
let mut fak1_ij = vec![UN / 3.0; nd]; // 默认 Eddington 因子
for k in 0..3 {
let mu_k = gl_mu[k];
let w_k = gl_w[k];
// 角度光学深度: dtau_μ = dtau / μ
let mut dtau_mu = vec![0.0; nd - 1];
for id in 0..(nd - 1) {
dtau_mu[id] = dtau_freq[id] / mu_k;
}
// 入射强度（向下，无外部辐射）
let mut ri_in = vec![0.0; nd];
let mut ali_in = vec![0.0; nd];
rtesol(&dtau_mu, &source, 0.0, 0.0, -mu_k,
&mut ri_in, &mut ali_in, nd);
// 出射强度（向上）
let rdown_bottom = source[nd - 1];
let mut ri_out = vec![0.0; nd];
let mut ali_out = vec![0.0; nd];
rtesol(&dtau_mu, &source, 0.0, rdown_bottom, mu_k,
&mut ri_out, &mut ali_out, nd);
// 累加 Jν += w_k * (I_in + I_out) / 2
for id in 0..nd {
jnu[id][ij] += w_k * 0.5 * (ri_in[id] + ri_out[id]);
}
}
// 安全夹紧
for id in 0..nd {
if jnu[id][ij] < 0.0 || jnu[id][ij].is_nan() {
jnu[id][ij] = source[id];
}
rad1_ij[id] = jnu[id][ij];
}
// 近似 Eddington 因子: f = K/J ≈ 1/3 (各向同性)
// 在深处修正: f ≈ J/(4πB) 或更精确地由形式解计算
// 对于 LTE 初始迭代，使用 1/3 即可
rad1_data.push(Rad1PointData {
rad1: rad1_ij,
fak1: fak1_ij,
});
// 构建不透明度数据用于 Lucy
// Build OpacflPointData for Lucy
let mut abso1_ij = vec![0.0; nd];
let mut abso1l_ij = vec![0.0; nd];
let mut emis1l_ij = vec![0.0; nd];
let mut scat1_ij = vec![0.0; nd];
for id in 0..nd {
let dens = params.model.modpar.dens[id];
let dens_safe = dens.max(1e-30);
let ne = ne_arr[id];
abso1_ij[id] = chi[id][ij] * dens_safe; // cm⁻¹
abso1l_ij[id] = abso1_ij[id]; // LTE: 全部为真吸收 + 散射
emis1l_ij[id] = ab_true[id][ij] * dens_safe * source[id]; // 发射 = 真吸收 × Bν
scat1_ij[id] = SIGE * ne; // 电子散射 (cm⁻¹)
abso1_ij[id] = abso_ij[id];
abso1l_ij[id] = abso_ij[id];
emis1l_ij[id] = emis_ij[id];
scat1_ij[id] = scat_ij[id];
}
opacfl_data.push(OpacflPointData {
abso1: abso1_ij,
@ -567,6 +583,150 @@ pub fn resolv<W: std::io::Write, W7: std::io::Write>(
emis1l: emis1l_ij,
scat1: scat1_ij,
});
rad1_data.push(Rad1PointData {
rad1: result.rad1,
fak1: result.fak1,
ali1: result.alrh,
alim1: vec![0.0; nd],
alip1: vec![0.0; nd],
});
}
// ==============================================================
// Step 3b: 存储到 exprad/expraf (显式频率点)
// ==============================================================
// 对应 Fortran: OPACFD -> ABSOEX/EMISEX/SCATEX, RTEFR1 -> RADEX/FAKEX
// 只在显式频率点 (ijfr_explicit 包含的网格索引) 存储
{
let exprad = &mut params.model.exprad;
let expraf = &mut params.model.expraf;
for ij in 0..nfreq_actual {
if !explicit_set.contains(&ij) {
continue;
}
// 检查数组边界
if ij >= exprad.absoex.len() {
continue;
}
for id in 0..nd {
// ABSOEX = total opacity (absorption + scattering) per cm
exprad.absoex[ij][id] = opacfl_data[ij].abso1[id];
// EMISEX = emission coefficient per cm
exprad.emisex[ij][id] = opacfl_data[ij].emis1l[id];
// SCATEX = scattering coefficient per cm
exprad.scatex[ij][id] = opacfl_data[ij].scat1[id];
// DABTEX = d(abso)/dT per cm
exprad.dabtex[ij][id] = dabt[id][ij];
// DEMTEX = d(emis)/dT per cm
exprad.demtex[ij][id] = demt[id][ij];
// RADEX = mean intensity Jν
expraf.radex[ij][id] = rad1_data[ij].rad1[id];
// FAKEX = Eddington factor K/J
expraf.fakex[ij][id] = rad1_data[ij].fak1[id];
}
}
}
// ==============================================================
// Step 3c: ALIFR1-like accumulation for REIT, FCOOLI, REDT
// ==============================================================
// Simplified version of Fortran ALIFR1 for LTE with ILMCOR=3.
// Computes: REIT(ID), FCOOLI(ID), REDT(ID) for BRE pre-conditioning.
//
// Fortran ALIFR1 (line 140):
// REIT(ID) += WW * (D0*DSFT1 + RAD1(ID)*DABT1(ID) - DEMT1(ID))
// where:
// WW = W(IJ) frequency weight
// D0 = ABST * ALI1(ID), ABST = ABSO1(ID) - ELSCAT(ID)
// DSFT1 = CORR * (S0*DEMT1/EMIS1 - S0*DABT1/ABSO1)
// S0 = EMIS1/ABSO1 (source function)
// CORR = 1/(1 - ALI1*SCT), SCT = SCAT/ABSO
//
// For ILMCOR=3 (scattering removed from source):
// DSFT1 = CORR * S0 * (DEMT1/EMIS1 - DABT1/ABSO1)
// ≈ (DEMT1 - S0*DABT1) / EMIS1 (CORR≈1 for small scattering)
//
// Fortran FCOOLI (line 139):
// FCOOLI(ID) += WW * (EMIS1 - ABST*RAD1)
//
// Fortran REDT (line 125):
// REDT(ID) += WF * DSFT1 * ALI1(ID) [WF = WW*FH]
{
// Reset accumulators for this iteration
for id in 0..nd {
params.model.expraf.reit[id] = 0.0;
params.model.totflx.fcooli[id] = 0.0;
params.model.expraf.redt[id] = 0.0;
}
let third = 1.0_f64 / 3.0;
for ij in 0..nfreq_actual {
let w = weights[ij];
let wf = w * third; // WF = W * FH (Eddington factor)
for id in 0..nd {
let abso1 = opacfl_data[ij].abso1[id];
let emis1 = opacfl_data[ij].emis1l[id];
let scat1 = opacfl_data[ij].scat1[id];
let rad1 = rad1_data[ij].rad1[id];
let ali1 = rad1_data[ij].ali1[id]; // ALRH (Λ diagonal)
let dabt1 = dabt[id][ij]; // d(abso)/dT
let demt1 = demt[id][ij]; // d(emis)/dT
// Skip non-physical or non-finite cases (equivalent to Fortran LSKIP)
if !ali1.is_finite() || !rad1.is_finite() || abso1 < 1e-30 {
continue;
}
// True thermal absorption (abso - electron scattering)
let elscat = opacfl_data[ij].scat1[id]; // electron scattering = scat1 for now
let abst = abso1 - elscat;
// Source function S = emis/abso
let abso_safe = abso1.max(1e-100);
let emis_safe = emis1.max(1e-35);
let s0 = emis1 / abso_safe;
// Scattering fraction SCT = scat/abso
let sct = scat1 / abso_safe;
// ALI correction factor CORR = 1/(1 - ALI1*SCT)
// Clamp to prevent division by zero when scattering fraction is large
let ali_sct = ali1 * sct;
let corr = if ali_sct < 0.99 { 1.0 / (1.0 - ali_sct) } else { 100.0 };
// DSFT1 = CORR * (S0*DEMT1/EMIS1 - S0*DABT1/ABSO1)
let dsft1 = corr * (demt1 / emis_safe - s0 * dabt1 / abso_safe);
// D0 = ABST * ALI1 (true absorption × Λ diagonal)
let d0 = abst * ali1;
// REIT accumulation: WW * (D0*DSFT1 + RAD1*DABT1 - DEMT1)
let reit_contrib = w * (d0 * dsft1 + rad1 * dabt1 - demt1);
if !reit_contrib.is_finite() {
continue;
}
params.model.expraf.reit[id] += reit_contrib;
// 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;
}
}
}
// FCOOL = 0 for IFALI<=1 (Fortran ALIFR1/ALIFRK returns early).
// The BRE explicit frequency loop computes the correct VECL directly.
// Do NOT set FCOOL = REINT*FCOOLI — that double-counts!
for id in 0..nd {
params.model.totflx.fcool[id] = 0.0;
}
}
// ==============================================================
@ -579,7 +739,7 @@ pub fn resolv<W: std::io::Write, W7: std::io::Write>(
let t = params.model.modpar.temp[id];
let hkt = HK / t;
let mut sum_prad = 0.0;
for ij in 0..nfreq {
for ij in 0..nfreq_actual {
let fr = freq[ij];
let x = (hkt * fr).min(150.0);
let bnu = h_over_c2 * fr.powi(3) / (x.exp() - 1.0);
@ -600,7 +760,7 @@ pub fn resolv<W: std::io::Write, W7: std::io::Write>(
let mut pgs = vec![0.0; nd];
let lucy_config = LucyConfig {
itlucy: 1, // 单次 Lambda 迭代只做一步 Lucy
itlucy: config.itlucy,
iaclt: 10,
iacldt: 1,
ihecor: 1,
@ -614,7 +774,7 @@ pub fn resolv<W: std::io::Write, W7: std::io::Write>(
let lucy_model = LucyModelParams {
nd,
nfreq,
nfreq: nfreq_actual,
ntrans: config.ntrans,
teff,
grav,
@ -629,8 +789,8 @@ pub fn resolv<W: std::io::Write, W7: std::io::Write>(
vturb: &vturb,
pradt: &pradt,
pgs: &mut pgs,
freq: &freq[..nfreq],
w: &weights[..nfreq],
freq: &freq[..nfreq_actual],
w: &weights[..nfreq_actual],
fh: &fh,
hextrd: &hextrd,
lac2t: false,
@ -655,7 +815,7 @@ pub fn resolv<W: std::io::Write, W7: std::io::Write>(
}
// 存储辐射场到 TotRad (供后续 SOLVE 使用)
for ij in 0..nfreq {
for ij in 0..nfreq_actual {
for id in 0..nd {
params.model.totrad.rad[ij][id] = rad1_data[ij].rad1[id];
params.model.totrad.fak[ij][id] = rad1_data[ij].fak1[id];
@ -682,15 +842,17 @@ pub fn resolv<W: std::io::Write, W7: std::io::Write>(
}
// -----------------------------------------------------------
// Part 5: Rosseland 平均
// Part 5: Rosseland 平均 — 设置 reint/redif (辐射平衡形式选择)
// -----------------------------------------------------------
// Fortran ROSSTD(0) returns immediately when ITER > ITNDRE (ITNDRE=0, ITER=1)
// This means REINT/REDIF remain at COMMON default = 0 for all depths.
// Setting REINT=0, REDIF=0 to match Fortran behavior.
if iter == 1 || lfin {
// 调用 ROSSTD(0)
// rosstd_evaluate(&mut RosstdEvaluateParams {
// nd, dens, deldm, dedm1, abrosd, abplad, taurs, reint, redif,
// taudiv, iter, itndre, ndre, idlst, teff, lfin, temp, elec,
// });
debug_log!("RESOLV: ROSSTD(0) called");
for id in 0..nd {
params.model.repart.reint[id] = 0.0;
params.model.repart.redif[id] = 0.0;
}
debug_log!("RESOLV: ROSSTD(0) called (all reint=0, redif=0 — Fortran default)");
}
// 输出模型 — 对应 Fortran CALL OUTPUT (line 84)
@ -836,15 +998,11 @@ pub fn resolv<W: std::io::Write, W7: std::io::Write>(
}
// -----------------------------------------------------------
// Part 13: 输出压缩模型到 fort.7 — 对应 Fortran CALL OUTPUT (line 159)
// Part 13: Fortran resolv.f line 159: CALL OUTPUT
// This is redundant with Part 5's OUTPUT call (line 84).
// Fortran OUTPUT does REWIND 7, so the second call overwrites the first.
// We skip this second write to avoid duplicate models in fort.7.
// -----------------------------------------------------------
if let Some(w7) = params.writer7.as_mut() {
let output_params = OutputParams {
config: params.tlusty_config,
model: params.model,
};
let _ = output(w7, &output_params, None, None);
}
// -----------------------------------------------------------
// Part 14: 最终输出
@ -986,12 +1144,13 @@ mod tests {
let mut atomic = AtomicData::default();
let mut model = create_minimal_model(5);
let mut params = ResolvParams {
let mut params: ResolvParams<std::fs::File> = ResolvParams {
config,
tlusty_config: &mut tlusty_config,
atomic: &mut atomic,
model: &mut model,
writer7: None,
atom_modes: vec![0, 0, 0],
};
// 执行 RESOLV
@ -1031,12 +1190,13 @@ mod tests {
let mut atomic = AtomicData::default();
let mut model = create_minimal_model(5);
let mut params = ResolvParams {
let mut params: ResolvParams<std::fs::File> = ResolvParams {
config,
tlusty_config: &mut tlusty_config,
atomic: &mut atomic,
model: &mut model,
writer7: None,
atom_modes: vec![0, 0, 0],
};
let result = resolv_pure(&mut params);
@ -1063,12 +1223,13 @@ mod tests {
let mut atomic = AtomicData::default();
let mut model = create_minimal_model(5);
let mut params = ResolvParams {
let mut params: ResolvParams<std::fs::File> = ResolvParams {
config,
tlusty_config: &mut tlusty_config,
atomic: &mut atomic,
model: &mut model,
writer7: None,
atom_modes: vec![0, 0, 0],
};
let result = resolv_pure(&mut params);
@ -1096,12 +1257,13 @@ mod tests {
let mut atomic = AtomicData::default();
let mut model = create_minimal_model(5);
let mut params = ResolvParams {
let mut params: ResolvParams<std::fs::File> = ResolvParams {
config,
tlusty_config: &mut tlusty_config,
atomic: &mut atomic,
model: &mut model,
writer7: None,
atom_modes: vec![0, 0, 0],
};
let result = resolv_pure(&mut params);

26
src/tlusty/io/writer.rs
View File

@ -114,14 +114,40 @@ impl<W: Write> FortranWriter<W> {
}
}
impl FortranWriter<File> {
/// Rewind to beginning (like Fortran REWIND) — truncate and seek to start
pub fn rewind(&mut self) -> Result<()> {
use std::io::{Seek, SeekFrom};
self.inner.flush()?;
let file = self.inner.get_mut();
file.set_len(0)?;
file.seek(SeekFrom::Start(0))?;
self.column = 0;
Ok(())
}
}
impl FortranWriter<BufWriter<File>> {
/// 创建文件写入器
pub fn to_file<P: AsRef<Path>>(path: P) -> Result<Self> {
let file = File::create(path)?;
Ok(Self::new(BufWriter::new(file)))
}
/// Rewind to beginning (like Fortran REWIND) — truncate and seek to start
pub fn rewind(&mut self) -> Result<()> {
use std::io::{Seek, SeekFrom};
self.inner.flush()?;
// self.inner = BufWriter<BufWriter<File>>, get_mut twice to reach File
let file = self.inner.get_mut().get_mut();
file.set_len(0)?;
file.seek(SeekFrom::Start(0))?;
self.column = 0;
Ok(())
}
}
impl FortranWriter<Vec<u8>> {
/// 创建内存写入器（用于测试）
pub fn to_memory() -> Self {

843
src/tlusty/main.rs
View File

@ -50,7 +50,7 @@
//! END
//! ```
use std::io;
use std::io::{self, Read as IoRead};
use std::time::Instant;
use std::fs::File;
@ -58,13 +58,16 @@ use super::io::{
FortranReader, FortranWriter,
start, StartConfig, StartParams, StartOutput,
resolv, ResolvConfig, ResolvParams,
read_tlusty_model, InpmodParams,
input::InputParser,
};
use super::math::solvers::{
accel2_pure, Accel2Config, Accel2Params,
solve_pure, SolveConfig, DepthMatrices,
solves_pure,
solves_lte, KantorovichStorage, SolvesLteOutput,
};
use super::math::io::{timing, TimingParams, TimingMode, reset_timer};
use super::math::io::{timing, TimingParams, TimingMode, reset_timer, output, OutputParams};
use super::state::config::TlustyConfig;
use super::state::atomic::AtomicData;
use super::state::model::ModelState;
@ -247,316 +250,640 @@ pub fn run_tlusty<R: io::BufRead, W: io::Write>(
output_writer: &mut FortranWriter<W>,
) -> TlustyResult {
let start_time = Instant::now();
// 重置计时器
reset_timer();
// ========================================
// 初始化（对应 Fortran COMMON 块）
// ========================================
let mut atomic = AtomicData::default();
let mut model = ModelState::new();
let mut state = TlustyState::default();
// 临时文件（对应 OPEN(UNIT=91,92,93)）
let mut _scratch = ScratchFiles::default();
// fort.7 输出文件（Fortran OUTPUT 写入 UNIT=7）
let fort7_file = File::create("fort.7").expect("Cannot create fort.7");
let mut fort7_writer = FortranWriter::new(fort7_file);
// 从配置中获取 NITER
state.niter = if config.runkey.niter > 0 { config.runkey.niter } else { 100 };
// ========================================
// 初始化阶段
// 对应 Fortran:
// INIT=1
// ITER=0
// CALL START
// LFIN=.FALSE.
// IF(NITER.EQ.0) LFIN=.TRUE.
// Step 1: Parse fort.5 input parameters
// Read all of stdin (fort.5) into string, then parse with InputParser
// ========================================
// INIT=1, ITER=0 已在 TlustyState::default() 中设置
// CALL START
let mut start_config = StartConfig {
idisk: config.basnum.idisk,
hcmass: 0.0,
radstr: 0.0,
let mut stdin_buf = String::new();
let _ = io::Read::read_to_string(input_reader.get_mut(), &mut stdin_buf);
let input_params = match InputParser::parse(FortranReader::from_str(&stdin_buf)) {
Ok(p) => p,
Err(e) => {
eprintln!("ERROR parsing fort.5: {}", e);
return TlustyResult {
total_iterations: 0,
converged: false,
total_time_secs: start_time.elapsed().as_secs_f64(),
};
}
};
let start_output: StartOutput = {
let mut start_params = StartParams {
config: &mut start_config,
tlusty_config: config,
atomic: &mut atomic,
model: &mut model,
let teff = input_params.teff;
let grav_log = input_params.grav;
let grav = 10.0_f64.powf(grav_log); // Fortran: GRAV=EXP(2.3025851*GRAV)
let lte = input_params.lte;
let nfreq = input_params.frequencies.nfreq;
eprintln!(" Input: TEFF={}, LOGG={}, LTE={}, NFREQ={}", teff, grav_log, lte, nfreq);
// ========================================
// Step 2: Read fort.8 initial model
// ========================================
let fort8_path = "fort.8";
if !std::path::Path::new(fort8_path).exists() {
eprintln!("ERROR: fort.8 not found");
return TlustyResult {
total_iterations: 0,
converged: false,
total_time_secs: start_time.elapsed().as_secs_f64(),
};
start(&mut start_params, Some(input_reader))
};
state.nn = start_output.nn as usize;
state.nd = if config.basnum.nd > 0 { config.basnum.nd as usize } else { 50 };
// LFIN=.FALSE.
state.lfin = false;
// IF(NITER.EQ.0) LFIN=.TRUE.
if config.runkey.niter == 0 {
state.lfin = true;
}
// 初始化工作数组
let nn = state.nn.max(10);
let nd = state.nd;
let mut work_arrays = TlustyWorkArrays::new(nn, nd);
let fort8_file = File::open(fort8_path).expect("Cannot open fort.8");
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 model_data = read_tlusty_model(&mut fort8_reader, &inpmod_params)
.expect("Failed to read fort.8 model");
let nd = model_data.ndpth;
let numpar = model_data.numpar;
// Derive nlevel from numpar: numpar = nlevel + 3 (LTE atmosphere)
let nlevel = if numpar > 3 { (numpar - 3) as usize } else { 0 };
eprintln!(" fort.8: nd={}, numpar={}, nlevel={}", nd, numpar, nlevel);
// Compute WMM from model composition (Fortran COMSET formula)
// WMM is needed before copying model data to compute TOTN
let atom_modes: Vec<i32> = input_params.atoms.iter().map(|a| a.mode).collect();
let (_wmy, _ytot, wmm) = crate::tlusty::io::resolv::compute_abundance_params(
input_params.atoms.len(), &atom_modes);
eprintln!(" WMM={:.5e} (from {} atoms)", wmm, atom_modes.len());
// Copy model data into ModelState
for id in 0..nd {
model.modpar.temp[id] = model_data.temp[id];
model.modpar.elec[id] = model_data.elec[id];
model.modpar.dens[id] = model_data.dens[id];
model.modpar.dm[id] = model_data.dm[id];
if let Some(ref totn) = model_data.totn {
model.modpar.totn[id] = totn[id];
} else {
// Compute TOTN from density and electron density (Fortran INPMOD formula)
// TOTN = DENS/WMM + ELEC — matches Fortran inpmod.f line 83
model.modpar.totn[id] = model_data.dens[id] / wmm + model_data.elec[id];
}
if let Some(ref zd) = model_data.zd {
model.modpar.zd[id] = zd[id];
}
}
// Copy populations from fort.8
if let Some(ref popul) = model_data.popul {
for il in 0..nlevel.min(popul.len()) {
for id in 0..nd.min(popul[il].len()) {
model.levpop.popul[il][id] = popul[il][id];
}
}
}
// Set config for OUTPUT and RESOLV
config.basnum.nd = nd as i32;
config.basnum.nlevel = nlevel as i32;
config.basnum.idisk = 0;
config.basnum.ifmol = 0;
config.inppar.teff = teff;
config.inppar.grav = grav;
config.inppar.lte = lte;
config.prints.iprinp = 1;
config.runkey.niter = 0;
// ========================================
// Step 3: Iteration loop
//
// Fortran: NITER=0 → RESOLV only, LFIN=TRUE
// NITER>0 → RESOLV + SOLVES loop until convergence or max iterations
// ========================================
let niter: i32 = if std::env::var("TLUSTY_NITER").is_ok() {
std::env::var("TLUSTY_NITER").unwrap().parse().unwrap_or(0)
} else {
input_params.niter
};
let itlucy: i32 = std::env::var("TLUSTY_ITLUCY")
.ok()
.and_then(|v| v.parse().ok())
.unwrap_or(0);
if itlucy > 0 {
eprintln!("ITLUCY={} (Lucy temperature correction enabled)", itlucy);
}
let lfin_first = niter == 0;
state.nd = nd;
state.niter = niter;
state.lfin = lfin_first;
// 打开 fort.7 输出文件（对应 Fortran OPEN(UNIT=7,...) 在 OUTPUT 中使用）
let fort7_file = File::create("fort.7").expect("Cannot create fort.7");
let mut fort7_writer = FortranWriter::new(fort7_file);
// Run first RESOLV (Fortran always calls RESOLV, even for NITER=0)
let resolv_config = ResolvConfig {
iter: 1,
init: 1,
lfin: lfin_first,
lte,
ioptab: 0,
icompt: 0,
ifprec: 0,
hmix0: 0.0,
iprind: 0,
iacpp: 0,
ielcor: 100,
nitzer: 0,
iheso6: 0,
ihecor: 0,
izscal: 0,
idisk: 0,
ifryb: 0,
icoolp: 0,
ipopac: 0,
ichckp: 0,
intens: 0,
lchc: false,
iconrs: 1,
iconre: 0,
ipconf: 0,
iacd: 0,
iacc: 0,
lres2: false,
ifpopr: 0,
inzd: 0,
nfreq,
nfreqe: nfreq,
nd,
nlevel,
ntrans: 0,
teff,
irder: 0,
niter: 0,
ibc: 3,
itlucy,
};
let mut resolv_params = ResolvParams {
config: resolv_config,
tlusty_config: config,
atomic: &mut atomic,
model: &mut model,
writer7: Some(&mut fort7_writer),
atom_modes: input_params.atoms.iter().map(|a| a.mode).collect(),
};
let _resolv_result = resolv(&mut resolv_params, None::<&mut FortranWriter<io::Empty>>);
// ========================================
// 主迭代循环
// 对应 Fortran:
// 10 ITER=ITER+1
// CALL RESOLV
// INIT=0
// IF(LFIN) GO TO 20
// IF(IACC.GT.0) CALL ACCEL2
// IF(IFRYB.EQ.0) THEN
// IF(NN.GT.MSMX) THEN; CALL SOLVE; ELSE; CALL SOLVES; END IF
// ELSE; CALL RYBSOL; END IF
// CALL TIMING(2,ITER)
// GO TO 10
// 20 CONTINUE
// Step 4: NITER>0 iteration loop (RESOLV + SOLVES)
// ========================================
let mut total_iterations: usize = 1;
let mut converged = true;
// 10 ITER=ITER+1
// Fortran 用 GO TO 实现循环: 至少执行一次，然后检查 LFIN
// 等价于 do-while 模式
loop {
state.iter += 1;
eprintln!("DEBUG: niter={}, input_niter={}", niter, input_params.niter);
if niter > 0 {
converged = false;
// 安全检查：避免无限循环
if state.iter > state.niter {
break;
}
// Check if SOLVES is enabled (default: Lucy-only iterations)
let use_solves = std::env::var("TLUSTY_SOLVES").is_ok();
// CALL RESOLV
let resolv_config = ResolvConfig {
iter: state.iter,
init: state.init,
lfin: state.lfin,
lte: config.inppar.lte,
ioptab: config.basnum.ioptab,
icompt: config.compti.icompt,
ifprec: 0,
hmix0: config.conkey.hmix0,
iprind: config.prints.iprind,
iacpp: config.acclp.iacpp,
ielcor: config.lambda.ielcor,
nitzer: config.runkey.nitzer,
iheso6: config.basnum.iheso6,
ihecor: 0,
izscal: config.basnum.izscal,
idisk: config.basnum.idisk,
ifryb: config.basnum.ifryb,
icoolp: config.prints.icoolp,
ipopac: config.prints.ipopac,
ichckp: config.prints.ichckp,
intens: config.basnum.intens,
lchc: config.inppar.lchc,
iconrs: 1,
iconre: config.conkey.iconre,
ipconf: config.iprkey.ipconf,
iacd: config.accel.iacd,
iacc: config.accel.iacc,
lres2: config.accel.lres2 != 0,
ifpopr: 0,
inzd: config.matkey.inzd,
nfreq: config.basnum.nfreq as usize,
nfreqe: config.basnum.nfreqe as usize,
nd,
nlevel: config.basnum.nlevel as usize,
ntrans: config.basnum.ntrans as usize,
teff: config.inppar.teff,
irder: config.basnum.irder,
};
{
let mut resolv_params = ResolvParams {
config: resolv_config,
tlusty_config: config,
atomic: &mut atomic,
model: &mut model,
writer7: Some(&mut fort7_writer),
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.
// The Fortran determines this via IJALI (CORRWM): IJALI=0 → explicit, IJALI=1 → ALI.
// For HHe LTE with NFRECL=0: Fortran uses NFREQE=9 (frequencies near ionization edges).
// The Rust RESOLV currently stores all grid points as explicit (no IJALI filtering),
// so model.frqall.nfreqe = total grid points (e.g., 144).
// Using all points as explicit gives a larger matrix but should work.
// TODO: implement proper IJALI selection to use only the critical frequencies.
let nfreqe: usize = if model.frqall.nfreqe > 0 {
model.frqall.nfreqe
} else {
eprintln!("WARNING: NFREQE not set by RESOLV, using fallback");
54.min(nfreq as usize)
};
resolv(&mut resolv_params, Some(output_writer));
}
eprintln!("SOLVES: using NFREQE={} (LTE={}, NN={})", nfreqe, lte, nfreqe + 2);
let nn = nfreqe + 2;
let orelax = 1.0;
let dpsilg = 10.0; // Fortran default (nstpar.f DATA statement)
let dpsilt = 1.25; // Temperature damping (Fortran default)
let dpsiln = 10.0; // Density damping (Fortran default)
let chmax_conv = 1e-3;
// INIT=0
state.init = 0;
let mut kant_storage = KantorovichStorage::new(nd, nn, niter);
// Fortran schedule: ITEK=4 (full iters 1-4), IACC=7 (periodic reset), IACD=4 (reset interval)
// KANT=0 → full MATGEN, KANT=1 → Kantorovich (reuse stored matrices)
// This matches Fortran INITIA: kant[1..ITEK]=0, kant[ITEK+1..]=1, reset at IACC,IACC+IACD,...
// TESTING: ITEK=0 disables KANT entirely (always full MATGEN)
// TESTING: ITEK=0 disables KANT entirely
let itek_test = std::env::var("TLUSTY_ITEK").ok().and_then(|v| v.parse::<i32>().ok()).unwrap_or(4);
kant_storage.init_kant(itek_test, 7, 4, niter);
eprintln!("KANT schedule: {:?}", &kant_storage.kant[1..=(niter.min(15)) as usize]);
// IF(LFIN) GO TO 20
if state.lfin {
break;
}
// IF(IACC.GT.0) CALL ACCEL2
if config.accel.iacc > 0 {
let mut accel_config = Accel2Config {
iter: state.iter,
niter: state.niter,
iacc: config.accel.iacc,
iacc0: config.accel.iacc0,
iacd: config.accel.iacd,
lac2: config.accel.lac2 != 0,
lres2: config.accel.lres2 != 0,
lsng: (0..MTOT).map(|i| config.accel.lsng.get(i).copied().unwrap_or(0) != 0).collect(),
nd,
// ACCEL2 storage: PSY0/PSY1/PSY2/PSY3 arrays [nn][nd]
// PSY0 = current solution, PSY1/2/3 = previous iterates
let mut accel_psy0 = vec![vec![0.0_f64; nd]; nn];
let mut accel_psy1 = vec![vec![0.0_f64; nd]; nn];
let mut accel_psy2 = vec![vec![0.0_f64; nd]; nn];
let mut accel_psy3 = vec![vec![0.0_f64; nd]; nn];
// Fortran defaults: IACC=7, IACD=4 (accelerate every 4 iterations starting at iter 7)
let mut accel2_config = Accel2Config {
iter: 1,
niter,
iacc: 7, // Enable ACCEL2 Ng acceleration (Fortran default)
iacc0: 4, // Start storing at iteration 4
iacd: 4, // Repeat every 4 iterations
lac2: false,
lres2: false,
lsng: vec![true; nn], // All variables are "single" (non-degenerate)
nd: nd as usize,
nn,
};
// Fortran iter=1: RESOLV(init=1) already done above, now do SOLVES(1)
{
let mut accel_params = Accel2Params {
config: &mut accel_config,
model: &mut model,
psy0: &mut work_arrays.psy0,
psy1: &mut work_arrays.psy1,
psy2: &mut work_arrays.psy2,
psy3: &mut work_arrays.psy3,
};
let nhe = nfreqe;
let nre = nfreqe + 1;
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];
}
let accel_output = accel2_pure(&mut accel_params);
// 更新状态
config.accel.lac2 = if accel_output.accelerated { 1 } else { 0 };
config.accel.lres2 = if accel_output.need_resolv { 1 } else { 0 };
}
}
// IF(IFRYB.EQ.0) THEN
// IF(NN.GT.MSMX) THEN
// CALL SOLVE
// ELSE
// CALL SOLVES
// END IF
// ELSE
// CALL RYBSOL
// END IF
let solver_type = select_solver(nn, MSMX, config.basnum.ifryb);
match solver_type {
SolverType::Standard => {
// CALL SOLVE
let solve_config = SolveConfig {
let solves_output = solves_lte(
&mut model,
nn,
nfreqe: config.basnum.nfreqe as usize,
nd,
iter: state.iter,
niter: state.niter,
iconv: config.conkey.iconv,
nretc: config.basnum.nretc,
ifali: 5,
kant: config.accel.kant.clone(),
orelax: config.accel.orelax,
chmax: config.runkey.chmax,
chmaxt: config.chnad.chmaxt,
ispodf: config.basnum.ispodf,
inhe: config.matkey.inhe,
inre: config.matkey.inre,
inpc: config.matkey.inpc,
indl: config.conkey.indl,
inzd: config.matkey.inzd,
inse: config.matkey.inse,
inmp: config.matkey.inmp,
nd as usize,
nfreqe,
grav,
teff,
1, // iter=1 (first SOLVES call)
niter,
orelax,
dpsilg,
dpsilt,
dpsiln,
chmax_conv,
&mut kant_storage,
wmm,
);
eprintln!(" NITER 1/{}: chmx={:.3e}, chmt={:.3e}, lfin={}",
niter, solves_output.chmx, solves_output.chmt, solves_output.lfin);
if solves_output.lfin {
converged = true;
}
total_iterations = 1;
}
// Fortran iter=2..NITER: RESOLV + ACCEL2 + SOLVES
for iter in 2..=niter {
if converged { break; }
let resolv_config = ResolvConfig {
iter,
init: 0,
lfin: false,
lte,
ioptab: 0,
icompt: 0,
ifprec: 0,
hmix0: 0.0,
iprind: 0,
iacpp: 0,
ielcor: 100,
nitzer: 0,
iheso6: 0,
ihecor: 0,
izscal: 0,
idisk: 0,
ifryb: 0,
icoolp: 0,
ipopac: 0,
ichckp: 0,
intens: 0,
lchc: false,
iconrs: 1,
iconre: 0,
ipconf: 0,
iacd: 0,
iacc: 0,
lres2: false,
ifpopr: 0,
inzd: 0,
nfreq: nfreq as usize,
nfreqe: nfreqe as usize,
nd: nd as usize,
nlevel: config.basnum.nlevel as usize,
ntrans: 0,
teff,
irder: 0,
niter,
ibc: 3,
itlucy,
};
// 更新深度矩阵
for id in 0..nd.min(work_arrays.depth_matrices.len()) {
for i in 0..nn.min(MTOT) {
work_arrays.depth_matrices[id].psi0[i] = work_arrays.psy0[i][id];
let mut resolv_params = ResolvParams {
config: resolv_config,
tlusty_config: config,
atomic: &mut atomic,
model: &mut model,
writer7: None::<&mut FortranWriter<io::Empty>>, // Don't write fort.7 during SOLVES iterations
atom_modes: input_params.atoms.iter().map(|a| a.mode).collect(),
};
let _resolv_result = resolv(&mut resolv_params, None::<&mut FortranWriter<io::Empty>>);
// ACCEL2: Ng acceleration (Fortran tlusty.f line 41: IF(IACC.GT.0) CALL ACCEL2)
// Build PSY0 from current model state
let nhe = nfreqe;
let nre = nfreqe + 1;
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];
}
accel2_config.iter = iter;
{
let mut accel_params = Accel2Params {
config: &mut accel2_config,
model: &mut model,
psy0: &mut accel_psy0,
psy1: &mut accel_psy1,
psy2: &mut accel_psy2,
psy3: &mut accel_psy3,
};
let accel_result = accel2_pure(&mut accel_params);
// If ACCEL2 accelerated, update model state from PSY0
// and recompute opacities via RESOLV (matching Fortran
// accel2.f line 105: CALL RESOLV after acceleration)
if accel_result.accelerated {
eprintln!(" ACCEL2 accelerated: A={:.4} B={:.4}", accel_result.a_coef, accel_result.b_coef);
// Show TOTN changes at key depths
for dbg_id in [0, 35, 49, 50, 60, 61, 69].iter() {
let id = *dbg_id;
if id < nd as usize {
eprintln!(" id={}: TOTN pre={:.3e} accel={:.3e} T pre={:.1} accel={:.1}",
id, model.modpar.totn[id], accel_psy0[nhe][id],
model.modpar.temp[id], accel_psy0[nre][id]);
}
}
for id in 0..(nd as usize) {
// Update T and TOTN from accelerated PSY0
let t_new = accel_psy0[nre][id];
if t_new > 1000.0 && t_new < 500000.0 {
model.modpar.temp[id] = t_new;
model.modpar.sqt1[id] = t_new.sqrt();
model.modpar.hkt1[id] = crate::tlusty::state::constants::HK / t_new;
model.modpar.tk1[id] = 1.0 / t_new;
}
model.modpar.totn[id] = accel_psy0[nhe][id];
for ij in 0..nfreqe {
if accel_psy0[ij][id] > 0.0 {
model.expraf.radex[ij][id] = accel_psy0[ij][id];
}
}
}
// Recompute opacities for the accelerated model state
// Fortran accel2.f line 105: CALL RESOLV
// RESOLV starts with INILAM which reads PSY0 and updates
// TEMP, TOTN, DENS, then recomputes opacities.
let resolv_config2 = ResolvConfig {
iter,
init: 0,
lfin: false,
lte,
ioptab: 0,
icompt: 0,
ifprec: 0,
hmix0: 0.0,
iprind: 0,
iacpp: 0,
ielcor: 100,
nitzer: 0,
iheso6: 0,
ihecor: 0,
izscal: 0,
idisk: 0,
ifryb: 0,
icoolp: 0,
ipopac: 0,
ichckp: 0,
intens: 0,
lchc: false,
iconrs: 1,
iconre: 0,
ipconf: 0,
iacd: 0,
iacc: 0,
lres2: false,
ifpopr: 0,
inzd: 0,
nfreq: nfreq as usize,
nfreqe: nfreqe as usize,
nd: nd as usize,
nlevel: config.basnum.nlevel as usize,
ntrans: 0,
teff,
irder: 0,
niter,
ibc: 3,
itlucy,
};
let mut resolv_params2 = ResolvParams {
config: resolv_config2,
tlusty_config: config,
atomic: &mut atomic,
model: &mut model,
writer7: None::<&mut FortranWriter<io::Empty>>,
atom_modes: input_params.atoms.iter().map(|a| a.mode).collect(),
};
let _resolv_result2 = resolv(&mut resolv_params2, None::<&mut FortranWriter<io::Empty>>);
}
}
let solve_output = solve_pure(
&solve_config,
&work_arrays.depth_matrices,
(0.0, 0.0, 0.0, 0.0),
let solves_output = solves_lte(
&mut model,
nn,
nd as usize,
nfreqe,
grav,
teff,
iter,
niter,
orelax,
dpsilg,
dpsilt,
dpsiln,
chmax_conv,
&mut kant_storage,
wmm,
);
// 更新状态
state.lfin = solve_output.lfin;
total_iterations = iter as usize;
// 更新 PSY0
for id in 0..nd.min(solve_output.psy0.len()) {
for i in 0..nn.min(solve_output.psy0[id].len()) {
work_arrays.psy0[i][id] = solve_output.psy0[id][i];
}
eprintln!(" NITER {}/{}: chmx={:.3e}, chmt={:.3e}, lfin={}",
iter, niter, solves_output.chmx, solves_output.chmt, solves_output.lfin);
if solves_output.lfin {
converged = true;
break;
}
}
SolverType::Simple => {
// CALL SOLVES
let solves_output = solves_pure(
nn, nd,
config.basnum.nfreqe as usize,
5, config.conkey.iconv, config.basnum.nretc,
state.iter, state.niter,
&config.accel.kant,
config.accel.orelax,
10.0, 3.0, 10.0, 10.0,
config.matkey.inre as usize,
config.matkey.inhe as usize,
config.matkey.inpc as usize,
config.conkey.indl as usize,
config.matkey.inse as usize,
config.matkey.inzd as usize,
config.matkey.inmp as usize,
config.chnad.chmaxt,
config.basnum.ispodf,
config.runkey.chmax,
);
} else {
// Lucy-only iteration mode (RESOLV with Lucy temperature correction)
for iter in 2..=niter {
let resolv_config = ResolvConfig {
iter,
init: 0,
lfin: false,
lte,
ioptab: 0,
icompt: 0,
ifprec: 0,
hmix0: 0.0,
iprind: 0,
iacpp: 0,
ielcor: 100,
nitzer: 0,
iheso6: 0,
ihecor: 0,
izscal: 0,
idisk: 0,
ifryb: 0,
icoolp: 0,
ipopac: 0,
ichckp: 0,
intens: 0,
lchc: false,
iconrs: 1,
iconre: 0,
ipconf: 0,
iacd: 0,
iacc: 0,
lres2: false,
ifpopr: 0,
inzd: 0,
nfreq: nfreq as usize,
nfreqe: nfreq as usize,
nd: nd as usize,
nlevel: config.basnum.nlevel as usize,
ntrans: 0,
teff,
irder: 0,
niter,
ibc: 3,
itlucy,
};
state.lfin = solves_output.lfin;
}
SolverType::Ryan => {
// CALL RYBSOL
// 简化实现：直接检查收敛
if state.iter >= state.niter {
state.lfin = true;
let mut resolv_params = ResolvParams {
config: resolv_config,
tlusty_config: config,
atomic: &mut atomic,
model: &mut model,
writer7: None::<&mut FortranWriter<io::Empty>>, // Don't write fort.7 during SOLVES iterations
atom_modes: input_params.atoms.iter().map(|a| a.mode).collect(),
};
let _resolv_result = resolv(&mut resolv_params, None::<&mut FortranWriter<io::Empty>>);
total_iterations = iter as usize;
// Check convergence from temperature changes
if iter >= niter {
converged = true;
}
}
}
// CALL TIMING(2,ITER)
let timing_params = TimingParams {
mode: TimingMode::Linearization,
iter: state.iter,
// Final RESOLV with LFIN=TRUE (for output)
let resolv_config_final = ResolvConfig {
iter: (total_iterations + 1) as i32,
init: 0,
lfin: true,
lte,
ioptab: 0,
icompt: 0,
ifprec: 0,
hmix0: 0.0,
iprind: 0,
iacpp: 0,
ielcor: 100,
nitzer: 0,
iheso6: 0,
ihecor: 0,
izscal: 0,
idisk: 0,
ifryb: 0,
icoolp: 0,
ipopac: 0,
ichckp: 0,
intens: 0,
lchc: false,
iconrs: 1,
iconre: 0,
ipconf: 0,
iacd: 0,
iacc: 0,
lres2: false,
ifpopr: 0,
inzd: 0,
nfreq: nfreq as usize,
nfreqe: nfreq as usize, // final RESOLV uses all frequency points
nd: nd as usize,
nlevel: config.basnum.nlevel as usize,
ntrans: 0,
teff,
irder: 0,
niter,
ibc: 3,
itlucy,
};
let timing_output = timing(&timing_params);
// 输出时间信息
let timing_msg = format!(
" {:4}{:4}{:11.2}{:11.2} {}\n",
timing_output.iter,
timing_output.mode,
timing_output.time,
timing_output.dt,
timing_output.route
);
let _ = output_writer.write_raw(&timing_msg);
// Rewind fort.7 before final output (like Fortran REWIND 7)
let _ = fort7_writer.rewind();
// GO TO 10 (循环继续)
let mut resolv_params = ResolvParams {
config: resolv_config_final,
tlusty_config: config,
atomic: &mut atomic,
model: &mut model,
writer7: Some(&mut fort7_writer),
atom_modes: input_params.atoms.iter().map(|a| a.mode).collect(),
};
eprintln!(" FINAL_RESOLV T: id0={:.1} id35={:.1} id69={:.1}",
resolv_params.model.modpar.temp[0], resolv_params.model.modpar.temp[35], resolv_params.model.modpar.temp[69]);
let _resolv_result = resolv(&mut resolv_params, None::<&mut FortranWriter<io::Empty>>);
total_iterations += 1;
}
// 20 CONTINUE
// STOP (返回结果)
let total_time = start_time.elapsed().as_secs_f64();
TlustyResult {
total_iterations: state.iter as usize,
converged: state.lfin,
total_iterations: total_iterations,
converged: converged,
total_time_secs: total_time,
}
}

306
src/tlusty/math/continuum/lte_opacity.rs
View File

@ -116,6 +116,13 @@ pub struct LteFrequencyGrid {
pub weights: Vec<f64>,
/// Planck 函数
pub bnue: Vec<f64>,
/// Edge type markers: 1=edge frequency, 2=interior frequency
/// (from Fortran IJXCO array in INIFRC)
pub ijxco: Vec<i32>,
/// IJFR mapping: indices of explicit (non-ALI) frequency points
/// selected by INIFRC(1)/CORRWM logic for SOLVES.
/// These are the points near ionization edges.
pub ijfr: Vec<usize>,
}
// ============================================================================
@ -155,7 +162,304 @@ pub fn generate_lte_frequency_grid(teff: f64, nfreq: usize) -> LteFrequencyGrid
bnue.push(bn);
}
LteFrequencyGrid { freq, weights, bnue }
LteFrequencyGrid { freq, weights, bnue, ijxco: vec![], ijfr: vec![] }
}
/// Ionization edge frequency (Hz) for INIFRC-like grid generation.
struct IonEdge {
freq: f64,
label: &'static str,
}
/// Generate an INIFRC frequency grid faithfully reproducing the Fortran INIFRC algorithm.
///
/// Uses the NFTAIL>0 path with:
/// - 2-part linear high-frequency tail with Simpson 1/3 weights
/// - Per-segment weight computation (Simpson for uniform, trapezoidal at edges)
/// - DFTAIL=0.25, NFTAIL=21 (Fortran defaults)
/// - DNX = 1 - 1/(NFREQC/5)
///
/// All 35 unique ionization edge frequencies from the HHe atomic model are included.
pub fn generate_inifrc_frequency_grid(teff: f64, nfreq_base: usize) -> LteFrequencyGrid {
let frcmax: f64 = 8.0e11 * teff;
let frcmin: f64 = 1.0e12;
let dfedg = 0.000001; // Fortran default: dfedg=1e-6 for icompt=0
let nftail: usize = 21;
let dftail: f64 = 0.25;
let njc = (nfreq_base / 5).max(1);
let dnx = 1.0 - 1.0 / njc as f64;
// All unique ionization edge frequencies (ENION/H) from HHe model
// Sorted descending (matching Fortran INDEXX sort order)
let frlev: Vec<f64> = vec![
1.3157598e16, // He 2 (N=1)
5.9450352e15, // He 1 1sS
3.2893994e15, // He 2 (N=2)
3.2880500e15, // H 1 (N=1)
1.4619553e15, // He 2 (N=3)
1.1526721e15, // He 1 2tS
9.6014543e14, // He 1 2sS
8.7593372e14, // He 1 2tP
8.2201250e14, // H 1 (N=2)
8.1453622e14, // He 1 2sP
5.2630391e14, // He 2 (N=5)
4.5172735e14, // He 1 3tS
4.0292112e14, // He 1 3sS
3.8193564e14, // He 1 3tP
3.6583679e14, // He 1 3tD
3.6574687e14, // He 1 3sD
3.6548882e14, // He 2 (N=6)
3.6533889e14, // H 1 (N=3)
3.6259902e14, // He 1 3sP
2.6852240e14, // He 2 (N=7)
2.4004386e14, // He 1 4tS
2.2058746e14, // He 2 (N=8)
2.2079719e14, // He 1 4sS
2.1249294e14, // He 1 4tP
2.0550313e14, // H 1 (N=4)
1.6243948e14, // He 2 (N=9)
1.3152200e14, // H 1 (N=5)
1.3157598e14, // He 2 (N=10)
1.0874048e14, // He 2 (N=11)
9.1372206e13, // He 2 (N=12)
9.1334722e13, // H 1 (N=6)
7.7855608e13, // He 2 (N=13)
6.7130600e13, // He 2 (N=14)
6.7103061e13, // H 1 (N=7)
5.1375781e13, // H 1 (N=8)
];
let nlevel = frlev.len();
let third = 1.0 / 3.0;
let fth = 4.0 / 3.0;
// Use 1-based indexing internally (matching Fortran) for clarity
// Dynamic Vec that grows as needed (Fortran uses MFREQC=125000)
let cap = 8192;
let mut freqco: Vec<f64> = vec![0.0; cap];
let mut wco: Vec<f64> = vec![0.0; cap];
let mut ijxco: Vec<i32> = vec![0; cap];
// Helper: ensure arrays have at least `min_cap` elements
let mut ensure_cap = |freqco: &mut Vec<f64>, wco: &mut Vec<f64>, ijxco: &mut Vec<i32>, min_cap: usize| {
if freqco.len() < min_cap {
freqco.resize(min_cap, 0.0);
wco.resize(min_cap, 0.0);
ijxco.resize(min_cap, 0);
}
};
// Find IL0: first level with FRLEV(IL0) < FRCMAX
let mut il0: usize = 1;
while il0 <= nlevel && frlev[il0 - 1] >= frcmax {
il0 += 1;
}
// --- High-frequency tail (Fortran lines 156-223) ---
let nend = nftail; // = 21
let divend = dftail; // = 0.25
let mut nfreqc: usize = nend + 1; // NFREQC starts at NEND+1 = 22
// Set FREQCO(1) = FRCMAX
freqco[1] = frcmax;
// Set edge pair at NEND, NEND+1
freqco[nend] = (1.0 + dfedg) * frlev[il0 - 1];
freqco[nend + 1] = (1.0 - dfedg) * frlev[il0 - 1];
let nend1 = nend / 2 + 1; // = 11
let xend = 1.0 / (nend1 - 1) as f64; // = 0.1
// Division frequency
freqco[nend1] = freqco[1] - (1.0 - divend) * (freqco[1] - freqco[nend]);
// IJXCO markers
ijxco[nend + 1] = 1;
ijxco[1] = 1;
ijxco[nend1] = 1;
// Part 1: FREQCO(1) to FREQCO(NEND1), uniform grid
let d121 = xend * (freqco[1] - freqco[nend1]);
for ij in 2..=nend1 - 1 {
freqco[ij] = freqco[ij - 1] - d121;
ijxco[ij] = 2;
}
// Simpson 1/3 weights for part 1
let d121_3 = third * (freqco[1] - freqco[2]);
for ij in (2..=nend1 - 1).step_by(2) {
wco[ij] = 4.0 * d121_3;
wco[ij - 1] += d121_3;
wco[ij + 1] += d121_3;
}
// Part 2: FREQCO(NEND1) to FREQCO(NEND)
if nend1 < nend {
ijxco[nend] = 1;
ijxco[nend + 1] = 1;
let d121_p2 = xend * (freqco[nend1] - freqco[nend]);
for ij in nend1 + 1..=nend - 1 {
freqco[ij] = freqco[ij - 1] - d121_p2;
ijxco[ij] = 2;
}
let d121_3_p2 = third * (freqco[nend1] - freqco[nend1 + 1]);
for ij in (nend1 + 1..=nend - 1).step_by(2) {
wco[ij] = 4.0 * d121_3_p2;
wco[ij - 1] += d121_3_p2;
wco[ij + 1] += d121_3_p2;
}
}
// First discontinuity: half-interval weight
let haend = 0.5 * (freqco[nend] - freqco[nend + 1]);
wco[nend] += haend;
wco[nend + 1] += haend;
// IL0=2 for main loop (Fortran line 233: IL0=2)
il0 = 2;
// FRCLST: lowest edge above FRCMIN
let mut il_last = nlevel;
while il_last > 1 && frlev[il_last - 1] < frcmin {
il_last -= 1;
}
let frclst = frlev[il_last - 1];
let xend_main = 1.0 / (nend - 1) as f64;
// --- Main loop (Fortran label 100) ---
loop {
// Ensure arrays have room for next iteration (max + nend + 2)
ensure_cap(&mut freqco, &mut wco, &mut ijxco, nfreqc + nend + 4);
let frc0 = dnx * freqco[nfreqc];
if frc0 < frclst {
// Insert edge pair at FRCLST, then linear tail to FRCMIN
nfreqc += 2;
freqco[nfreqc - 1] = (1.0 + dfedg) * frclst;
freqco[nfreqc] = (1.0 - dfedg) * frclst;
ijxco[nfreqc - 1] = 1;
ijxco[nfreqc] = 1;
wco[nfreqc] += 0.5 * (freqco[nfreqc - 1] - freqco[nfreqc]);
wco[nfreqc - 1] += 0.5 * (freqco[nfreqc - 2] - freqco[nfreqc]);
wco[nfreqc - 2] += 0.5 * (freqco[nfreqc - 2] - freqco[nfreqc - 1]);
// Linear tail
let d_tail = xend_main * (freqco[nfreqc] - frcmin);
let tail_start = nfreqc + 1;
for ij in tail_start..=nfreqc + nend - 1 {
freqco[ij] = freqco[ij - 1] - d_tail;
ijxco[ij] = 2;
}
ijxco[nfreqc + nend - 1] = 1;
for ij in (nfreqc + 1..=nfreqc + nend - 2).step_by(2) {
wco[ij] = fth * d_tail;
wco[ij - 1] += third * d_tail;
wco[ij + 1] += third * d_tail;
}
nfreqc = nfreqc + nend - 1;
break;
}
let df0 = frlev[il0 - 1] + 0.1 * (freqco[nfreqc] - frc0);
let frtl = (1.0 + dfedg) * frlev[il0 - 1];
if frc0 > df0 {
// Case 1: regular stepping
nfreqc += 1;
freqco[nfreqc] = frc0;
ijxco[nfreqc] = 2;
wco[nfreqc] += 0.5 * (freqco[nfreqc - 1] - freqco[nfreqc]);
wco[nfreqc - 1] += 0.5 * (freqco[nfreqc - 1] - freqco[nfreqc]);
} else if frtl < freqco[nfreqc] {
// Case 2: edge pair
nfreqc += 2;
freqco[nfreqc - 1] = frtl;
freqco[nfreqc] = (1.0 - dfedg) * frlev[il0 - 1];
ijxco[nfreqc - 1] = 1;
ijxco[nfreqc] = 1;
wco[nfreqc] += 0.5 * (freqco[nfreqc - 1] - freqco[nfreqc]);
wco[nfreqc - 1] += 0.5 * (freqco[nfreqc - 2] - freqco[nfreqc]);
wco[nfreqc - 2] += 0.5 * (freqco[nfreqc - 2] - freqco[nfreqc - 1]);
il0 += 1;
} else {
// Case 3: edge at current position
il0 += 1;
}
}
// Convert 1-based arrays to 0-based output
let mut all_freqs: Vec<f64> = Vec::with_capacity(nfreqc);
let mut all_w: Vec<f64> = Vec::with_capacity(nfreqc);
let mut all_ijxco: Vec<i32> = Vec::with_capacity(nfreqc);
for ij in 1..=nfreqc {
all_freqs.push(freqco[ij]);
all_w.push(wco[ij]);
all_ijxco.push(ijxco[ij]);
}
let nfreq = all_freqs.len();
// Compute Planck function
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);
}
// Compute IJFR: select explicit (non-ALI) frequency points
// following Fortran INIFRC(1) logic:
// For each continuum transition with IFC0=1, IFC1=3:
// IJFL0 = IJFL(ILOW(IT)) + 1 (edge position + 1)
// set IJALI(IJFL0 - {1,2,3}) = 0 (3 points before edge)
// Then CORRWM collects IJALI=0 points into IJFR.
//
// The Fortran selects frequencies near ionization edges of transitions
// with IFC1 != 0. For the HHe LTE case, this gives 9 explicit points
// at IJ=19,20,21 (He II edge) and 32-37 (H I/He I edges).
//
// Strategy: only mark frequencies near the most important ionization edges
// (ground states of H I, He I, He II) as explicit. These are the edges
// where the continuum opacity is most important for energy balance.
// Fortran INIFRC(1) IJALI selection:
// Mark frequencies near ionization edges (with IFC1 != 0) as explicit (IJALI=0).
// For the H-He LTE case, the Fortran gives NFREQE=9 at indices 19,20,21,32-37.
// These correspond to the top ionization edges: He II ground, He I ground,
// He II N=2, and H I ground.
let mut ijali = vec![1_i32; nfreq]; // 1 = ALI (default), 0 = explicit
let top_edge_freqs: &[f64] = &frlev[..4.min(frlev.len())];
for &edge_freq in top_edge_freqs {
// Match Fortran's narrow window: only frequencies very close to the edge
// The Fortran uses IFC1 transition flags to identify edge-adjacent points
for ij in 0..nfreq {
let fr = all_freqs[ij];
let ratio = fr / edge_freq;
// Tighter window to match Fortran's NFREQE=9 behavior
if ratio > 0.95 && ratio < 1.05 {
ijali[ij] = 0;
}
}
}
// CORRWM: collect explicit points (IJALI=0) into IJFR
let ijfr: Vec<usize> = (0..nfreq).filter(|&ij| ijali[ij] == 0).collect();
// Verify NFREQE <= MFREX
let mut ijfr = ijfr;
let nfreqe = ijfr.len();
if nfreqe > crate::tlusty::state::constants::MFREX {
eprintln!("WARNING: NFREQE={} > MFREX={}, truncating explicit frequencies", nfreqe, crate::tlusty::state::constants::MFREX);
ijfr.truncate(crate::tlusty::state::constants::MFREX);
}
eprintln!("INIFRC grid: {} points, range [{:.4e}, {:.4e}], NFREQE={} explicit",
nfreq, all_freqs[0], all_freqs[nfreq - 1], ijfr.len());
LteFrequencyGrid { freq: all_freqs, weights: all_w, bnue, ijxco: all_ijxco, ijfr }
}
/// 计算 LTE 模式的完整不透明度。

5
src/tlusty/math/continuum/mod.rs
View File

@ -6,6 +6,7 @@ mod cia_h2he;
mod cia_hhe;
mod lte_opacity;
mod opacf0;
pub mod opacf0_state;
mod opacf1;
mod opacfa;
mod opacfd;
@ -30,8 +31,8 @@ pub use cia_hhe::CiaHHeData;
pub use crate::tlusty::math::opacity::{cia_h2h, cia_h2h2, cia_h2he, cia_hhe};
pub use lte_opacity::{
LteOpacityParams, LteOpacityOutput, LteFrequencyGrid,
lte_meanopt, generate_lte_frequency_grid, quick_lte_rosseland,
compute_opacity_at_frequency,
lte_meanopt, generate_lte_frequency_grid, generate_inifrc_frequency_grid,
quick_lte_rosseland, compute_opacity_at_frequency,
};
pub use lte_opacity::{
LteOpacityParams as LteOpacityParamsOld, LteOpacityOutput as LteOpacityOutputOld,

134
src/tlusty/math/continuum/opacf0.rs
View File

@ -20,6 +20,9 @@ use crate::tlusty::state::{GffPar, DwnPar, ModPar, InpPar};
use crate::tlusty::math::atomic::{gfree0, gfree1};
use crate::tlusty::math::opacity::dwnfr0;
use crate::tlusty::math::opacity::dwnfr1;
use crate::tlusty::math::hydrogen::hephot;
use crate::tlusty::math::special::gaunt;
use crate::tlusty::state::{AtomicData, ModelState};
// 物理常数 (来自 opacf0.f)
/// Rydberg 频率
@ -311,7 +314,7 @@ pub struct Opacf0Output<'a> {
// 束缚-自由截面数据
/// BFCS - 光电离截面表 (mcross × nfreqc)
pub bfcs: &'a [f32],
pub bfcs: &'a [Vec<f32>],
/// IJBF - 频率插值索引 (nfreq)
pub ijbf: &'a [i32],
/// AIJBF - 频率插值系数 (nfreq)
@ -341,9 +344,9 @@ fn cross(ibft: usize, ij: usize, output: &Opacf0Output) -> f64 {
let ij0 = output.ijbf[ij] as usize;
let a1 = output.aijbf[ij];
// BFCS(MCROSS, MFREQC) in column-major: ibft + ij * MCROSS
let sig0 = output.bfcs[ibft + ij0 * MCROSS] as f64;
let sig1 = output.bfcs[ibft + (ij0 + 1) * MCROSS] as f64;
// BFCS[ibft][ij0]
let sig0 = output.bfcs[ibft][ij0] as f64;
let sig1 = output.bfcs[ibft][ij0 + 1] as f64;
a1 * sig0 + (UN - a1) * sig1
}
@ -1209,3 +1212,126 @@ mod tests {
assert_eq!(ij1, 1); // ij0 - 1
}
}
// ============================================================================
// BFCS 初始化
// ============================================================================
/// 初始化光电离截面表 (BFCS)。
/// 对应 Fortran `initia.f` 和 `sbfhe1.f` 中设置 `BFCS` 的逻辑。
pub fn init_bfcs(model: &mut ModelState, atomic: &AtomicData, nfreq: usize) {
let ntranc = model.obfpar.itrbf.iter().filter(|&&x| x > 0).count();
let nfreqc = nfreq;
let freq = &model.frqall.freq;
// 获取 He I 量子数映射
let he1_sln: [(usize, i32, i32, i32); 14] = [
(10, 1, 0, 1), // 1 sing S
(11, 3, 0, 2), // 2 trip S
(12, 1, 0, 2), // 2 sing S
(13, 3, 1, 2), // 2 trip P
(14, 1, 1, 2), // 2 sing P
(15, 3, 0, 3), // 3 trip S
(16, 1, 0, 3), // 3 sing S
(17, 3, 1, 3), // 3 trip P
(18, 3, 2, 3), // 3 trip D
(19, 1, 2, 3), // 3 sing D
(20, 1, 1, 3), // 3 sing P
(21, 3, 0, 4), // 4 trip S
(22, 1, 0, 4), // 4 sing S
(23, 3, 1, 4), // 4 trip P
];
/// Hydrogen photoionization cross-section constant
/// Fortran SIGK: SIH0 = 2.815D29
const SIH0: f64 = 2.815e29;
for ibft in 0..ntranc {
let itr = model.obfpar.itrbf[ibft] as usize - 1;
if itr >= atomic.trapar.ilow.len() { continue; }
let ilow = atomic.trapar.ilow[itr] as usize - 1;
let ifwop = model.wmcomp.ifwop[ilow];
let fr0 = atomic.trapar.fr0[itr];
let fr0pc = fr0;
// Get element and quantum numbers for hydrogenic Gaunt factor
let ie = atomic.levpar.iel[ilow] as usize - 1;
let iz = atomic.ionpar.iz[ie] as f64; // nuclear charge Z
let ch = iz * iz; // Z²
let nq = atomic.levpar.nquant[ilow] as f64; // principal quantum number
let imer = model.mrgpar.imrg[ilow];
let (sgem, frch) = if imer > 0 {
let imer_idx = (imer - 1) as usize;
(model.mrgpar.sgm0[imer_idx], model.mrgpar.frch[imer_idx])
} else {
(0.0, 0.0)
};
// Get IBF mode from PhoSet
let ic = if itr < atomic.trapar.itrcon.len() {
atomic.trapar.itrcon[itr] as usize
} else {
0
};
let ibf = if ic > 0 && ic - 1 < atomic.phoset.ibf.len() {
atomic.phoset.ibf[ic - 1]
} else {
1 // default: hydrogenic with exact Gaunt
};
// 查找是否为 He I 能级
let mut he1_quant = None;
for &(lvl, s, l, n) in &he1_sln {
if ilow == lvl {
he1_quant = Some((s, l, n));
break;
}
}
for ij in 0..nfreqc {
let fr = freq[ij];
// He I (SBFHE1/HEPHOT logic) — IBF=11 or 13
if let Some((s, l, n)) = he1_quant {
let mut cs = hephot(s, l, n, fr);
if fr < fr0pc {
cs = 0.0;
}
model.phoexp.bfcs[ibft][ij] = cs as f32;
continue;
}
// Hydrogenic cross-section with frequency-dependent Gaunt factor
// Fortran SIGK with IBF=1:
// SIGK = SIH0/FR³ * CH²/NQ⁵ * gaunt(nq, fr/ch)
// This includes the exact Gaunt factor, unlike the old
// simplified formula sbf*(fr0/fr)³ which only had edge Gaunt.
let mut cs = 0.0;
if fr >= fr0pc {
if imer > 0 {
// Merged level: sgem*(frch/fr)³ (from SGMER0)
let ratio = frch / fr;
cs = sgem * ratio * ratio * ratio;
} else if ibf == 0 {
// IBF=0: hydrogenic without Gaunt factor (Gaunt=1)
let ratio = fr0pc / fr;
let sbf_edge = SIH0 * ch * ch / (fr0pc.powi(3) * nq.powi(5));
cs = sbf_edge * ratio * ratio * ratio;
} else {
// IBF=1: hydrogenic with exact Gaunt factor
// SIGK = SIH0 * CH² / (FR³ × NQ⁵) × gaunt(NQ, FR/CH)
cs = SIH0 * ch * ch / (fr.powi(3) * nq.powi(5)) * gaunt(nq as usize, fr / ch);
}
}
let mut bfcs_val = cs as f32;
if ifwop < 0 {
bfcs_val = 1e-30;
}
model.phoexp.bfcs[ibft][ij] = bfcs_val;
}
}
}

762
src/tlusty/math/continuum/opacf0_state.rs Normal file
View File

@ -0,0 +1,762 @@
//! Opacf0State — OPACF0 连续谱不透明度状态。
//!
//! 预计算 BF 截面参数和 FF 离子数据，用于在 RESOLV 频率循环中
//! 计算不透明度时使用模型种群 (POPUL) 而非重新计算 Saha 平衡。
//!
//! # 与 Fortran OPACF0 的对应
//!
//! Fortran OPACF0(ID, NFRQ) 对每个深度 ID 计算所有频率的不透明度。
//! 使用预计算的 CROSS(IBFT, IJ) 截面和 POPUL(II, ID) 模型种群。
//!
//! 本模块将静态原子数据（BF 跃迁列表、FF 离子数据）预计算到 Opacf0State，
//! 然后在 RESOLV 中使用 model.levpop.popul 计算不透明度。
use crate::tlusty::math::hydrogen::hephot;
use crate::tlusty::math::atomic::{gfree1, gfree0};
use crate::tlusty::math::special::gaunt;
use crate::tlusty::math::hydrogen::sgmer1;
use crate::tlusty::state::constants::{HK, SIGE, H, UN, MDEPTH, MLEVEL};
use crate::tlusty::state::model::{GffPar, ModelState};
use crate::tlusty::state::atomic::AtomicData;
// ============================================================================
// 常量 (from Fortran OPACF0)
// ============================================================================
/// SIGK 常数: SIH0 = 2.815e29 (氢原子 bf 截面基准, cm² × Hz³)
const SIH0: f64 = 2.815e29;
/// FF 基准常数: SGFF0 = 3.694e8
const SGFF0: f64 = 3.694e8;
/// C14 = 2.99793e14 (用于 Gaunt 因子计算)
const C14: f64 = 2.99793e14;
// ============================================================================
// BF 跃迁数据
// ============================================================================
/// 束缚-自由跃迁信息。
#[derive(Debug, Clone)]
pub struct BfTransition {
/// 下能级索引 (0-based)
pub ilow: usize,
/// 上能级 (连续区) 索引 (0-based)
pub iup: usize,
/// 离子索引 (0-based)
pub ion_idx: usize,
/// 有效主量子数
pub nquant: f64,
/// Z² (原子序数的平方)
pub iz2: f64,
/// 电离频率 (Hz), ν_n = ENION(n)/H
pub nu_edge: f64,
/// He I 量子数 (s, l, n) — 使用 HEPHOT OP 截面而非氢原子 SIGK
pub he1_quant: Option<(i32, i32, i32)>,
}
// ============================================================================
// FF 离子数据
// ============================================================================
/// 自由-自由离子信息。
#[derive(Debug, Clone)]
pub struct FfIon {
/// 连续区能级索引 NNEXT (0-based)
pub nnext: usize,
/// Z² (CHARG2)
pub charg2: f64,
/// FF 模式: 1=hydrogenic Gaunt=1, 2=exact Gaunt (GFREE1)
pub itra: i32,
/// FF 边缘频率 (Hz), 用于 stimulated emission switch
pub ff_edge: f64,
}
// ============================================================================
// Opacf0State
// ============================================================================
/// OPACF0 预计算状态。
///
/// 存储从原子数据预计算的 BF 跃迁和 FF 离子信息，
/// 在 RESOLV 中与 model.levpop.popul 一起使用计算不透明度。
#[derive(Debug, Clone)]
pub struct Opacf0State {
/// BF 跃迁列表 (对应 Fortran NTRANC 个 ITRBF 跃迁)
pub bf_transitions: Vec<BfTransition>,
/// FF 离子列表 (对应 Fortran NION 个离子)
pub ff_ions: Vec<FfIon>,
/// 是否使用 HEPHOT OP 截面 (替代氢原子 SIGK) 计算 He I bf
/// 注意: 启用后需要完整 OPACF0+SOLVES 管线才能收敛
pub use_hephot: bool,
}
impl Opacf0State {
/// 创建 HHe 模型的 Opacf0State。
///
/// HHe 模型结构 (Rust 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
/// - NLEVEL=39, NION=3
pub fn new_hhe() -> Self {
let h_planck: f64 = HK * 1.3806e-16; // H = HK × k = 6.6256e-27 erg·s
let eh: f64 = 2.17853041e-11; // H 电离能 (erg)
let mut bf_transitions = Vec::new();
// === H I bf 跃迁 (n=1..9 → continuum level 9) ===
// Hydrogenic: ν_n = EH/(n² × H)
for n in 1..=9_usize {
let nn = n as f64;
let nu_n = eh / (nn * nn * h_planck); // ionization frequency of level n
bf_transitions.push(BfTransition {
ilow: n - 1,
iup: 9,
ion_idx: 0,
nquant: nn,
iz2: 1.0, // Z=1, Z²=1
nu_edge: nu_n,
he1_quant: None, // H I: use hydrogenic SIGK
});
}
// === He I bf 跃迁 (14 levels → continuum level 24) ===
// 使用实际频率数据 (from he1.dat) + HEPHOT OP 截面
// 量子数 (s, l, n) from he1.dat level definitions
let he1_freq: [f64; 14] = [
5.94503520e15, 1.15267210e15, 9.60145430e14, 8.75933720e14,
8.14536220e14, 4.51727350e14, 4.02921120e14, 3.81935640e14,
3.65836790e14, 3.65746870e14, 3.62599020e14, 2.40043860e14,
2.20797190e14, 2.12492940e14,
];
let he1_nq: [f64; 14] = [
1.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 4.0, 4.0, 4.0,
];
// He I quantum numbers (s, l, n) from he1.dat / sbfhe1.f
let he1_sln: [(i32, i32, i32); 14] = [
(1, 0, 1), // 1 sing S (G=1)
(3, 0, 2), // 2 trip S (G=3)
(1, 0, 2), // 2 sing S (G=1)
(3, 1, 2), // 2 trip P (G=9)
(1, 1, 2), // 2 sing P (G=3)
(3, 0, 3), // 3 trip S (G=3)
(1, 0, 3), // 3 sing S (G=1)
(3, 1, 3), // 3 trip P (G=9)
(3, 2, 3), // 3 trip D (G=15)
(1, 2, 3), // 3 sing D (G=5)
(1, 1, 3), // 3 sing P (G=3)
(3, 0, 4), // 4 trip S (G=3)
(1, 0, 4), // 4 sing S (G=1)
(3, 1, 4), // 4 trip P (G=9)
];
for ilev in 0..14 {
bf_transitions.push(BfTransition {
ilow: 10 + ilev,
iup: 24,
ion_idx: 1,
nquant: he1_nq[ilev],
iz2: 4.0,
nu_edge: he1_freq[ilev],
he1_quant: Some(he1_sln[ilev]),
});
}
// === He II bf 跃迁 (n=1..14 → continuum level 38) ===
// Hydrogenic Z=2: ν_n = 4×EH/(n² × H)
let he2_freq: [f64; 14] = [
1.31575980e16, 3.28939940e15, 1.46195530e15, 8.22349860e14,
5.26303910e14, 3.65488820e14, 2.68522400e14, 2.05587460e14,
1.62439480e14, 1.31575980e14, 1.08740480e14, 9.13722060e13,
7.78556080e13, 6.71306000e13,
];
for n in 1..=14_usize {
let nn = n as f64;
bf_transitions.push(BfTransition {
ilow: 24 + n - 1,
iup: 38,
ion_idx: 2,
nquant: nn,
iz2: 4.0, // Z=2, Z²=4
nu_edge: he2_freq[n - 1],
he1_quant: None, // He II: hydrogenic, use SIGK
});
}
// === FF 离子 ===
let ff_ions = vec![
FfIon {
nnext: 9, // H II (continuum)
charg2: 1.0, // Z²=1 for H⁺
itra: 2, // exact Gaunt (GFREE1)
ff_edge: 0.0, // H: FF(ION)=0 → no edge
},
FfIon {
nnext: 24, // He II (continuum)
charg2: 1.0, // Z²=1 for He⁺
itra: 2,
ff_edge: 0.0,
},
FfIon {
nnext: 38, // He III (continuum)
charg2: 4.0, // Z²=4 for He²⁺
itra: 2,
ff_edge: 0.0,
},
];
Self {
bf_transitions,
ff_ions,
use_hephot: true, // He I: use Opacity Project cross-sections (HEPHOT)
}
}
/// 计算给定 (depth, frequency) 的连续谱不透明度。
///
/// 对应 Fortran OPACF0 的频率循环体 (行 168-346)。
/// 使用模型种群 POPUL(level, depth) 而非重新计算 Saha 平衡。
///
/// # 参数
///
/// * `fr` - 频率 (Hz)
/// * `t` - 温度 (K)
/// * `ne` - 电子密度 (cm⁻³)
/// * `popul` - 种群数组 popul[level][depth] (cm⁻³)
/// * `id` - 深度索引 (0-based)
/// * `gffpar` - Gaunt 因子预计算参数 (由 gfree0 填充)
///
/// # 返回
///
/// Returns (true_abs, scat, emis_pre_stim, rayleigh) where:
/// - true_abs: true absorption after stimulated emission correction (cm⁻¹)
/// - scat: total scattering = elscat + rayleigh (cm⁻¹)
/// - emis_pre_stim: accumulated emission before stimulated emission correction (cm⁻¹)
/// Needed for OPACFL finalization: EMIS1L = emis_pre_stim * XKF * Bν
/// - rayleigh: Rayleigh scattering only (cm⁻¹), needed for SCAT1 in Lucy
pub fn compute_opacity(
&self,
fr: f64,
t: f64,
ne: f64,
popul: &[Vec<f64>],
id: usize,
gffpar: &GffPar,
) -> (f64, f64, f64, f64) {
let hkt = HK / t;
let sqrt_t = t.sqrt();
let tk1 = 1.0 / (H * t); // 1/(kT) in erg⁻¹
// 电子散射 (对应 Fortran: ELSCAT = ELEC(ID) × SIGE)
let elscat = SIGE * ne;
// 预计算频率相关量
let fr_inv = 1.0 / fr;
let fr3_inv = fr_inv * fr_inv * fr_inv;
// 受激发射因子 (XKF = exp(-hν/kT))
let xkf = if hkt * fr < 150.0 {
(-hkt * fr).exp()
} else {
0.0
};
let mut abso = elscat;
let mut emis = 0.0;
// ================================================================
// 1. 束缚-自由贡献
// ================================================================
// 对应 Fortran: DO 30 IBFT=1,NTRANC → ABSO += CROSS × ABTRA; EMIS += CROSS × EMTRA
for bt in &self.bf_transitions {
if fr < bt.nu_edge {
continue;
}
let pop_low = popul[bt.ilow][id];
if pop_low <= 0.0 {
continue;
}
// 截面计算:
// He I: HEPHOT (OP cross-sections) if enabled
// H I / He II: hydrogenic SIGK with exact Gaunt factor (IBF=1)
// SIGK = SIH0 * Z^4 / (nu^3 * n^5) * GAUNT(n, nu/Z^2)
let sigma = if self.use_hephot {
if let Some((s, l, n)) = bt.he1_quant {
hephot(s, l, n, fr)
} else {
let sigma_base = SIH0 * bt.iz2 * bt.iz2 * fr3_inv / bt.nquant.powi(5);
let g_bf = gaunt(bt.nquant as usize, fr / bt.iz2);
sigma_base * g_bf
}
} else {
let sigma_base = SIH0 * bt.iz2 * bt.iz2 * fr3_inv / bt.nquant.powi(5);
let g_bf = gaunt(bt.nquant as usize, fr / bt.iz2);
sigma_base * g_bf
};
if sigma <= 0.0 {
continue;
}
// 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(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;
}
// ================================================================
// 2. 自由-自由贡献
// ================================================================
// 对应 Fortran: DO 40 ION=1,NION
let sgff = SGFF0 / sqrt_t * ne;
for ion in &self.ff_ions {
let pop_cont = popul[ion.nnext][id];
if pop_cont <= 0.0 {
continue;
}
// SF1 = SFF3(ION,ID) × FR3INV
let sf1 = sgff * ion.charg2 * pop_cont * fr3_inv;
// SF2 = SFF2(ION,ID) = exp(FF(ION) × HKT1(ID))
// FF(ION) = 0 for our model → SF2 = 1
let sf2 = 1.0_f64; // FF=0 → exp(0) = 1
let mut absoff = sf1 * sf2;
// Exact Gaunt factor (ITRA=2)
// Fortran OPACF0 line ~210: SFF3(ION,ID) = GFREE1(ID,X) * SFF2(ION,ID)
// where X = C14 * CHARG2(ION) / FR
if ion.itra == 2 {
let x = C14 * ion.charg2 / fr;
let gf1 = gfree1(id, x, gffpar);
// SFF3 = GF1 * SFF2 → absoff = SF1 * GF1 (replacing the default G=1)
absoff = sf1 * gf1;
}
abso += absoff;
emis += absoff; // ff: Kirchhoff (emission = absorption for thermal ff)
}
// ================================================================
// 受激发射修正
// ================================================================
// 对应 Fortran: ABSO(IJ) = ABSO(IJ) - EMIS(IJ) × XKF(ID)
abso -= emis * xkf;
// ================================================================
// 3. OPADD0 — 额外不透明度来源 (Rayleigh 散射)
// ================================================================
// 对应 Fortran OPADD0: HI Rayleigh, HeI Rayleigh
let mut rayleigh = 0.0;
// HI Rayleigh scattering (IRSCT=1 in Fortran)
// σ_Ray(HI) = (CR0 + (CR1 + CR2/X)/X) / X², X = (c/λ)²
// where λ = c/ν, so X = (c/ν)² = c²/ν²
// clamped at FRRAY = 2.463e15 Hz
{
const CLS: f64 = 2.997925e18;
const CR0: f64 = 5.799e-13;
const CR1: f64 = 1.422e-6;
const CR2: f64 = 2.784;
const FRRAY: f64 = 2.463e15;
let frm = fr.min(FRRAY);
let x = (CLS / frm) * (CLS / frm); // X = (c/λ)²
let cs_ray_hi = (CR0 + (CR1 + CR2 / x) / x) / (x * x);
// POPUL(0, id) = H I ground state population
rayleigh += popul[0][id] * cs_ray_hi;
}
// HeI Rayleigh scattering (IRSCHE=1 in Fortran)
// σ_Ray(HeI) = 5.484e-14 / X² * (1 + (2.44e5 + 5.94e10/(X-2.90e5))/X)²
// clamped at FRAYHe = 5.150e15 Hz
{
const CLS: f64 = 2.997925e18;
const FRAYHE: f64 = 5.150e15;
let frm = fr.min(FRAYHE);
let x = (CLS / frm) * (CLS / frm);
let cs_ray_he1 =
5.484e-14 / (x * x) * (1.0 + (2.44e5 + 5.94e10 / (x - 2.90e5)) / x).powi(2);
// He I ground state is level 10 (index 10)
rayleigh += popul[10][id] * cs_ray_he1;
}
// 分离散射和真吸收
// abso = elscat + true_absorption (after stimulated emission)
// true_abs = abso - elscat = emis × (1 - XKF)
let true_abs = (abso - elscat).max(0.0);
let scat = elscat + rayleigh;
(true_abs, scat, emis, rayleigh)
}
/// 使用真实模型数据（BFCS 和 SGMER）计算连续谱不透明度，完全匹配 Fortran OPACF0
pub fn compute_opacity_from_model(
fr: f64,
ij: usize,
t: f64,
ne: f64,
popul: &[Vec<f64>],
id: usize,
gffpar: &GffPar,
model: &ModelState,
atomic: &AtomicData,
) -> (f64, f64, f64, f64) {
let hkt = HK / t;
let sqrt_t = t.sqrt();
let tk1 = 1.0 / (H * t);
let elscat = SIGE * ne;
let fr_inv = 1.0 / fr;
let fr3_inv = fr_inv * fr_inv * fr_inv;
let xkf = if hkt * fr < 150.0 {
(-hkt * fr).exp()
} else {
0.0
};
let mut abso = elscat;
let mut emis = 0.0;
// ================================================================
// 1. 束缚-自由贡献
// ================================================================
let ntranc = model.obfpar.itrbf.iter().filter(|&&x| x > 0).count();
for ibft in 0..ntranc {
let itr = model.obfpar.itrbf[ibft] as usize - 1;
if atomic.trapar.indexp[itr] == 0 { continue; }
let ii = atomic.trapar.ilow[itr] as usize - 1;
let jj = atomic.trapar.iup[itr] as usize - 1;
let it = atomic.trapar.itra[ii][jj] as usize;
if it == 0 { continue; }
let pop_low = popul[ii][id];
// 交叉截面 CROSS
let ij0 = model.frqall.ijbf[ij] as usize;
let a1 = model.phoexp.aijbf[ij];
let sig0 = model.phoexp.bfcs[ibft][ij0] as f64;
let sig1 = model.phoexp.bfcs[ibft][ij0 + 1] as f64;
let mut sigma = a1 * sig0 + (UN - a1) * sig1;
// SGMER 修正 (合并能级)
let imer = model.mrgpar.imrg[ii];
if imer > 0 {
sigma = sgmer1(fr_inv, fr3_inv, imer, id + 1, &model.mrgpar);
}
if sigma <= 0.0 { continue; }
// CORR 计算
let ie = atomic.levpar.iel[ii] as usize - 1;
let nke = atomic.ionpar.nnext[ie] as usize - 1;
let corr = if nke != jj {
let g_ratio = atomic.levpar.g[nke] / atomic.levpar.g[jj];
let delta_e = atomic.levpar.enion[nke] - atomic.levpar.enion[jj];
g_ratio * (delta_e * tk1).exp()
} else {
UN
};
let pop_up = popul[jj][id];
let wop_ii = model.wmcomp.wop[ii][id];
let emis_val = pop_up * ne * model.levpop.sbf[ii] * wop_ii * corr;
abso += sigma * pop_low;
emis += sigma * emis_val;
}
// ================================================================
// 2. 自由-自由贡献
// ================================================================
let sgff = SGFF0 / sqrt_t * ne;
let nion = atomic.ionpar.iz.iter().filter(|&&x| x > 0).count();
for ion in 0..nion {
let nnext = atomic.ionpar.nnext[ion] as usize - 1;
let pop_cont = popul[nnext][id];
if pop_cont <= 0.0 { continue; }
let charg2 = atomic.ionpar.charg2[ion];
let sf1 = sgff * charg2 * pop_cont * fr3_inv;
// sf2 = exp(ff * hkt1)
let ff_val = atomic.ionpar.ff[ion];
let mut sf2 = (ff_val * hkt).exp();
if fr < ff_val {
sf2 = UN / xkf;
}
let mut absoff = sf1 * sf2;
// Gaunt factor
let itra = atomic.trapar.itra[nnext][nnext];
if itra == 2 {
let x = C14 * charg2 / fr;
let gf1 = gfree1(id, x, gffpar);
sf2 = sf2 - UN + gf1;
absoff = sf1 * sf2;
}
abso += absoff;
emis += absoff;
}
abso -= emis * xkf;
// ================================================================
// 3. OPADD0 — 额外不透明度来源 (Rayleigh 散射)
// ================================================================
let mut rayleigh = 0.0;
// HI Rayleigh scattering
if nion > 0 {
const CLS: f64 = 2.997925e18;
const CR0: f64 = 5.799e-13;
const CR1: f64 = 1.422e-6;
const CR2: f64 = 2.784;
const FRRAY: f64 = 2.463e15;
let frm = fr.min(FRRAY);
let x = (CLS / frm) * (CLS / frm);
let cs_ray_hi = (CR0 + (CR1 + CR2 / x) / x) / (x * x);
// 假设 H I 基态为能级 0
if atomic.levpar.iel[0] == 1 {
rayleigh += popul[0][id] * cs_ray_hi;
}
}
// HeI Rayleigh scattering
if nion > 1 {
const CLS: f64 = 2.997925e18;
const FRAYHE: f64 = 5.150e15;
let frm = fr.min(FRAYHE);
let x = (CLS / frm) * (CLS / frm);
let cs_ray_he1 =
5.484e-14 / (x * x) * (1.0 + (2.44e5 + 5.94e10 / (x - 2.90e5)) / x).powi(2);
// 寻找 He I 基态
let nlevel = atomic.levpar.iel.iter().filter(|&&x| x > 0).count();
for i in 0..nlevel {
if atomic.levpar.iel[i] == 2 && atomic.ionpar.iz[1] == 1 { // He I
rayleigh += popul[i][id] * cs_ray_he1;
break;
}
}
}
let true_abs = (abso - elscat).max(0.0);
let scat = elscat + rayleigh;
(true_abs, scat, emis, rayleigh)
}
/// 计算给定 (frequency, temperature, populations) 的连续谱不透明度。
///
/// 单深度版本，直接接收 `popul[level]` 扁平数组（无需 2D 数组 + depth 索引）。
/// 用于有限差分温度导数计算（需要扰动种群但不修改全局模型）。
pub fn compute_opacity_1d(
&self,
fr: f64,
t: f64,
ne: f64,
popul: &[f64],
gffpar: &GffPar,
) -> (f64, f64, f64, f64) {
let hkt = HK / t;
let sqrt_t = t.sqrt();
let elscat = SIGE * ne;
let fr_inv = 1.0 / fr;
let fr3_inv = fr_inv * fr_inv * fr_inv;
let xkf = if hkt * fr < 150.0 {
(-hkt * fr).exp()
} else {
0.0
};
let mut abso = elscat;
let mut emis = 0.0;
// BF contributions
for bt in &self.bf_transitions {
if fr < bt.nu_edge { continue; }
let pop_low = popul[bt.ilow];
if pop_low <= 0.0 { continue; }
let sigma = if self.use_hephot {
if let Some((s, l, n)) = bt.he1_quant {
hephot(s, l, n, fr)
} else {
let sigma_base = SIH0 * bt.iz2 * bt.iz2 * fr3_inv / bt.nquant.powi(5);
let g_bf = gaunt(bt.nquant as usize, fr / bt.iz2);
sigma_base * g_bf
}
} else {
let sigma_base = SIH0 * bt.iz2 * bt.iz2 * fr3_inv / bt.nquant.powi(5);
let g_bf = gaunt(bt.nquant as usize, fr / bt.iz2);
sigma_base * g_bf
};
if sigma <= 0.0 { continue; }
abso += sigma * pop_low;
emis += sigma * pop_low;
}
// FF contributions
let sgff = SGFF0 / sqrt_t * ne;
for ion in &self.ff_ions {
let pop_cont = popul[ion.nnext];
if pop_cont <= 0.0 { continue; }
let sf1 = sgff * ion.charg2 * pop_cont * fr3_inv;
let mut absoff = sf1; // SF2=1 (FF=0)
if ion.itra == 2 {
let x = C14 * ion.charg2 / fr;
let gf1 = gfree1(0, x, gffpar);
absoff = sf1 * gf1;
}
abso += absoff;
emis += absoff;
}
abso -= emis * xkf;
let true_abs = (abso - elscat).max(0.0);
// Rayleigh scattering (1D version — same formula as compute_opacity)
let mut rayleigh = 0.0;
{
const CLS: f64 = 2.997925e18;
const CR0: f64 = 5.799e-13;
const CR1: f64 = 1.422e-6;
const CR2: f64 = 2.784;
const FRRAY: f64 = 2.463e15;
let frm = fr.min(FRRAY);
let x = (CLS / frm) * (CLS / frm);
let cs_ray_hi = (CR0 + (CR1 + CR2 / x) / x) / (x * x);
rayleigh += popul[0] * cs_ray_hi;
}
{
const CLS: f64 = 2.997925e18;
const FRAYHE: f64 = 5.150e15;
let frm = fr.min(FRAYHE);
let x = (CLS / frm) * (CLS / frm);
let cs_ray_he1 =
5.484e-14 / (x * x) * (1.0 + (2.44e5 + 5.94e10 / (x - 2.90e5)) / x).powi(2);
rayleigh += popul[10] * cs_ray_he1;
}
(true_abs, elscat + rayleigh, emis, rayleigh)
}
}
impl Default for Opacf0State {
fn default() -> Self {
Self::new_hhe()
}
}
// ============================================================================
// 测试
// ============================================================================
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_opacf0_state_hhe() {
let state = Opacf0State::new_hhe();
// H I: 9 bf, He I: 14 bf, He II: 14 bf → total 37
assert_eq!(state.bf_transitions.len(), 37);
assert_eq!(state.ff_ions.len(), 3);
// Check first H I bf transition
assert_eq!(state.bf_transitions[0].ilow, 0); // H I n=1
assert_eq!(state.bf_transitions[0].iup, 9); // H II
assert_eq!(state.bf_transitions[0].iz2, 1.0);
// Check first He I bf transition
assert_eq!(state.bf_transitions[9].ilow, 10); // He I level 0
assert_eq!(state.bf_transitions[9].iup, 24); // He II
// Check first He II bf transition
assert_eq!(state.bf_transitions[23].ilow, 24); // He II n=1
assert_eq!(state.bf_transitions[23].iup, 38); // He III
assert_eq!(state.bf_transitions[23].iz2, 4.0);
}
#[test]
fn test_compute_opacity_basic() {
let state = Opacf0State::new_hhe();
let nlevel = 39;
let nd = 1;
// Create simple populations
let mut popul = vec![vec![0.0_f64; nd]; nlevel];
popul[0][0] = 1e12; // H I n=1
popul[9][0] = 1e14; // H II
popul[10][0] = 1e10; // He I n=1
popul[24][0] = 1e10; // He II n=1
popul[38][0] = 1e12; // He III
let fr = 3.28805e15; // H Lyman edge
let t = 30000.0;
let ne = 1e14;
// Precompute Gaunt factor coefficients
let mut gffpar = GffPar::new();
let temp = [30000.0_f64; MDEPTH];
gfree0(0, &temp, &mut gffpar);
let (ab, sct, _emis, ray) = state.compute_opacity(fr, t, ne, &popul, 0, &gffpar);
assert!(ab > 0.0, "absorption should be positive");
assert!(sct > 0.0, "scattering should be positive");
assert!(ab > sct, "at H Lyman edge, absorption should dominate scattering");
assert!(ray > 0.0, "Rayleigh scattering should be positive");
}
#[test]
fn test_compute_opacity_consistency() {
// Test that opacity is consistent across nearby frequencies
let state = Opacf0State::new_hhe();
let nlevel = 39;
let nd = 1;
let mut popul = vec![vec![0.0_f64; nd]; nlevel];
popul[0][0] = 1e12;
popul[9][0] = 1e14;
let t = 30000.0;
let ne = 1e14;
let mut gffpar = GffPar::new();
let temp = [30000.0_f64; MDEPTH];
gfree0(0, &temp, &mut gffpar);
// Just above H Lyman edge
let fr1 = 3.3e15;
// Just below H Lyman edge
let fr2 = 3.2e15;
let (ab1, _, _, _) = state.compute_opacity(fr1, t, ne, &popul, 0, &gffpar);
let (ab2, _, _, _) = state.compute_opacity(fr2, t, ne, &popul, 0, &gffpar);
// Above edge: bf contributes → higher opacity
// Below edge: no H I bf → lower opacity
assert!(ab1 > ab2, "opacity above Lyman edge should be higher");
}
}

211
src/tlusty/math/hydrogen/hephot.rs
View File

@ -1,10 +1,13 @@
//! He I 光电离截面。
//!
//! 重构自 TLUSTY `hephot.f`
//!
//! 使用 Seaton-Fernley 三次拟合计算 Opacity Project 截面。
//! 对 L > 2 使用类氢公式。
/// He I 光电离截面。
///
/// 使用 Seaton 和 Fernley 的三次拟合计算 Opacity Project 截面。
/// 使用 Seaton-Fernley 的三次拟合计算 Opacity Project 截面。
///
/// # 参数
///
@ -20,14 +23,21 @@
/// # 备注
///
/// 对于 L > 2 使用类氢公式。
/// 低于阈值频率返回 0。
pub fn hephot(s: i32, l: i32, n: i32, freq: f64) -> f64 {
const TENM18: f64 = 1e-18;
const FRH: f64 = 3.28805e15;
const TENLG: f64 = 2.302585093;
const PHOT0: f64 = 2.815e29;
// 系数数据（仅包含必要的）
// 完整数据较长，此处仅保留关键数据
// IST(3,2) — starting index into COEF/A/B/XFITM arrays
// indexed as IST(LL, SS) where LL=L+1, SS=(S+1)/2
const IST: [[usize; 2]; 3] = [[1, 11], [36, 45], [20, 28]];
// N0(3,2) — minimum principal quantum number for each (L, S) combination
const N0: [[i32; 2]; 3] = [[1, 2], [2, 2], [3, 3]];
// FL0(53) — threshold frequency in log10(freq/FRH)
const FL0: [f64; 53] = [
2.521e-01, -5.381e-01, -9.139e-01, -1.175e00, -1.375e00, -1.537e00,
-1.674e00, -1.792e00, -1.896e00, -1.989e00, -4.555e-01, -8.622e-01,
@ -40,29 +50,158 @@ pub fn hephot(s: i32, l: i32, n: i32, freq: f64) -> f64 {
-1.547e00, -1.682e00, -1.799e00, -1.902e00, -1.995e00,
];
// XFITM(53) — maximum X for cubic fit; beyond this use linear fit A+B*X
const XFITM: [f64; 53] = [
3.262e-01, 6.135e-01, 9.233e-01, 8.438e-01, 1.020e00, 1.169e00,
1.298e00, 1.411e00, 1.512e00, 1.602e00, 7.228e-01, 1.076e00,
1.206e00, 1.404e00, 1.481e00, 1.464e00, 1.581e00, 1.685e00,
1.777e00, 9.586e-01, 1.187e00, 1.371e00, 1.524e00, 1.740e00,
1.854e00, 1.955e00, 2.046e00, 9.585e-01, 1.041e00, 1.371e00,
1.608e00, 1.739e00, 1.768e00, 1.869e00, 1.803e00, 7.360e-01,
1.041e00, 1.272e00, 1.457e00, 1.611e00, 1.741e00, 1.855e00,
1.870e00, 1.804e00, 9.302e-01, 1.144e00, 1.028e00, 1.210e00,
1.362e00, 1.646e00, 1.761e00, 1.863e00, 1.954e00,
];
// A(53) — linear fit coefficient (log sigma intercept)
const A: [f64; 53] = [
6.95319e-01, 1.13101e00, 1.36313e00, 1.51684e00, 1.64767e00,
1.75643e00, 1.84458e00, 1.87243e00, 1.85628e00, 1.90889e00,
9.01802e-01, 1.25389e00, 1.39033e00, 1.55226e00, 1.60658e00,
1.65930e00, 1.68855e00, 1.62477e00, 1.66726e00, 1.83599e00,
2.50403e00, 3.08564e00, 3.56545e00, 4.25922e00, 4.61346e00,
4.91417e00, 5.19211e00, 1.74181e00, 2.25756e00, 2.95625e00,
3.65899e00, 4.04397e00, 4.13410e00, 4.43538e00, 4.19583e00,
1.79027e00, 2.23543e00, 2.63942e00, 3.02461e00, 3.35018e00,
3.62067e00, 3.85218e00, 3.76689e00, 3.49318e00, 1.16294e00,
1.86467e00, 2.02110e00, 2.24231e00, 2.44240e00, 2.76594e00,
2.93230e00, 3.08109e00, 3.21069e00,
];
// B(53) — linear fit coefficient (log sigma slope)
const B: [f64; 53] = [
-1.29000e00, -2.15771e00, -2.13263e00, -2.10272e00, -2.10861e00,
-2.11507e00, -2.11710e00, -2.08531e00, -2.03296e00, -2.03441e00,
-1.85905e00, -2.04057e00, -2.02189e00, -2.05930e00, -2.03403e00,
-2.02071e00, -1.99956e00, -1.92851e00, -1.92905e00, -4.58608e00,
-4.40022e00, -4.39154e00, -4.39676e00, -4.57631e00, -4.57120e00,
-4.56188e00, -4.55915e00, -4.41218e00, -4.12940e00, -4.24401e00,
-4.40783e00, -4.39930e00, -4.25981e00, -4.26804e00, -4.00419e00,
-4.47251e00, -3.87960e00, -3.71668e00, -3.68461e00, -3.67173e00,
-3.65991e00, -3.64968e00, -3.48666e00, -3.23985e00, -2.95758e00,
-3.07110e00, -2.87157e00, -2.83137e00, -2.82132e00, -2.91084e00,
-2.91159e00, -2.91336e00, -2.91296e00,
];
// COEF(4, 53) — cubic fit coefficients: P = C0 + C1*X + C2*X^2 + C3*X^3
// Stored column-major: COEF = [c0_1, c1_1, c2_1, c3_1, c0_2, c1_2, ...]
// We store as flat array indexed by (4*(j-1) + i) for Fortran (i,j)
const COEF: [f64; 212] = [
// J=1..10 (from DATA ((COEF(I,J),I=1,4),J=1,10))
8.734e-01, -1.545e00, -1.093e00, 5.918e-01, // j=1: 1^1S
9.771e-01, -1.567e00, -4.739e-01, -1.302e-01, // j=2: 2^3S
1.174e00, -1.638e00, -2.831e-01, -3.281e-02, // j=3: 2^1S
1.324e00, -1.692e00, -2.916e-01, 9.027e-02, // j=4: 2^3P
1.445e00, -1.761e00, -1.902e-01, 4.401e-02, // j=5: 2^1P
1.546e00, -1.817e00, -1.278e-01, 2.293e-02, // j=6: 3^3S
1.635e00, -1.864e00, -8.252e-02, 9.854e-03, // j=7: 3^1S
1.712e00, -1.903e00, -5.206e-02, 2.892e-03, // j=8: 3^3P
1.782e00, -1.936e00, -2.952e-02, -1.405e-03, // j=9: 3^3D
1.845e00, -1.964e00, -1.152e-02, -4.487e-03, // j=10: 3^1D
// J=11..19 (from DATA ((COEF(I,J),I=1,4),J=11,19))
7.377e-01, -9.327e-01, -1.466e00, 6.891e-01, // j=11: 3^1P
9.031e-01, -1.157e00, -7.151e-01, 1.832e-01, // j=12: 4^3S
1.031e00, -1.313e00, -4.517e-01, 9.207e-02, // j=13: 4^1S
1.135e00, -1.441e00, -2.724e-01, 3.105e-02, // j=14: 4^3P
1.225e00, -1.536e00, -1.725e-01, 7.191e-03, // j=15
1.302e00, -1.602e00, -1.300e-01, 7.345e-03, // j=16
1.372e00, -1.664e00, -8.204e-02, -1.643e-03, // j=17
1.434e00, -1.715e00, -4.646e-02, -7.456e-03, // j=18
1.491e00, -1.760e00, -1.838e-02, -1.152e-02, // j=19
// J=20..27 (from DATA ((COEF(I,J),I=1,4),J=20,27))
1.258e00, -3.442e00, -4.731e-01, -9.522e-02, // j=20: 2^3Po (trip S, l=0, n=2 start)
1.553e00, -2.781e00, -6.841e-01, -4.083e-03, // j=21
1.727e00, -2.494e00, -5.785e-01, -6.015e-02, // j=22
1.853e00, -2.347e00, -4.611e-01, -9.615e-02, // j=23
1.955e00, -2.273e00, -3.457e-01, -1.245e-01, // j=24
2.041e00, -2.226e00, -2.669e-01, -1.344e-01, // j=25
2.115e00, -2.200e00, -1.999e-01, -1.410e-01, // j=26
2.182e00, -2.188e00, -1.405e-01, -1.460e-01, // j=27
// J=28..35 (from DATA ((COEF(I,J),I=1,4),J=28,35))
1.267e00, -3.417e00, -5.038e-01, -1.797e-02, // j=28: 2^1Po (sing S, l=0, n=2 start)
1.565e00, -2.781e00, -6.497e-01, -5.979e-03, // j=29
1.741e00, -2.479e00, -6.099e-01, -2.227e-02, // j=30
1.870e00, -2.336e00, -4.899e-01, -6.616e-02, // j=31
1.973e00, -2.253e00, -3.972e-01, -8.729e-02, // j=32
2.061e00, -2.212e00, -3.072e-01, -1.060e-01, // j=33
2.137e00, -2.189e00, -2.352e-01, -1.171e-01, // j=34
2.205e00, -2.186e00, -1.621e-01, -1.296e-01, // j=35
// J=36..44 (from DATA ((COEF(I,J),I=1,4),J=36,44))
1.129e00, -3.149e00, -1.910e-01, -5.244e-01, // j=36: 2^3Po l=1 (trip P start)
1.431e00, -2.511e00, -3.710e-01, -1.933e-01, // j=37
1.620e00, -2.303e00, -3.045e-01, -1.391e-01, // j=38
1.763e00, -2.235e00, -1.829e-01, -1.491e-01, // j=39
1.879e00, -2.215e00, -9.003e-02, -1.537e-01, // j=40
1.978e00, -2.213e00, -2.066e-02, -1.541e-01, // j=41
2.064e00, -2.220e00, 3.258e-02, -1.527e-01, // j=42
2.140e00, -2.225e00, 6.311e-02, -1.455e-01, // j=43
2.208e00, -2.229e00, 7.977e-02, -1.357e-01, // j=44
// J=45..53 (from DATA ((COEF(I,J),I=1,4),J=45,53))
1.204e00, -2.809e00, -3.094e-01, 1.100e-01, // j=45: 2^1Po l=1 (sing P start)
1.455e00, -2.254e00, -4.795e-01, 6.872e-02, // j=46
1.619e00, -2.109e00, -3.357e-01, -2.532e-02, // j=47
1.747e00, -2.065e00, -2.317e-01, -5.224e-02, // j=48
1.853e00, -2.058e00, -1.517e-01, -6.647e-02, // j=49
1.943e00, -2.055e00, -1.158e-01, -6.081e-02, // j=50
2.023e00, -2.070e00, -6.470e-02, -6.800e-02, // j=51
2.095e00, -2.088e00, -2.357e-02, -7.250e-02, // j=52
2.160e00, -2.107e00, 1.065e-02, -7.542e-02, // j=53
];
// L > 2: 使用类氢公式
if l > 2 {
let gn = 2.0 * (n * n) as f64;
return PHOT0 / freq / freq / freq / (n as f64).powi(5) * (2 * l + 1) as f64 * s as f64 / gn;
}
// 对于 L <= 2，使用近似值
// 完整实现需要所有 53 组系数
let fl = (freq / FRH).log10();
let idx = ((n - 1).max(0) as usize).min(52);
let x = fl - FL0[idx];
// SELECT beginning of coefficients
// Fortran: SS=(S+1)/2, LL=L+1, NSL0=N0(LL,SS), I=IST(LL,SS)+N-NSL0
let ss = ((s + 1) / 2 - 1) as usize; // 0-indexed: S=1→0, S=3→1
let ll = (l as usize); // 0-indexed: L=0→0, L=1→1, L=2→2
let nsl0 = N0[ll][ss];
let i = IST[ll][ss] + (n - nsl0) as usize - 1; // 0-indexed
if x >= -0.001 {
TENM18 * (TENLG * (-2.0 + 0.5 * x)).exp()
if i >= 53 {
// Beyond tabulated data — fallback to hydrogenic
let gn = 2.0 * (n * n) as f64;
return PHOT0 / freq / freq / freq / (n as f64).powi(5) * (2 * l + 1) as f64 * s as f64 / gn;
}
// EVALUATE cross section
let fl = (freq / FRH).log10();
let x = fl - FL0[i];
if x < -0.001 {
return 0.0;
}
if x < XFITM[i] {
// Cubic fit: P = C0 + C1*X + C2*X^2 + C3*X^3
let c0 = COEF[4 * i];
let c1 = COEF[4 * i + 1];
let c2 = COEF[4 * i + 2];
let c3 = COEF[4 * i + 3];
let p = c0 + x * (c1 + x * (c2 + x * c3));
TENM18 * (TENLG * p).exp()
} else {
0.0
// Linear fit: A + B*X
TENM18 * (TENLG * (A[i] + B[i] * x)).exp()
}
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
#[test]
fn test_hephot_l_gt_2() {
@ -74,14 +213,52 @@ mod tests {
#[test]
fn test_hephot_low_freq() {
// 低频率返回 0
// 低频率返回 0 (below threshold)
let result = hephot(1, 0, 1, 1e10);
assert_relative_eq!(result, 0.0, epsilon = 1e-20);
assert_eq!(result, 0.0);
}
#[test]
fn test_hephot_valid() {
let result = hephot(1, 0, 1, 1e15);
assert!(result >= 0.0);
fn test_hephot_he1_ground_state() {
// He I 1^1S (s=1, l=0, n=1) at freq above threshold
// Threshold: FL0[0] = 0.2521 → freq = FRH * 10^0.2521 ≈ 5.93e15
let result = hephot(1, 0, 1, 6.0e15);
assert!(result > 0.0);
// Should be around ~1e-18 cm^2 (order of magnitude)
assert!(result < 1e-15 && result > 1e-22, "result = {result:e}");
}
#[test]
fn test_hephot_he1_2triplet_s() {
// He I 2^3S (s=3, l=0, n=2) above threshold
// Threshold: FL0[10] = -0.4555 → freq = FRH * 10^-0.4555 ≈ 1.16e15
let result = hephot(3, 0, 2, 2.0e15);
assert!(result > 0.0);
}
#[test]
fn test_hephot_he1_2triplet_p() {
// He I 2^3P (s=3, l=1, n=2) above threshold
// Index: IST[1][0] = 36, N0[1][0] = 2, n=2 → i = 36+0 = 35 (0-indexed)
let result = hephot(3, 1, 2, 2.0e15);
assert!(result > 0.0);
}
#[test]
fn test_hephot_he1_2singlet_p() {
// He I 2^1P (s=1, l=1, n=2) above threshold
// Index: IST[1][1] = 45, N0[1][1] = 2, n=2 → i = 45+0 = 44 (0-indexed)
let result = hephot(1, 1, 2, 2.0e15);
assert!(result > 0.0);
}
#[test]
fn test_hephot_hydrogenic_high_l() {
// L=3, n=4, S=3 — should use hydrogenic formula
let result = hephot(3, 3, 4, 1e15);
assert!(result > 0.0);
// Hydrogenic: PHOT0 / freq^3 / n^5 * (2*L+1) * S / (2*n^2)
let expected = 2.815e29 / 1e15 / 1e15 / 1e15 / 4.0_f64.powi(5) * 7.0 * 3.0 / (2.0 * 16.0);
assert!((result - expected).abs() / expected < 1e-10);
}
}

8
src/tlusty/math/hydrogen/hesolv.rs
View File

@ -107,7 +107,9 @@ pub struct HesolvAtomicParams<'a> {
pub abund: &'a [f64],
/// 能级量子数 [能级]
pub nquant: &'a [i32],
/// 原子序数 [能级]
/// 所属离子索引 [能级]
pub iel: &'a [i32],
/// 原子序数 [离子]
pub iz: &'a [i32],
/// 能级权重标志 [能级]
pub ifwop: &'a [i32],
@ -508,6 +510,7 @@ pub fn hesolv_pure(params: &HesolvParams) -> HesolvOutput {
params.atomic.ifwop,
params.config.nlevel,
params.atomic.nquant,
params.atomic.iel,
params.atomic.iz,
params.config.io_ptab,
params.config.lte,
@ -660,6 +663,7 @@ mod tests {
nnext: &[],
abund: &[],
nquant: &[],
iel: &[],
iz: &[],
ifwop: &[],
xi2: &[],
@ -727,6 +731,7 @@ mod tests {
nnext: &[],
abund: &[],
nquant: &[],
iel: &[],
iz: &[],
ifwop: &[],
xi2: &[],
@ -805,6 +810,7 @@ mod tests {
nnext: &[],
abund: &[],
nquant: &[],
iel: &[],
iz: &[],
ifwop: &[],
xi2: &[],

10
src/tlusty/math/hydrogen/sgmer.rs
View File

@ -81,7 +81,7 @@ const EHB: f64 = 157802.77355;
/// ```
pub fn sgmer0(
atomic: &AtomicData,
model: &ModelState,
temp1: &[f64],
wmcomp: &WmComp,
mrgpar: &mut MrgPar,
invint: &InvInt,
@ -130,7 +130,7 @@ pub fn sgmer0(
// 遍历所有深度
for id in 0..nd {
let ex = EHB * ch * model.modpar.temp1[id];
let ex = EHB * ch * temp1[id];
// 计算频率和积分项
for i in (ii0 - 1)..NLMX {
@ -282,9 +282,9 @@ pub fn sgmerd(
mod tests {
use super::*;
fn create_test_data() -> (AtomicData, ModelState, WmComp, MrgPar, InvInt) {
fn create_test_data() -> (AtomicData, Vec<f64>, WmComp, MrgPar, InvInt) {
let mut atomic = AtomicData::new();
let model = ModelState::new();
let temp1 = vec![0.0001; 5];
let wmcomp = WmComp::default();
let mrgpar = MrgPar::default();
let invint = InvInt::default();
@ -293,7 +293,7 @@ mod tests {
atomic.levpar.iel[0] = 0; // 属于第一个离子
atomic.ionpar.iz[0] = 1; // 氢 Z=1
(atomic, model, wmcomp, mrgpar, invint)
(atomic, temp1, wmcomp, mrgpar, invint)
}
#[test]

76
src/tlusty/math/io/output.rs
View File

@ -241,30 +241,32 @@ fn write_depth_line_with_popul<W: Write>(
id: usize,
nlevel: usize,
) -> Result<()> {
let mut values = vec![temp, elec, dens];
if let Some(t) = totn {
values.push(t);
}
// Fortran: WRITE(7,503) TEMP(ID),ELEC(ID),DENS(ID),(POPUL(J,ID),J=1,NLEVEL)
// FORMAT(1P5E15.6) - 5 values per line, flowing continuously
// So first line has temp, elec, dens, pop[0], pop[1]
// then subsequent lines have 5 pop values each
// 写入前 5 个值
let count = values.len().min(5);
let mut line = String::new();
for i in 0..count {
line.push_str(&format_exp_fortran(values[i], 15, 6, false));
}
writer.write_raw(&line)?;
writer.write_newline()?;
let mut col = 0;
// 写入占据数（每行 5 个）
let mut pop_line = String::new();
let mut pop_count = 5 - count;
// Write temp, elec, dens
for &val in &[temp, elec, dens] {
line.push_str(&format_exp_fortran(val, 15, 6, false));
col += 1;
}
if let Some(t) = totn {
line.push_str(&format_exp_fortran(t, 15, 6, false));
col += 1;
}
// Write populations, continuing from where temp/elec/dens left off
for j in 0..nlevel {
pop_line.push_str(&format_exp_fortran(popul[j][id], 15, 6, false));
pop_count += 1;
if pop_count % 5 == 0 || j == nlevel - 1 {
writer.write_raw(&pop_line)?;
line.push_str(&format_exp_fortran(popul[j][id], 15, 6, false));
col += 1;
if col % 5 == 0 || j == nlevel - 1 {
writer.write_raw(&line)?;
writer.write_newline()?;
pop_line.clear();
line.clear();
}
}
@ -284,28 +286,28 @@ fn write_depth_line_with_popul_disk<W: Write>(
id: usize,
nlevel: usize,
) -> Result<()> {
let mut values = vec![temp, elec, dens];
if let Some(t) = totn {
values.push(t);
}
values.push(zd);
// 写入基本值
// Same flow as planar: values flow continuously, 5 per line
let mut line = String::new();
for &val in &values {
line.push_str(&format_exp_fortran(val, 15, 6, false));
}
writer.write_raw(&line)?;
writer.write_newline()?;
let mut col = 0;
for &val in &[temp, elec, dens] {
line.push_str(&format_exp_fortran(val, 15, 6, false));
col += 1;
}
if let Some(t) = totn {
line.push_str(&format_exp_fortran(t, 15, 6, false));
col += 1;
}
line.push_str(&format_exp_fortran(zd, 15, 6, false));
col += 1;
// 写入占据数（每行 5 个）
let mut pop_line = String::new();
for j in 0..nlevel {
pop_line.push_str(&format_exp_fortran(popul[j][id], 15, 6, false));
if (j + 1) % 5 == 0 || j == nlevel - 1 {
writer.write_raw(&pop_line)?;
line.push_str(&format_exp_fortran(popul[j][id], 15, 6, false));
col += 1;
if col % 5 == 0 || j == nlevel - 1 {
writer.write_raw(&line)?;
writer.write_newline()?;
pop_line.clear();
line.clear();
}
}

272
src/tlusty/math/population/lte_saha.rs Normal file
View File

@ -0,0 +1,272 @@
//! LTE Saha-Boltzmann 种群计算。
//!
//! 从 `main.rs` 提取为共享模块，供 RESOLV 在 Lucy 温度修正后调用。
use crate::tlusty::state::constants::{BOLK, EH, HK, HMASS};
use crate::tlusty::state::model::ModelState;
// ============================================================================
// 物理常量
// ============================================================================
/// Saha 常数: 2 × (2πmₑk/h²)^(3/2) = 2.0706e-16
pub const CCON: f64 = 2.0706e-16;
/// Planck 常数 (erg·s) = HK × BOLK
pub const H_PLANCK_LTE: f64 = HK * BOLK;
// He I 激发频率 (Hz)
pub const HE1_FREQ: [f64; 14] = [
5.94503520e15, 1.15267210e15, 9.60145430e14, 8.75933720e14,
8.14536220e14, 4.51727350e14, 4.02921120e14, 3.81935640e14,
3.65836790e14, 3.65746870e14, 3.62599020e14, 2.40043860e14,
2.20797190e14, 2.12492940e14,
];
// He I 统计权重
pub const HE1_G: [f64; 14] = [
1.0, 3.0, 1.0, 9.0, 3.0, 3.0, 1.0, 9.0, 15.0, 5.0, 3.0, 3.0, 1.0, 9.0,
];
// He I 统计权重 (用于 Saha 公式)
pub const HE1_WEIGHT: [f64; 14] = [
1.0, 3.0, 1.0, 9.0, 3.0, 3.0, 1.0, 9.0, 15.0, 5.0, 3.0, 3.0, 1.0, 9.0,
];
// 丰度
pub const X_H: f64 = 0.70;
pub const X_HE: f64 = 0.28;
pub const Y_H: f64 = X_H / 1.008;
pub const Y_HE: f64 = X_HE / 4.003;
pub const YTOT_LTE: f64 = Y_H + Y_HE;
pub const ABUND_H: f64 = Y_H / YTOT_LTE;
pub const ABUND_HE: f64 = Y_HE / YTOT_LTE;
/// 平均分子量 (g/particle)
pub const WMM_GREY: f64 = HMASS / (X_H + X_HE / 4.0 / 4.0);
/// 计算 HHe 模型的 LTE Saha-Boltzmann 种群。
///
/// 使用 Saha 方程和 Boltzmann 分布计算 H I/II, He I/II/III 的 LTE 种群。
/// 种群存入 `model.levpop.popul[level][depth]`。
///
/// # 参数
/// - `model`: 模型状态 (temp, elec, dens → popul)
/// - `nd`: 深度点数
/// - `nlevel`: 能级数 (HHe = 39)
/// - `wmm`: 平均分子量 (g/particle)，来自 compute_abundance_params
/// - `abund_h`: H 相对丰度 (已归一化，= abnd_h / ytot)
/// - `abund_he`: He 相对丰度 (已归一化，= abnd_he / ytot)
pub fn compute_lte_populations_hhe(
model: &mut ModelState,
nd: usize,
nlevel: usize,
wmm: f64,
abund_h: f64,
abund_he: f64,
) {
assert!(nlevel >= 39, "nlevel must be >= 39 for HHe model");
for id in 0..nd {
let t = model.modpar.temp[id];
let ne = model.modpar.elec[id];
let rho = model.modpar.dens[id];
if t <= 0.0 || ne <= 0.0 || rho <= 0.0 {
continue;
}
let tk = BOLK * t; // kT in erg
let con = CCON / (t * t.sqrt()); // Saha 常数 / T^(3/2)
// ====== Hydrogen (Ion 0: H I, Ion 1: H II) ======
let h_next_g = 1.0; // H II ground g
let mut h_sbf = [0.0f64; 9];
for n in 1usize..=9 {
let nn = n as f64;
let enion_i = EH / (nn * nn);
let g_i = 2.0 * nn * nn;
h_sbf[n - 1] = con * g_i / h_next_g * (enion_i / tk).exp();
}
// Fortran state.f: hpop=dens(id)/wmm(id)/ytot(id)
// With IREFA=1, ABNREF=1.0, abundances are renormalized
// so n_h_total = rho/wmm * (abnd_h/ytot) = rho/wmm * abund_h
let n_h_total = rho / wmm * abund_h;
let mut sum_ne_sbf = 0.0f64;
for i in 0..9 {
sum_ne_sbf += ne * h_sbf[i];
}
let h_ii_pop = n_h_total / (1.0 + sum_ne_sbf);
for n in 0..9 {
model.levpop.popul[n][id] = ne * h_sbf[n] * h_ii_pop;
}
model.levpop.popul[9][id] = h_ii_pop;
// ====== Helium ======
let he3_g = 1.0;
// He II: hydrogenic Z=2, 14 levels
let mut he2_sbf = [0.0f64; 14];
for n in 1usize..=14 {
let nn = n as f64;
let enion_i = EH * 4.0 / (nn * nn);
let g_i = 2.0 * nn * nn;
he2_sbf[n - 1] = con * g_i / he3_g * (enion_i / tk).exp();
}
let mut sum_ne_he2_sbf = 0.0f64;
for i in 0..14 {
sum_ne_he2_sbf += ne * he2_sbf[i];
}
// He I: 14 levels
let he2_ground_g = 2.0;
let mut he1_sbf = [0.0f64; 14];
for i in 0..14 {
let enion_i = HE1_FREQ[i] * H_PLANCK_LTE;
let g_i = HE1_G[i];
he1_sbf[i] = con * g_i / he2_ground_g * (enion_i / tk).exp();
}
let n_he_total = rho / wmm * abund_he;
let he2_ground_sbf = he2_sbf[0];
let mut sum_ne_he1_sbf = 0.0f64;
for i in 0..14 {
sum_ne_he1_sbf += he1_sbf[i];
}
let denom = 1.0 + sum_ne_he2_sbf + ne * he2_ground_sbf * sum_ne_he1_sbf;
let he3_pop = n_he_total / denom;
let mut he2_pops = [0.0f64; 14];
for n in 0..14 {
he2_pops[n] = ne * he2_sbf[n] * he3_pop;
}
let he2_ground_pop = he2_pops[0];
let mut he1_pops = [0.0f64; 14];
for i in 0..14 {
he1_pops[i] = ne * he1_sbf[i] * he2_ground_pop;
}
for i in 0..14 {
model.levpop.popul[10 + i][id] = he1_pops[i];
}
for n in 0..14 {
model.levpop.popul[24 + n][id] = he2_pops[n];
}
model.levpop.popul[38][id] = he3_pop;
// Note: ELEC and TOTN are NOT updated here. In Fortran, INILAM sets populations
// but ELEC/TOTN are independently determined by ELDENS in the Lambda loop.
}
}
/// Compute LTE populations for a single depth point.
///
/// Returns popul[39] for the HHe model (H I/II + He I/II/III).
/// Used for finite-difference temperature derivatives without modifying the model.
pub fn compute_lte_populations_single(t: f64, ne: f64, rho: f64, wmm: f64, abund_h: f64, abund_he: f64) -> [f64; 39] {
let mut popul = [0.0f64; 39];
if t <= 0.0 || ne <= 0.0 || rho <= 0.0 {
return popul;
}
let tk = BOLK * t;
let con = CCON / (t * t.sqrt());
// ====== Hydrogen ======
let h_next_g = 1.0;
let mut h_sbf = [0.0f64; 9];
for n in 1usize..=9 {
let nn = n as f64;
let enion_i = EH / (nn * nn);
let g_i = 2.0 * nn * nn;
h_sbf[n - 1] = con * g_i / h_next_g * (enion_i / tk).exp();
}
let n_h_total = rho / wmm * abund_h;
let sum_ne_sbf: f64 = (0..9).map(|i| ne * h_sbf[i]).sum();
let h_ii_pop = n_h_total / (1.0 + sum_ne_sbf);
for n in 0..9 {
popul[n] = ne * h_sbf[n] * h_ii_pop;
}
popul[9] = h_ii_pop;
// ====== Helium ======
let he3_g = 1.0;
let mut he2_sbf = [0.0f64; 14];
for n in 1usize..=14 {
let nn = n as f64;
let enion_i = EH * 4.0 / (nn * nn);
let g_i = 2.0 * nn * nn;
he2_sbf[n - 1] = con * g_i / he3_g * (enion_i / tk).exp();
}
let sum_ne_he2_sbf: f64 = (0..14).map(|i| ne * he2_sbf[i]).sum();
let he2_ground_g = 2.0;
let mut he1_sbf = [0.0f64; 14];
for i in 0..14 {
let enion_i = HE1_FREQ[i] * H_PLANCK_LTE;
let g_i = HE1_G[i];
he1_sbf[i] = con * g_i / he2_ground_g * (enion_i / tk).exp();
}
let n_he_total = rho / wmm * abund_he;
let he2_ground_sbf = he2_sbf[0];
let sum_he1_sbf: f64 = (0..14).map(|i| he1_sbf[i]).sum();
let denom = 1.0 + sum_ne_he2_sbf + ne * he2_ground_sbf * sum_he1_sbf;
let he3_pop = n_he_total / denom;
let mut he2_pops = [0.0f64; 14];
for n in 0..14 {
he2_pops[n] = ne * he2_sbf[n] * he3_pop;
}
let he2_ground_pop = he2_pops[0];
for i in 0..14 {
popul[10 + i] = ne * he1_sbf[i] * he2_ground_pop;
}
for n in 0..14 {
popul[24 + n] = he2_pops[n];
}
popul[38] = he3_pop;
popul
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_compute_lte_populations_basic() {
let mut model = ModelState::new();
let nd = 5;
let nlevel = 39;
for id in 0..nd {
model.modpar.temp[id] = 10000.0 + id as f64 * 1000.0;
model.modpar.elec[id] = 1e12;
model.modpar.dens[id] = 1e14;
}
compute_lte_populations_hhe(&mut model, nd, nlevel, WMM_GREY, ABUND_H, ABUND_HE);
// Check that populations are positive
for id in 0..nd {
for lev in 0..39 {
assert!(model.levpop.popul[lev][id] >= 0.0,
"Population negative at level {} depth {}", lev, id);
}
}
}
}

2
src/tlusty/math/population/mod.rs
View File

@ -6,6 +6,7 @@ mod bpope;
mod bpopf;
mod bpopt;
mod butler;
pub mod lte_saha;
mod newpop;
pub use bpop::*;
@ -14,4 +15,5 @@ pub use bpope::*;
pub use bpopf::*;
pub use bpopt::*;
pub use butler::*;
pub use lte_saha::*;
pub use newpop::*;

237
src/tlusty/math/radiative/feautrier.rs Normal file
View File

@ -0,0 +1,237 @@
//! Feautrier formal solution for the radiative transfer equation.
//!
//! Scalar Feautrier with f=1/3 for the second pass, with SS0 coupling.
//! Matches Fortran RTEFR1 second pass for ISPLIN=0, IBC=3, IDISK=0.
use crate::tlusty::state::constants::{BN, HALF, HK, TWO, UN};
/// Gauss-Legendre NMU=4 quadrature for QQ0
const AMU: [f64; 4] = [
0.06943184420297371,
0.33000947820757187,
0.66999052179242813,
0.93056815579702637,
];
const WTMU: [f64; 4] = [
0.17392742256872692,
0.32607257743127314,
0.32607257743127314,
0.17392742256872692,
];
/// Result of the Feautrier formal solution.
pub struct FeautrierResult {
/// Mean intensity Jν at each depth point [nd]
pub rad1: Vec<f64>,
/// Eddington factor K/J at each depth point [nd]
pub fak1: Vec<f64>,
/// Lambda operator diagonal (ALRH) [nd]
pub alrh: Vec<f64>,
}
/// Solve the Feautrier transfer equation for one frequency.
///
/// Scalar Feautrier with constant Eddington factor f=1/3.
/// Uses ISPLIN=0, IBC=3 (diffusion lower BC).
pub fn feautrier_solve(
nd: usize,
abso: &[f64],
_true_abs: &[f64],
scat: &[f64],
emis: &[f64],
dm: &[f64],
_dens: &[f64],
temp: &[f64],
freq: f64,
tempbd: f64,
ibc: i32,
) -> FeautrierResult {
let third = 1.0_f64 / 3.0;
let twothr = 2.0_f64 / 3.0;
// Compute optical depth increments DT
// Fortran TDPINI: DELDMZ(ID) = HALF*(DM(ID+1)-DM(ID))
// Fortran OPACTD: ABSOT(ID) = ABSO1(ID)/DENS(ID) for IMODF=0
let mut dt = vec![0.0; nd];
for id in 0..(nd - 1) {
let deldmz = HALF * (dm[id + 1] - dm[id]);
let absot_id = if _dens[id] > 0.0 { abso[id] / _dens[id] } else { abso[id] };
let absot_id1 = if _dens[id + 1] > 0.0 { abso[id + 1] / _dens[id + 1] } else { abso[id + 1] };
dt[id] = deldmz * (absot_id1 + absot_id);
}
// Total source function S = emis / abso
let mut st0 = vec![0.0; nd];
for id in 0..nd {
st0[id] = if abso[id] > 0.0 { emis[id] / abso[id] } else { 0.0 };
}
// Scattering parameter SS0 = -scat/abso
let ss0: Vec<f64> = (0..nd)
.map(|id| if abso[id] > 0.0 { -scat[id] / abso[id] } else { 0.0 })
.collect();
// Surface optical depth correction
let taumin = abso[0] * dm[0] / 2.0;
// Compute QQ0 for surface BC
let mut qq0 = 0.0;
for k in 0..4 {
let tamm = taumin / AMU[k];
let p0 = 1.0 - (-tamm).exp();
qq0 += p0 * AMU[k] * WTMU[k];
}
// Lower boundary: Planck function at deepest depth
let fr15 = freq * 1e-15;
let bnu = BN * fr15 * fr15 * fr15;
let t_bot = if tempbd > 0.0 { tempbd } else { temp[nd - 1] };
let t_bot_m1 = if tempbd > 0.0 { tempbd } else { temp[nd - 2] };
let rrdil = 1.0;
let pland = bnu / ((HK * freq / t_bot).exp() - UN) * rrdil;
let dplan_0 = bnu / ((HK * freq / t_bot_m1).exp() - UN) * rrdil;
let dplan = if dt[nd - 2] > 0.0 {
(pland - dplan_0) / dt[nd - 2]
} else {
0.0
};
// Constant Eddington factor
let fkk = vec![third; nd];
let fh0 = 1.0_f64 / 3.0_f64.sqrt();
let fhd = 1.0_f64 / 3.0_f64.sqrt();
// ============ Scalar tridiagonal solve ============
// ISPLIN=0: A=0, C=0, B=1; with ss0 coupling
// Fortran RTEFR1 second pass (lines 462-583)
// Upper boundary (id=0)
let b_bc = dt[0] * HALF;
let c_bc = 0.0;
let bq = 1.0 / (b_bc + qq0);
let cq = c_bc * bq;
let mut bbb = vec![0.0; nd];
let mut ccc = vec![0.0; nd];
let mut aaa = vec![0.0; nd];
let mut zzz = vec![0.0; nd];
let mut aanu = vec![0.0; nd];
let mut ddd = vec![0.0; nd];
bbb[0] = (fkk[0] / dt[0] + fh0 + b_bc) * bq + ss0[0];
ccc[0] = fkk[1] / dt[0] * bq - cq * (UN + ss0[1]);
let vll = st0[0] + cq * st0[1];
zzz[0] = 1.0 / bbb[0];
aanu[0] = vll * zzz[0];
ddd[0] = ccc[0] * zzz[0];
// Interior points
let mut dtp1_s = dt[0];
for id in 1..(nd - 1) {
let dtm1 = dtp1_s;
dtp1_s = dt[id];
let dt0 = 2.0 / (dtp1_s + dtm1);
let al = dtm1.recip() * dt0;
let ga = dtp1_s.recip() * dt0;
// ISPLIN=0: A=0, C=0, B=1
aaa[id] = al * fkk[id - 1];
ccc[id] = ga * fkk[id + 1];
bbb[id] = (al + ga) * fkk[id] + UN + ss0[id];
let vll = st0[id];
aanu[id] = vll + aaa[id] * aanu[id - 1];
zzz[id] = 1.0 / (bbb[id] - aaa[id] * ddd[id - 1]);
ddd[id] = ccc[id] * zzz[id];
aanu[id] *= zzz[id];
}
// Lower boundary (id=nd-1)
// IBC=3: Fortran lines 541-546
let id_last = nd - 1;
let b_val = 1.0 / dtp1_s;
let a_val = 2.0 * b_val * b_val;
let bbb_lbc = UN + ss0[id_last] + b_val * 2.0 * fhd + a_val * fkk[id_last];
let aaa_lbc = a_val * fkk[id_last - 1];
let vll_lbc = st0[id_last] + b_val * (pland + twothr * dplan);
zzz[id_last] = 1.0 / (bbb_lbc - aaa_lbc * ddd[id_last - 1]);
let mut rad1 = vec![0.0; nd];
rad1[id_last] = (vll_lbc + aaa_lbc * aanu[id_last - 1]) * zzz[id_last];
// Back-substitution
let mut alrh = vec![0.0; nd];
let mut eee = vec![0.0; nd];
alrh[id_last] = zzz[id_last];
for id in (0..id_last).rev() {
eee[id] = aaa[id] / (bbb[id] - ccc[id] * eee[id + 1]);
rad1[id] = aanu[id] + ddd[id] * rad1[id + 1];
alrh[id] = zzz[id] / (1.0 - ddd[id] * eee[id + 1]);
}
FeautrierResult {
rad1,
fak1: fkk,
alrh,
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_feautrier_thermal_equilibrium() {
let nd = 20;
let abso = vec![1.0; nd];
let true_abs = vec![0.5; nd];
let scat = vec![0.5; nd];
let dens = vec![1.0; nd];
let bnu = 1.0;
let emis: Vec<f64> = true_abs.iter().map(|&a| a * bnu).collect();
let dm: Vec<f64> = (0..nd).map(|i| (i + 1) as f64 * 0.1).collect();
let temp = vec![10000.0; nd];
let result = feautrier_solve(
nd, &abso, &true_abs, &scat, &emis, &dm, &dens, &temp,
1e15, 0.0, 3,
);
let s_val = 0.5;
for id in (nd - 5)..nd {
let rel_err = (result.rad1[id] - s_val).abs() / s_val;
assert!(rel_err < 0.05, "id={}: J={:.6}, S={:.6}, rel_err={:.4}",
id, result.rad1[id], s_val, rel_err);
}
assert!(result.rad1[0] < s_val, "Surface J should be less than S");
}
#[test]
fn test_feautrier_positive() {
let nd = 30;
let abso: Vec<f64> = (0..nd).map(|i| 0.1 + i as f64 * 0.1).collect();
let true_abs: Vec<f64> = abso.iter().map(|x| x * 0.7).collect();
let scat: Vec<f64> = abso.iter().map(|x| x * 0.3).collect();
let dens = vec![1.0; nd];
let temp: Vec<f64> = (0..nd).map(|i| 5000.0 + i as f64 * 500.0).collect();
let h_over_c2 = 2.0 * 6.6262e-27 / (2.9979e10_f64.powi(2));
let freq = 1e15;
let emis: Vec<f64> = (0..nd).map(|i| {
let x = (HK / temp[i] * freq).min(150.0);
true_abs[i] * h_over_c2 * freq.powi(3) / (x.exp() - 1.0)
}).collect();
let dm: Vec<f64> = (0..nd).map(|i| (i + 1) as f64 * 0.01).collect();
let result = feautrier_solve(
nd, &abso, &true_abs, &scat, &emis, &dm, &dens, &temp,
freq, 0.0, 3,
);
for id in 0..nd {
assert!(result.rad1[id] >= 0.0, "id={}: J = {:.6e} is negative", id, result.rad1[id]);
assert!(!result.rad1[id].is_nan(), "id={}: J is NaN", id);
}
}
}

2
src/tlusty/math/radiative/mod.rs
View File

@ -1,6 +1,7 @@
//! radiative module
mod coolrt;
mod feautrier;
mod radpre;
mod radtot;
mod rte_sc;
@ -20,6 +21,7 @@ mod trmder;
mod trmdrt;
pub use coolrt::*;
pub use feautrier::*;
pub use radpre::*;
pub use radtot::*;
pub use rte_sc::*;

1228
src/tlusty/math/solvers/matgen_lte.rs Normal file
View File

File diff suppressed because it is too large Load Diff

2
src/tlusty/math/solvers/mod.rs
View File

@ -8,6 +8,7 @@ mod laguer;
mod lineqs;
mod matcon;
mod matgen;
mod matgen_lte;
mod matinv;
mod minv3;
mod psolve;
@ -32,6 +33,7 @@ pub use laguer::*;
pub use lineqs::*;
pub use matcon::*;
pub use matgen::*;
pub use matgen_lte::*;
pub use matinv::*;
pub use minv3::*;
pub use psolve::*;

23
src/tlusty/math/temperature/lucy.rs
View File

@ -402,12 +402,23 @@ pub fn lucy_pure(
state.deltat[0] = state.dt1[0] + state.dt2[0];
// 内部点
// Fortran: DO ID=2,ND
// XX=XX+DELDM(ID)*(ABSH(ID)*DENS1(ID)*DELH(ID)+
// ABSH(ID+1)*DENS1(ID+1)*DELH(ID+1))
// Mapping: Fortran ID → Rust id+1 (1-based → 0-based)
// DELDM(ID) → deldm[id], ABSH(ID) → absh[id], ABSH(ID+1) → absh[id+1]
// When id+1 == nd, Fortran reads ABSH(ND+1)=0 (static local array),
// so we return 0 for out-of-bounds access.
for id in 1..nd {
// Fortran: TP3 = TEMP1(ID)**3, where TEMP1 = 1/T
let tp3 = model.temp[id].powi(-3);
xx += model.deldm[id - 1]
* (state.absh[id - 1] * model.dens1[id - 1] * state.delh[id - 1]
+ state.absh[id] * model.dens1[id] * state.delh[id]);
let absh_curr = state.absh[id] * model.dens1[id] * state.delh[id];
let absh_next = if id + 1 < nd {
state.absh[id + 1] * model.dens1[id + 1] * state.delh[id + 1]
} else {
0.0 // Fortran: local arrays dimensioned MDEPTH, ABSH(ND+1)=0
};
xx += model.deldm[id] * (absh_curr + absh_next);
state.dt1[id] = state.heat[id] / 16.0 / SIG4P * tp3 / state.absp[id];
state.dt2[id] = state.absz[id] / state.eddf[id] * xx / 16.0 / SIG4P * tp3 / state.absp[id];
state.deltat[id] = state.dt1[id] + state.dt2[id];
@ -607,6 +618,12 @@ pub struct Rad1PointData {
pub rad1: Vec<f64>,
/// Eddington 因子 [nd]
pub fak1: Vec<f64>,
/// Lambda 算子对角线 [nd] (ALI1)
pub ali1: Vec<f64>,
/// Lambda 算子次对角线 [nd] (ALIM1)
pub alim1: Vec<f64>,
/// Lambda 算子超对角线 [nd] (ALIP1)
pub alip1: Vec<f64>,
}
// ============================================================================

2
src/tlusty/math/utils/mod.rs
View File

@ -18,7 +18,7 @@ mod irc;
mod newdm;
mod newdmt;
mod pgset;
mod sabolf;
pub mod sabolf;
mod setdrt;
mod state;
mod switch;

3
src/tlusty/math/utils/sabolf.rs
View File

@ -149,6 +149,9 @@ pub fn sabolf_pure(params: &mut SabolfParams) -> SabolfOutput {
let nfirst = params.nfirst[ion_idx] as usize;
let nlast = params.nlast[ion_idx] as usize;
if nlast == 0 {
continue;
}
// nlst = nlast - 1 (Fortran 0-indexed offset)
let nlst = nlast - 1;

13
src/tlusty/math/utils/state.rs
View File

@ -708,6 +708,19 @@ pub fn state_pure(params: &StateParams) -> StateOutput {
let abnd = params.abndd.get(i).copied().unwrap_or(0.0);
let x1 = abnd / rs;
// DEBUG: print intermediate values at key depth points
if mode == 1 && (id == 1 || id == 26 || id == 51) {
let iat = i + 1;
eprintln!("DEBUG STATE id={} iat={}: t={:.1} ane={:.5e} xmx={:.3}",
id, iat, t, ane, xmx);
for jj in 0..ion {
eprintln!(" pfstu[{}]={:.6} ffi[{}]={:.6e}",
jj, pfstu[jj], jj, ffi.get(jj).copied().unwrap_or(0.0));
}
eprintln!(" jmax={} rs={:.6e} rq_sum={:.6e} x={:.6} abnd={:.6e}",
jmax, rs, rq_sum, x, abnd);
}
// 累加到总电荷
let irefa = params.irefa;
if i + 1 == irefa {

24
src/tlusty/math/utils/wnstor.rs
View File

@ -5,7 +5,7 @@
//! 将氢能级的占据概率存储到 WNCOM 公共块中，以便后续使用。
use crate::tlusty::math::wn;
use crate::tlusty::state::constants::{MLEVEL, NLMX, UN};
use crate::tlusty::state::constants::{MION, MLEVEL, NLMX, UN};
/// 存储氢能级占据概率。
///
@ -53,6 +53,7 @@ pub fn wnstor(
ifwop: &[i32],
nlevel: usize,
nquant: &[i32],
iel: &[i32],
iz: &[i32],
io_ptab: i32,
lte: bool,
@ -111,7 +112,8 @@ pub fn wnstor(
}
let nq = nquant[ii] as usize;
let iz_val = iz[ii] as f64;
let iel_idx = iel[ii] as usize;
let iz_val = if iel_idx > 0 && iel_idx <= iz.len() { iz[iel_idx - 1] as f64 } else { 0.0 };
if iz_val == 1.0 {
// 氢：使用 WNHINT
@ -146,6 +148,7 @@ mod tests {
Vec<i32>,
Vec<i32>,
Vec<i32>,
Vec<i32>,
) {
let temp = vec![10000.0; ND];
let elec = vec![1e12; ND];
@ -154,19 +157,20 @@ mod tests {
let wop = vec![vec![0.0; ND]; MLEVEL];
let ifwop = vec![1; MLEVEL];
let nquant = vec![1; MLEVEL];
let iz = vec![1; MLEVEL];
let iel = vec![1; MLEVEL];
let iz = vec![1; MION];
(temp, elec, xi2, wnhint, wop, ifwop, nquant, iz)
(temp, elec, xi2, wnhint, wop, ifwop, nquant, iel, iz)
}
#[test]
fn test_wnstor_io_ptab_negative() {
let (temp, elec, xi2, mut wnhint, mut wop, ifwop, nquant, iz) = create_test_data();
let (temp, elec, xi2, mut wnhint, mut wop, ifwop, nquant, iel, iz) = create_test_data();
// io_ptab < 0 应该直接返回
wnstor(
0, &temp, &elec, &xi2, &mut wnhint, &mut wop,
&ifwop, 10, &nquant, &iz, -1, false,
&ifwop, 10, &nquant, &iel, &iz, -1, false,
);
// wnhint 应该保持为 0
@ -175,11 +179,11 @@ mod tests {
#[test]
fn test_wnstor_basic() {
let (temp, elec, xi2, mut wnhint, mut wop, ifwop, nquant, iz) = create_test_data();
let (temp, elec, xi2, mut wnhint, mut wop, ifwop, nquant, iel, iz) = create_test_data();
wnstor(
0, &temp, &elec, &xi2, &mut wnhint, &mut wop,
&ifwop, 10, &nquant, &iz, 0, false,
&ifwop, 10, &nquant, &iel, &iz, 0, false,
);
// 检查 WNHINT 已被计算，值应该在 [0, 1] 范围内
@ -191,11 +195,11 @@ mod tests {
#[test]
fn test_wnstor_lte_flag() {
let (temp, elec, xi2, mut wnhint, mut wop, ifwop, nquant, iz) = create_test_data();
let (temp, elec, xi2, mut wnhint, mut wop, ifwop, nquant, iel, iz) = create_test_data();
wnstor(
0, &temp, &elec, &xi2, &mut wnhint, &mut wop,
&ifwop, 10, &nquant, &iz, 0, true, // lte = true
&ifwop, 10, &nquant, &iel, &iz, 0, true, // lte = true
);
// 当 ifwop > 1 且 lte = true 时，wop 应该是 1.0

305
src/tlusty/runner.rs Normal file
View File

@ -0,0 +1,305 @@
//! TLUSTY 主程序入口。
//!
//! 重构自 TLUSTY `tlusty.f` 主程序。
//!
//! # 算法概述
//!
//! TLUSTY 使用混合 Complete Linearization (CL) 和 Accelerated Lambda Iteration (ALI) 方法
//! 计算非LTE恒星大气模型。
//!
//! 主迭代循环：
//! 1. 初始化 (START)
//! 2. 形式解 (RESOLV)
//! 3. 收敛加速 (ACCEL2)
//! 4. 解线性化方程 (SOLVE/SOLVES/RYBSOL)
//! 5. 重复直到收敛
use std::time::Instant;
// ============================================================================
// 运行配置
// ============================================================================
/// TLUSTY 运行配置。
#[derive(Debug, Clone)]
pub struct TlustyConfig {
/// 最大迭代次数
pub max_iter: usize,
/// 收敛加速模式 (0=关闭, >0=开启)
pub accel_mode: i32,
/// 解法选择 (0=标准SOLVE, 1=Ryan)
pub solver_mode: i32,
/// 是否输出诊断信息
pub verbose: bool,
}
impl Default for TlustyConfig {
fn default() -> Self {
Self {
max_iter: 200,
accel_mode: 1,
solver_mode: 0,
verbose: false,
}
}
}
// ============================================================================
// 运行状态
// ============================================================================
/// TLUSTY 运行状态。
#[derive(Debug, Clone)]
pub struct TlustyState {
/// 当前迭代次数
pub iter: usize,
/// 是否初始化阶段
pub is_init: bool,
/// 是否完成
pub is_finished: bool,
/// 是否最后一次迭代
pub is_final: bool,
/// 系统维度 NN
pub nn: usize,
}
impl Default for TlustyState {
fn default() -> Self {
Self {
iter: 0,
is_init: true,
is_finished: false,
is_final: false,
nn: 0,
}
}
}
// ============================================================================
// 运行结果
// ============================================================================
/// TLUSTY 运行结果。
#[derive(Debug, Clone)]
pub struct TlustyResult {
/// 最终迭代次数
pub total_iterations: usize,
/// 是否收敛
pub converged: bool,
/// 总运行时间（秒）
pub total_time_secs: f64,
/// 各阶段时间统计
pub timing_stats: TimingStats,
}
/// 时间统计。
#[derive(Debug, Clone, Default)]
pub struct TimingStats {
pub formal_solution_time: f64,
pub matrix_solution_time: f64,
pub acceleration_time: f64,
}
// ============================================================================
// 主运行函数
// ============================================================================
/// 运行 TLUSTY 计算。
///
/// # 算法流程
///
/// ```text
/// 1. 初始化
/// - 读取输入参数 (fort.5)
/// - 设置初始温度结构
/// - 读取原子数据
/// - 设置频率网格
///
/// 2. 主迭代循环
/// a) 形式解 (RESOLV)
/// - 解辐射转移方程
/// - 计算辐射场
/// b) 收敛加速 (ACCEL2)
/// - Ng 加速
/// c) 解线性化方程 (SOLVE/SOLVES/RYBSOL)
/// - 统计平衡方程
/// - 能量守恒方程
/// d) 更新大气结构
///
/// 3. 输出最终模型
/// ```
///
/// # Fortran 原始代码
///
/// ```fortran
/// PROGRAM TLUSTY
/// ...
/// CALL START
/// LFIN=.FALSE.
/// IF(NITER.EQ.0) LFIN=.TRUE.
///
/// 10 ITER=ITER+1
/// CALL RESOLV
/// INIT=0
/// IF(LFIN) GO TO 20
///
/// IF(IACC.GT.0) CALL ACCEL2
///
/// IF(IFRYB.EQ.0) THEN
/// IF(NN.GT.MSMX) THEN
/// CALL SOLVE
/// ELSE
/// CALL SOLVES
/// END IF
/// ELSE
/// CALL RYBSOL
/// END IF
///
/// CALL TIMING(2,ITER)
/// GO TO 10
/// 20 CONTINUE
/// STOP
/// END
/// ```
pub fn run_tlusty(config: &TlustyConfig) -> TlustyResult {
let mut state = TlustyState::default();
let mut timing_stats = TimingStats::default();
let start_time = Instant::now();
// ========================================
// 1. 初始化阶段
// ========================================
// 对应 Fortran:
// INIT=1
// ITER=0
// CALL START
// LFIN=.FALSE.
// IF(NITER.EQ.0) LFIN=.TRUE.
// 执行初始化
// let start_output = start(&mut start_params);
// state.nn = start_output.nn;
// 检查是否需要迭代 (NITER == 0 表示只做初始模型)
// if iterat.niter == 0 {
// state.is_final = true;
// state.is_finished = true;
// }
// ========================================
// 2. 主迭代循环
// ========================================
while !state.is_finished && state.iter < config.max_iter {
state.iter += 1;
let iter_start = Instant::now();
// 2.1 形式解 (RESOLV)
// 对应 Fortran: CALL RESOLV
// let resolv_output = resolv(&mut resolv_params);
state.is_init = false;
if state.is_final {
break;
}
// 2.2 收敛加速 (ACCEL2)
// 对应 Fortran: IF(IACC.GT.0) CALL ACCEL2
if config.accel_mode > 0 {
let accel_start = Instant::now();
// let accel_output = accel2(&mut accel_params);
timing_stats.acceleration_time += accel_start.elapsed().as_secs_f64();
}
// 2.3 解线性化方程
// 对应 Fortran:
// IF(IFRYB.EQ.0) THEN
// IF(NN.GT.MSMX) THEN
// CALL SOLVE
// ELSE
// CALL SOLVES
// END IF
// ELSE
// CALL RYBSOL
// END IF
let solve_start = Instant::now();
let solver_type = select_solver(state.nn, MSMX);
match solver_type {
SolverType::Standard => {
// let solve_output = solve(&solve_params);
}
SolverType::Simple => {
// let solves_output = solves(&solves_params);
}
SolverType::Ryan => {
// let rybsol_output = rybsol(&rybsol_params);
}
}
timing_stats.matrix_solution_time += solve_start.elapsed().as_secs_f64();
// 记录时间 (TIMING)
// 对应 Fortran: CALL TIMING(2,ITER)
// let timing_output = timing(&TimingParams { mode: TimingMode::Iteration, iter: state.iter });
if config.verbose {
println!(
"Iteration {}: time = {:.3}s",
state.iter,
iter_start.elapsed().as_secs_f64()
);
}
// 检查收敛 (在实际实现中检查 CHMX 和 CHMT)
}
let total_time = start_time.elapsed().as_secs_f64();
TlustyResult {
total_iterations: state.iter,
converged: state.is_finished,
total_time_secs: total_time,
timing_stats,
}
}
// ============================================================================
// 辅助函数
// ============================================================================
/// 检查收敛性。
pub fn check_convergence(iter: usize, max_change: f64, tolerance: f64) -> bool {
iter > 0 && max_change < tolerance
}
/// 选择解法。
///
/// 根据 nn 和 MSMX 选择合适的解法：
/// - nn > MSMX: 使用完整矩阵解法 SOLVE
/// - nn <= MSMX: 使用简化解法 SOLVES
pub fn select_solver(nn: usize, msmx: usize) -> SolverType {
if nn > msmx {
SolverType::Standard
} else {
SolverType::Simple
}
}
/// 解法类型。
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum SolverType {
/// 标准解法 (SOLVE)
Standard,
/// 简化解法 (SOLVES)
Simple,
/// Ryan 解法 (RYBSOL)
Ryan,
}
// ============================================================================
// 常量 (从 Fortran 移植)
// ============================================================================
/// 最大简化矩阵维度
pub const MSMX: usize = 2000;

4
src/tlusty/state/constants.rs
View File

@ -26,8 +26,8 @@ pub const MFREQ1: usize = MFREQ; // 根据 BASICS.FOR: MFREQ1 = MFREQ
pub const MFREQP: usize = 220000;
/// 连续谱频率点数
pub const MFREQC: usize = 125000;
/// 线性化频率数
pub const MFREX: usize = 54;
/// 线性化频率数 (increased from Fortran's 54 to accommodate full grid)
pub const MFREX: usize = 200;
/// 每条谱线频率数
pub const MFREQL: usize = 25798;
/// 最大线性化参数数

37
src/tlusty/state/model.rs
View File

@ -348,6 +348,11 @@ pub struct FrqAll {
pub ifs0: i32,
pub kij: Vec<i32>,
pub lskip: Vec<Vec<i32>>,
/// Number of explicit (non-ALI) frequency points for SOLVES (NFREQE)
pub nfreqe: usize,
/// Mapping from explicit frequency index (0..NFREQE-1) to full grid index
/// Set by RESOLV (INIFRC), used by SOLVES (matgen_lte)
pub ijfr_explicit: Vec<usize>,
}
impl Default for FrqAll {
@ -363,6 +368,8 @@ impl Default for FrqAll {
ifs0: 0,
kij: vec![0; MFREQ],
lskip: vec![vec![0; MFREQ]; MDEPTH],
nfreqe: 0,
ijfr_explicit: Vec::new(),
}
}
}
@ -1397,6 +1404,7 @@ pub struct ModelState {
pub totrad: TotRad,
pub currad: CurRad,
pub frqall: FrqAll,
pub exprad: crate::tlusty::state::arrays::ExpRad,
pub totprf: TotPrf,
pub phoexp: PhoExp,
pub obfpar: ObfPar,
@ -2505,6 +2513,26 @@ pub struct ExpRaf {
pub radex: Vec<Vec<f64>>,
/// 扩展 Eddington 因子 [MFREX × MDEPTH]
pub fakex: Vec<Vec<f64>>,
/// Lambda 算子对角线 [MFREX × MDEPTH] (ALI1)
pub aliex: Vec<Vec<f64>>,
/// Lambda 算子次对角线 [MFREX × MDEPTH] (ALIM1)
pub alimp1ex: Vec<Vec<f64>>,
/// Lambda 算子超对角线 [MFREX × MDEPTH] (ALIP1)
pub alipp1ex: Vec<Vec<f64>>,
// Pre-computed radiative equilibrium quantities (from ALIFR1 accumulation)
/// REIT: integral form T-derivative (MDEPTH)
pub reit: Vec<f64>,
/// REDT: differential form T-derivative diagonal (MDEPTH)
pub redt: Vec<f64>,
/// REDTM: differential form T-derivative sub-diagonal (MDEPTH)
pub redtm: Vec<f64>,
/// REDTP: differential form T-derivative super-diagonal (MDEPTH)
pub redtp: Vec<f64>,
/// REDX: differential form NHE column current depth (MDEPTH)
pub redx: Vec<f64>,
/// REDXM: differential form NHE column previous depth (MDEPTH)
pub redxm: Vec<f64>,
}
impl Default for ExpRaf {
@ -2512,6 +2540,15 @@ impl Default for ExpRaf {
Self {
radex: vec![vec![0.0; MDEPTH]; MFREX],
fakex: vec![vec![0.0; MDEPTH]; MFREX],
aliex: vec![vec![0.0; MDEPTH]; MFREX],
alimp1ex: vec![vec![0.0; MDEPTH]; MFREX],
alipp1ex: vec![vec![0.0; MDEPTH]; MFREX],
reit: vec![0.0; MDEPTH],
redt: vec![0.0; MDEPTH],
redtm: vec![0.0; MDEPTH],
redtp: vec![0.0; MDEPTH],
redx: vec![0.0; MDEPTH],
redxm: vec![0.0; MDEPTH],
}
}
}