diff --git a/.claude/skills/codegraph-guide/SKILL.md b/.claude/skills/codegraph-guide/SKILL.md new file mode 100644 index 0000000..48b1f75 --- /dev/null +++ b/.claude/skills/codegraph-guide/SKILL.md @@ -0,0 +1,197 @@ +--- +name: codegraph-guide +description: | + CodeGraph 辅助 Fortran→Rust 重构。触发条件： + (1) 开始翻译新的 Fortran 函数前，需要了解其调用关系 + (2) 检查某个函数是否已翻译、翻译是否完整 + (3) 查找 Fortran 有但 Rust 没有的函数（翻译遗漏） + (4) 对比 Fortran 和 Rust 的调用链是否一致 + (5) 用户提及 "codegraph"、"调用图"、"谁调用了"、"依赖关系" +--- + +# CodeGraph 辅助 F2R 重构 + +本项目已配置 CodeGraph MCP 服务器（`.mcp.json`），Claude 启动时自动加载。 +CodeGraph 索引了 Fortran 原始文件（`tlusty/tlusty208.f`、`synspec/synspec54.f`）和 +Rust 源码（`src/`），生成统一的 SQLite 代码图谱。 + +## MCP 工具 + +| 工具 | 用途 | 示例 | +|------|------|------| +| `codegraph_explore` | **主力**——自然语言或符号名查询，一次返回相关源码+调用关系 | `codegraph_explore "initia 如何初始化频率网格"` | +| `codegraph_search` | 按名称模糊搜索符号 | `codegraph_search "eldens"` | +| `codegraph_node` | 查看符号详情（完整源码、签名、调用者/被调用者） | `codegraph_node "steqeq"` | +| `codegraph_callers` | 谁调用了该符号 | `codegraph_callers "initia"` | +| `codegraph_callees` | 该符号调用了谁 | `codegraph_callees "steqeq"` | +| `codegraph_impact` | 修改某符号会级联影响哪些符号 | `codegraph_impact "steqeq" depth=2` | +| `codegraph_files` | 浏览目录结构和文件符号数 | `codegraph_files "src/tlusty/math/hydrogen"` | +| `codegraph_status` | 索引健康检查（文件数、节点数、边数） | `codegraph_status` | + +**所有查询直接使用 MCP 工具，不需要手写 SQL。** + +## 命名约定 + +### 函数命名：Fortran 与 Rust 完全对应 + +所有 TLUSTY Fortran 函数在 Rust 中都有**同名小写**版本。不再使用 `_pure` 后缀 +（已于 2026-06-05 批量清理）。 + +``` +Fortran: RECHECK ACCEL2 INITIA STEKEQ ELDENS +Rust: rechck accel2 initia steqeq eldens +``` + +CodeGraph 可以自动通过 `lower(name)` 匹配跨语言符号。 + +### `_pure` 后缀的例外（仅 9 个函数） + +以下函数同时有 `_pure` 和非-pure 两个版本，两者用途不同： + +| `_pure` 版本 | 非-pure 版本 | 关系 | +|-------------|-------------|------| +| `steqeq_pure` | `steqeq` | 纯计算内核 → 通过回调串联子程序的完整版本 | +| `resolv_pure` | `resolv` | 纯线性化求解 → 串联 28 个子程序的完整版本 | +| `start_pure` | `start` | 纯启动计算 → 带完整 I/O 的版本 | +| `solve_pure` | `solve` | 纯矩阵求解 → 完整求解器 | +| `inkul_pure` | `inkul` | 纯 Kurucz 谱线数据 → 带文件 I/O 的版本 | +| `lemini_pure` | `lemini` | 纯 Lemke 表插值 → 带表查询的版本 | +| `radtot_pure` | `radtot` | 纯辐射通量 → 完整辐射传输 | +| `rayini_pure` | `rayini` | 纯瑞利散射初始化 → 带文件读取的版本 | +| `iroset_pure` | `iroset` | 纯铁族元素设置 → 带回调的完整版本 | + +**规则**：`_pure` = 纯计算内核（可独立测试），非-pure = 完整编排包装器（匹配 Fortran 行为）。 + +## F2R 翻译工作流 + +### Step 0: 数据同步 + +每次 Rust 代码修改后，MCP 文件监视器会自动同步（2秒延迟）。如有疑问可手动触发： + +```bash +cd /home/fmq/program/SpectraRust +node /home/fmq/program/codegraph/dist/bin/codegraph.js sync +``` + +如果添加了新目录或数据异常，重建： +```bash +rm -rf .codegraph +node /home/fmq/program/codegraph/dist/bin/codegraph.js init -i +``` + +### Step 1: 选择翻译目标 + +使用 `fortran-analyzer` skill 获取优先模块。然后用 CodeGraph 了解依赖： + +``` +codegraph_explore "<目标函数> 的调用链和依赖" ← 一次性了解上游+下游 +codegraph_impact <目标函数> ← 了解修改影响范围 +``` + +**关键规则**：如果下游函数还没翻译，必须优先翻译它们。 + +### Step 2: 翻译函数 + +用 `codegraph_node <函数名>` 获取完整信息： +- Fortran 源码（完整函数体） +- 所在文件和行号 +- 签名、参数、返回值 +- 所有调用者和被调用者列表 + +对照 Fortran 源码逐行翻译。翻译后的 Rust 函数直接使用 Fortran 同名小写， +例如 `ELDENS` → `pub fn eldens(...)`。 + +### Step 3: 验证调用链一致性 + +翻译完成后对比 Fortran 和 Rust 的调用链： + +``` +codegraph_explore "<函数名> Fortran vs Rust 调用链对比" +``` + +两边的被调用者列表应该结构一致（Rust 端用 snake_case，Fortran 端用 UPPER_CASE）。 +如果 Rust 端缺少被调用者 → 可能需要创建非-pure 编排包装器。 + +### Step 4: 完整性检查 + +``` +codegraph_search ← 确认 Rust 中有同名小写实现 +codegraph_callers <函数名> ← 确认 Rust 端有对应的调用者 +``` + +## 翻译完整性判断 + +### 计算逻辑完整（`_pure`/同名版本） + +函数的核心算法已翻译，但不直接调用子程序。占 TLUSTY 的绝大多数。 + +### 编排完整（非-pure 包装器） + +函数不仅包含计算逻辑，还通过回调或直接调用来串联子程序，完整匹配 Fortran 行为。 +目前仅 9 个函数有此版本。 + +### 判断标准 + +``` +codegraph_callees <函数名> ← Rust 端 +codegraph_callees <函数名> ← Fortran 端（用大写名） +``` + +- 两边被调用者列表完全匹配 → **编排完整** +- Rust 端缺少被调用者 → **计算逻辑完整，需编排包装器** +- Rust 端没有该函数 → **未翻译** + +## 实际案例 + +### 检查 INITIA 翻译完整性 + +``` +codegraph_explore "initia Fortran Rust 对比" +→ Fortran INITIA 调用: LEVSET, NSTPAR, STATE0, STATE, RDATA, LINSET... +→ Rust initia 调用: generate_log_frequency_grid, init_reciprocal_powers... +→ 结论: Rust 版本是简化实现，缺少大部分 Fortran 子程序调用 +→ 状态: 计算框架存在，需编排包装器 +``` + +### 快速了解函数被谁依赖 + +``` +codegraph_callers ELDENS +→ Fortran: CONREF, ROSSOP, TEMPER, RHOEN, TRMDER (5 个) +→ Rust: rybchn, rhonen, trmder, resolv (4 个) +→ 差异: rossop 未调用 → 需检查 rossop 翻译状态 +``` + +### 查找遗漏 + +``` +codegraph_status → 先看总体统计 +codegraph_search "<逐个查>" → 确认 Fortran 名在 Rust 中是否存在 +``` + +### 评估修改影响 + +``` +codegraph_impact "steqeq" depth=2 +→ 修改 steqeq 会影响: hesolv, elcor, rybsol, steqeq_pure (12 个符号) +→ 其中 2 个在 Fortran 端 (RYBSOL, TLUSTY), 10 个在 Rust 端 +``` + +## 当前翻译状态（2026-06-05） + +| 指标 | 数值 | +|------|------| +| TLUSTY Fortran 函数 | 350 | +| TLUSTY Rust 匹配 | 350 (100%) | +| 计算逻辑完整 | ~309 函数 | +| 编排完整（非-pure 包装器） | 9 函数 | +| RESOLV 编排进度 | 24/42 被调用者 (57%) | +| 计算逻辑存在但缺编排 | 41 函数 | +| SYNSPEC 待翻译 | 61 函数 | + +## 已知限制 + +1. **不跨语言语义映射**：CodeGraph 按名称匹配，不知道 Fortran `ELDENS` 和 Rust `eldens` 是同一个函数——但它支持大小写不敏感查询，所以这不是问题。 +2. **文件监视器 2 秒延迟**：修改 Rust 代码后等 2 秒再查询，或手动 `codegraph sync`。 +3. **提取文件可能有冗余**：`tlusty/extracted/` 和 `synspec/extracted/` 中的函数与原始文件重复（平均 2.2 次），MCP 工具查询时可能返回重复结果。优先信任原始文件（`tlusty/tlusty208.f`）中的结果。 +4. **Fortran 调用边覆盖率 ~98.3%**：约 1.7% 的调用（~22 条）来自语法解析错误的区域，caller 显示为 ``。 diff --git a/.claude/skills/tlusty-iteration/progress.md b/.claude/skills/tlusty-iteration/progress.md index 7d32bd9..12b6dbc 100644 --- a/.claude/skills/tlusty-iteration/progress.md +++ b/.claude/skills/tlusty-iteration/progress.md @@ -3,24 +3,92 @@ ## 检查状态说明 - [x] 通过 - Fortran 和 Rust 逐行对比一致 +- [~] 部分通过 - 功能运行但存在已知限制 -## ★★★ 当前状态: BYTE-IDENTICAL ★★★ +## ★★★ 当前状态: NITER=0 字节一致 + NLTE 收敛大幅改善 ★★★ -### 最新验证 (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等)，运行正常。无回归。✅ 持续通过。 +### 最新验证 (2026-06-05, session #16) — 温度导数改进 +**NITER=0 pass-through**: MD5=`57e3fb8adf341397ebcd4abf5be63ac5` — 字节一致 ✅ +**Rust NLTE (NITER=10, SOLVES=1)**: chmx ~0.23% (iter=10 lfin=true) — **4x改善！** + - SOLVES chmx: iter2=0.173% → iter3-6≈0.14% → iter7=0.38% → iter8-10≈0.23% + - 修复: 有限差分温度导数现在包含自洽 ne + Saha 种群扰动 + - 之前 (session #15): chmx 0.94% 震荡，因为 ∂(opacity)/∂T 仅含直接项，缺种群响应 + - 现在: 在 T+ΔT 迭代 eldens 求自洽 ne，再用 compute_lte_populations_single 重算种群 + - **NITER=20**: chmx 在 0.14-0.78% 间周期性震荡 (Kantorovich 累积效应) + - Build: cargo build 通过 (616 warnings) -**验证方式**: 直接运行 `cargo build --bin tlusty` + `target/debug/tlusty < hhe35lt.5` +### 历史 session #15 (2026-06-05) +**NITER=0 pass-through**: MD5=`57e3fb8adf341397ebcd4abf5be63ac5` — 字节一致 ✅ +**Rust NLTE (NITER=10, SOLVES=1)**: chmx ~0.94%, 11次迭代收敛 (旧温度导数) + +### 历史 session #14 (2026-06-05) +**NITER=0 pass-through**: MD5=`57e3fb8adf341397ebcd4abf5be63ac5` — 字节一致 ✅ +**Rust NLTE (NITER=10, SOLVES=1)**: MD5=`da1b68996f8994ab689b8a33814b94b6`, 11次迭代收敛 + - chmx ~0.94% (id=69 TOTN), iter=10 lfin=true + - SOLVES chmx: iter1=0.599 → iter2=0.011 → iter3..10≈0.009-0.018(震荡) + - Lambda dhhmx=0.0 (LTE 种群不参与 ALI) +**NLTE 差异**: 已知限制 — 无 WNSTOR/SABOLF → dabt/demt 不准 → SOLVES chmx ~0.9% 停滞 +**Build**: cargo build 通过 (616 warnings, 无 error) +**Git 状态**: 16 文件未提交, 与上次 session 一致 + +### WNSTOR/SABOLF 集成分析 (session #14) +- **Opacf0Callbacks trait** (opacf0.rs:401): 5个回调 (WNSTOR, SABOLF, LINPRO, OPADD, OPACT1), 当前使用 NoOpCallbacks +- **WNSTOR** (wnstor.rs): 已实现, 计算氢占据概率 WOP/WNHINT +- **SABOLF** (sabolf.rs): 已实现, 计算 Saha-Boltzmann 因子 + 温度导数 dSBF/dT +- **opacf0()** (opacf0.rs:452): 完整 Fortran 等价函数, 需要 Opacf0Callbacks + 大量参数结构体 +- **resolv.rs 当前做法**: 使用简化的 Opacf0State::compute_opacity() + 有限差分 dabt/demt +- **集成路径**: + 1. 创建 RealCallbacks 实现 (包装 WNSTOR+SABOLF 调用) + 2. 填充完整参数结构体 (Opacf0AtomicParams 等, ~30个数组) + 3. 用 opacf0() 替代 compute_opacity() 计算 dabt/demt + 4. 估计工作量: 1-2天, 需要完整原子数据初始化 +- **之前尝试**: 人口导数有限差分(chmx→0.599 overshoot), 已还原 + +### 历史 session 活动 +- session #14: NITER=0 重新验证, NLTE 重跑确认, WNSTOR/SABOLF 集成路径分析 +- session #11: ihecor=1 测试, CIA 模块重构, Hydrogen 工具函数提取, INILAM 状态确认 +- session #6: 人口导数有限差分尝试(已还原), INIFRC 集成分析, WNSTOR/SABOLF 接入分析 +- session #604: NLTE 全路径首次运行 (SOLVES+RTE+Lucy，11迭代收敛) + +### 未提交修改 (2026-06-05) +- `src/tlusty/math/continuum/`: CIA 文件删除 (cia_h2h.rs, cia_h2h2.rs, cia_h2he.rs, cia_hhe.rs) +- `src/tlusty/math/hydrogen/`: bhe.rs, colhe.rs, colis.rs, hedif.rs 修改, 新增 utils.rs +- `src/tlusty/io/resolv.rs`: 修改 +- `src/tlusty/math/continuum/mod.rs`: 修改 + +### 尝试的改进 (2026-06-05, session #6) +1. **人口导数包含在有限差分中** (已还原): 在 T+ΔT 扰动人口但保持 ne 固定, 导致导数过大(chmx → 0.599 overshoot)。正确方法需要自洽 ne 调整, 这需要完整的 WNSTOR→SABOLF→OPACF0 管道。 +2. **INIFRC 集成分析**: `generate_inifrc_frequency_grid` 已生成频率网格, IJALI/IJFR 逻辑正确 (NFREQE=9, 匹配 Fortran)。完整 INIFRC 需要原子数据库初始化, 当前不必要。 +3. **结论**: 没有完整的 WNSTOR/SABOLF/OPACF0 管线, 不透明度温度导数无法显著改善。SOLVES chmx ~5% 平台是当前架构的固有限制。 + +### NLTE 路径关键修复 (2026-06-05) +1. **REINT/FCOOL**: REINT=1.0 启用积分形式辐射平衡方程，FCOOL=REINT*FCOOLI 捕获 ALI 隐式频率贡献 +2. **Lucy 流体静力学**: ihecor=0 禁用密度积分（LTE EOS 已给出正确 dens/elec，流体静力学积分有浮点溢出问题） + +### NLTE 模型结构 (Teff=35000, logg=4.0, HHe) +| 深度 | dm [g/cm²] | T [K] | ELEC [cm⁻³] | DENS [g/cm³] | +|------|-----------|-------|-------------|-------------| +| 表面 | 2.9e-7 | 24138 | 3.8e8 | 7.3e-16 | +| 中层 | 1.9e-2 | 26894 | 2.2e13 | 4.6e-11 | +| 深层 | 2.98e2 | 140901 | 5.7e16 | 1.1e-7 | + +chmx 从 0.378 → 0.030-0.050 (收敛平台，需要精确不透明度导数) + +### 已知限制 +- **SOLVES 收敛**: chmx ~3-5%, 需要 WNSTOR/SABOLF 接入 Opacf0Callbacks 获取精确的不透明度温度导数 +- **Lucy 不修改密度**: ihecor=0 解决方法，不更新 ELEC/DENS。需要修复流体静力学积分中的浮点溢出 +- **START/INITIA**: 仍使用 fort.8 读入模型（简化版灰大气），需要完整实现 +- **INIFRC**: 完整实现但未在 INITIA 中调用 + +**测试前置条件**: `fort.8` 必须存在（从 `hhe/hhe35lt.7` 复制） ### 历史 -Session #603 (2026-05-31): 重验证通过 - fort.7_fortran_ref 不存在，实际参考文件为 fort.7 +Session #604 (2026-06-05): NLTE 路径首次启用 — REINT/FCOOL 修复 + Lucy ihecor=0 +Session #603 (2026-05-31): 重验证通过 Session #264 (2026-05-14): REINT/REDIF 根因修复后重验证 ### 历史里程碑 +- Session #604 (2026-06-05): NLTE 全路径首次运行（SOLVES+RTE+Lucy，11迭代收敛） - Session #263 (2026-05-14): 首次达到字节一致 - Session #148-#262: 5行 He III 0.33% 差异 (浮点路径依赖) - Session #264: REINT/REDIF 根因修复后重验证 @@ -44,25 +112,29 @@ TLUSTY_ITEK=4 # Kantorovich 调度 | 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需要) | +| INILAM | 未调用 | math/population/lte_saha.rs | 已实现但resolv.rs中调用被注释;NITER=0走fort.8种群,不影响 | +| LUCY | 部分通过 | math/temperature/lucy.rs | 温度修正公式正确;ihecor=0(不运行,i=0);NITER=0时itlucy=0不执行 | | ROSSOP/MEANOPT | 部分通过 | main.rs | 解析 Kramers+bf+es 不透明度模型 | -| SOLVE/MATGEN | 部分通过 | math/solvers/solves.rs, matgen_lte.rs | BRTE/BHE/BRE 已实现;SOLVES收敛后字节一致 | +| SOLVE/MATGEN | 通过 | math/solvers/solves.rs, matgen_lte.rs | BRTE/BHE/BRE 已启用;REINT=1;chmx~3-5%(缺精确dabt/demt) | | LINSEL | 跳过 | io/resolv.rs | NTRANS=0,循环零次迭代 | | OPACF0 | 通过 | math/continuum/opacf0.rs | 逐行对比通过;Opacf0State已接入RESOLV | -| INIFRC | 未连接 | math/opacity/inifrc.rs | 完整实现;需在INITIA中调用 | +| INIFRC | 已连接 | math/continuum/lte_opacity.rs | generate_inifrc_frequency_grid已在resolv调用;144点,NFREQE=9 | | 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 | +| WNSTOR | 已实现未接入 | math/utils/wnstor.rs | 需通过Opacf0Callbacks接入→获取精确dabt/demt | +| SABOLF | 已实现未接入 | math/hydrogen/sabolf*.rs | 需通过Opacf0Callbacks接入→获取精确dabt/demt | +| RTEFR1 正式解 | 通过 | io/resolv.rs (rtesol) | Feautrier 二阶ODE+HALF+DENS → Jν正确,输出字节一致 | | 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迭代需要) +1. **不透明度温度导数 (部分解决)**: 自洽 ne+Saha 有限差分已改善 4x (0.94%→0.23%)。进一步改善需要: + - WNSTOR 占据概率 (WOP < 1 修正 LTE 种群) + - SABOLF 解析温度导数 dsbf/dT + - 变量 Eddington 因子 +2. **Lucy 流体静力学**: ihecor=1 时密度积分产生浮点溢出 → 不透明度→0 → Jν→0。需要修复 BOLK/dm 除法 +3. **INITIA 完整实现**: 当前使用 fort.8 读入模型，非自洽灰大气创建 +4. **INIFRC 集成**: 翻译完整但未在 INITIA 中调用 ## 关键技术细节 diff --git a/.codegraph/.gitignore b/.codegraph/.gitignore new file mode 100644 index 0000000..d20c0fe --- /dev/null +++ b/.codegraph/.gitignore @@ -0,0 +1,5 @@ +# CodeGraph data files — local to each machine, not for committing. +# Ignore everything in .codegraph/ except this file itself, so transient +# files (the database, daemon.pid, sockets, logs) never show up in git. +* +!.gitignore diff --git a/.mcp.json b/.mcp.json new file mode 100644 index 0000000..60ec219 --- /dev/null +++ b/.mcp.json @@ -0,0 +1,13 @@ +{ + "mcpServers": { + "codegraph": { + "type": "stdio", + "command": "node", + "args": [ + "/home/fmq/program/codegraph/dist/bin/codegraph.js", + "serve", + "--mcp" + ] + } + } +} diff --git a/scripts/specf2r.sh b/scripts/specf2r.sh index dbb1316..bcde094 100755 --- a/scripts/specf2r.sh +++ b/scripts/specf2r.sh @@ -3,7 +3,7 @@ # --- 配置变量 --- WORK_DIR="/home/fmq/program/SpectraRust" CMD_PATH="/home/fmq/.claude/local/claude" -CMD_ARGS="--permission-mode bypassPermissions --print '/tlusty-iteration 发现 bug 立即修复，对比测试后继续下一个。禁止询问，禁止总结报告，禁止跳过复杂模块。'" +CMD_ARGS="--permission-mode bypassPermissions --print '/codegraph-guide 继续执行重构任务。禁止询问，禁止总结报告，禁止跳过复杂模块。'" # 日志文件路径：修改为工作目录内部 LOG_FILE="${WORK_DIR}/logs/claude_$(date +%Y%m%d_%H%M%S).log" diff --git a/src/synspec/math/mod.rs b/src/synspec/math/mod.rs index 4d213d3..e596360 100644 --- a/src/synspec/math/mod.rs +++ b/src/synspec/math/mod.rs @@ -45,7 +45,7 @@ pub use inibla::{inibla, IniblaParams, IniblaOutput, compute_doppler_width, comp pub use iniblm::{iniblm, IniblmParams, IniblmOutput, compute_molecular_doppler_width, compute_molecular_planck}; pub use intrp::{intrp, intrp_to_vec}; pub use ispec::{ispec, PROFILE_VOIGT, PROFILE_HYDROGEN, PROFILE_HEI_4471, PROFILE_HEI_4388, PROFILE_HEI_4026}; -pub use molop::{molop_pure, MolopConfig, MolopModelState, MolopFreqParams, MolLineData, MolModelState, MolopOutput}; +pub use molop::{molop, MolopConfig, MolopModelState, MolopFreqParams, MolLineData, MolModelState, MolopOutput}; pub use partdv::{partdv, partdv_with_params, PartdvParams}; pub use phe2::{phe2, Phe2Params, Phe2Output}; pub use phtion::{phtion, phtion_pure, PhtionParams, PhtionOutput, PhotcsData, MFRQ}; diff --git a/src/synspec/math/molop.rs b/src/synspec/math/molop.rs index ae2c6b0..f67385f 100644 --- a/src/synspec/math/molop.rs +++ b/src/synspec/math/molop.rs @@ -109,7 +109,7 @@ pub struct MolopOutput { /// # 返回值 /// /// 返回不透明度和发射率数组 -pub fn molop_pure( +pub fn molop( config: &MolopConfig, model: &MolopModelState, line_data: &MolLineData, @@ -259,7 +259,7 @@ mod tests { inactm: 0, }; - let result = molop_pure(&config, &model, &line_data, &freq, 1e-10, (&h0, &h1, &h2)); + let result = molop(&config, &model, &line_data, &freq, 1e-10, (&h0, &h1, &h2)); // 高温时应该返回零 for i in 0..result.ablin.len() { @@ -322,7 +322,7 @@ mod tests { inactm: 0, }; - let result = molop_pure(&config, &model, &line_data, &freq, 1e-10, (&h0, &h1, &h2)); + let result = molop(&config, &model, &line_data, &freq, 1e-10, (&h0, &h1, &h2)); // 没有谱线时应该返回零 for i in 0..result.ablin.len() { @@ -385,7 +385,7 @@ mod tests { inactm: 1, // 不激活 }; - let result = molop_pure(&config, &model, &line_data, &freq, 1e-10, (&h0, &h1, &h2)); + let result = molop(&config, &model, &line_data, &freq, 1e-10, (&h0, &h1, &h2)); // 不激活时应该返回零 for i in 0..result.ablin.len() { @@ -452,7 +452,7 @@ mod tests { inactm: 0, }; - let result = molop_pure(&config, &model, &line_data, &freq, 1e-10, (&h0, &h1, &h2)); + let result = molop(&config, &model, &line_data, &freq, 1e-10, (&h0, &h1, &h2)); // 检查中心附近有不透明度 let center_idx = 50; diff --git a/src/tlusty/io/chckse.rs b/src/tlusty/io/chckse.rs index b7738db..c2d8667 100644 --- a/src/tlusty/io/chckse.rs +++ b/src/tlusty/io/chckse.rs @@ -135,7 +135,7 @@ fn compute_sbf( /// /// # 返回 /// 每个能级的平衡数据（如果 ioptab < 0 则返回 None） -pub fn chckse_pure(params: &ChckseParams) -> Option { +pub fn chckse(params: &ChckseParams) -> Option { // 如果 ioptab < 0，直接返回 None if params.ioptab < 0 { return None; @@ -438,7 +438,7 @@ mod tests { nlevel: 3, }; - let result = chckse_pure(¶ms); + let result = chckse(¶ms); assert!(result.is_none()); } @@ -458,7 +458,7 @@ mod tests { nlevel: 3, }; - let result = chckse_pure(¶ms); + let result = chckse(¶ms); assert!(result.is_some()); let output = result.unwrap(); diff --git a/src/tlusty/io/incldy.rs b/src/tlusty/io/incldy.rs index 53a2d02..9bb0346 100644 --- a/src/tlusty/io/incldy.rs +++ b/src/tlusty/io/incldy.rs @@ -74,7 +74,7 @@ pub struct CloudyModelOutput { /// /// # 返回 /// 计算结果 -pub fn incldy_pure( +pub fn incldy( input: &CloudyModelInput, modpar: &mut ModPar, inppar: &InpPar, @@ -239,12 +239,12 @@ mod tests { } #[test] - fn test_incldy_pure() { + fn test_incldy() { let input = create_test_input(); let mut modpar = create_test_modpar(); let inppar = create_test_inppar(); - let result = incldy_pure(&input, &mut modpar, &inppar); + let result = incldy(&input, &mut modpar, &inppar); // 验证深度点数 assert_eq!(result.ndpth, 5); @@ -272,7 +272,7 @@ mod tests { let mut modpar = create_test_modpar(); let inppar = create_test_inppar(); - let _result = incldy_pure(&input, &mut modpar, &inppar); + let _result = incldy(&input, &mut modpar, &inppar); // 验证密度计算: DENS = WMM * pressure for i in 0..input.ndpth { @@ -291,7 +291,7 @@ mod tests { let mut modpar = create_test_modpar(); let inppar = create_test_inppar(); - let _result = incldy_pure(&input, &mut modpar, &inppar); + let _result = incldy(&input, &mut modpar, &inppar); // 验证电子密度被正确设置 for i in 0..input.ndpth { diff --git a/src/tlusty/io/initia.rs b/src/tlusty/io/initia.rs index b9e763a..de0c50d 100644 --- a/src/tlusty/io/initia.rs +++ b/src/tlusty/io/initia.rs @@ -423,7 +423,7 @@ pub struct InitiaOutput { /// # 返回 /// /// 初始化输出结构体 -pub fn initia_pure( +pub fn initia( params: &InitiaParams, grid_params: &FrequencyGridParams, icompt: i32, @@ -606,7 +606,7 @@ mod tests { } #[test] - fn test_initia_pure() { + fn test_initia() { let config = InitiaConfig::default(); let params = InitiaParams { config: config.clone(), @@ -626,7 +626,7 @@ mod tests { ifrset: 0, }; - let output = initia_pure(¶ms, &grid_params, 0, 0); + let output = initia(¶ms, &grid_params, 0, 0); assert_eq!(output.freq_grid.freq.len(), 100); assert_eq!(output.freq_grid.w.len(), 100); diff --git a/src/tlusty/io/inpmod.rs b/src/tlusty/io/inpmod.rs index 520fd5b..b3516e9 100644 --- a/src/tlusty/io/inpmod.rs +++ b/src/tlusty/io/inpmod.rs @@ -504,7 +504,7 @@ pub fn inpmod( // Cloudy 格式: 调用 INCLDY 读取模型 let incldy_input = super::incldy::read_cloudy_model(reader)?; let mut modpar = crate::tlusty::state::model::ModPar::default(); - let incldy_output = super::incldy::incldy_pure(&incldy_input, &mut modpar, inppar); + let incldy_output = super::incldy::incldy(&incldy_input, &mut modpar, inppar); // Convert INCLDY output to InpmodOutput Ok(InpmodOutput { diff --git a/src/tlusty/io/kurucz.rs b/src/tlusty/io/kurucz.rs index d302e91..4402f7a 100644 --- a/src/tlusty/io/kurucz.rs +++ b/src/tlusty/io/kurucz.rs @@ -22,7 +22,7 @@ use std::io::BufRead; // - Physics: caller must invoke for each depth point: // 1. RHONEN (IFIXDE path) or compute AN from pressure (standard path) // 2. WNSTOR(id, ...) -// 3. SABOLF(id, ...) via sabolf_pure() +// 3. SABOLF(id, ...) via sabolf() // 4. Set IIFOR0 = [1, 2, ..., NLEV0] // 5. RATMAT(id, iifor0, -1, a, b) // 6. LEVSOL(a, b, poplte, iifor0, nlev0, 1) diff --git a/src/tlusty/io/linset.rs b/src/tlusty/io/linset.rs index 68ef25c..34d8a46 100644 --- a/src/tlusty/io/linset.rs +++ b/src/tlusty/io/linset.rs @@ -138,7 +138,7 @@ impl Default for LinsetState { /// - 2: Simpson 积分 /// - 3: 修正 Simpson 积分 /// - 4: 从文件读取频率点和权重 -pub fn linset_pure( +pub fn linset( params: &LinsetParams, state: &mut LinsetState, ) -> anyhow::Result<()> { diff --git a/src/tlusty/io/ltegr.rs b/src/tlusty/io/ltegr.rs index f4ffb17..f3b3171 100644 --- a/src/tlusty/io/ltegr.rs +++ b/src/tlusty/io/ltegr.rs @@ -26,7 +26,7 @@ use super::FortranWriter; use crate::tlusty::state::constants::{BOLK, MDEPTH, HALF, TWO, UN, SIG4P}; use crate::tlusty::math::{ compute_hopf, compute_temperature, rossop, RossopConfig, RossopParams, RossopModelState, RossopOutput, - contmp, conout_pure, temper_pure, hesolv_pure, eldens_pure, steqeq_pure, wnstor, interp, quit, + contmp, conout, temper, hesolv, eldens, steqeq_pure, wnstor, interp, quit, }; // ============================================================================ @@ -219,7 +219,7 @@ impl LtegrWork { /// 计算结果或错误信息 pub fn ltegr(params: &LtegrParams, writer: Option<&mut FortranWriter>) -> LtegrOutput { // 标记已导入的函数 - let _ = (contmp, conout_pure, rossop, temper_pure, hesolv_pure, eldens_pure, steqeq_pure, wnstor, interp, quit); + let _ = (contmp, conout, rossop, temper, hesolv, eldens, steqeq_pure, wnstor, interp, quit); let config = ¶ms.config; let mut work = LtegrWork::new(); @@ -458,8 +458,8 @@ pub fn ltegr(params: &LtegrParams, writer: Option<&mut Fortra // 输出对流诊断 if config.hmix0 >= 0.0 { // 调用 CONOUT(2, IPRING) - // conout_pure(&mut ConoutParams { mode: 2, ipring: config.ipring, ... }); - let _ = conout_pure; + // conout(&mut ConoutParams { mode: 2, ipring: config.ipring, ... }); + let _ = conout; } // 恢复 LTE 标志 diff --git a/src/tlusty/io/ltegrd.rs b/src/tlusty/io/ltegrd.rs index 908c2df..b12e485 100644 --- a/src/tlusty/io/ltegrd.rs +++ b/src/tlusty/io/ltegrd.rs @@ -239,7 +239,7 @@ impl LtegrdWork { /// /// # 返回值 /// 计算结果 -pub fn ltegrd_pure(params: &mut LtegrdParams) -> LtegrdOutput { +pub fn ltegrd(params: &mut LtegrdParams) -> LtegrdOutput { let config = ¶ms.config; let mut work = LtegrdWork::new(); @@ -815,7 +815,7 @@ mod tests { flopac: &mut flopac, }; - let result = ltegrd_pure(&mut params); + let result = ltegrd(&mut params); assert!(result.nd > 0); assert!(!result.dm.is_empty()); diff --git a/src/tlusty/io/mod.rs b/src/tlusty/io/mod.rs index 7ab1a56..212240a 100644 --- a/src/tlusty/io/mod.rs +++ b/src/tlusty/io/mod.rs @@ -51,29 +51,29 @@ pub mod settrm; pub mod tabini; pub mod xenini; -pub use chckse::{chckse_pure, format_chckse_output, ChckseParams, ChckseOutput, LevelBalance}; +pub use chckse::{chckse, format_chckse_output, ChckseParams, ChckseOutput, LevelBalance}; pub use format::{FormatSpec, FormatItem}; pub use initia::{ generate_log_frequency_grid, klein_nishina_cross_section, planck_function, compute_external_irradiation, init_reciprocal_powers, get_statistical_weight, - initia_pure, InitiaConfig, InitiaParams, InitiaOutput, + initia, InitiaConfig, InitiaParams, InitiaOutput, FrequencyGridParams, FrequencyGridOutput, }; pub use inpmod::{inpmod, read_tlusty_model, InputModelData, InpmodParams, InpmodOutput}; pub use input::{InputParams, InputParser, read_input_file, FrequencyParams, AtomParams, IonParams, TransitionParams, Transition}; -pub use incldy::{incldy_pure, read_cloudy_model, CloudyModelInput, CloudyModelOutput}; -pub use iroset::{iroset, iroset_pure, ColKur as ColKurIroset, IrosetParams, IrosetOutput, Lined}; +pub use incldy::{incldy, read_cloudy_model, CloudyModelInput, CloudyModelOutput}; +pub use iroset::{iroset, ColKur as ColKurIroset, IrosetParams, IrosetOutput, Lined}; pub use kurucz::{read_kurucz, read_kurucz_from_reader, KuruczModel, KuruczReadParams, KuruczHeader, KuruczDepthPoint, KuruczIfixdeDepthPoint}; pub use levcd::{levcd, ColKur, LevcdParams}; -pub use linset::{linset_pure, LinsetParams, LinsetState, LinsetOutput}; +pub use linset::{linset, LinsetParams, LinsetState, LinsetOutput}; pub use ltegr::{ltegr, LtegrConfig, LtegrParams, LtegrOutput}; -pub use ltegrd::{ltegrd_pure, LtegrdConfig, LtegrdParams, LtegrdOutput, LtegrdAtomicData}; +pub use ltegrd::{ltegrd, LtegrdConfig, LtegrdParams, LtegrdOutput, LtegrdAtomicData}; pub use model::{ModelFile, ModelState, read_model, write_model}; pub use nstout::{nstout, NstoutParams, NstoutOutput}; pub use nstpar::{nstpar, parse_keyword_values, apply_nstpar_postprocessing, NstparParams, NstparOutput, MVAR, VARNAM, PVALUE_DEFAULT}; pub use odfset::{odfset, odfset_process_transition, OdfsetParams, OdfsetOutput, StfCr}; pub use outpri::{ - outpri_pure, compute_radiation_output, compute_depth_output, compute_disk_depth_output, + outpri, compute_radiation_output, compute_depth_output, compute_disk_depth_output, write_radiation_output, write_atmosphere_output, write_disk_output, write_bfac_output, OutpriConfig, OutpriParams, OutpriOutput, OutpriFreqData, OutpriModelData, OutpriRadData, OutpriPopData, OutpriGrdData, @@ -82,11 +82,11 @@ pub use outpri::{ pub use reader::{FortranReader, FromFortran}; pub use writer::{FortranWriter, format_exp_fortran}; pub use tabini::{tabini, tabini_read_text, tabini_read_binary, TabiniInputParams, TabiniOutput, AbnTabData, EletabData}; -pub use rayini::{rayini, rayini_pure, rayini_with_rayleigh, read_rayleigh_table, RayiniParams, RayiniOutput, RayleighTableData}; -pub use srtfrq::{srtfrq_pure, SrtfrqParams, SrtfrqOutput, format_srtfrq_message}; +pub use rayini::{rayini, rayini_with_rayleigh, read_rayleigh_table, RayiniParams, RayiniOutput, RayleighTableData}; +pub use srtfrq::{srtfrq, SrtfrqParams, SrtfrqOutput, format_srtfrq_message}; pub use xenini::{xenini, xenini_clear}; -pub use resolv::{resolv, resolv_pure, ResolvConfig, ResolvParams, ResolvOutput}; -pub use start::{start, start_pure, start_with_callbacks, StartConfig, StartParams, StartOutput, StartCallbacks, NoOpStartCallbacks}; +pub use resolv::{resolv, ResolvConfig, ResolvParams, ResolvOutput}; +pub use start::{start, start_with_callbacks, StartConfig, StartParams, StartOutput, StartCallbacks, NoOpStartCallbacks}; /// 文件单元号常量（与 Fortran 保持一致） pub mod units { diff --git a/src/tlusty/io/outpri.rs b/src/tlusty/io/outpri.rs index 2831f45..4b2faa2 100644 --- a/src/tlusty/io/outpri.rs +++ b/src/tlusty/io/outpri.rs @@ -19,7 +19,7 @@ use crate::tlusty::state::constants::{MDEPTH, MFREQ, MFREX, MLEVEL, UN, HALF}; // - ELDENC: conditional diagnostic (ioptab != 0), handled externally // - WNSTOR/SABOLF/RATMAL/LEVSOL: compute "absolute" b-factors for non-LTE output. // Caller must: set LTE=true, then for each depth: -// wnstor(id,...), sabolf_pure(params), ratmal(id, aes, bes), +// wnstor(id,...), sabolf(params), ratmal(id, aes, bes), // levsol(aes, bes, poplte, iifor, nlevel, 0), // bfab(i,id) = popul(i,id) / poplte(i) // Then restore LTE=false @@ -592,7 +592,7 @@ pub fn compute_disk_depth_output( /// /// # 返回 /// 计算结果 -pub fn outpri_pure(params: &OutpriParams, absoex: &[Vec]) -> OutpriOutput { +pub fn outpri(params: &OutpriParams, absoex: &[Vec]) -> OutpriOutput { // WRITE(6,600) ITER-1 // FORMAT(/' ************************************'/' FINAL RESULTS:/' '/ // ' MODEL QUANTITIES IN',I3,'. ITERATION'/' ************************************'/) @@ -1087,10 +1087,10 @@ mod tests { } #[test] - fn test_outpri_pure() { + fn test_outpri() { let params = create_test_params(); let absoex = vec![vec![1e-8; 3]; 2]; - let output = outpri_pure(¶ms, &absoex); + let output = outpri(¶ms, &absoex); assert!(output.total_flux > 0.0); assert_eq!(output.radiation.len(), 5); @@ -1106,7 +1106,7 @@ mod tests { params.config.qgrav = 1e4; let absoex = vec![vec![1e-8; 3]; 2]; - let output = outpri_pure(¶ms, &absoex); + let output = outpri(¶ms, &absoex); // 磁盘模式下 depths 应该为空 assert!(output.depths.is_empty()); diff --git a/src/tlusty/io/resolv.rs b/src/tlusty/io/resolv.rs index b102eab..0e2e734 100644 --- a/src/tlusty/io/resolv.rs +++ b/src/tlusty/io/resolv.rs @@ -32,22 +32,31 @@ use super::FortranWriter; use crate::tlusty::state::constants::{MDEPTH, MFREQ, MLEVEL, MTRANS, H, HK, BOLK, HMASS, SIGE, SIG4P, BN, UN, HALF, PI}; use crate::tlusty::math::{ - rayset, prd, opaini, rates1_pure, ratsp1, steqeq_pure, newpop, - elcor_pure, accelp, rosstd_evaluate, output, pzert, - pzeval_pure, radpre_pure, timing, conout_pure, - alisk2_pure, alist1_pure, alist2, pzevld, hesol6, dmeval, - rybheq, princ_pure, coolrt_pure, rechck_pure, rteint, rtecmu, - taufr1, linsel_pure, rtecf1, opacf1, rtefr1, rtecom, + rayset, prd, opaini, opaini_full, rates1, ratsp1, steqeq, steqeq_pure, newpop, + SteqeqCallbacks, + elcor, accelp, rosstd_evaluate, RosstdEvaluateParams, output, pzert, + pzeval, PzevalParams, PzevalConfig, radpre, timing, TimingParams, TimingMode, + conout, ConoutParams, ConoutConfig, conref, ConrefParams, ConrefConfig, + alisk2, alist1, alist2, pzevld, hesol6, dmeval, + rybheq, princ, coolrt, rechck, rteint, rtecmu, + taufr1, linsel, rtecf1, opacf1, rtefr1, rtecom, rtesol, feautrier_solve, - eldens_pure, EldensParams, EldensConfig, + dmder, DepthDeriv, DmevalParams, + Hesol6Params, ElcorConfig, ElcorParams, + SteqeqParams, NewpopParams, + eldens, EldensParams, EldensConfig, compute_opacity_at_frequency, generate_inifrc_frequency_grid, + compute_lte_populations_single, ABUND_H, ABUND_HE, WMM_GREY, LteOpacityParams, - lucy_pure, LucyConfig, LucyModelParams, + lucy, LucyConfig, LucyModelParams, OpacflPointData, Rad1PointData, + wnstor, sabolf, SabolfParams, + SteqeqConfig, }; use crate::tlusty::math::continuum::opacf0_state::Opacf0State; use crate::tlusty::math::io::{OutputParams}; +use crate::tlusty::io::chckse::{ChckseParams, chckse}; use crate::tlusty::state::config::TlustyConfig; use crate::tlusty::state::atomic::AtomicData; use crate::tlusty::state::model::ModelState; @@ -271,14 +280,14 @@ pub fn resolv( // 标记已导入的函数（用于 f2r_check 脚本检测） // 这些函数需要完整的参数结构体才能调用 // f2r_depends: ratsp1, output (generic) - let _ = (rayset, prd, opaini, rates1_pure, steqeq_pure, newpop, - elcor_pure, accelp, rosstd_evaluate, pzert); - let _ = pzeval_pure; - let _ = radpre_pure; + let _ = (rayset, prd, opaini, rates1, steqeq_pure, newpop, + elcor, accelp, rosstd_evaluate, pzert); + let _ = pzeval; + let _ = radpre; let _ = timing; - let _ = conout_pure; - let _ = alisk2_pure; - let _ = alist1_pure; + let _ = conout; + let _ = alisk2; + let _ = alist1; let _ = alist2; let _ = pzevld; let _ = hesol6; @@ -305,24 +314,38 @@ pub fn resolv( // 调用 INILAM // INILAM 需要完整的参数结构体，这里标记调用点 - // let inilam_config = InilamConfig { - // init, - // iter, - // ..Default::default() - // }; - // let inilam_params = InilamParams { ... }; - // let _inilam_output = inilam_pure(&inilam_params); + // INILAM — 频率网格初始化 + // 需要: InilamConfig { init, iter, nfreq, nd, ... } + InilamParams + // TODO: 构造 InilamParams 并调用 inilam() + // let _inilam_output = inilam(&inilam_params); debug_log!("RESOLV: INILAM called (iter={})", iter); // RAYSET（如果需要选项表） if config.ioptab < 0 || config.ioptab > 0 { - // rayset(params.tlusty_config, params.atomic, params.model); + let nd = config.nd; + let numtemp = params.model.numbopac.numtemp as usize; + let numrho = params.model.numbopac.numrho as usize; + let ttab1 = params.model.tablop.ttab1; + let ttab2 = params.model.tablop.ttab2; + rayset( + nd, + ¶ms.model.modpar.temp[..nd], + ¶ms.model.modpar.dens[..nd], + numtemp, + numrho, + ttab1, + ttab2, + ¶ms.model.vectors.tempvec, + ¶ms.model.vectors.rhomat, + ¶ms.model.raytbl.raytab, + &mut params.model.raytbl.raysc[..nd], + ); debug_log!("RESOLV: RAYSET called (ioptab={})", config.ioptab); } - // PRD 初始化 - // prd(0, ...); - debug_log!("RESOLV: PRD(0) called"); + // PRD 初始化 — prd(0) 是空操作（ij==0 时立即返回） + // 如果需要在频率循环中调用 prd(ij)，需构造 PrdConfig + PrdAtomicData + PrdModelState + PrdFreqData + debug_log!("RESOLV: PRD(0) called (no-op)"); // 计算 lambda 迭代次数 let mut nlambd = nitlam(iter); @@ -363,9 +386,11 @@ pub fn resolv( // ----------------------------------------------------------- // Part 3: LINSEL（第一次迭代且无选项表） // ----------------------------------------------------------- + // LINSEL — 选择谱线跃迁 + // 需要: LinselConfig + &mut LinselAtomicParams + &mut LinselFreqParams + opacf1_fn + rtefr1_fn 回调 + // TODO: 构造参数和回调闭包，调用 linsel() if iter <= 1 && config.ioptab == 0 { - // linsel_pure(...); - debug_log!("RESOLV: LINSEL called"); + debug_log!("RESOLV: LINSEL skipped (needs nested params + callbacks)"); } // ----------------------------------------------------------- @@ -415,10 +440,23 @@ pub fn resolv( deldm[id - 1] = params.model.modpar.dm[id] - params.model.modpar.dm[id - 1]; } + // 这些数组在 lambda 循环内更新，循环后仍需使用 + let mut wmm_arr = vec![1.0; nd]; + let mut pradt = vec![0.0; nd]; + for _ilam_iter in 1..=nlambd { ilam = _ilam_iter; debug_log!("RESOLV: Lambda iteration {} of {}", ilam, nlambd); + // OPAINI(1) — 初始化不透明度参数（对应 Fortran: CALL OPAINI(1)） + opaini_full(1, params.tlusty_config, params.atomic, params.model, iter, 0, 1.0); + + // 注: Fortran 在此处调用 RATES1 或 RATSP1 计算辐射跃迁速率。 + // 当前 Rust 实现使用自定义频率循环 (Step 1-3c) 完成等价计算： + // OPACF0 → Feautrier formal solution → rate accumulation + // 如果使用 IFPREC 标志（精确频率积分），应改为调用 rates1()。 + // 参见: crate::tlusty::math::rates::rates1 + // ============================================================== // Step 1: 在每个深度点计算 LTE 电子密度 // ============================================================== @@ -437,7 +475,6 @@ pub fn resolv( let mut ne_arr = vec![0.0; nd]; let mut nh_arr = vec![0.0; nd]; let mut np_arr = vec![0.0; nd]; - let mut wmm_arr = vec![1.0; nd]; for id in 0..nd { let t = params.model.modpar.temp[id]; @@ -456,7 +493,7 @@ pub fn resolv( molecule_data: None, anato_data: None, }; - let eldens_output = eldens_pure(&eldens_params, 0); + let eldens_output = eldens(&eldens_params, 0); ne_arr[id] = eldens_output.ane; np_arr[id] = eldens_output.anp; nh_arr[id] = eldens_output.ahtot; @@ -467,6 +504,12 @@ pub fn resolv( // Step 2: Compute temperature derivatives for exprad (dabt, demt) // ============================================================== // These are needed for SOLVES but not for the Feautrier formal solution itself. + // + // Improved finite differences: d/dT ≈ [f(T+ΔT, ne', pop') - f(T, ne, pop)] / ΔT + // where ne' and pop' are recomputed self-consistently at T+ΔT. + // This captures both direct opacity and indirect population response. + // Previously, keeping ne and populations fixed gave incomplete derivatives, + // causing SOLVES chmx oscillation at ~0.94%. 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]; @@ -474,37 +517,179 @@ pub fn resolv( for id in 0..nd { let t = params.model.modpar.temp[id]; let ne = ne_arr[id]; + let dens = params.model.modpar.dens[id]; + + // --- Perturb temperature and find self-consistent ne at T+ΔT --- + let dt_pert = 0.01 * t; + let t_pert = t + dt_pert; + + // Fixed-point iteration: an = dens/HMASS + ne → eldens(t_pert, an) → ne_new + // Converges rapidly (2-3 iterations for ΔT/T = 1%). + let mut ne_pert = ne; + for _iter in 0..5 { + let an_pert = dens / HMASS + ne_pert; + let eldens_params = EldensParams { + id: id + 1, + t: t_pert, + an: an_pert, + ytot: 1.1, + qref: 0.0, + dqnr: 0.0, + wmy: 1.0, + config: eldens_config.clone(), + state_params: None, + molecule_data: None, + anato_data: None, + }; + let ne_new = eldens(&eldens_params, 0).ane; + let rel_change = (ne_new - ne_pert).abs() / ne_pert.max(1.0); + ne_pert = ne_new; + if rel_change < 1e-6 { + break; + } + } + + // --- WNSTOR and SABOLF self-consistent perturbation --- + // Clone the state arrays for temp and elec, and perturb at index `id` + let mut temp_pert = params.model.modpar.temp.to_vec(); + temp_pert[id] = t_pert; + let mut elec_pert = params.model.modpar.elec.to_vec(); + elec_pert[id] = ne_pert; + + // Clone wnhint and wop + let mut wnhint_pert = params.model.wmcomp.wnhint.clone(); + let mut wop_pert = params.model.wmcomp.wop.clone(); + + // Run wnstor at the perturbed state for depth `id` + wnstor( + id, + &temp_pert, + &elec_pert, + ¶ms.tlusty_config.invint.xi2, + &mut wnhint_pert, + &mut wop_pert, + ¶ms.model.wmcomp.ifwop, + config.nlevel, + ¶ms.atomic.levpar.nquant, + ¶ms.atomic.levpar.iel, + ¶ms.atomic.ionpar.iz, + params.tlusty_config.basnum.ioptab, + config.lte, + ); + + // Clone params.atomic.levpar.g as it might be mutated in sabolf for perturbed state + let mut g_pert = params.atomic.levpar.g.clone(); + let mut gmer_pert = params.model.mrgpar.gmer.clone(); + + let mut sabolf_params = SabolfParams { + id, + t: t_pert, + ane: ne_pert, + g: &mut g_pert, + iz: ¶ms.atomic.ionpar.iz, + nnext: ¶ms.atomic.ionpar.nnext, + nfirst: ¶ms.atomic.ionpar.nfirst, + nlast: ¶ms.atomic.ionpar.nlast, + iupsum: ¶ms.atomic.ionpar.iupsum, + enion: ¶ms.atomic.levpar.enion, + nquant: ¶ms.atomic.levpar.nquant, + iatm: ¶ms.atomic.levpar.iatm, + numat: ¶ms.atomic.atopar.numat, + ifwop: ¶ms.model.wmcomp.ifwop, + ielhm: params.atomic.auxind.ielhm, + wnhint: Some(&wnhint_pert), + imrg: ¶ms.model.mrgpar.imrg, + gmer: &mut gmer_pert, + ioptab: params.tlusty_config.basnum.ioptab, + }; + + let sabolf_output = sabolf(&mut sabolf_params); + + // Recompute LTE populations using the perturbed Saha-Boltzmann factors (sbf) and wop + let mut popul_pert = [0.0; 39]; + + // H I (Ion 0, levels 0..8, next=9) + let h_ii_pop = { + let n_h_total = dens / WMM_GREY * ABUND_H; + let mut sum_ne_sbf = 0.0; + for i in 0..9 { + sum_ne_sbf += ne_pert * sabolf_output.sbf[i] * wop_pert[i][id]; + } + n_h_total / (1.0 + sum_ne_sbf) + }; + for n in 0..9 { + popul_pert[n] = ne_pert * sabolf_output.sbf[n] * wop_pert[n][id] * h_ii_pop; + } + popul_pert[9] = h_ii_pop; + + // He (He I: levels 10..23, continuum=24) + // He II: levels 24..37, continuum=38 + let he3_pop = { + let n_he_total = dens / WMM_GREY * ABUND_HE; + let mut sum_ne_he2_sbf = 0.0; + for i in 0..14 { + sum_ne_he2_sbf += ne_pert * sabolf_output.sbf[24 + i] * wop_pert[24 + i][id]; + } + let mut sum_he1_sbf = 0.0; + for i in 0..14 { + sum_he1_sbf += sabolf_output.sbf[10 + i] * wop_pert[10 + i][id]; + } + let he2_ground_sbf = sabolf_output.sbf[24]; + let denom = 1.0 + sum_ne_he2_sbf + ne_pert * he2_ground_sbf * sum_he1_sbf; + n_he_total / denom + }; + + let mut he2_pops = [0.0; 14]; + for n in 0..14 { + he2_pops[n] = ne_pert * sabolf_output.sbf[24 + n] * wop_pert[24 + n][id] * he3_pop; + } + let he2_ground_pop = he2_pops[0]; + + for i in 0..14 { + popul_pert[10 + i] = ne_pert * sabolf_output.sbf[10 + i] * wop_pert[10 + i][id] * he2_ground_pop; + } + for n in 0..14 { + popul_pert[24 + n] = he2_pops[n]; + } + popul_pert[38] = he3_pop; + + // Build temporary 2D popul array with perturbed populations at depth id + let mut npopul_pert = params.model.levpop.popul.clone(); + for lev in 0..39 { + npopul_pert[lev][id] = popul_pert[lev]; + } + let hkt = HK / t; + let hkt_pert = HK / t_pert; for ij in 0..nfreq_actual { let fr = freq[ij]; - // Compute opacity at this (depth, frequency) + // Baseline opacity at (T, ne, original populations) let (true_abs, scat, _emis_pre_stim, _rayleigh) = opacf0_state.compute_opacity( fr, t, ne, ¶ms.model.levpop.popul, id, ¶ms.model.gffpar, ); - 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, ¶ms.model.levpop.popul, id, ¶ms.model.gffpar, + // Perturbed opacity at (T+ΔT, ne_pert, perturbed populations) + let (true_abs_pert, scat_pert, _emis_pert, _ray_pert) = opacf0_state.compute_opacity( + fr, t_pert, ne_pert, &npopul_pert, id, ¶ms.model.gffpar, ); - let abso_cm_pert = true_abs_pert + scat; + let abso_cm_pert = true_abs_pert + scat_pert; + + // d(abso)/dT = [abso(T+ΔT) - abso(T)] / ΔT dabt[id][ij] = (abso_cm_pert - abso_cm) / dt_pert; - // d(emis)/dT: emis = true_abs * Bν(T) per cm - let x = hkt * fr; + // d(emis)/dT: emis = true_abs * Bν(T) + // Uses full Bν at each temperature for consistency + let x = (hkt * fr).min(150.0); 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; + let x_pert = (hkt_pert * fr).min(150.0); + let bnu_pert = h_over_c2 * fr.powi(3) / (x_pert.exp() - 1.0).max(1e-100); + + let emis_base = true_abs * bnu; + let emis_pert = true_abs_pert * bnu_pert; + demt[id][ij] = (emis_pert - emis_base) / dt_pert; } } @@ -721,12 +906,11 @@ pub fn resolv( } } - // 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; - } + // FCOOL is set below in Part 5 (REINT/REDIF block) as: + // FCOOL(ID) = REINT(ID)*FCOOLI(ID) + // This captures the ALI (implicit) frequency contribution. + // bre_lte subtracts the explicit frequency part and adds + // the Jν coupling matrix elements. } // ============================================================== @@ -734,7 +918,7 @@ pub fn resolv( // ============================================================== // dP_rad/dm = (4π/c) Σ_ν κ_ν^true × (J_ν - B_ν) × w_ν // Lucy 用 pradt 来修正流体静力学平衡 - let mut pradt = vec![0.0; nd]; + pradt.fill(0.0); for id in 0..nd { let t = params.model.modpar.temp[id]; let hkt = HK / t; @@ -763,7 +947,7 @@ pub fn resolv( itlucy: config.itlucy, iaclt: 10, iacldt: 1, - ihecor: 1, + ihecor: 0, // Disable hydrostatic integration for LTE (EOS already gives correct dens/elec) lte: config.lte, lchc: config.lchc, ielcor: config.ielcor, @@ -796,7 +980,7 @@ pub fn resolv( lac2t: false, }; - let lucy_output = lucy_pure(&lucy_config, &lucy_model, &opacfl_data, &rad1_data); + let lucy_output = lucy(&lucy_config, &lucy_model, &opacfl_data, &rad1_data); // ============================================================== // Step 5: 更新模型状态 @@ -842,19 +1026,58 @@ pub fn resolv( } // ----------------------------------------------------------- - // Part 5: Rosseland 平均 — 设置 reint/redif (辐射平衡形式选择) + // Part 5: Rosseland 平均 — ROSSTD(0) // ----------------------------------------------------------- - // 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 { + // Fortran: IF(ITER.EQ.1.OR.LFIN) CALL ROSSTD(0) + // ROSSTD 计算 Rosseland 平均不透明度并设置 REINT/REDIF。 + { + let mut abrosd_mut = params.model.opmean.abrosd.clone(); + let mut taurs_mut = params.model.taurss.tross.clone(); + let mut reint_mut = params.model.repart.reint.clone(); + let mut redif_mut = params.model.repart.redif.clone(); + let pld_dummy = vec![0.0; nd]; + + let mut rosstd_params = RosstdEvaluateParams { + nd, + dens: ¶ms.model.modpar.dens[..nd], + deldm: &deldm, + dedm1: 0.0, + abrosd: &mut abrosd_mut, + abplad: ¶ms.model.opmean.abplad[..nd], + taurs: &mut taurs_mut, + reint: &mut reint_mut, + redif: &mut redif_mut, + taudiv: 10.0, + iter, + itndre: 0, + ndre: params.tlusty_config.matkey.ndre, + idlst: nd - 1, + teff, + lfin, + temp: ¶ms.model.modpar.temp[..nd], + elec: ¶ms.model.modpar.elec[..nd], + dm: ¶ms.model.modpar.dm[..nd], + pld: &pld_dummy, + }; + + let _rosstd_output = rosstd_evaluate(&mut rosstd_params); + + // 写回修改的数组 for id in 0..nd { - params.model.repart.reint[id] = 0.0; - params.model.repart.redif[id] = 0.0; + params.model.opmean.abrosd[id] = abrosd_mut[id]; + params.model.taurss.tross[id] = taurs_mut[id]; + params.model.repart.reint[id] = reint_mut[id]; + params.model.repart.redif[id] = redif_mut[id]; } - debug_log!("RESOLV: ROSSTD(0) called (all reint=0, redif=0 — Fortran default)"); } + // FCOOL(ID) = REINT(ID)*FCOOLI(ID) + for id in 0..nd { + params.model.totflx.fcool[id] = params.model.repart.reint[id] + * params.model.totflx.fcooli[id]; + } + debug_log!("RESOLV: REINT set to 1.0, FCOOL = REINT*FCOOLI"); + // 输出模型 — 对应 Fortran CALL OUTPUT (line 84) if let Some(w7) = params.writer7.as_mut() { let output_params = OutputParams { @@ -868,23 +1091,69 @@ pub fn resolv( // Part 6: 压力评估 // ----------------------------------------------------------- if iter <= config.nitzer { - // pzert(params.tlusty_config, params.atomic, params.model); + pzert(params.tlusty_config, params.atomic, params.model); debug_log!("RESOLV: PZERT called (iter <= nitzer={})", config.nitzer); } if (config.iheso6 != 0 || config.hmix0 > 0.0) && init == 1 { - // pzeval_pure(&mut PzevalParams { ... }); + // PZEVAL — 压力评估 + 对流梯度 + let vturb = vec![0.0f64; nd]; + let mut ptotal_mut = params.model.pressr.ptotal.clone(); + let mut pgs_mut = params.model.pressr.pgs.clone(); + let mut delta_mut = params.model.modcon.delta.clone(); + + let mut pzeval_params = PzevalParams { + config: PzevalConfig { + nd, + teff, + grav, + idisk: config.idisk, + hmix0: config.hmix0, + iconv: 0, + indl: 0, + ipnzev: 0, + iter, + lfin, + iconrs: config.iconrs, + iconre: config.iconre, + ipconf: config.ipconf, + }, + temp: ¶ms.model.modpar.temp[..nd], + dens: ¶ms.model.modpar.dens[..nd], + wmm: &wmm_arr, + elec: ¶ms.model.modpar.elec[..nd], + dm: ¶ms.model.modpar.dm[..nd], + vturb: &vturb, + abrosd: ¶ms.model.opmean.abrosd[..nd], + pradt: &pradt, + prd0, + ptotal: &mut ptotal_mut, + pgs: &mut pgs_mut, + delta: &mut delta_mut, + }; + + let _pzeval_output = pzeval(&mut pzeval_params); + + // 写回修改的数组 + for id in 0..nd { + params.model.pressr.ptotal[id] = ptotal_mut[id]; + params.model.pressr.pgs[id] = pgs_mut[id]; + params.model.modcon.delta[id] = delta_mut[id]; + } debug_log!("RESOLV: PZEVAL called"); } // ----------------------------------------------------------- // Part 7: 辐射压力 // ----------------------------------------------------------- - // radpre_pure(&RadpreParams { ... }); + // RADPRE — 计算辐射压力 + // 需要: RadpreConfig + RadpreModelState + &mut RadpreFreqParamsMut + &mut RadpreAliParamsMut + &mut RadpreOutputStateMut + // Radpre 内部遍历所有频率点，调用 opacf1/rtefr1 + // TODO: 构造嵌套参数结构体并调用 radpre() debug_log!("RESOLV: RADPRE called"); // 计时 - // timing(&TimingParams { iter_type: 1, iter }); + timing(&TimingParams { mode: TimingMode::FormalSolution, iter }); debug_log!("RESOLV: TIMING(1, {}) called", iter); // ----------------------------------------------------------- @@ -897,32 +1166,138 @@ pub fn resolv( }; if !(ipng == 0 && iter >= config.iacc && config.lres2) { - // 输出对流信息 + let zd = vec![0.0f64; nd]; + let vturb = vec![0.0f64; nd]; + + // --- CONREF (先执行，使用 clone/writeback) --- + if config.hmix0 > 0.0 && config.iconre > 0 && iter <= config.iconre && iter >= config.iconrs { + let mut temp_mut = params.model.modpar.temp.clone(); + let mut elec_mut = params.model.modpar.elec.clone(); + let mut flxc_mut = params.model.modcon.flxc.clone(); + let mut delta_mut = params.model.modcon.delta.clone(); + + let mut conref_params = ConrefParams { + nd, + teff, + config: ConrefConfig { + hmix0: config.hmix0, + iconv: 0, + indl: 0, + crflim: 0.1, + ideepc: 0, + ndcgap: 2, + idconz: 25, + ilgder: 0, + idisk: config.idisk, + ioptab: params.tlusty_config.basnum.ioptab, + icbegp: 0, + icendp: 0, + ifryb: config.ifryb, + iter, + imucon: 9999, + ipconf: config.ipconf, + grav, + qgrav: 0.0, + aconml: 1.0, + bconml: 1.0, + cconml: 1.0, + }, + temp: &mut temp_mut, + ptotal: ¶ms.model.pressr.ptotal[..nd], + pgs: ¶ms.model.pressr.pgs[..nd], + dens: ¶ms.model.modpar.dens[..nd], + vturb: &vturb, + abrosd: ¶ms.model.opmean.abrosd[..nd], + flrd: ¶ms.model.totflx.flrd[..nd], + flxc: &mut flxc_mut, + delta: &mut delta_mut, + wmm: &wmm_arr, + elec: &mut elec_mut, + taurs: ¶ms.model.taurss.tross[..nd], + thetav: &zd, + zd: &zd, + }; + let _conref_output = conref(&mut conref_params); + for id in 0..nd { + params.model.modpar.temp[id] = temp_mut[id]; + params.model.modpar.elec[id] = elec_mut[id]; + params.model.modcon.flxc[id] = flxc_mut[id]; + params.model.modcon.delta[id] = delta_mut[id]; + } + } + + // --- CONOUT (conref 之后，使用不可变借用 + 可变 clone) --- + let mut abrosd_mut = params.model.opmean.abrosd.clone(); + let mut flxc_mut = params.model.modcon.flxc.clone(); + let mut delta_mut = params.model.modcon.delta.clone(); + let mut reint_mut = params.model.repart.reint.clone(); + let mut redif_mut = params.model.repart.redif.clone(); + + let mut conout_params = ConoutParams { + imod: 1, + iprin: 0, + nd, + teff, + config: ConoutConfig { + hmix0: config.hmix0, + iconv: 0, + idisk: config.idisk, + ioptab: params.tlusty_config.basnum.ioptab, + ilgder: 0, + grav, + aconml: 1.0, + bconml: 1.0, + cconml: 1.0, + }, + temp: ¶ms.model.modpar.temp[..nd], + elec: ¶ms.model.modpar.elec[..nd], + dens: ¶ms.model.modpar.dens[..nd], + wmm: &wmm_arr, + dm: ¶ms.model.modpar.dm[..nd], + zd: &zd, + ptotal: ¶ms.model.pressr.ptotal[..nd], + pgs: ¶ms.model.pressr.pgs[..nd], + vturb: &vturb, + abrosd: &mut abrosd_mut, + flrd: ¶ms.model.totflx.flrd[..nd], + flxc: &mut flxc_mut, + delta: &mut delta_mut, + reint: &mut reint_mut, + redif: &mut redif_mut, + thetav: &zd, + qgrav: 0.0, + pradt: &pradt, + }; + if config.hmix0 == 0.0 { if let Some(w) = writer.as_mut() { - // WRITE(6,611) iter-1 let _ = w.write_raw(&format!("** CONVECTIVE FLUX: RESOLV; GLOBAL ITERATION ={:-3}\n", iter - 1)); - // call conout(1, ipconf) } + let _conout_output = conout(&mut conout_params); } else if config.hmix0 > 0.0 { - if config.iconre > 0 && iter <= config.iconre && iter >= config.iconrs { - // conref_pure(&mut ConrefParams { ... }); - } if config.ipconf > 0 || (config.ipconf == 0 && lfin) { if let Some(w) = writer.as_mut() { - // WRITE(6,611) iter-1 let _ = w.write_raw(&format!("** CONVECTIVE FLUX: RESOLV; GLOBAL ITERATION ={:-3}\n", iter - 1)); - // conout_pure(&mut ConoutParams { ... }); } + let _conout_output = conout(&mut conout_params); } } + + // 写回 conout 修改的数组 + for id in 0..nd { + params.model.opmean.abrosd[id] = abrosd_mut[id]; + params.model.modcon.flxc[id] = flxc_mut[id]; + params.model.modcon.delta[id] = delta_mut[id]; + params.model.repart.reint[id] = reint_mut[id]; + params.model.repart.redif[id] = redif_mut[id]; + } } // ----------------------------------------------------------- // Part 9: ALI 参数评估 // ----------------------------------------------------------- - // OPAINI(0) - // opaini(&OpainiParams { mode: 0, ... }); + // OPAINI(0) — 更新不透明度初始化参数 + opaini_full(0, params.tlusty_config, params.atomic, params.model, iter, 0, 1.0); debug_log!("RESOLV: OPAINI(0) called"); if config.icompt != 0 && ilam > 1 { @@ -936,17 +1311,22 @@ pub fn resolv( let use_kant = false; // kant(iter) == 1 || lfin if use_kant || lfin { - // ALISK2 - // alisk2_pure(...); + // ALISK2 — Kantorovich ALI 算法 + // 需要: Alisk2Config + Alisk2FreqParams + Alisk2AtomicParams + Alisk2ModelState + &mut Alisk2OutputState + // Alisk2OutputState 有 ~30 个 &mut [f64] 字段，来自 ModelState 各子结构 + // TODO: 构造 Alisk2OutputState（从 model 中借用 ~30 个 &mut 切片）并调用 alisk2() debug_log!("RESOLV: ALISK2 called (use_kant={} lfin={})", use_kant, lfin); } else { if config.irder == 0 { - // ALIST1 - // alist1_pure(...); + // ALIST1 — 标准 ALI（无导数） + // 需要: Alist1Config + Alist1FreqParams + Alist1AtomicParams + Alist1ModelState + &mut Alist1OutputState + // Alist1OutputState 有 ~40 个 &mut 字段，需从 ModelState 的 expraf/totflx/levpop/opmean/pressr 等借用 + // TODO: 构造 Alist1OutputState（~40 个 &mut 切片）并调用 alist1() debug_log!("RESOLV: ALIST1 called (irder=0)"); } else { - // ALIST2 - // alist2(...); + // ALIST2 — ALI with derivatives + // 需要: 类似 ALIST1 但增加温度导数相关字段 + // TODO: 构造 Alist2OutputState 并调用 alist2() debug_log!("RESOLV: ALIST2 called (irder={})", config.irder); } } @@ -956,10 +1336,147 @@ pub fn resolv( // ----------------------------------------------------------- if config.ifpopr == 2 { debug_log!("RESOLV: IFPOPR=2, updating populations"); + let natom = params.tlusty_config.basnum.natom as usize; + let nlevel = config.nlevel; + let nion = params.tlusty_config.basnum.nion as usize; + for id in 0..config.nd { - // steqeq_pure(&SteqeqParams { ... }, 1); + // 提取当前深度点的丰度和占据数 + let abund_id: Vec = (0..natom) + .map(|iat| params.atomic.atopar.abund[iat][id]) + .collect(); + let popul_id: Vec = (0..nlevel) + .map(|il| params.model.levpop.popul[il][id]) + .collect(); + + let t = params.model.modpar.temp[id]; + let elec = params.model.modpar.dens[id]; + let dens = params.model.modpar.dens[id]; + let wmm = params.tlusty_config.inppar.wmm[id]; + let an = dens / wmm + params.model.modpar.elec[id]; + + // STEQEQ — 统计平衡方程求解 + // 构造 SteqeqParams 和空操作的 callbacks + // 注: 空 callbacks 意味着 sabolf/ratmat/levsol 不被调用 + // 结果: pop1 直接从 popul 拷贝（无变化） + // TODO: 实现真正的 SteqeqCallbacks，调用 sabolf/ratmat/levsol + let matrix_a: Vec> = vec![vec![0.0; nlevel]; nlevel]; + let vector_b: Vec = vec![0.0; nlevel]; + let mut pop0: Vec = popul_id.clone(); + + // 提取 2D 数组在深度 id 处的列向量 + let wop_id: Vec = (0..nlevel).map(|i| params.model.wmcomp.wop[i][id]).collect(); + let sbpsi_id: Vec = (0..nlevel).map(|i| params.model.levref.sbpsi[i][id]).collect(); + let iltref_id: Vec = (0..nlevel).map(|i| params.model.levref.iltref[i][id]).collect(); + let ipzero_id: Vec = (0..nlevel).map(|i| params.model.popzr0.ipzero[i][id]).collect(); + let nrefs_id: Vec = (0..natom).map(|iat| params.atomic.atopar.nrefs[iat][id]).collect(); + + let steqeq_config = SteqeqConfig::default(); + let steqeq_params = SteqeqParams { + id: id + 1, + temp: t, + elec: params.model.modpar.elec[id], + dens, + wmm, + ytot: 1.1, + abund: &abund_id, + popul: &popul_id, + sbf: ¶ms.model.levpop.sbf, + wop: &wop_id, + sbpsi: &sbpsi_id, + iifor: ¶ms.atomic.levpar.iifor, + iltref: &iltref_id, + imodl: ¶ms.atomic.levpar.imodl, + iatm: ¶ms.atomic.levpar.iatm, + iifix: ¶ms.atomic.atopar.iifix, + ipzero: &ipzero_id, + n0a: ¶ms.atomic.atopar.n0a, + nka: ¶ms.atomic.atopar.nka, + nrefs: &nrefs_id, + ilk: ¶ms.atomic.levpar.ilk, + nfirst: ¶ms.atomic.ionpar.nfirst, + nlast: ¶ms.atomic.ionpar.nlast, + nnext: ¶ms.atomic.ionpar.nnext, + natom, + nlevel, + nion, + matrix_a: &matrix_a, + vector_b: &vector_b, + pop0: &mut pop0, + config: steqeq_config, + kant: ¶ms.tlusty_config.accel.kant, + }; + + // 使用空操作的 callbacks（sabolf/ratmat/levsol 不会被调用） + struct NoopCallbacks; + impl SteqeqCallbacks for NoopCallbacks {} + let steqeq_output = steqeq(&steqeq_params, 1, NoopCallbacks); + debug_log!("RESOLV: STEQEQ at depth {}", id); + + // NEWPOP — 使用 steqeq 输出的 pop1 更新占据数 + { + let mut newpop_params = NewpopParams { + config: params.tlusty_config, + atomic: params.atomic, + model: params.model, + id, + pop1: &steqeq_output.pop1, + }; + let _newpop_result = newpop(&mut newpop_params); + } + debug_log!("RESOLV: NEWPOP at depth {}", id); + + // ELCOR — 电子密度修正 if !config.lchc && iter < config.ielcor { - // elcor_pure(&ElcorParams { ... }); + let t = params.model.modpar.temp[id]; + let elec = params.model.modpar.elec[id]; + let dens = params.model.modpar.dens[id]; + let wmm = params.tlusty_config.inppar.wmm[id]; + + let elcor_config = ElcorConfig { + ifixde: 0, + ifmol: 0, + tmolim: 1e10, + ioptab: params.tlusty_config.basnum.ioptab, + iatref: 1, + max_iter: 10, + tolerance: 1e-6, + }; + let steqeq_config = SteqeqConfig::default(); + + let elcor_params = ElcorParams { + config: elcor_config, + id: id + 1, // 1-indexed + temp: t, + elec, + dens, + wmm, + ytot: 1.1, + abund: &abund_id, + qfix: 0.0, + natom, + iifix: ¶ms.atomic.atopar.iifix, + n0a: ¶ms.atomic.atopar.n0a, + nka: ¶ms.atomic.atopar.nka, + ilk: ¶ms.atomic.levpar.ilk, + iel: ¶ms.atomic.levpar.iel, + iz: ¶ms.atomic.ionpar.iz, + usum: ¶ms.model.levpop.usum, + imodl: ¶ms.atomic.levpar.imodl, + popul: &popul_id, + qadd: 0.0, + steqeq_config, + state_params: None, + moleq_params: None, + steqeq_params: None, + }; + let elcor_result = elcor(&elcor_params); + + // 写回修正后的电子密度和密度 + params.model.modpar.elec[id] = elcor_result.elec; + params.model.modpar.dens[id] = elcor_result.dens; + + debug_log!("RESOLV: ELCOR at depth {}, relane={:.4e}", id, elcor_result.relane); } } } @@ -976,25 +1493,126 @@ pub fn resolv( if config.ihecor >= -2 && config.izscal == 0 { if config.inzd > 0 || (config.idisk == 1 && config.ifryb > 0) { if config.iheso6 == 0 { - // PZEVLD - // pzevld(...); + // PZEVLD — 需要 ComputeArrays 和 DepthDeriv + // DepthDeriv 可以从 dm 数组计算 + let deriv = dmder(¶ms.model.modpar.dm, nd); + // ComputeArrays 需要从外部传入或创建临时实例 + // TODO: 将 ComputeArrays 添加到 ResolvParams 或从调用方传入 + // 目前使用默认的空 ComputeArrays（pzevld 仅使用 exprad.absoex） + let mut tmp_arrays = crate::tlusty::state::arrays::ComputeArrays::new(); + pzevld( + params.tlusty_config, + params.atomic, + params.model, + &mut tmp_arrays, + &deriv, + ); debug_log!("RESOLV: PZEVLD called"); } else { - // HESOL6 - // hesol6(&mut Hesol6Params { ... }); + // HESOL6 — clone/writeback 模式 + let nd = config.nd; + let nfreqe = config.nfreqe; + let teff = config.teff; + let iter = config.iter as usize; + let qgrav = params.tlusty_config.inppar.qgrav; + let bolk = BOLK; + let sig4p = SIG4P; + let pck = 4.19168946e-10; // PCK 常数 (4 * PI * K / C) + + let abrosd = params.model.opmean.abrosd[..nd].to_vec(); + let fprd = params.model.totflx.fprd[..nd].to_vec(); + let wmm = params.tlusty_config.inppar.wmm[..nd].to_vec(); + let pradt = params.model.pressr.pradt[..nd].to_vec(); + + // 可变字段需要 clone + writeback + let mut temp = params.model.modpar.temp[..nd].to_vec(); + let mut dm = params.model.modpar.dm[..nd].to_vec(); + let mut zd = params.model.modpar.zd[..nd].to_vec(); + let mut dens = params.model.modpar.dens[..nd].to_vec(); + let mut elec = params.model.modpar.elec[..nd].to_vec(); + let mut pgs = params.model.pressr.pgs[..nd].to_vec(); + let mut ptotal = params.model.pressr.ptotal[..nd].to_vec(); + + // 可选辐射场数据（浅拷贝/空） + let ijfr: Vec = params.model.freaux.ijfr[..nfreqe] + .iter().map(|&x| x as usize).collect(); + let lskip_flat: Vec = Vec::new(); // TODO: 需要展平的 lskip + let w_freq = ¶ms.model.frqall.w[..config.nfreq]; + let fh = ¶ms.model.surfac.fh[..config.nfreq]; + let radex: Vec = Vec::new(); // TODO: 需要 radex 数组 + let hextrd = ¶ms.model.totrad.hextrd[..config.nfreq]; + let absoe1: Vec = Vec::new(); // TODO: 需要 absoe1 数组 + + let grd = params.model.grdpra.grd[0]; + + let mut hesol6_params = Hesol6Params { + nd, + iter, + nfreqe, + teff, + sig4p, + pck, + qgrav, + bolk, + abrosd: &abrosd, + fprd: &fprd, + grd, + temp: &mut temp, + dm: &mut dm, + zd: &mut zd, + dens: &mut dens, + elec: &mut elec, + wmm: &wmm, + pradt: &pradt, + pgs: &mut pgs, + ptotal: &mut ptotal, + ijfr: &ijfr, + lskip: &lskip_flat, + w: w_freq, + fh, + radex: &radex, + hextrd, + absoe1: &absoe1, + }; + let _result = hesol6(&mut hesol6_params); + + // Writeback + for id in 0..nd { + params.model.modpar.temp[id] = temp[id]; + params.model.modpar.dm[id] = dm[id]; + params.model.modpar.zd[id] = zd[id]; + params.model.modpar.dens[id] = dens[id]; + params.model.modpar.elec[id] = elec[id]; + params.model.pressr.pgs[id] = pgs[id]; + params.model.pressr.ptotal[id] = ptotal[id]; + } debug_log!("RESOLV: HESOL6 called"); } } } if config.izscal == 1 { - // dmeval(&mut DmevalParams { ... }); + // DMEVAL — 需要 ComputeArrays 和 IterControl + // TODO: 需要从调用方传入 ComputeArrays + let mut tmp_arrays = crate::tlusty::state::arrays::ComputeArrays::new(); + let tmp_iter = crate::tlusty::state::iterat::IterControl::default(); + let mut dmeval_params = DmevalParams { + config: params.tlusty_config, + atomic: params.atomic, + model: params.model, + arrays: &mut tmp_arrays, + iter: &tmp_iter, + }; + // dmeval 内部修改 model 的 dm/dens/pressr + let _result = dmeval(&mut dmeval_params); debug_log!("RESOLV: DMEVAL called"); } if config.ifryb > 0 { - // rybheq(&RybheqParams { ... }); - debug_log!("RESOLV: RYBHEQ called"); + // RYBHEQ — 需要频率维度的辐射场数据 (rad1_all, fak1_all, abso1_all) + // 这些数据在 lambda 循环中计算，此处可能不可用 + // TODO: 需要从 lambda 循环中保存辐射场数据后才能调用 + debug_log!("RESOLV: RYBHEQ skipped (needs per-frequency radiation data)"); } // ----------------------------------------------------------- @@ -1043,44 +1661,56 @@ fn final_output( let config = ¶ms.config; if !config.lte { - // PRINC - 主输出 - // princ_pure(&PrincParams { ... }); + // PRINC - 主输出（跃迁速率、占据数等诊断信息） + // 需要: freq, rad, popul, rru, rrd, abso1, emis1, scat1 等频率维度数据 + // TODO: 这些数据需要从 lambda 循环中累积后传入 + debug_log!("RESOLV final_output: PRINC skipped (needs per-frequency data)"); } // OUTPRI - 输出模型 - // outpri_pure(&OutpriParams { ... }); + // TODO: outpri 需要 InpMod 格式写入 // COOLRT - 冷却速率输出 if config.icoolp != 0 || config.ipopac != 0 { - // coolrt_pure(&CoolrtParams { ... }); + // COOLRT 需要频率循环中的不透明度/辐射场数据 + // (abso1, emis1, rad1, absoti, emisti 等每个频率点的值) + // TODO: 需要从 lambda 循环中累积不透明度数据后传入 + debug_log!("RESOLV final_output: COOLRT skipped (needs per-frequency opacity data)"); } // RECHCK - 检查电荷守恒 - // rechck_pure(&RechckParams { ... }); + // 需要: w, abso1, rad1, emis1 等频率维度数据 + // TODO: 需要从 lambda 循环中累积后传入 // CHCKSE - 检查谱线平衡 if config.ichckp != 0 { - // chckse_pure(&ChckseParams { ... }); + let chckse_params = ChckseParams { + model: params.model, + atomic: params.atomic, + ioptab: params.tlusty_config.basnum.ioptab, + natom: params.tlusty_config.basnum.natom as usize, + n0a: ¶ms.atomic.atopar.n0a, + nka: ¶ms.atomic.atopar.nka, + ntrans: config.ntrans, + nlevel: config.nlevel, + }; + let _ = chckse(&chckse_params); } - // RTEINT - 强度计算 + // RTEINT - 强度计算（沿指定角度输出辐射强度） + // 需要: RteIntConfig + RteIntModelState + RteIntFreqParams + RteIntPhysics + + // &mut RteIntOpacity + RteIntFlux + &mut RteIntOutput + opacf1 回调 + // TODO: 构造嵌套参数结构体和 opacf1 闭包，调用 rteint() if config.intens > 0 { - // rteint(...); + debug_log!("RESOLV final_output: RTEINT skipped (needs nested params + callback)"); } // 康普顿散射最终处理 if config.icompt > 0 { - // RTECMU - // rtecmu(...); - - // OPAINI(0) - // opaini(&OpainiParams { mode: 0, ... }); - - // 循环所有频率点 - // for ij in 0..config.nfreq { - // opacf1(ij, ...); - // taufr1(&Taufr1Params { ... }); - // } + // RTECMU — 康普顿散射修正 + // 需要: 频率循环 + opacf1 + taufr1 + // TODO: 循环所有频率点调用 rtecmu + opacf1 + taufr1 + debug_log!("RESOLV final_output: RTECMU+OPAINI+freq loop skipped (icompt={})", config.icompt); } ResolvOutput { diff --git a/src/tlusty/io/srtfrq.rs b/src/tlusty/io/srtfrq.rs index 331d9f3..5dc9fcf 100644 --- a/src/tlusty/io/srtfrq.rs +++ b/src/tlusty/io/srtfrq.rs @@ -98,7 +98,7 @@ pub struct SrtfrqParams<'a> { /// # 返回值 /// /// 输出信息 -pub fn srtfrq_pure(params: &mut SrtfrqParams) -> SrtfrqOutput { +pub fn srtfrq(params: &mut SrtfrqParams) -> SrtfrqOutput { // 标记已导入的函数 let _ = (indexx, quit); diff --git a/src/tlusty/main.rs b/src/tlusty/main.rs index 195daa9..ad0c4b1 100644 --- a/src/tlusty/main.rs +++ b/src/tlusty/main.rs @@ -62,9 +62,9 @@ use super::io::{ input::InputParser, }; use super::math::solvers::{ - accel2_pure, Accel2Config, Accel2Params, - solve_pure, SolveConfig, DepthMatrices, - solves_pure, + accel2, Accel2Config, Accel2Params, + SolveConfig, DepthMatrices, + solves, solves_lte, KantorovichStorage, SolvesLteOutput, }; use super::math::io::{timing, TimingParams, TimingMode, reset_timer, output, OutputParams}; @@ -631,7 +631,7 @@ pub fn run_tlusty( psy2: &mut accel_psy2, psy3: &mut accel_psy3, }; - let accel_result = accel2_pure(&mut accel_params); + let accel_result = accel2(&mut accel_params); // If ACCEL2 accelerated, update model state from PSY0 // and recompute opacities via RESOLV (matching Fortran diff --git a/src/tlusty/math/ali/alisk1.rs b/src/tlusty/math/ali/alisk1.rs index 6007ea3..e08bbb2 100644 --- a/src/tlusty/math/ali/alisk1.rs +++ b/src/tlusty/math/ali/alisk1.rs @@ -236,7 +236,7 @@ pub struct Alisk1Output { /// /// 此函数是一个框架实现，实际调用 OPACF1、RTEFR1、ALIFRK、ROSSTD /// 需要在完整系统中实现。当前版本主要用于结构验证。 -pub fn alisk1_pure( +pub fn alisk1( config: &Alisk1Config, freq_params: &Alisk1FreqParams, atomic_params: &Alisk1AtomicParams, @@ -732,7 +732,7 @@ mod tests { rad1: &mut rad1, }; - let output = alisk1_pure(&config, &freq_params, &atomic_params, &model_state, &mut output_state); + let output = alisk1(&config, &freq_params, &atomic_params, &model_state, &mut output_state); assert!(output.computed); assert!(output.lross); // 因为 iter=1 且 ndre=0 @@ -869,7 +869,7 @@ mod tests { rad1: &mut rad1, }; - let output = alisk1_pure(&_config, &freq_params, &atomic_params, &model_state, &mut output_state); + let output = alisk1(&_config, &freq_params, &atomic_params, &model_state, &mut output_state); assert!(output.computed); } diff --git a/src/tlusty/math/ali/alisk2.rs b/src/tlusty/math/ali/alisk2.rs index b712342..b88bc38 100644 --- a/src/tlusty/math/ali/alisk2.rs +++ b/src/tlusty/math/ali/alisk2.rs @@ -235,7 +235,7 @@ pub struct Alisk2Output { /// # 返回值 /// /// 返回 `Alisk2Output`，包含计算结果信息。 -pub fn alisk2_pure( +pub fn alisk2( config: &Alisk2Config, freq_params: &Alisk2FreqParams, atomic_params: &Alisk2AtomicParams, @@ -822,7 +822,7 @@ mod tests { fcool: &mut fcool, }; - let output = alisk2_pure(&config, &freq_params, &atomic_params, &model_state, &mut output_state); + let output = alisk2(&config, &freq_params, &atomic_params, &model_state, &mut output_state); assert!(output.computed); assert!(output.lross); // 因为 iter=1 且 ndre=0 @@ -956,7 +956,7 @@ mod tests { fcool: &mut fcool, }; - let _ = alisk2_pure(&config, &freq_params, &atomic_params, &model_state, &mut output_state); + let _ = alisk2(&config, &freq_params, &atomic_params, &model_state, &mut output_state); // FLEXP 应该被初始化为零 for id in 0..nd { diff --git a/src/tlusty/math/ali/alist1.rs b/src/tlusty/math/ali/alist1.rs index e55a895..3ff57fb 100644 --- a/src/tlusty/math/ali/alist1.rs +++ b/src/tlusty/math/ali/alist1.rs @@ -306,7 +306,7 @@ pub struct Alist1Output { /// /// # 返回 /// - `Alist1Output`: 包含 prd0 和 prdx -pub fn alist1_pure( +pub fn alist1( config: &Alist1Config, freq_params: &Alist1FreqParams, atomic: &Alist1AtomicParams, diff --git a/src/tlusty/math/continuum/cia_h2h.rs b/src/tlusty/math/continuum/cia_h2h.rs deleted file mode 100644 index 439e115..0000000 --- a/src/tlusty/math/continuum/cia_h2h.rs +++ /dev/null @@ -1,136 +0,0 @@ -//! CIA H2-H 不透明度。 -//! -//! 重构自 TLUSTY `cia_h2h.f` -//! -//! # 功能 -//! -//! 计算 H2-H 碰撞诱导吸收 (CIA) 不透明度。 -//! 数据来源：TURBOSPEC - -/// Amagat 单位转换常数 (cm^-3) -const AMAGAT: f64 = 2.6867774e19; -/// 不透明度转换因子 -const FAC: f64 = 1.0 / (AMAGAT * AMAGAT); -/// 光速 (cm/s) -const CAS: f64 = 2.997925e10; -/// 温度表点数 -const NTEMP: usize = 4; -/// 频率/波长表行数 -const NLINES: usize = 67; - -/// CIA H2-H 温度表 -static TEMP_TABLE: [f64; NTEMP] = [1000.0, 1500.0, 2000.0, 2500.0]; - -/// CIA H2-H 不透明度表数据。 -pub struct CiaH2HData { - /// 频率数组 - freq: Vec, - /// 对数不透明度 alpha(freq, temp) - alpha: Vec>, - /// 是否已初始化 - loaded: bool, -} - -impl Default for CiaH2HData { - fn default() -> Self { - Self { - freq: Vec::new(), - alpha: Vec::new(), - loaded: false, - } - } -} - -impl CiaH2HData { - /// 加载 CIA 数据表。 - pub fn load(&mut self) { - if self.loaded { - return; - } - println!("Reading in H2-H CIA opacity tables..."); - - // 在完整实现中，这里会从 "./data/CIA_H2H.dat" 读取数据 - self.freq = vec![0.0; NLINES]; - self.alpha = vec![vec![0.0; NTEMP]; NLINES]; - self.loaded = true; - } - - /// 计算 CIA H2-H 不透明度。 - /// - /// # 参数 - /// - /// * `t` - 温度 (K) - /// * `ah2` - H2 数密度 - /// * `ah` - H 数密度 - /// * `ff` - 频率 (Hz) - pub fn opacity(&mut self, t: f64, ah2: f64, ah: f64, ff: f64) -> f64 { - // H2-H 数据仅适用于 T <= 2500 K - if t > 2500.0 { - return 0.0; - } - - self.load(); - - let f = ff / CAS; - - let j = locate(&TEMP_TABLE, t); - - if j == 0 { - println!(); - println!( - "Warning: requested temperature is below{:6.0} K", - TEMP_TABLE[0] - ); - println!("CIA H2-H CIA opacity set to zero"); - println!(); - return 0.0; - } - - let i = locate(&self.freq, f); - - let alp = if j == NTEMP { - let y1 = self.alpha[i][j - 1]; - let y2 = self.alpha[i + 1][j - 1]; - let tt = (f - self.freq[i]) / (self.freq[i + 1] - self.freq[i]); - (1.0 - tt) * y1 + tt * y2 - } else if i == 0 || i == NLINES { - -50.0 - } else { - let y1 = self.alpha[i][j - 1]; - let y2 = self.alpha[i + 1][j - 1]; - let y3 = self.alpha[i + 1][j]; - let y4 = self.alpha[i][j]; - - let tt = (f - self.freq[i]) / (self.freq[i + 1] - self.freq[i]); - let uu = (t - TEMP_TABLE[j - 1]) / (TEMP_TABLE[j] - TEMP_TABLE[j - 1]); - - (1.0 - tt) * (1.0 - uu) * y1 - + tt * (1.0 - uu) * y2 - + tt * uu * y3 - + (1.0 - tt) * uu * y4 - }; - - let alp = alp.exp(); - - FAC * ah2 * ah * alp - } -} - -/// 便捷函数：计算 CIA H2-H 不透明度（匹配 Fortran SUBROUTINE CIA_H2H 签名）。 -pub fn cia_h2h(t: f64, ah2: f64, ah: f64, ff: f64, opac: &mut f64) { - let mut data = CiaH2HData::default(); - *opac = data.opacity(t, ah2, ah, ff); -} - -/// 在有序数组中定位 x 的位置。 -fn locate(arr: &[f64], x: f64) -> usize { - if x < arr[0] { - return 0; - } - for i in 1..arr.len() { - if x < arr[i] { - return i; - } - } - arr.len() -} diff --git a/src/tlusty/math/continuum/cia_h2h2.rs b/src/tlusty/math/continuum/cia_h2h2.rs deleted file mode 100644 index 4d75148..0000000 --- a/src/tlusty/math/continuum/cia_h2h2.rs +++ /dev/null @@ -1,147 +0,0 @@ -//! CIA H2-H2 不透明度。 -//! -//! 重构自 TLUSTY `cia_h2h2.f` -//! -//! # 功能 -//! -//! 计算 H2-H2 碰撞诱导吸收 (CIA) 不透明度。 -//! 数据来源：Borysow A., Jorgensen U.G., Fu Y. 2001, JQSRT 68, 235 - -/// Amagat 单位转换常数 (cm^-3) -const AMAGAT: f64 = 2.6867774e19; -/// 不透明度转换因子 -const FAC: f64 = 1.0 / (AMAGAT * AMAGAT); -/// 光速 (cm/s) -const CAS: f64 = 2.997925e10; -/// 温度表点数 -const NTEMP: usize = 7; -/// 频率/波长表行数 -const NLINES: usize = 1000; - -/// CIA H2-H2 温度表 -static TEMP_TABLE: [f64; NTEMP] = [1000.0, 2000.0, 3000.0, 4000.0, 5000.0, 6000.0, 7000.0]; - -/// CIA H2-H2 不透明度表数据。 -/// -/// 在首次调用时从文件加载，之后缓存。 -pub struct CiaH2H2Data { - /// 频率数组 - freq: Vec, - /// 对数不透明度 alpha(freq, temp) - alpha: Vec>, - /// 是否已初始化 - loaded: bool, -} - -impl Default for CiaH2H2Data { - fn default() -> Self { - Self { - freq: Vec::new(), - alpha: Vec::new(), - loaded: false, - } - } -} - -impl CiaH2H2Data { - /// 加载 CIA 数据表。 - pub fn load(&mut self) { - if self.loaded { - return; - } - println!("Reading in H2-H2 CIA opacity tables..."); - - // 在完整实现中，这里会从 "./data/CIA_H2H2.dat" 读取数据 - // 暂时初始化为空表 - self.freq = vec![0.0; NLINES]; - self.alpha = vec![vec![0.0; NTEMP]; NLINES]; - self.loaded = true; - } - - /// 计算 CIA H2-H2 不透明度。 - /// - /// # 参数 - /// - /// * `t` - 温度 (K) - /// * `ah2` - H2 数密度 - /// * `ff` - 频率 (Hz) - /// - /// # 返回值 - /// - /// CIA 不透明度 - pub fn opacity(&mut self, t: f64, ah2: f64, ff: f64) -> f64 { - self.load(); - - // 输入频率为 Hz，但需要波数 (cm^-1) - let f = ff / CAS; - - // 在温度数组中定位 - let j = locate(&TEMP_TABLE, t); - - if j == 0 { - println!(); - println!( - "Warning: requested temperature is below{:6.0} K", - TEMP_TABLE[0] - ); - println!("CIA H2-H2 opacity set to 0"); - println!(); - return 0.0; - } - - // 在频率数组中定位 - let i = locate(&self.freq, f); - - let alp = if j == NTEMP { - // 高温端外推：保持恒定 - let y1 = self.alpha[i][j - 1]; - let y2 = self.alpha[i + 1][j - 1]; - let tt = (f - self.freq[i]) / (self.freq[i + 1] - self.freq[i]); - (1.0 - tt) * y1 + tt * y2 - } else if i == 0 || i == NLINES { - // 频率表外：设置非常小的值 - -50.0 - } else { - // 在表内双线性插值 - let y1 = self.alpha[i][j - 1]; - let y2 = self.alpha[i + 1][j - 1]; - let y3 = self.alpha[i + 1][j]; - let y4 = self.alpha[i][j]; - - let tt = (f - self.freq[i]) / (self.freq[i + 1] - self.freq[i]); - let uu = (t - TEMP_TABLE[j - 1]) / (TEMP_TABLE[j] - TEMP_TABLE[j - 1]); - - (1.0 - tt) * (1.0 - uu) * y1 - + tt * (1.0 - uu) * y2 - + tt * uu * y3 - + (1.0 - tt) * uu * y4 - }; - - let alp = alp.exp(); - - // 最终不透明度 - FAC * ah2 * ah2 * alp - } -} - -/// 便捷函数：计算 CIA H2-H2 不透明度（匹配 Fortran SUBROUTINE CIA_H2H2 签名）。 -pub fn cia_h2h2(t: f64, ah2: f64, ff: f64, opac: &mut f64) { - let mut data = CiaH2H2Data::default(); - *opac = data.opacity(t, ah2, ff); -} - -/// 在有序数组中定位 x 的位置。 -/// -/// 返回索引 j，使得 arr[j-1] <= x < arr[j]。 -/// 如果 x < arr[0]，返回 0。 -fn locate(arr: &[f64], x: f64) -> usize { - if x < arr[0] { - return 0; - } - for i in 1..arr.len() { - if x < arr[i] { - return i; - } - } - arr.len() -} diff --git a/src/tlusty/math/continuum/cia_h2he.rs b/src/tlusty/math/continuum/cia_h2he.rs deleted file mode 100644 index 2efd60a..0000000 --- a/src/tlusty/math/continuum/cia_h2he.rs +++ /dev/null @@ -1,132 +0,0 @@ -//! CIA H2-He 不透明度。 -//! -//! 重构自 TLUSTY `cia_h2he.f` -//! -//! # 功能 -//! -//! 计算 H2-He 碰撞诱导吸收 (CIA) 不透明度。 -//! 数据来源：Jorgensen U.G., Hammer D., Borysow A., Falkesgaard J., 2000, -//! Astronomy & Astrophysics 361, 283 - -/// Amagat 单位转换常数 (cm^-3) -const AMAGAT: f64 = 2.6867774e19; -/// 不透明度转换因子 -const FAC: f64 = 1.0 / (AMAGAT * AMAGAT); -/// 光速 (cm/s) -const CAS: f64 = 2.997925e10; -/// 温度表点数 -const NTEMP: usize = 7; -/// 频率/波长表行数 -const NLINES: usize = 242; - -/// CIA H2-He 温度表 -static TEMP_TABLE: [f64; NTEMP] = [1000.0, 2000.0, 3000.0, 4000.0, 5000.0, 6000.0, 7000.0]; - -/// CIA H2-He 不透明度表数据。 -pub struct CiaH2HeData { - /// 频率数组 - freq: Vec, - /// 对数不透明度 alpha(freq, temp) - alpha: Vec>, - /// 是否已初始化 - loaded: bool, -} - -impl Default for CiaH2HeData { - fn default() -> Self { - Self { - freq: Vec::new(), - alpha: Vec::new(), - loaded: false, - } - } -} - -impl CiaH2HeData { - /// 加载 CIA 数据表。 - pub fn load(&mut self) { - if self.loaded { - return; - } - println!("Reading in H2-He CIA opacity tables..."); - - // 在完整实现中，这里会从 "./data/CIA_H2He.dat" 读取数据 - self.freq = vec![0.0; NLINES]; - self.alpha = vec![vec![0.0; NTEMP]; NLINES]; - self.loaded = true; - } - - /// 计算 CIA H2-He 不透明度。 - /// - /// # 参数 - /// - /// * `t` - 温度 (K) - /// * `ah2` - H2 数密度 - /// * `ahe` - He 数密度 - /// * `ff` - 频率 (Hz) - pub fn opacity(&mut self, t: f64, ah2: f64, ahe: f64, ff: f64) -> f64 { - self.load(); - - let f = ff / CAS; - - let j = locate(&TEMP_TABLE, t); - - if j == 0 { - println!(); - println!( - "Warning: requested temperature is below{:6.0} K", - TEMP_TABLE[0] - ); - println!("CIA H2-He opacity set to 0"); - println!(); - return 0.0; - } - - let i = locate(&self.freq, f); - - let alp = if j == NTEMP { - let y1 = self.alpha[i][j - 1]; - let y2 = self.alpha[i + 1][j - 1]; - let tt = (f - self.freq[i]) / (self.freq[i + 1] - self.freq[i]); - (1.0 - tt) * y1 + tt * y2 - } else if i == 0 || i == NLINES { - -50.0 - } else { - let y1 = self.alpha[i][j - 1]; - let y2 = self.alpha[i + 1][j - 1]; - let y3 = self.alpha[i + 1][j]; - let y4 = self.alpha[i][j]; - - let tt = (f - self.freq[i]) / (self.freq[i + 1] - self.freq[i]); - let uu = (t - TEMP_TABLE[j - 1]) / (TEMP_TABLE[j] - TEMP_TABLE[j - 1]); - - (1.0 - tt) * (1.0 - uu) * y1 - + tt * (1.0 - uu) * y2 - + tt * uu * y3 - + (1.0 - tt) * uu * y4 - }; - - let alp = alp.exp(); - - FAC * ah2 * ahe * alp - } -} - -/// 便捷函数：计算 CIA H2-He 不透明度（匹配 Fortran SUBROUTINE CIA_H2HE 签名）。 -pub fn cia_h2he(t: f64, ah2: f64, ahe: f64, ff: f64, opac: &mut f64) { - let mut data = CiaH2HeData::default(); - *opac = data.opacity(t, ah2, ahe, ff); -} - -/// 在有序数组中定位 x 的位置。 -fn locate(arr: &[f64], x: f64) -> usize { - if x < arr[0] { - return 0; - } - for i in 1..arr.len() { - if x < arr[i] { - return i; - } - } - arr.len() -} diff --git a/src/tlusty/math/continuum/cia_hhe.rs b/src/tlusty/math/continuum/cia_hhe.rs deleted file mode 100644 index 8adc426..0000000 --- a/src/tlusty/math/continuum/cia_hhe.rs +++ /dev/null @@ -1,133 +0,0 @@ -//! CIA H-He 不透明度。 -//! -//! 重构自 TLUSTY `cia_hhe.f` -//! -//! # 功能 -//! -//! 计算 H-He 碰撞诱导吸收 (CIA) 不透明度。 -//! 数据来源：Gustafsson M., Frommhold, L. 2001, ApJ 546, 1168 - -/// Amagat 单位转换常数 (cm^-3) -const AMAGAT: f64 = 2.6867774e19; -/// 不透明度转换因子 -const FAC: f64 = 1.0 / (AMAGAT * AMAGAT); -/// 光速 (cm/s) -const CAS: f64 = 2.997925e10; -/// 温度表点数 -const NTEMP: usize = 11; -/// 频率/波长表行数 -const NLINES: usize = 43; - -/// CIA H-He 温度表 -static TEMP_TABLE: [f64; NTEMP] = [ - 1000.0, 1500.0, 2250.0, 3000.0, 4000.0, 5000.0, 6000.0, 7000.0, 8000.0, 9000.0, 10000.0, -]; - -/// CIA H-He 不透明度表数据。 -pub struct CiaHHeData { - /// 频率数组 - freq: Vec, - /// 对数不透明度 alpha(freq, temp) - alpha: Vec>, - /// 是否已初始化 - loaded: bool, -} - -impl Default for CiaHHeData { - fn default() -> Self { - Self { - freq: Vec::new(), - alpha: Vec::new(), - loaded: false, - } - } -} - -impl CiaHHeData { - /// 加载 CIA 数据表。 - pub fn load(&mut self) { - if self.loaded { - return; - } - println!("Reading in H-He CIA opacity tables..."); - - // 在完整实现中，这里会从 "./data/CIA_HHe.dat" 读取数据 - self.freq = vec![0.0; NLINES]; - self.alpha = vec![vec![0.0; NTEMP]; NLINES]; - self.loaded = true; - } - - /// 计算 CIA H-He 不透明度。 - /// - /// # 参数 - /// - /// * `t` - 温度 (K) - /// * `ah` - H 数密度 - /// * `ahe` - He 数密度 - /// * `ff` - 频率 (Hz) - pub fn opacity(&mut self, t: f64, ah: f64, ahe: f64, ff: f64) -> f64 { - self.load(); - - let f = ff / CAS; - - let j = locate(&TEMP_TABLE, t); - - if j == 0 { - println!(); - println!( - "Warning: requested temperature is below{:6.0} K", - TEMP_TABLE[0] - ); - println!("CIA H-He opacity set to 0"); - println!(); - return 0.0; - } - - let i = locate(&self.freq, f); - - let alp = if j == NTEMP { - let y1 = self.alpha[i][j - 1]; - let y2 = self.alpha[i + 1][j - 1]; - let tt = (f - self.freq[i]) / (self.freq[i + 1] - self.freq[i]); - (1.0 - tt) * y1 + tt * y2 - } else if i == 0 || i == NLINES { - -50.0 - } else { - let y1 = self.alpha[i][j - 1]; - let y2 = self.alpha[i + 1][j - 1]; - let y3 = self.alpha[i + 1][j]; - let y4 = self.alpha[i][j]; - - let tt = (f - self.freq[i]) / (self.freq[i + 1] - self.freq[i]); - let uu = (t - TEMP_TABLE[j - 1]) / (TEMP_TABLE[j] - TEMP_TABLE[j - 1]); - - (1.0 - tt) * (1.0 - uu) * y1 - + tt * (1.0 - uu) * y2 - + tt * uu * y3 - + (1.0 - tt) * uu * y4 - }; - - let alp = alp.exp(); - - FAC * ah * ahe * alp - } -} - -/// 便捷函数：计算 CIA H-He 不透明度（匹配 Fortran SUBROUTINE CIA_HHE 签名）。 -pub fn cia_hhe(t: f64, ah: f64, ahe: f64, ff: f64, opac: &mut f64) { - let mut data = CiaHHeData::default(); - *opac = data.opacity(t, ah, ahe, ff); -} - -/// 在有序数组中定位 x 的位置。 -fn locate(arr: &[f64], x: f64) -> usize { - if x < arr[0] { - return 0; - } - for i in 1..arr.len() { - if x < arr[i] { - return i; - } - } - arr.len() -} diff --git a/src/tlusty/math/continuum/mod.rs b/src/tlusty/math/continuum/mod.rs index 36c31f9..00a5e8c 100644 --- a/src/tlusty/math/continuum/mod.rs +++ b/src/tlusty/math/continuum/mod.rs @@ -1,9 +1,5 @@ //! continuum module -mod cia_h2h; -mod cia_h2h2; -mod cia_h2he; -mod cia_hhe; mod lte_opacity; mod opacf0; pub mod opacf0_state; @@ -23,12 +19,16 @@ mod opdata; mod opfrac; mod opacity_table; -pub use cia_h2h::CiaH2HData; -pub use cia_h2h2::CiaH2H2Data; -pub use cia_h2he::CiaH2HeData; -pub use cia_hhe::CiaHHeData; -// Re-export free functions from opacity module (the authoritative implementations) -pub use crate::tlusty::math::opacity::{cia_h2h, cia_h2h2, cia_h2he, cia_hhe}; +// CIA 不透明度：使用 opacity 模块中的权威实现，并提供向后兼容别名 +pub use crate::tlusty::math::opacity::{ + cia_h2h, cia_h2h2, cia_h2he, cia_hhe, + CiaH2hData, CiaH2h2Data, CiaH2heData, CiaHheData, +}; +// 大写后缀的旧版别名（向后兼容） +pub use crate::tlusty::math::opacity::CiaH2hData as CiaH2HData; +pub use crate::tlusty::math::opacity::CiaH2h2Data as CiaH2H2Data; +pub use crate::tlusty::math::opacity::CiaH2heData as CiaH2HeData; +pub use crate::tlusty::math::opacity::CiaHheData as CiaHHeData; pub use lte_opacity::{ LteOpacityParams, LteOpacityOutput, LteFrequencyGrid, lte_meanopt, generate_lte_frequency_grid, generate_inifrc_frequency_grid, diff --git a/src/tlusty/math/continuum/opactr.rs b/src/tlusty/math/continuum/opactr.rs index 47b4633..a2f1cc3 100644 --- a/src/tlusty/math/continuum/opactr.rs +++ b/src/tlusty/math/continuum/opactr.rs @@ -264,7 +264,7 @@ pub type PgsetFn = fn(usize); /// ... /// END /// ``` -pub fn opactr_pure( +pub fn opactr( ij: usize, config: &OpactrConfig, model_state: &mut OpactrModelState, diff --git a/src/tlusty/math/continuum/opfrac.rs b/src/tlusty/math/continuum/opfrac.rs index 20e9a17..855f3ec 100644 --- a/src/tlusty/math/continuum/opfrac.rs +++ b/src/tlusty/math/continuum/opfrac.rs @@ -142,7 +142,7 @@ pub struct OpfracOutput { /// /// # 返回 /// 配分函数和电离分数 -pub fn opfrac_pure(params: &OpfracParams, pfoptb: &PfOptB) -> OpfracOutput { +pub fn opfrac(params: &OpfracParams, pfoptb: &PfOptB) -> OpfracOutput { let iat = params.iat; // 如果 IAT == 0，返回默认值 @@ -345,7 +345,7 @@ mod tests { } #[test] - fn test_opfrac_pure_zero_iat() { + fn test_opfrac_zero_iat() { let pfoptb = PfOptB::new(); let params = OpfracParams { @@ -355,13 +355,13 @@ mod tests { ane: 1.0e12, }; - let result = opfrac_pure(¶ms, &pfoptb); + let result = opfrac(¶ms, &pfoptb); assert!((result.pf - 1.0).abs() < 1e-10); assert!((result.fra - 1.0).abs() < 1e-10); } #[test] - fn test_opfrac_pure_with_data() { + fn test_opfrac_with_data() { let mut pfoptb = PfOptB::new(); // 设置一些测试数据 @@ -386,7 +386,7 @@ mod tests { ane: 1.0e12, }; - let result = opfrac_pure(¶ms, &pfoptb); + let result = opfrac(¶ms, &pfoptb); assert!(result.pf > 0.0, "PF should be positive"); } diff --git a/src/tlusty/math/convection/concor.rs b/src/tlusty/math/convection/concor.rs index c71ee84..b000c9e 100644 --- a/src/tlusty/math/convection/concor.rs +++ b/src/tlusty/math/convection/concor.rs @@ -130,7 +130,7 @@ pub struct ConcorOutput { /// ... /// END /// ``` -pub fn concor_pure(params: &mut ConcorParams) -> ConcorOutput { +pub fn concor(params: &mut ConcorParams) -> ConcorOutput { // 检查是否需要执行 if params.config.indl == 0 { return ConcorOutput { @@ -313,7 +313,7 @@ mod tests { prd0: 0.0, }; - let output = concor_pure(&mut params); + let output = concor(&mut params); assert!(!output.computed); assert!(!output.temp_modified); @@ -348,7 +348,7 @@ mod tests { prd0: 0.0, }; - let output = concor_pure(&mut params); + let output = concor(&mut params); assert!(output.computed); assert!(output.temp_modified); @@ -385,7 +385,7 @@ mod tests { prd0: 0.0, }; - let output = concor_pure(&mut params); + let output = concor(&mut params); assert!(output.computed); // 在盘模式下，ptotal 会被重新计算 diff --git a/src/tlusty/math/convection/conout.rs b/src/tlusty/math/convection/conout.rs index b75d5f6..8e7041d 100644 --- a/src/tlusty/math/convection/conout.rs +++ b/src/tlusty/math/convection/conout.rs @@ -191,7 +191,7 @@ pub struct CubconData { /// ... /// END /// ``` -pub fn conout_pure(params: &mut ConoutParams) -> ConoutOutput { +pub fn conout(params: &mut ConoutParams) -> ConoutOutput { // f2r_depends: convec, meanop, meanopt, opacf0 (generic) let _ = (convec, meanop, meanopt); @@ -639,7 +639,7 @@ mod tests { #[test] fn test_conout_basic() { let mut params = TestParamsBuilder::new(50).build(); - let output = conout_pure(&mut params); + let output = conout(&mut params); // 验证基本输出 assert_eq!(output.depth_results.len(), 50); @@ -671,7 +671,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = conout_pure(&mut params); + let output = conout(&mut params); // 禁用对流时不应该有对流区 assert_eq!(output.icbeg, 0); @@ -686,7 +686,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = conout_pure(&mut params); + let output = conout(&mut params); // 验证基本功能 assert_eq!(output.depth_results.len(), 50); @@ -700,7 +700,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = conout_pure(&mut params); + let output = conout(&mut params); // 验证基本功能 assert_eq!(output.depth_results.len(), 50); @@ -714,7 +714,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = conout_pure(&mut params); + let output = conout(&mut params); // 盘模式应该正常工作 assert_eq!(output.depth_results.len(), 50); diff --git a/src/tlusty/math/convection/conref.rs b/src/tlusty/math/convection/conref.rs index 9d68abe..51f8e6c 100644 --- a/src/tlusty/math/convection/conref.rs +++ b/src/tlusty/math/convection/conref.rs @@ -340,7 +340,7 @@ fn call_convc1( /// ... /// END /// ``` -pub fn conref_pure(params: &mut ConrefParams) -> ConrefOutput { +pub fn conref(params: &mut ConrefParams) -> ConrefOutput { let nd = params.nd; let config = ¶ms.config.clone(); @@ -978,7 +978,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = conref_pure(&mut params); + let output = conref(&mut params); assert_eq!(output.icbeg, 0); assert_eq!(output.icend, 0); @@ -993,7 +993,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = conref_pure(&mut params); + let output = conref(&mut params); assert_eq!(output.icbeg, 0); } @@ -1006,7 +1006,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = conref_pure(&mut params); + let output = conref(&mut params); assert!(output.icbeg <= 50); assert!(output.icend <= 50); @@ -1059,7 +1059,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = conref_pure(&mut params); + let output = conref(&mut params); assert!(output.icbeg <= 50); } @@ -1073,7 +1073,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = conref_pure(&mut params); + let output = conref(&mut params); assert!(output.icbeg <= 50); } diff --git a/src/tlusty/math/convection/contmd.rs b/src/tlusty/math/convection/contmd.rs index fe4caab..9ea7cab 100644 --- a/src/tlusty/math/convection/contmd.rs +++ b/src/tlusty/math/convection/contmd.rs @@ -171,7 +171,7 @@ pub struct ContmdOutput { /// 计算盘模型对流层的温度 (CONTMD) - 无回调版本。 /// /// 用于纯计算模式，不更新电子密度和不透明度。 -pub fn contmd_pure(params: &mut ContmdParams) -> ContmdOutput { +pub fn contmd(params: &mut ContmdParams) -> ContmdOutput { contmd_impl(params, &mut NoOpCallbacks) } @@ -613,7 +613,7 @@ mod tests { #[test] fn test_contmd_basic() { let mut params = create_test_params(); - let output = contmd_pure(&mut params); + let output = contmd(&mut params); // 验证迭代次数在合理范围内 assert!(output.iconit <= params.config.nconit); @@ -631,7 +631,7 @@ mod tests { // 禁用对流 params.config.hmix0 = -1.0; - let output = contmd_pure(&mut params); + let output = contmd(&mut params); // 禁用对流时不应该有对流区 for &iconv in &output.iconv { @@ -646,7 +646,7 @@ mod tests { // 保存原始温度 let orig_temp = params.temp.to_vec(); - let _output = contmd_pure(&mut params); + let _output = contmd(&mut params); // 温度可能被更新 // 检查温度仍然是有限值 @@ -731,7 +731,7 @@ mod tests { hr1: Box::leak(Box::new(hr1)), }; - let output = contmd_pure(&mut params); + let output = contmd(&mut params); assert_eq!(output.iconv.len(), nd); assert!(output.iconit > 0); diff --git a/src/tlusty/math/eos/eldenc.rs b/src/tlusty/math/eos/eldenc.rs index 6d19bbb..bab5e1e 100644 --- a/src/tlusty/math/eos/eldenc.rs +++ b/src/tlusty/math/eos/eldenc.rs @@ -7,9 +7,9 @@ //! - 分析各元素对电子密度的贡献 //! - 输出电子供体信息 -use crate::tlusty::math::{moleq_pure, MoleqParams}; -use crate::tlusty::math::{rhonen_pure, RhonenParams}; -use crate::tlusty::math::{state_pure, StateParams}; +use crate::tlusty::math::{moleq, MoleqParams}; +use crate::tlusty::math::{rhonen, RhonenParams}; +use crate::tlusty::math::{state, StateParams}; use crate::tlusty::state::constants::{MDEPTH, MLEVEL}; // f2r_depends: RHONEN, STATE, MOLEQ @@ -130,7 +130,7 @@ pub struct EldencOutput { /// /// # 返回值 /// 包含插值电子密度和贡献者信息的输出 -pub fn eldenc_pure(params: &EldencParams) -> EldencOutput { +pub fn eldenc(params: &EldencParams) -> EldencOutput { let nd = params.nd; let config = ¶ms.config; @@ -168,7 +168,7 @@ pub fn eldenc_pure(params: &EldencParams) -> EldencOutput { // 计算 LTE 电子密度 if let Some(rhonen_params) = ¶ms.rhonen_params { - let rhonen_out = rhonen_pure(rhonen_params); + let rhonen_out = rhonen(rhonen_params); ane_lte[id] = rhonen_out.ane; } @@ -190,12 +190,12 @@ pub fn eldenc_pure(params: &EldencParams) -> EldencOutput { let rho = params.dens[id]; if let Some(rhonen_params) = ¶ms.rhonen_params { - let rhonen_out = rhonen_pure(rhonen_params); + let rhonen_out = rhonen(rhonen_params); let aein = rhonen_out.ane; // 调用 MOLEQ 计算分子平衡 if let Some(moleq_params) = ¶ms.moleq_params { - let moleq_out = moleq_pure(moleq_params); + let moleq_out = moleq(moleq_params); // Fortran: elcon(ia,id)=anion(ia,id)/elec(id) // moleq_out.anio0 包含离子数密度 for ia in 0..30.min(moleq_out.anio0.len()) { @@ -214,7 +214,7 @@ pub fn eldenc_pure(params: &EldencParams) -> EldencOutput { } else { // 使用 STATE if let Some(state_params) = ¶ms.state_params { - let _state_out = state_pure(state_params); + let _state_out = state(state_params); for ia in 0..30 { let iat = params.iatex[ia]; @@ -402,7 +402,7 @@ mod tests { moleq_params: None, }; - let output = eldenc_pure(¶ms); + let output = eldenc(¶ms); // 检查输出维度 assert_eq!(output.elecg.len(), 5); diff --git a/src/tlusty/math/eos/eldens.rs b/src/tlusty/math/eos/eldens.rs index f08c3f6..c819e87 100644 --- a/src/tlusty/math/eos/eldens.rs +++ b/src/tlusty/math/eos/eldens.rs @@ -9,9 +9,9 @@ //! - 计算内能和熵 use crate::tlusty::math::lineqs; -use crate::tlusty::math::{moleq_pure, MoleqParams, MoleculeEqData}; +use crate::tlusty::math::{moleq, MoleqParams, MoleculeEqData}; use crate::tlusty::math::{mpartf, MpartfResult}; -use crate::tlusty::math::{state_pure, StateParams}; +use crate::tlusty::math::{state, StateParams}; use crate::tlusty::math::{entene, EnteneParams, EnteneOutput}; use crate::tlusty::state::constants::{BOLK, HMASS, UN, TWO, HALF}; @@ -138,7 +138,7 @@ pub struct EldensOutput { /// /// # 返回值 /// 包含电子密度、内能、熵等的输出结构体 -pub fn eldens_pure(params: &EldensParams, ipri: i32) -> EldensOutput { +pub fn eldens(params: &EldensParams, ipri: i32) -> EldensOutput { // 检查 ioptab 标志 if params.config.ioptab < -1 { return EldensOutput { @@ -183,7 +183,7 @@ pub fn eldens_pure(params: &EldensParams, ipri: i32) -> EldensOutput { // 如果包含分子且温度低于分子温度上限，调用 MOLEQ if params.config.ifmol > 0 && t < params.config.tmolim { let aein = an * anerel; - // 调用 moleq_pure 计算分子平衡 + // 调用 moleq 计算分子平衡 let default_mol_data: MoleculeEqData = Default::default(); let moleq_params = MoleqParams { id: params.id, @@ -204,7 +204,7 @@ pub fn eldens_pure(params: &EldensParams, ipri: i32) -> EldensOutput { ifmol: params.config.ifmol, moltab: 0, }; - let moleq_output = moleq_pure(&moleq_params); + let moleq_output = moleq(&moleq_params); let ane = moleq_output.ane; anerel = ane / an; @@ -292,7 +292,7 @@ pub fn eldens_pure(params: &EldensParams, ipri: i32) -> EldensOutput { ifoppf: state_params.ifoppf, lrm: state_params.lrm, }; - let state_output = state_pure(&updated_params); + let state_output = state(&updated_params); q = state_output.q; dqn = state_output.dqn; } @@ -620,7 +620,7 @@ mod tests { anato_data: None, }; - let output = eldens_pure(¶ms, 0); + let output = eldens(¶ms, 0); // 验证输出 assert!(output.ane > 0.0); @@ -648,7 +648,7 @@ mod tests { anato_data: None, }; - let output = eldens_pure(¶ms, 0); + let output = eldens(¶ms, 0); // 低温时电子密度应该较低 assert!(output.anerel < 0.5); diff --git a/src/tlusty/math/eos/moleq.rs b/src/tlusty/math/eos/moleq.rs index 63c5716..4245705 100644 --- a/src/tlusty/math/eos/moleq.rs +++ b/src/tlusty/math/eos/moleq.rs @@ -193,7 +193,7 @@ pub fn parse_molecule_data(content: &str, exclude_large_carbon: bool) -> Molecul /// /// # 返回值 /// 包含电子密度、熵、内能、密度等的输出结构体 -pub fn moleq_pure(params: &MoleqParams) -> MoleqOutput { +pub fn moleq(params: &MoleqParams) -> MoleqOutput { // 如果不包含分子，直接返回空结果 if params.ifmol == 0 { return MoleqOutput { @@ -516,7 +516,7 @@ mod tests { moltab: 1, }; - let output = moleq_pure(¶ms); + let output = moleq(¶ms); // 当 ifmol=0 时，应返回初始估计 assert!((output.ane - params.aein).abs() < 1e-10); diff --git a/src/tlusty/math/eos/rhonen.rs b/src/tlusty/math/eos/rhonen.rs index ee6ee55..5ecd9df 100644 --- a/src/tlusty/math/eos/rhonen.rs +++ b/src/tlusty/math/eos/rhonen.rs @@ -6,7 +6,7 @@ //! - 从给定的温度和质量密度迭代求解总粒子密度和电子密度 //! - 使用 eldens 计算电子密度 -use crate::tlusty::math::{eldens_pure, EldensConfig, EldensOutput, EldensParams}; +use crate::tlusty::math::{eldens, EldensConfig, EldensOutput, EldensParams}; use crate::tlusty::state::constants::{HMASS, UN}; /// RHONEN 输入参数 @@ -58,7 +58,7 @@ pub struct RhonenOutput { /// /// # 返回值 /// 包含粒子密度、电子密度等的输出结构体 -pub fn rhonen_pure(params: &RhonenParams) -> RhonenOutput { +pub fn rhonen(params: &RhonenParams) -> RhonenOutput { let id = params.id; let t = params.t; let rho = params.rho; @@ -125,7 +125,7 @@ pub fn rhonen_pure(params: &RhonenParams) -> RhonenOutput { anato_data: None, }; - let eldens_output = eldens_pure(&eldens_params, 0); + let eldens_output = eldens(&eldens_params, 0); ane = eldens_output.ane; enrgi = eldens_output.energ; @@ -202,7 +202,7 @@ mod tests { eldens_wmy: 1.0, }; - let result = rhonen_pure(¶ms); + let result = rhonen(¶ms); // 验证基本物理约束 assert!(result.an > 0.0, "粒子密度应为正"); @@ -230,7 +230,7 @@ mod tests { eldens_wmy: 1.0, }; - let result = rhonen_pure(¶ms); + let result = rhonen(¶ms); // 验证基本物理约束 assert!(result.an > 0.0, "粒子密度应为正"); @@ -258,7 +258,7 @@ mod tests { eldens_wmy: 1.0, }; - let result = rhonen_pure(¶ms); + let result = rhonen(¶ms); assert!(result.an > 0.0); assert!(result.ane > 0.0); @@ -287,7 +287,7 @@ mod tests { eldens_wmy: 1.0, }; - let result = rhonen_pure(¶ms); + let result = rhonen(¶ms); assert!(result.ane >= 0.0, "温度 {} K 时电子密度应为非负", t); } } diff --git a/src/tlusty/math/hydrogen/bhe.rs b/src/tlusty/math/hydrogen/bhe.rs index a92352d..b67a452 100644 --- a/src/tlusty/math/hydrogen/bhe.rs +++ b/src/tlusty/math/hydrogen/bhe.rs @@ -5,6 +5,7 @@ //! 这些函数填充矩阵 A, B, C 的流体静力学平衡行 //! (NFREQE+INHE)-th row。 +use crate::tlusty::math::special::erfcx; use crate::tlusty::state::constants::{BOLK, HALF, MDEPTH, MFREQ, MLEVEL, MTOT, TWO, UN}; // ============================================================================ @@ -275,35 +276,6 @@ const PCK: f64 = 4.0 * std::f64::consts::PI / 3.0; // 辅助函数 // ============================================================================ -/// 计算 erfcx(x) = exp(x²) * erfc(x) -/// 使用 Fortran 代码中的系数 -fn erfcx(x: f64) -> f64 { - if x < 3.0 { - // 使用近似公式 - 8.86226925e-1 * (x * x).exp() * erfc_approx(x) - } else { - // 渐近展开 - HALF * (UN - HALF / (x * x)) / x - } -} - -/// erfc 近似 -fn erfc_approx(x: f64) -> f64 { - // 简单近似，精度足够 - if x < 0.0 { - 2.0 - erfc_approx(-x) - } else if x < 6.0 { - let t = 1.0 / (1.0 + 0.5 * x); - let poly = t * (0.17087277 + t - * (-0.82215223 + t - * (1.48851587 + t - * (-1.13520398 + t - * (0.27886807 + t * (-0.18628806 + t * 0.4206475)))))); - t * (-x * x).exp() * poly - } else { - 0.0 - } -} // ============================================================================ // BHE 主函数 diff --git a/src/tlusty/math/hydrogen/colhe.rs b/src/tlusty/math/hydrogen/colhe.rs index f527bd3..098e4d9 100644 --- a/src/tlusty/math/hydrogen/colhe.rs +++ b/src/tlusty/math/hydrogen/colhe.rs @@ -11,7 +11,7 @@ use crate::tlusty::math::collhe; use crate::tlusty::math::cspec; use crate::tlusty::math::irc; use crate::tlusty::state::constants::{EH, HK, MLEVEL, UN}; -// f2r_depends: COLLHE, CSPEC, IRC, CHEAV +use super::expi_approx; // ============================================================================ // 常量和数据 @@ -46,21 +46,7 @@ static A: [[f64; 10]; 6] = [ -0.72724497, 1.0879648, 5.6239786, 9.5323009, 16.150818], ]; -// ============================================================================ -// 辅助函数 -// ============================================================================ - -/// 计算指数积分 E1(x) 的近似值（Abramowitz-Stegun 公式）。 -fn expi_approx(u0: f64) -> f64 { - if u0 <= UN { - // 小参数展开 - -u0.ln() + (-0.57721566 + u0 * (0.99999193 + u0 * (-0.24991055 + u0 * (0.05519968 + u0 * (-0.00976004 + u0 * 0.00107857))))) - } else { - // 大参数渐近展开 - (-u0).exp() * ((0.2677734343 + u0 * (8.6347608925 + u0 * (18.059016973 + u0 * 8.5733287401))) - / (3.9584969228 + u0 * (21.0996530827 + u0 * (25.6329561486 + u0 * 9.5733223454)))) / u0 - } -} +// f2r_depends: COLLHE, CSPEC, IRC, CHEAV // ============================================================================ // 输入/输出结构体 diff --git a/src/tlusty/math/hydrogen/colis.rs b/src/tlusty/math/hydrogen/colis.rs index 9218903..bc4b93c 100644 --- a/src/tlusty/math/hydrogen/colis.rs +++ b/src/tlusty/math/hydrogen/colis.rs @@ -16,30 +16,14 @@ use crate::tlusty::math::cspec; use crate::tlusty::math::irc; use crate::tlusty::math::ylintp; use crate::tlusty::state::constants::{EH, HK, MLEVEL, TWO, UN}; +use super::expi_approx; + // f2r_depends: COLH, COLHE, CSPEC, IRC, CION, YLININTP // ============================================================================ // 常量 // ============================================================================ -/// 指数积分展开系数 -const EXPIA1: f64 = -0.57721566; -const EXPIA2: f64 = 0.99999193; -const EXPIA3: f64 = -0.24991055; -const EXPIA4: f64 = 0.05519968; -const EXPIA5: f64 = -0.00976004; -const EXPIA6: f64 = 0.00107857; - -const EXPIB1: f64 = 0.2677734343; -const EXPIB2: f64 = 8.6347608925; -const EXPIB3: f64 = 18.059016973; -const EXPIB4: f64 = 8.5733287401; - -const EXPIC1: f64 = 3.9584969228; -const EXPIC2: f64 = 21.0996530827; -const EXPIC3: f64 = 25.6329561486; -const EXPIC4: f64 = 9.5733223454; - /// 最大碰撞类型数 pub const MXTCOL: usize = 3; /// 最大碰撞拟合点数 @@ -79,16 +63,6 @@ fn cupsx(x: f64, u: f64, a: f64) -> f64 { 8.631e-6 / x * (-u).exp() * a } -/// 指数积分 E1 近似 -fn expi_approx(u0: f64) -> f64 { - if u0 <= UN { - -u0.ln() + EXPIA1 + u0 * (EXPIA2 + u0 * (EXPIA3 + u0 * (EXPIA4 + u0 * (EXPIA5 + u0 * EXPIA6)))) - } else { - (-u0).exp() * ((EXPIB1 + u0 * (EXPIB2 + u0 * (EXPIB3 + u0 * EXPIB4))) - / (EXPIC1 + u0 * (EXPIC2 + u0 * (EXPIC3 + u0 * EXPIC4)))) / u0 - } -} - // ============================================================================ // 输入/输出结构体 // ============================================================================ diff --git a/src/tlusty/math/hydrogen/hedif.rs b/src/tlusty/math/hydrogen/hedif.rs index 014e2e2..2511912 100644 --- a/src/tlusty/math/hydrogen/hedif.rs +++ b/src/tlusty/math/hydrogen/hedif.rs @@ -4,6 +4,7 @@ //! 计算深度相关的氦丰度剖面。 use crate::tlusty::io::{FortranWriter, Result}; +use crate::tlusty::math::solvers::raph; use crate::tlusty::state::atomic::AtomicData; use crate::tlusty::state::config::TlustyConfig; use crate::tlusty::state::constants::MDEPTH; @@ -19,15 +20,7 @@ const A2: f64 = 4.0; // 氦原子量 const BIGG: f64 = 6.6732e-8; // 引力常数 const PI: f64 = std::f64::consts::PI; -/// 计算扩散因子 RAPH。 -/// RAPH = diffusion factor for He/H separation -fn raph(gam: f64, z1: f64, z2: f64, a1: f64, a2: f64) -> f64 { - // 简化的扩散因子计算 - // 实际公式更复杂，这里使用近似 - let dz = z2 - z1; - let da = a2 - a1; - gam * (dz / da).abs().min(1.0) -} + /// HEDIF 参数结构体 pub struct HedifParams<'a> { diff --git a/src/tlusty/math/hydrogen/hesolv.rs b/src/tlusty/math/hydrogen/hesolv.rs index 5e73e7a..a24e2d3 100644 --- a/src/tlusty/math/hydrogen/hesolv.rs +++ b/src/tlusty/math/hydrogen/hesolv.rs @@ -9,7 +9,7 @@ use crate::tlusty::math::erfcx; use crate::tlusty::math::matinv; -use crate::tlusty::math::{rhonen_pure, RhonenParams}; +use crate::tlusty::math::{rhonen, RhonenParams}; use crate::tlusty::math::{steqeq_pure, SteqeqConfig, SteqeqParams}; use crate::tlusty::math::wnstor; use crate::tlusty::state::constants::{HALF, MDEPTH, TWO, UN}; @@ -216,7 +216,7 @@ pub struct HesolvOutput { /// /// # 返回值 /// 包含更新后的压力、密度、深度等的结果结构体 -pub fn hesolv_pure(params: &HesolvParams) -> HesolvOutput { +pub fn hesolv(params: &HesolvParams) -> HesolvOutput { let nd = params.model.nd; let mut p = vec![0.0; nd]; let mut vsnd2 = params.aux.vsnd2.clone(); @@ -492,7 +492,7 @@ pub fn hesolv_pure(params: &HesolvParams) -> HesolvOutput { eldens_wmy: params.config.eldens_wmy, }; - let rhonen_output = rhonen_pure(&rhonen_params); + let rhonen_output = rhonen(&rhonen_params); elec[id] = rhonen_output.ane; anerel = rhonen_output.anerel; @@ -685,7 +685,7 @@ mod tests { config, }; - let output = hesolv_pure(¶ms); + let output = hesolv(¶ms); // 验证输出 assert_eq!(output.ptotal.len(), nd); @@ -753,7 +753,7 @@ mod tests { config, }; - let output = hesolv_pure(¶ms); + let output = hesolv(¶ms); // 验证压力随深度增加 for i in 1..nd { @@ -832,7 +832,7 @@ mod tests { config, }; - let output = hesolv_pure(¶ms); + let output = hesolv(¶ms); // 所有密度应为正 for (i, &d) in output.dens.iter().enumerate() { diff --git a/src/tlusty/math/hydrogen/mod.rs b/src/tlusty/math/hydrogen/mod.rs index cc223c1..8bf8892 100644 --- a/src/tlusty/math/hydrogen/mod.rs +++ b/src/tlusty/math/hydrogen/mod.rs @@ -1,5 +1,8 @@ //! hydrogen module +mod utils; +pub use utils::expi_approx; + mod bhe; mod bre; mod brez; diff --git a/src/tlusty/math/hydrogen/sigave.rs b/src/tlusty/math/hydrogen/sigave.rs index e8ea0dc..282f2ea 100644 --- a/src/tlusty/math/hydrogen/sigave.rs +++ b/src/tlusty/math/hydrogen/sigave.rs @@ -212,7 +212,7 @@ fn interpolate_cross_section( /// SIGAVE 纯计算函数 /// /// 从读取的截面数据计算指定频率点的截面值 -pub fn sigave_pure( +pub fn sigave( transitions: &[(usize, Vec)], nfreqb: usize, freq: &[f64], @@ -355,7 +355,7 @@ fn sigave_impl(params: &SigaveParams, mut reader: FortranReader) }; // 计算截面 - sigave_pure( + sigave( &transitions, nfreqb, params.freq, @@ -513,7 +513,7 @@ mod tests { } #[test] - fn test_sigave_pure() { + fn test_sigave() { let transitions = vec![( 2, vec![ @@ -532,7 +532,7 @@ mod tests { let ifreqb: Vec = (0..50).map(|i| i as i32).collect(); let mut bfcs = vec![vec![0.0_f32; 100]; 10]; - sigave_pure(&transitions, 100, &freq, &ifreqb, 0, &mut bfcs); + sigave(&transitions, 100, &freq, &ifreqb, 0, &mut bfcs); // 检查截面是否被计算 assert!(bfcs[1].iter().any(|&x| x > 0.0)); diff --git a/src/tlusty/math/hydrogen/utils.rs b/src/tlusty/math/hydrogen/utils.rs new file mode 100644 index 0000000..2edbc89 --- /dev/null +++ b/src/tlusty/math/hydrogen/utils.rs @@ -0,0 +1,66 @@ +//! hydrogen 公共辅助函数 +//! +//! 提取自 colhe.rs 和 colis.rs 中的重复实现。 + +use crate::tlusty::state::constants::UN; + +// 指数积分展开系数（来自 Abramowitz-Stegun） +const EXPIA1: f64 = -0.57721566; +const EXPIA2: f64 = 0.99999193; +const EXPIA3: f64 = -0.24991055; +const EXPIA4: f64 = 0.05519968; +const EXPIA5: f64 = -0.00976004; +const EXPIA6: f64 = 0.00107857; + +const EXPIB1: f64 = 0.2677734343; +const EXPIB2: f64 = 8.6347608925; +const EXPIB3: f64 = 18.059016973; +const EXPIB4: f64 = 8.5733287401; + +const EXPIC1: f64 = 3.9584969228; +const EXPIC2: f64 = 21.0996530827; +const EXPIC3: f64 = 25.6329561486; +const EXPIC4: f64 = 9.5733223454; + +/// 指数积分 E1(x) 的近似值（Abramowitz-Stegun 公式 5.1.53 & 5.1.56）。 +/// +/// 精度约为 2×10⁻⁷。 +#[inline] +pub fn expi_approx(u0: f64) -> f64 { + if u0 <= UN { + // 小参数：级数展开 + -u0.ln() + + EXPIA1 + + u0 * (EXPIA2 + u0 * (EXPIA3 + u0 * (EXPIA4 + u0 * (EXPIA5 + u0 * EXPIA6)))) + } else { + // 大参数：有理近似（渐近展开） + (-u0).exp() + * ((EXPIB1 + u0 * (EXPIB2 + u0 * (EXPIB3 + u0 * EXPIB4))) + / (EXPIC1 + u0 * (EXPIC2 + u0 * (EXPIC3 + u0 * EXPIC4)))) + / u0 + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_expi_small() { + let result = expi_approx(0.5); + assert_relative_eq!(result, 0.5598, epsilon = 0.01); + } + + #[test] + fn test_expi_large() { + let result = expi_approx(5.0); + assert_relative_eq!(result, 0.001148, epsilon = 1e-4); + } + + #[test] + fn test_expi_unity() { + let result = expi_approx(1.0); + assert_relative_eq!(result, 0.2194, epsilon = 0.01); + } +} diff --git a/src/tlusty/math/io/princ.rs b/src/tlusty/math/io/princ.rs index 6627e72..58ad8e7 100644 --- a/src/tlusty/math/io/princ.rs +++ b/src/tlusty/math/io/princ.rs @@ -21,7 +21,7 @@ use crate::tlusty::state::constants::{UN, HALF, HK, BOLK, BN, MDEPTH, MFREQ, MLE const NPTR: usize = 30; use crate::tlusty::state::atomic::AtomicData; use crate::tlusty::state::model::ModelState; -use crate::tlusty::math::{sabolf_pure, SabolfParams, SabolfOutput}; +use crate::tlusty::math::{sabolf, SabolfParams, SabolfOutput}; use crate::tlusty::math::{linpro, LinproParams, LinproOutput}; use crate::tlusty::math::dwnfr; use crate::tlusty::math::cross; @@ -173,7 +173,7 @@ fn get_sbf(id: usize, model: &ModelState, atomic: &AtomicData) -> SabolfOutput { ioptab: 0, }; - sabolf_pure(&mut params) + sabolf(&mut params) } /// 获取谱线轮廓在指定频率点的值 @@ -235,7 +235,7 @@ fn get_bf_cross_section( /// /// # 返回 /// 各跃迁的详细信息 -pub fn princ_pure(params: &PrincParams) -> PrincOutput { +pub fn princ(params: &PrincParams) -> PrincOutput { let nd = params.model.modpar.temp.len(); let nct = params.nct; diff --git a/src/tlusty/math/io/prnt.rs b/src/tlusty/math/io/prnt.rs index dfb4e5e..026858b 100644 --- a/src/tlusty/math/io/prnt.rs +++ b/src/tlusty/math/io/prnt.rs @@ -12,7 +12,7 @@ use crate::tlusty::state::config::InpPar; use crate::tlusty::state::constants::HK; use crate::tlusty::state::model::{CraTes, LevPop, ModPar, MrgPar, RrRates, WmComp}; -use crate::tlusty::math::{sabolf_pure, SabolfParams}; +use crate::tlusty::math::{sabolf, SabolfParams}; // ============================================================================ // 输出结构体 @@ -77,7 +77,7 @@ pub struct PrntParams<'a> { /// /// # 返回 /// 各能级的速率平衡结果 -pub fn prnt_pure(params: &PrntParams) -> PrntOutput { +pub fn prnt(params: &PrntParams) -> PrntOutput { let mut balances = Vec::new(); let nd = params.modpar.temp.len(); @@ -118,7 +118,7 @@ pub fn prnt_pure(params: &PrntParams) -> PrntOutput { gmer: &mut gmer_work, ioptab: 0, }; - let sabolf_result = sabolf_pure(&mut sabolf_params); + let sabolf_result = sabolf(&mut sabolf_params); let sbf = &sabolf_result.sbf; let usum = &sabolf_result.usum; @@ -556,7 +556,7 @@ mod tests { ipop: &ipop, }; - let result = prnt_pure(¶ms); + let result = prnt(¶ms); // 由于测试数据是空的，结果应该为空或只有有限的结果 println!("Number of balances: {}", result.balances.len()); @@ -593,7 +593,7 @@ mod tests { ipop: &ipop, }; - let result = prnt_pure(¶ms); + let result = prnt(¶ms); // 验证结果 for balance in &result.balances { diff --git a/src/tlusty/math/io/pzeval.rs b/src/tlusty/math/io/pzeval.rs index 5e34360..ad700e3 100644 --- a/src/tlusty/math/io/pzeval.rs +++ b/src/tlusty/math/io/pzeval.rs @@ -164,7 +164,7 @@ pub struct PzevalOutput { /// ... /// END /// ``` -pub fn pzeval_pure(params: &mut PzevalParams) -> PzevalOutput { +pub fn pzeval(params: &mut PzevalParams) -> PzevalOutput { let nd = params.config.nd; let mut depth_results = Vec::with_capacity(nd); let mut conref_called = false; @@ -359,7 +359,7 @@ mod tests { #[test] fn test_pzeval_basic() { let mut params = TestParamsBuilder::new(50).build(); - let output = pzeval_pure(&mut params); + let output = pzeval(&mut params); assert_eq!(output.depth_results.len(), 50); assert!(!output.conref_called); @@ -378,7 +378,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = pzeval_pure(&mut params); + let output = pzeval(&mut params); assert_eq!(output.depth_results.len(), 50); } @@ -393,7 +393,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = pzeval_pure(&mut params); + let output = pzeval(&mut params); // iter=5 在 iconrs=3 和 iconre=10 之间，应该触发 CONREF assert!(output.conref_called); @@ -409,7 +409,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = pzeval_pure(&mut params); + let output = pzeval(&mut params); // ipnzev=0 AND lfin=true → CONOUT(1,1) assert!(output.conout_1); diff --git a/src/tlusty/math/io/rdata.rs b/src/tlusty/math/io/rdata.rs index 65d7338..fd492e0 100644 --- a/src/tlusty/math/io/rdata.rs +++ b/src/tlusty/math/io/rdata.rs @@ -771,7 +771,7 @@ pub fn compute_hydrogen_oscillator_strength( /// /// # 返回值 /// 处理后的能级数据数组 -pub fn process_levels_pure( +pub fn process_levels( params: &RdataParams, level_inputs: &[LevelInputData], ) -> Vec { @@ -800,7 +800,7 @@ pub fn process_levels_pure( } /// 处理连续跃迁数据（纯计算部分）。 -pub fn process_continua_pure( +pub fn process_continua( params: &RdataParams, inputs: &[ContinuumInputData], enion: &[f64], diff --git a/src/tlusty/math/io/rdatax.rs b/src/tlusty/math/io/rdatax.rs index bb54879..d69aa9b 100644 --- a/src/tlusty/math/io/rdatax.rs +++ b/src/tlusty/math/io/rdatax.rs @@ -367,7 +367,7 @@ pub fn compute_cross_sections( /// /// 这是一个简化版本，只处理纯计算部分。 /// 完整的 I/O 操作需要在外部处理。 -pub fn rdatax_pure( +pub fn rdatax( transitions: &[TransitionData], freq: &[f64], nfreq: usize, @@ -441,23 +441,23 @@ mod tests { } #[test] - fn test_rdatax_pure_basic() { + fn test_rdatax_basic() { let transitions = vec![create_test_transition()]; let freq = vec![1e15, 2e15, 3e15]; - let output = rdatax_pure(&transitions, &freq, 3, 0, 0); + let output = rdatax(&transitions, &freq, 3, 0, 0); assert_eq!(output.ntrx, 1); assert_eq!(output.bfcs.len(), 1); } #[test] - fn test_rdatax_pure_with_ibfint() { + fn test_rdatax_with_ibfint() { let transitions = vec![create_test_transition()]; let freq = vec![1e15, 2e15, 3e15]; // ibfint > 0 时使用 nfreqc - let output = rdatax_pure(&transitions, &freq, 3, 2, 1); + let output = rdatax(&transitions, &freq, 3, 2, 1); assert_eq!(output.bfcs[0].len(), 2); // 使用 nfreqc=2 } @@ -467,7 +467,7 @@ mod tests { let transitions: Vec = vec![]; let freq = vec![1e15, 2e15]; - let output = rdatax_pure(&transitions, &freq, 2, 0, 0); + let output = rdatax(&transitions, &freq, 2, 0, 0); assert_eq!(output.ntrx, 0); assert!(output.bfcs.is_empty()); diff --git a/src/tlusty/math/io/rechck.rs b/src/tlusty/math/io/rechck.rs index edeb073..049ede0 100644 --- a/src/tlusty/math/io/rechck.rs +++ b/src/tlusty/math/io/rechck.rs @@ -100,7 +100,7 @@ pub struct RechckOutput { /// ... /// END /// ``` -pub fn rechck_pure(params: &RechckParams) -> RechckOutput { +pub fn rechck(params: &RechckParams) -> RechckOutput { let nd = params.nd; let nfreq = params.nfreq; @@ -215,7 +215,7 @@ mod tests { emis1: &emis1, }; - let output = rechck_pure(¶ms); + let output = rechck(¶ms); assert_eq!(output.depth_results.len(), 3); @@ -251,7 +251,7 @@ mod tests { emis1: &emis1, }; - let output = rechck_pure(¶ms); + let output = rechck(¶ms); assert_eq!(output.depth_results.len(), 1); let result = &output.depth_results[0]; @@ -281,7 +281,7 @@ mod tests { emis1: &emis1, }; - let output = rechck_pure(¶ms); + let output = rechck(¶ms); // 当发射为零时，相对误差应为 0 assert_relative_eq!(output.depth_results[0].re, 0.0, epsilon = 1e-15); diff --git a/src/tlusty/math/opacity/inilam.rs b/src/tlusty/math/opacity/inilam.rs index c44ebaf..ea35b61 100644 --- a/src/tlusty/math/opacity/inilam.rs +++ b/src/tlusty/math/opacity/inilam.rs @@ -346,7 +346,7 @@ fn set_2d(arr: &mut [T], i: usize, j: usize, nrows: usize, val: T) { /// # 注意 /// 此函数仅实现核心计算逻辑，不包含 I/O 操作。 /// 外部函数调用（如 TDPINI、WNSTOR 等）需要通过回调实现。 -pub fn inilam_pure( +pub fn inilam( config: &InilamConfig, model: &mut InilamModelState, atomic: &InilamAtomicParams, @@ -825,7 +825,7 @@ mod tests { hkt1: &hkt1, }; - let output = inilam_pure(&config, &mut model, &atomic, &freq); + let output = inilam(&config, &mut model, &atomic, &freq); // INIT=1 分支应该返回 assert_relative_eq!(output.anerel, 0.5); // TEFF >= 8000 diff --git a/src/tlusty/math/opacity/inpdis.rs b/src/tlusty/math/opacity/inpdis.rs index c4b8900..b608f34 100644 --- a/src/tlusty/math/opacity/inpdis.rs +++ b/src/tlusty/math/opacity/inpdis.rs @@ -114,7 +114,7 @@ impl Default for InpDisParams { /// /// # 返回值 /// 返回计算结果 InpDisResult -pub fn inpdis_pure(params: &mut InpDisParams, cnu1: f64) -> InpDisResult { +pub fn inpdis(params: &mut InpDisParams, cnu1: f64) -> InpDisResult { let un = 1.0; // 处理 fractv 和 dmvisc @@ -340,7 +340,7 @@ pub fn inpdis_io( writer.write_newline()?; // 调用纯计算函数 - let result = inpdis_pure(params, cnu1); + let result = inpdis(params, cnu1); // 写入输出结果 writer.write_raw(&format!("TEFF ={:#10.0}", result.teff))?; @@ -411,7 +411,7 @@ mod tests { reynum: 0.0, }; - let result = inpdis_pure(&mut params, 1e15); + let result = inpdis(&mut params, 1e15); // 验证基本物理量 assert!(result.teff > 0.0, "TEFF should be positive"); @@ -445,7 +445,7 @@ mod tests { reynum: 0.0, }; - let result = inpdis_pure(&mut params, 1e15); + let result = inpdis(&mut params, 1e15); // GR 修正应该生效 assert!(result.teff > 0.0); @@ -470,7 +470,7 @@ mod tests { reynum: 0.0, }; - let result = inpdis_pure(&mut params, 1e15); + let result = inpdis(&mut params, 1e15); assert_relative_eq!(result.teff, 35000.0, epsilon = 1e-6); assert_relative_eq!(result.qgrav, 1e4, epsilon = 1e-6); @@ -493,7 +493,7 @@ mod tests { reynum: 0.0, }; - let result = inpdis_pure(&mut params, 0.0); + let result = inpdis(&mut params, 0.0); // fractv 应该被重新计算 assert!(params.fractv > 0.0); @@ -516,7 +516,7 @@ mod tests { }; let cnu1 = 2.0e15; - let result = inpdis_pure(&mut params, cnu1); + let result = inpdis(&mut params, cnu1); // FRLMAX = 1e11 * cnu1 * teff let expected_frlmax = 1e11 * cnu1 * 30000.0; diff --git a/src/tlusty/math/opacity/linsel.rs b/src/tlusty/math/opacity/linsel.rs index ce7bca1..3c5d9c5 100644 --- a/src/tlusty/math/opacity/linsel.rs +++ b/src/tlusty/math/opacity/linsel.rs @@ -209,7 +209,7 @@ pub struct LinselDebugLine { /// END /// ``` #[allow(clippy::too_many_arguments)] -pub fn linsel_pure( +pub fn linsel( config: &LinselConfig, atomic: &mut LinselAtomicParams, freq: &mut LinselFreqParams, @@ -942,7 +942,7 @@ mod tests { vec![1.0; 10] } - let output = linsel_pure(&config, &mut atomic, &mut freq_params, mock_opacf1, mock_rtefr1); + let output = linsel(&config, &mut atomic, &mut freq_params, mock_opacf1, mock_rtefr1); // ODF 模式：nppx = nfreq assert_eq!(output.stats.nppx, nfreq as i32); diff --git a/src/tlusty/math/opacity/profsp.rs b/src/tlusty/math/opacity/profsp.rs index 99ae7f6..db4704d 100644 --- a/src/tlusty/math/opacity/profsp.rs +++ b/src/tlusty/math/opacity/profsp.rs @@ -8,7 +8,7 @@ //! - 基于 Klaus Werner 公式 A.3.4 //! - 归一化到单位 -use crate::tlusty::math::{sabolf_pure, SabolfParams, SabolfOutput}; +use crate::tlusty::math::{sabolf, SabolfParams, SabolfOutput}; use crate::tlusty::math::ubeta; use crate::tlusty::math::voigt; use crate::tlusty::state::atomic::AtomicData; @@ -287,7 +287,7 @@ fn compute_zmikro(params: &ProfspParams, id: usize) -> f64 { ioptab: 0, }; - let sabolf_result = sabolf_pure(&mut sabolf_params); + let sabolf_result = sabolf(&mut sabolf_params); // 添加离子贡献 let nion = params.atomic.ionpar.iz.len(); diff --git a/src/tlusty/math/partition/partf.rs b/src/tlusty/math/partition/partf.rs index 348e246..ba947ef 100644 --- a/src/tlusty/math/partition/partf.rs +++ b/src/tlusty/math/partition/partf.rs @@ -126,7 +126,7 @@ fn init_arrays() { /// /// # 返回 /// 配分函数及其导数 -pub fn partf_pure(params: &PartfParams) -> PartfOutput { +pub fn partf(params: &PartfParams) -> PartfOutput { // 初始化辅助数组 init_arrays(); @@ -433,7 +433,7 @@ fn partf_user_defined(params: &PartfParams) -> PartfOutput { /// Opacity Project 配分函数 fn partf_opacity_project(params: &PartfParams) -> PartfOutput { // 调用 OPFRAC - use crate::tlusty::math::{opfrac_pure, OpfracParams, PfOptB}; + use crate::tlusty::math::{opfrac, OpfracParams, PfOptB}; let opfrac_params = OpfracParams { iat: params.iat, @@ -443,7 +443,7 @@ fn partf_opacity_project(params: &PartfParams) -> PartfOutput { }; let pfoptb = PfOptB::new(); - let result = opfrac_pure(&opfrac_params, &pfoptb); + let result = opfrac(&opfrac_params, &pfoptb); PartfOutput { u: result.pf, @@ -473,7 +473,7 @@ fn partf_ni(params: &PartfParams) -> PartfOutput { /// 重元素配分函数 (Z > 30) fn partf_heavy(params: &PartfParams) -> PartfOutput { // 调用 PFHEAV - use crate::tlusty::math::{pfheav_pure, PfheavParams}; + use crate::tlusty::math::{pfheav, PfheavParams}; let pfheav_params = PfheavParams { iiz: params.iat, @@ -483,7 +483,7 @@ fn partf_heavy(params: &PartfParams) -> PartfOutput { ane: params.ane, }; - let result = pfheav_pure(&pfheav_params); + let result = pfheav(&pfheav_params); PartfOutput { u: result.u, dut: 0.0, @@ -524,7 +524,7 @@ mod tests { iirwin: 0, }; - let result = partf_pure(¶ms); + let result = partf(¶ms); assert!(result.u > 0.0, "Partition function should be positive"); } @@ -540,7 +540,7 @@ mod tests { iirwin: 0, }; - let result = partf_pure(¶ms); + let result = partf(¶ms); assert!(result.u > 0.0, "Partition function should be positive"); } @@ -556,7 +556,7 @@ mod tests { iirwin: 0, }; - let result = partf_pure(¶ms); + let result = partf(¶ms); assert!(result.u > 0.0, "Iron partition function should be positive"); } @@ -572,7 +572,7 @@ mod tests { iirwin: 0, }; - let result = partf_pure(¶ms); + let result = partf(¶ms); assert!(result.u > 0.0, "Nickel partition function should be positive"); } @@ -589,7 +589,7 @@ mod tests { iirwin: 0, }; - let result = partf_pure(¶ms); + let result = partf(¶ms); assert!(result.u > 0.0, "High ionization partition function should be positive"); } @@ -605,7 +605,7 @@ mod tests { iirwin: 0, }; - let result = partf_pure(¶ms); + let result = partf(¶ms); assert!(result.u > 0.0, "Opacity Project partition function should be positive"); } diff --git a/src/tlusty/math/partition/pfheav.rs b/src/tlusty/math/partition/pfheav.rs index 1c32760..a00681a 100644 --- a/src/tlusty/math/partition/pfheav.rs +++ b/src/tlusty/math/partition/pfheav.rs @@ -125,7 +125,7 @@ pub struct PfheavOutput { /// /// # 返回 /// 配分函数值 -pub fn pfheav_pure(params: &PfheavParams) -> PfheavOutput { +pub fn pfheav(params: &PfheavParams) -> PfheavOutput { let iiz = params.iiz; let jnion = params.jnion; let mode = params.mode; @@ -272,7 +272,7 @@ mod tests { ane: 1.0e12, }; - let result = pfheav_pure(¶ms); + let result = pfheav(¶ms); assert!((result.u - 1.0).abs() < 1e-10); } @@ -287,7 +287,7 @@ mod tests { ane: 1.0e12, }; - let result = pfheav_pure(¶ms); + let result = pfheav(¶ms); assert!(result.u > 0.0, "Partition function should be positive"); } @@ -302,7 +302,7 @@ mod tests { ane: 1.0e12, }; - let result = pfheav_pure(¶ms); + let result = pfheav(¶ms); assert!(result.u > 0.0, "Partition function should be positive"); } @@ -317,7 +317,7 @@ mod tests { ane: 1.0e15, }; - let result = pfheav_pure(¶ms); + let result = pfheav(¶ms); assert!(result.u > 0.0, "Partition function should be positive"); } @@ -332,7 +332,7 @@ mod tests { ane: 1.0e10, }; - let result = pfheav_pure(¶ms); + let result = pfheav(¶ms); assert!(result.u > 0.0, "Partition function should be positive"); } diff --git a/src/tlusty/math/population/bpopc.rs b/src/tlusty/math/population/bpopc.rs index c34e732..9a1744a 100644 --- a/src/tlusty/math/population/bpopc.rs +++ b/src/tlusty/math/population/bpopc.rs @@ -15,7 +15,7 @@ //! - APM = (占据数向量) * (dAJ/dN) use crate::tlusty::state::constants::*; -use crate::tlusty::math::{state_pure, StateParams, StateOutput}; +use crate::tlusty::math::{state, StateParams, StateOutput}; // ============================================================================ // 输入参数结构体 @@ -170,7 +170,7 @@ pub struct BpopcOutput { /// /// # 返回 /// BPOPC 输出结构 -pub fn bpopc_pure(params: &BpopcParams) -> Option { +pub fn bpopc(params: &BpopcParams) -> Option { let id = params.id; // 计算索引 @@ -212,7 +212,7 @@ pub fn bpopc_pure(params: &BpopcParams) -> Option { lrm: params.state_lrm, }; - let state_out = state_pure(&state_params); + let state_out = state(&state_params); let qq = if params.ioptab > 0 { state_out.q / params.ytot @@ -411,7 +411,7 @@ mod tests { state_natoms: 2, }; - let result = bpopc_pure(¶ms); + let result = bpopc(¶ms); assert!(result.is_none(), "Should return None when INPC=0"); } @@ -468,7 +468,7 @@ mod tests { state_natoms: 2, }; - let result = bpopc_pure(¶ms); + let result = bpopc(¶ms); assert!(result.is_some(), "Should return Some when INPC>0"); let output = result.unwrap(); @@ -529,7 +529,7 @@ mod tests { state_natoms: 1, }; - let result = bpopc_pure(¶ms); + let result = bpopc(¶ms); assert!(result.is_some()); let output = result.unwrap(); diff --git a/src/tlusty/math/radiative/coolrt.rs b/src/tlusty/math/radiative/coolrt.rs index ea9bf14..c733e57 100644 --- a/src/tlusty/math/radiative/coolrt.rs +++ b/src/tlusty/math/radiative/coolrt.rs @@ -135,7 +135,7 @@ impl Default for CoolrtOutput { /// 2. 累积各离子的加热率: HTRAT += W * ABSOTI * RAD1 /// 3. 累积净辐射冷却率: CLHT2 += W * (EM - ABSO1 * RAD1) /// 4. 累积纯发射冷却率: CLHT3 += W * EMIS1 -pub fn coolrt_pure(params: &CoolrtParams) -> CoolrtOutput { +pub fn coolrt(params: &CoolrtParams) -> CoolrtOutput { let nd = params.nd; let nion = params.nion; let nfreq = params.nfreq; @@ -408,9 +408,9 @@ mod tests { } #[test] - fn test_coolrt_pure_basic() { + fn test_coolrt_basic() { let params = create_test_params(); - let result = coolrt_pure(¶ms); + let result = coolrt(¶ms); // 验证输出数组大小 assert_eq!(result.clht1.len(), params.nd); @@ -429,9 +429,9 @@ mod tests { } #[test] - fn test_coolrt_pure_clht1_sum() { + fn test_coolrt_clht1_sum() { let params = create_test_params(); - let result = coolrt_pure(¶ms); + let result = coolrt(¶ms); // CLHT1 应该是所有离子 (CLRAT - HTRAT) 的和 for id in 0..params.nd { diff --git a/src/tlusty/math/radiative/radpre.rs b/src/tlusty/math/radiative/radpre.rs index db64d2c..fbfc08c 100644 --- a/src/tlusty/math/radiative/radpre.rs +++ b/src/tlusty/math/radiative/radpre.rs @@ -199,7 +199,7 @@ pub struct RadpreOutput { /// # 返回值 /// RADPRE 输出结果 #[allow(clippy::too_many_arguments)] -pub fn radpre_pure( +pub fn radpre( config: &RadpreConfig, model: &RadpreModelState, freq: &mut RadpreFreqParamsMut, diff --git a/src/tlusty/math/radiative/trmder.rs b/src/tlusty/math/radiative/trmder.rs index 5971038..0549ade 100644 --- a/src/tlusty/math/radiative/trmder.rs +++ b/src/tlusty/math/radiative/trmder.rs @@ -15,7 +15,7 @@ //! 2. 调用 ELDENS 计算密度、能量和熵 //! 3. 使用有限差分计算导数 -use crate::tlusty::math::{eldens_pure, EldensConfig, EldensParams}; +use crate::tlusty::math::{eldens, EldensConfig, EldensParams}; use crate::tlusty::math::StateParams; use crate::tlusty::state::constants::{BOLK, HMASS, UN}; @@ -202,7 +202,7 @@ pub fn trmder(params: &TrmderParams) -> TrmderOutput { anato_data: None, }; - let result = eldens_pure(&eldens_params, 0); + let result = eldens(&eldens_params, 0); rhon[i] = result.rhoter; // 能量密度 = (1.5*P + 内能 + 3*Prad*(T'/T)^4) / rho @@ -377,7 +377,7 @@ mod tests { let result = trmder(¶ms); // 熵模式下应该也能正常工作 - // 注意：由于 eldens_pure 的简化实现，熵值可能为零或很小 + // 注意：由于 eldens 的简化实现，熵值可能为零或很小 // 这可能导致 heatcp 和 grdadb 不正确，所以我们只验证基本有限性 assert!(result.rho > 0.0, "密度应该为正"); assert!(result.heatcp.is_finite(), "定压比热应该是有限的"); diff --git a/src/tlusty/math/rates/rates1.rs b/src/tlusty/math/rates/rates1.rs index ee8c18c..9b9b0ab 100644 --- a/src/tlusty/math/rates/rates1.rs +++ b/src/tlusty/math/rates/rates1.rs @@ -241,7 +241,7 @@ pub struct Opacf1Result { /// /// # 返回值 /// 包含所有计算结果的结构体 -pub fn rates1_pure(params: &mut Rates1Params) -> Rates1Output { +pub fn rates1(params: &mut Rates1Params) -> Rates1Output { let nd = params.config.nd; let nfreq = params.config.nfreq; let ntrans = params.config.ntrans; @@ -665,7 +665,7 @@ mod tests { } #[test] - fn test_rates1_pure_basic() { + fn test_rates1_basic() { // 创建基本测试参数 let config = Rates1Config { imor: 0, @@ -780,7 +780,7 @@ mod tests { rosstd_contribute: None, }; - let output = rates1_pure(&mut params); + let output = rates1(&mut params); // 验证输出维度 assert_eq!(output.rru.len(), 1); diff --git a/src/tlusty/math/solvers/accel2.rs b/src/tlusty/math/solvers/accel2.rs index f25b0e6..de030dd 100644 --- a/src/tlusty/math/solvers/accel2.rs +++ b/src/tlusty/math/solvers/accel2.rs @@ -22,7 +22,7 @@ use crate::tlusty::state::constants::{MDEPTH, MFREQ, MLEVEL}; use crate::tlusty::state::model::ModelState; use crate::tlusty::state::config::TlustyConfig; use crate::tlusty::state::atomic::AtomicData; -use crate::tlusty::io::resolv_pure; +use crate::tlusty::io::resolv; // ============================================================================ // 配置参数 @@ -117,7 +117,7 @@ pub struct Accel2Output { /// /// # 注意 /// 当 `need_resolv` 为 true 时，调用者需要执行 RESOLV 重新计算 -pub fn accel2_pure(params: &mut Accel2Params) -> Accel2Output { +pub fn accel2(params: &mut Accel2Params) -> Accel2Output { // 标记 RESOLV 依赖（Fortran 在加速后调用 RESOLV） // 标记依赖（确保编译器知道 accel2 依赖 resolv_pure） let _ = "accel2 depends on resolv_pure"; @@ -287,7 +287,7 @@ pub fn accel2_io( params: &mut Accel2Params, writer: &mut W, ) -> Accel2Output { - let result = accel2_pure(params); + let result = accel2(params); if result.accelerated { let _ = writeln!(writer, " **** ACCEL2, ITER={}", params.config.iter); @@ -343,7 +343,7 @@ mod tests { psy3: &mut psy3, }; - let result = accel2_pure(&mut params); + let result = accel2(&mut params); assert!(!result.accelerated); assert!(!result.need_resolv); @@ -380,7 +380,7 @@ mod tests { psy3: &mut psy3, }; - let result = accel2_pure(&mut params); + let result = accel2(&mut params); // 应该存储到 PSY3，但不执行加速 assert!(!result.accelerated); @@ -426,7 +426,7 @@ mod tests { psy3: &mut psy3, }; - let result = accel2_pure(&mut params); + let result = accel2(&mut params); // 应该执行加速 assert!(result.accelerated); @@ -470,7 +470,7 @@ mod tests { psy3: &mut psy3, }; - let result = accel2_pure(&mut params); + let result = accel2(&mut params); // ipng != 0，应该只做移位，不执行加速 // iter - iacc = 5 - 3 = 2, iacd = 2, ipng = 2 % 2 = 0 @@ -514,7 +514,7 @@ mod tests { psy3: &mut psy3, }; - let result = accel2_pure(&mut params); + let result = accel2(&mut params); assert!(result.accelerated); diff --git a/src/tlusty/math/solvers/rybchn.rs b/src/tlusty/math/solvers/rybchn.rs index 802e562..ae21982 100644 --- a/src/tlusty/math/solvers/rybchn.rs +++ b/src/tlusty/math/solvers/rybchn.rs @@ -8,7 +8,7 @@ //! - 使用 ELDENS 计算电子密度 //! - 使用 PGSET 迭代计算气体压力 -use crate::tlusty::math::{eldens_pure, EldensConfig, EldensOutput, EldensParams}; +use crate::tlusty::math::{eldens, EldensConfig, EldensOutput, EldensParams}; use crate::tlusty::math::{pgset, PgsetParams, PgsetOutput}; use crate::tlusty::state::constants::{BOLK, HALF, TWO, UN}; @@ -142,7 +142,7 @@ pub struct RybchnOutput { /// /// # 返回值 /// 包含更新后的温度、密度、压力等的输出结构体 -pub fn rybchn_pure(params: &RybchnParams) -> RybchnOutput { +pub fn rybchn(params: &RybchnParams) -> RybchnOutput { let nd = params.nd; let config = ¶ms.config; @@ -400,7 +400,7 @@ fn compute_eldens(params: &RybchnParams, id: usize, t: f64, an: f64) -> EldensOu anato_data: None, }; - eldens_pure(&eldens_params, 1) + eldens(&eldens_params, 1) } /// ERFCX 近似函数（缩放互补误差函数） @@ -463,7 +463,7 @@ mod tests { #[test] fn test_rybchn_basic() { let params = create_test_params(); - let output = rybchn_pure(¶ms); + let output = rybchn(¶ms); // 检查输出维度 assert_eq!(output.temp.len(), params.nd); @@ -488,7 +488,7 @@ mod tests { // 设置较大的温度变化 params.changt = vec![0.5, 0.5, 0.5, 0.5, 0.5]; - let output = rybchn_pure(¶ms); + let output = rybchn(¶ms); // 检查温度变化被限制 // 由于 dplp = dpsilt - 1 = 2 - 1 = 1, dplm = 1/dpsilt - 1 = 0.5 - 1 = -0.5 diff --git a/src/tlusty/math/solvers/rybsol.rs b/src/tlusty/math/solvers/rybsol.rs index ee7ffdc..9275bff 100644 --- a/src/tlusty/math/solvers/rybsol.rs +++ b/src/tlusty/math/solvers/rybsol.rs @@ -263,7 +263,7 @@ pub struct RybsolOutput { /// ... /// END /// ``` -pub fn rybsol_pure( +pub fn rybsol( params: &mut RybsolParams, rybmtx: &mut RybmtxWork, // 回调函数 @@ -711,7 +711,7 @@ mod tests { let mut rybmtx = RybmtxWork::new(); // 使用空回调函数 - let output = rybsol_pure( + let output = rybsol( &mut params, &mut rybmtx, |_ij| {}, // opactr diff --git a/src/tlusty/math/solvers/solves.rs b/src/tlusty/math/solvers/solves.rs index 5bdefa4..ab69bf0 100644 --- a/src/tlusty/math/solvers/solves.rs +++ b/src/tlusty/math/solvers/solves.rs @@ -360,7 +360,7 @@ pub fn back_solution_step( /// 求解线性系统。 /// /// 完整版本需要集成 matgen、matinv、wnstor 等。 -pub fn solves_pure( +pub fn solves( nn: usize, nd: usize, nfreqe: usize, @@ -515,10 +515,10 @@ mod tests { } #[test] - fn test_solves_pure_basic() { + fn test_solves_basic() { let kant = vec![0; 200]; - let output = solves_pure( + let output = solves( 10, // nn 5, // nd 3, // nfreqe diff --git a/src/tlusty/math/temperature/elcor.rs b/src/tlusty/math/temperature/elcor.rs index 2e21504..ab1b4ff 100644 --- a/src/tlusty/math/temperature/elcor.rs +++ b/src/tlusty/math/temperature/elcor.rs @@ -7,8 +7,8 @@ //! - 使用迭代方法求解电荷守恒方程 //! - 与统计平衡方程耦合 -use crate::tlusty::math::{moleq_pure, MoleqOutput, MoleqParams}; -use crate::tlusty::math::{state_pure, StateOutput, StateParams}; +use crate::tlusty::math::{moleq, MoleqOutput, MoleqParams}; +use crate::tlusty::math::{state, StateOutput, StateParams}; use crate::tlusty::math::{steqeq_pure, SteqeqConfig, SteqeqParams}; use crate::tlusty::math::wnstor; use crate::tlusty::state::constants::{HALF, MDEPTH, MLEVEL, UN}; @@ -125,7 +125,7 @@ pub struct ElcorOutput { /// /// # 返回值 /// 包含更新后的电子密度和密度 -pub fn elcor_pure(params: &ElcorParams) -> ElcorOutput { +pub fn elcor(params: &ElcorParams) -> ElcorOutput { let config = ¶ms.config; // 检查 ioptab 标志 @@ -262,7 +262,7 @@ fn compute_qq(params: &ElcorParams, _id: usize, t: f64, an: f64, ane: f64, dens: if config.ifmol == 0 || t >= config.tmolim { // 使用 STATE 函数计算电荷 if let Some(state_params) = ¶ms.state_params { - let state_out = state_pure(state_params); + let state_out = state(state_params); let q = state_out.q; let qq = if config.ioptab > 0 { @@ -277,7 +277,7 @@ fn compute_qq(params: &ElcorParams, _id: usize, t: f64, an: f64, ane: f64, dens: // 使用 MOLEQ 分子平衡计算 // f2r_depends: MOLEQ if let Some(moleq_params) = ¶ms.moleq_params { - let _moleq_out = moleq_pure(moleq_params); + let _moleq_out = moleq(moleq_params); } params.qadd } @@ -348,7 +348,7 @@ mod tests { #[test] fn test_elcor_basic() { let params = create_test_params(); - let output = elcor_pure(¶ms); + let output = elcor(¶ms); // 检查输出 println!("ELEC: {}", output.elec); @@ -366,7 +366,7 @@ mod tests { let mut params = create_test_params(); params.config.ioptab = 1; // 应该跳过计算 - let output = elcor_pure(¶ms); + let output = elcor(¶ms); // 检查跳过后的输出等于输入 assert!((output.elec - params.elec).abs() < 1e-10); diff --git a/src/tlusty/math/temperature/greyd.rs b/src/tlusty/math/temperature/greyd.rs index 58f3ff0..0009b8a 100644 --- a/src/tlusty/math/temperature/greyd.rs +++ b/src/tlusty/math/temperature/greyd.rs @@ -250,7 +250,7 @@ where /// /// # 返回 /// 计算结果 -pub fn greyd_pure( +pub fn greyd( config: &GreydConfig, state: &mut GreydState, rhonen_fn: &mut F_Rhonen, @@ -385,7 +385,7 @@ mod tests { } #[test] - fn test_greyd_pure() { + fn test_greyd() { let config = GreydConfig { teff: 35000.0, qgrav: 1e4, @@ -420,7 +420,7 @@ mod tests { (1e-5, 1e-3) // 返回简化的不透明度 }; - let output = greyd_pure( + let output = greyd( &config, &mut state, &mut rhonen_fn, diff --git a/src/tlusty/math/temperature/lucy.rs b/src/tlusty/math/temperature/lucy.rs index 359fc7c..1f5354f 100644 --- a/src/tlusty/math/temperature/lucy.rs +++ b/src/tlusty/math/temperature/lucy.rs @@ -262,7 +262,7 @@ impl LucyState { /// ! 积分流体静力学平衡 /// END /// ``` -pub fn lucy_pure( +pub fn lucy( config: &LucyConfig, model: &LucyModelParams, opacfl_data: &[OpacflPointData], @@ -720,7 +720,7 @@ mod tests { ..config.clone() }; - let output = lucy_pure(&config_zero, &model, &opacfl_data, &rad1_data); + let output = lucy(&config_zero, &model, &opacfl_data, &rad1_data); assert_eq!(output.ilucy, 0); } diff --git a/src/tlusty/math/temperature/temcor.rs b/src/tlusty/math/temperature/temcor.rs index 25c6e3a..afa7052 100644 --- a/src/tlusty/math/temperature/temcor.rs +++ b/src/tlusty/math/temperature/temcor.rs @@ -164,7 +164,7 @@ pub struct CubconData { /// ... /// END /// ``` -pub fn temcor_pure(params: &mut TemcorParams) -> TemcorOutput { +pub fn temcor(params: &mut TemcorParams) -> TemcorOutput { let nd = params.nd; let mut depth_results = Vec::new(); let mut any_correction = false; @@ -526,7 +526,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = temcor_pure(&mut params); + let output = temcor(&mut params); // iconv <= 0 且 indl == 0 时应该跳过 assert_eq!(output.depth_results.len(), 0); @@ -541,7 +541,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = temcor_pure(&mut params); + let output = temcor(&mut params); // 应该处理所有深度点 assert_eq!(output.depth_results.len(), 49); // id from 2 to nd @@ -555,7 +555,7 @@ mod tests { ..Default::default() }; let mut params = TestParamsBuilder::new(50).config(config).build(); - let output = temcor_pure(&mut params); + let output = temcor(&mut params); // 盘模式应该正常工作 assert_eq!(output.depth_results.len(), 49); diff --git a/src/tlusty/math/temperature/temper.rs b/src/tlusty/math/temperature/temper.rs index f6ee484..fdceaaa 100644 --- a/src/tlusty/math/temperature/temper.rs +++ b/src/tlusty/math/temperature/temper.rs @@ -234,7 +234,7 @@ impl FactrsData { /// ... /// END /// ``` -pub fn temper_pure( +pub fn temper( params: &mut TemperParams, prsaux: &mut PrsauxData, flxaux: &mut FlxauxData, @@ -588,7 +588,7 @@ mod tests { flxaux.t4 = 35000.0_f64.powi(4); let factrs = FactrsData::new(50); - let output = temper_pure(&mut params, &mut prsaux, &mut flxaux, &factrs); + let output = temper(&mut params, &mut prsaux, &mut flxaux, &factrs); // 验证基本输出 assert!(output.t > 0.0); @@ -609,7 +609,7 @@ mod tests { let mut flxaux = FlxauxData::default(); let factrs = FactrsData::new(50); - let output = temper_pure(&mut params, &mut prsaux, &mut flxaux, &factrs); + let output = temper(&mut params, &mut prsaux, &mut flxaux, &factrs); // itgr > 1 且 itgmax > 0 时应该只计算电子密度 assert!(output.t > 0.0); @@ -626,7 +626,7 @@ mod tests { flxaux.t4 = 35000.0_f64.powi(4); let factrs = FactrsData::new(50); - let output = temper_pure(&mut params, &mut prsaux, &mut flxaux, &factrs); + let output = temper(&mut params, &mut prsaux, &mut flxaux, &factrs); // 第一个深度点应该设置 HG1, HR1, RR1 assert!(prsaux.hg1 >= 0.0); diff --git a/src/tlusty/math/utils/change.rs b/src/tlusty/math/utils/change.rs index ed16ed9..6a6e028 100644 --- a/src/tlusty/math/utils/change.rs +++ b/src/tlusty/math/utils/change.rs @@ -133,7 +133,7 @@ pub struct ChangeOutput { /// /// # 返回值 /// 包含更新后的粒子数 -pub fn change_pure(params: &ChangeParams) -> ChangeOutput { +pub fn change(params: &ChangeParams) -> ChangeOutput { let nd = params.nd; let nlevel = params.nlevel; let config = ¶ms.config; @@ -399,7 +399,7 @@ mod tests { steqeq_fn: None, }; - let output = change_pure(¶ms); + let output = change(¶ms); // 检查输出维度 assert_eq!(output.popul.len(), nlevel); diff --git a/src/tlusty/math/utils/newdm.rs b/src/tlusty/math/utils/newdm.rs index 33bd799..84ecf7d 100644 --- a/src/tlusty/math/utils/newdm.rs +++ b/src/tlusty/math/utils/newdm.rs @@ -204,7 +204,7 @@ impl Default for NewdmModelState { /// # 注意 /// /// 完整实现需要调用 TEMPER 和 HESOLV，这里只实现网格计算部分。 -pub fn newdm_pure( +pub fn newdm( config: &NewdmConfig, state: &mut NewdmModelState, _prs_aux: &mut PrsAux, @@ -545,7 +545,7 @@ mod tests { let mut prs_aux = PrsAux::default(); let mut factrs = Factrs::default(); - newdm_pure(&config, &mut state, &mut prs_aux, &mut factrs); + newdm(&config, &mut state, &mut prs_aux, &mut factrs); // 验证基本属性：dm[0] 应该是正数（因为它是从原始网格插值来的） // 注意：由于我们使用简化的测试网格，tauros 可能很小 @@ -584,7 +584,7 @@ mod tests { let mut prs_aux = PrsAux::default(); let mut factrs = Factrs::default(); - newdm_pure(&config, &mut state, &mut prs_aux, &mut factrs); + newdm(&config, &mut state, &mut prs_aux, &mut factrs); // 验证 theta 在有效范围内 // 注意：theta 的范围取决于 fractv 参数 diff --git a/src/tlusty/math/utils/newdmt.rs b/src/tlusty/math/utils/newdmt.rs index b342099..ee68906 100644 --- a/src/tlusty/math/utils/newdmt.rs +++ b/src/tlusty/math/utils/newdmt.rs @@ -177,7 +177,7 @@ impl Default for NewdmtModelState { /// # 注意 /// /// 完整实现需要调用 TEMPER 和 HESOLV，这里只实现网格计算部分。 -pub fn newdmt_pure( +pub fn newdmt( config: &NewdmtConfig, state: &mut NewdmtModelState, prs_aux: &mut NewdmtPrsAux, @@ -527,7 +527,7 @@ mod tests { let mut prs_aux = NewdmtPrsAux::default(); let mut factrs = NewdmtFactrs::default(); - newdmt_pure(&config, &mut state, &mut prs_aux, &mut factrs); + newdmt(&config, &mut state, &mut prs_aux, &mut factrs); // 验证基本属性 assert!(state.dm[0] > 0.0, "dm[0] should be positive"); @@ -568,7 +568,7 @@ mod tests { let mut prs_aux = NewdmtPrsAux::default(); let mut factrs = NewdmtFactrs::default(); - newdmt_pure(&config, &mut state, &mut prs_aux, &mut factrs); + newdmt(&config, &mut state, &mut prs_aux, &mut factrs); // 验证 theta 在有效范围内 for i in 0..config.nd { @@ -614,7 +614,7 @@ mod tests { let mut prs_aux = NewdmtPrsAux::default(); let mut factrs = NewdmtFactrs::default(); - newdmt_pure(&config, &mut state, &mut prs_aux, &mut factrs); + newdmt(&config, &mut state, &mut prs_aux, &mut factrs); // 验证 tauros 是单调递增的 for i in 1..config.nd { diff --git a/src/tlusty/math/utils/sabolf.rs b/src/tlusty/math/utils/sabolf.rs index 0162cf2..8ea6e0b 100644 --- a/src/tlusty/math/utils/sabolf.rs +++ b/src/tlusty/math/utils/sabolf.rs @@ -97,7 +97,7 @@ pub struct SabolfOutput { /// /// # 返回 /// Saha-Boltzmann 因子和上能级求和 -pub fn sabolf_pure(params: &mut SabolfParams) -> SabolfOutput { +pub fn sabolf(params: &mut SabolfParams) -> SabolfOutput { // 如果 ioptab < 0，返回空数组 if params.ioptab < 0 { let nlevels = params.enion.len(); @@ -269,7 +269,7 @@ pub fn sabolf_pure(params: &mut SabolfParams) -> SabolfOutput { }; let xmx = xmax * qz.sqrt(); - use crate::tlusty::math::partition::{partf_pure, PartfParams, PartfMode}; + use crate::tlusty::math::partition::{partf, PartfParams, PartfMode}; let partf_params = PartfParams { iat, izi: params.iz[ion_idx], @@ -279,7 +279,7 @@ pub fn sabolf_pure(params: &mut SabolfParams) -> SabolfOutput { mode: PartfMode::Standard, iirwin: 0, }; - let partf_result = partf_pure(&partf_params); + let partf_result = partf(&partf_params); let u_pf = partf_result.u; let dut_pf = partf_result.dut; let dun_pf = partf_result.dun; @@ -445,7 +445,7 @@ mod tests { #[test] fn test_sabolf_basic() { let mut params = create_test_params(); - let result = sabolf_pure(&mut params); + let result = sabolf(&mut params); assert_eq!(result.sbf.len(), 9); assert_eq!(result.usum.len(), 3); @@ -458,7 +458,7 @@ mod tests { fn test_sabolf_ioptab_negative() { let mut params = create_test_params(); params.ioptab = -1; - let result = sabolf_pure(&mut params); + let result = sabolf(&mut params); for &sb in &result.sbf { assert!((sb - 0.0).abs() < 1e-10); @@ -470,7 +470,7 @@ mod tests { let mut params = create_test_params(); params.t = 50000.0; params.ane = 1.0e15; - let result = sabolf_pure(&mut params); + let result = sabolf(&mut params); for &sb in &result.sbf { assert!(sb >= 0.0, "SBF should be non-negative"); diff --git a/src/tlusty/math/utils/state.rs b/src/tlusty/math/utils/state.rs index e9d5283..e4a7bed 100644 --- a/src/tlusty/math/utils/state.rs +++ b/src/tlusty/math/utils/state.rs @@ -10,7 +10,7 @@ //! - MODE=3: 类似 MODE=2，但计算导数 use crate::tlusty::state::constants::*; -use crate::tlusty::math::{partf_pure, PartfParams, PartfMode}; +use crate::tlusty::math::{partf, PartfParams, PartfMode}; // ============================================================================ // 常量 @@ -474,7 +474,7 @@ pub struct StateParams<'a> { /// /// # 返回 /// STATE 输出结构 -pub fn state_pure(params: &StateParams) -> StateOutput { +pub fn state(params: &StateParams) -> StateOutput { let mode = params.mode; let id = params.id; let t = params.t; @@ -497,7 +497,7 @@ pub fn state_pure(params: &StateParams) -> StateOutput { // 初始化 Opacity Project 离子分数（如果需要） if params.ifoppf > 0 { - use crate::tlusty::math::{opfrac_pure, OpfracParams, PfOptB}; + use crate::tlusty::math::{opfrac, OpfracParams, PfOptB}; let pfoptb = PfOptB::new(); let opfrac_params = OpfracParams { iat: 0, @@ -505,7 +505,7 @@ pub fn state_pure(params: &StateParams) -> StateOutput { t: params.t, ane: params.ane, }; - let _ = opfrac_pure(&opfrac_params, &pfoptb); + let _ = opfrac(&opfrac_params, &pfoptb); } return output; } @@ -570,7 +570,7 @@ pub fn state_pure(params: &StateParams) -> StateOutput { iirwin: 0, }; - let pf_result = partf_pure(&partf_params); + let pf_result = partf(&partf_params); let um = pf_result.u; let dutm = pf_result.dut; let dunm = pf_result.dun; @@ -605,7 +605,7 @@ pub fn state_pure(params: &StateParams) -> StateOutput { iirwin: 0, }; - let pf_result_j = partf_pure(&partf_params_j); + let pf_result_j = partf(&partf_params_j); let u = pf_result_j.u; let dut = pf_result_j.dut; let dun = pf_result_j.dun; @@ -862,7 +862,7 @@ mod tests { } #[test] - fn test_state_pure_basic() { + fn test_state_basic() { // 创建测试参数 let abndd = vec![1.0; MATOM]; // 丰度 let ioniz = vec![2; MATOM]; // 电离级 @@ -887,7 +887,7 @@ mod tests { ifoppf: 0, }; - let result = state_pure(¶ms); + let result = state(¶ms); // 验证结果 assert!(result.q >= 0.0, "Total charge should be non-negative"); @@ -924,7 +924,7 @@ mod tests { ifoppf: 0, }; - let result = state_pure(¶ms); + let result = state(¶ms); assert!((result.qm - qm_expected).abs() < 1e-20, "H- charge mismatch"); }