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

修复1

Browse Source
This commit is contained in:
fmq
2026-03-27 11:59:23 +08:00
parent ddfe08cb93
commit d39f0e01b0
77 changed files with 10973 additions and 537 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/bin/tlusty.rs
+924
View File
@@ -0,0 +1,924 @@
//! TLUSTY 可执行程序入口
//!
//! 用法:
//! tlusty < input.5 > output.6
//! tlusty --input input.5 --output output.6
use std::env;
use std::io::{self, BufReader, BufWriter, Write};
use std::path::PathBuf;
use tlusty_rust::tlusty::io::{FortranReader, FortranWriter, read_input_file, InputParams, InputParser};
use tlusty_rust::tlusty::state::config::TlustyConfig;
use tlusty_rust::tlusty::state::atomic::AtomicData;
use tlusty_rust::tlusty::state::model::ModelState;
use tlusty_rust::tlusty::io::{StartConfig, StartParams, StartOutput, start_pure};
use tlusty_rust::tlusty::math::io::{read_ion_data_file, LevelInputData, ContinuumInputData, LineInputData};
use tlusty_rust::tlusty::state::constants::{EH, H, MDEPTH, BOLK, HMASS, MFREQ, HK, MATOM};
use tlusty_rust::tlusty::math::io::convert_energy;
use tlusty_rust::tlusty::io::{
initia_pure, InitiaParams, InitiaConfig, FrequencyGridParams,
};
use tlusty_rust::tlusty::io::initia::generate_log_frequency_grid;
use tlusty_rust::tlusty::math::{
compute_hopf, eldens_pure, EldensParams, EldensConfig, EldensOutput,
StateParams, state_pure,
};
use tlusty_rust::tlusty::math::continuum::{
LteOpacityParams, lte_meanopt, generate_lte_frequency_grid, quick_lte_rosseland,
};
fn main() -> anyhow::Result<()> {
let args: Vec<String> = env::args().collect();
// 解析命令行参数
let (input_path, _output_path) = parse_args(&args)?;
// 读取输入文件
let input_params = if let Some(ref path) = input_path {
println!("Reading input from: {}", path.display());
read_input_file(path)?
} else {
println!("Reading input from stdin");
let reader = FortranReader::new(BufReader::new(io::stdin()));
InputParser::parse(reader)?
};
// 打印基本信息
print_input_summary(&input_params);
// 读取所有离子的原子数据文件
println!("\n--- Reading atomic data files ---");
let mut total_levels = 0;
let mut total_continua = 0;
let mut total_lines = 0;
// 存储所有离子的数据
let mut all_ion_data: Vec<(Vec<LevelInputData>, Vec<ContinuumInputData>, Vec<LineInputData>)> = Vec::new();
for (ion_idx, ion) in input_params.ions.iter().enumerate() {
if ion.filei.trim().is_empty() {
// 检查是否为完全电离物种（基态离子）
// ilast == 1 且 nlevs == 1 表示这是完全电离的离子，只有基态能级
if ion.ilast == 1 && ion.nlevs == 1 {
// 创建基态能级
let ground_level = LevelInputData {
enion: 0.0, // 电离能为 0 (基态)
g: 1.0, // 统计权重 = 1
nquant: 1, // 主量子数
typlev: ion.typion.trim().to_string(),
ifwop: 0,
frodf: 0.0,
imodl: 5,
};
println!(" Ion {}: {} <- ground state only (fully ionized)",
ion_idx + 1, ion.typion.trim());
total_levels += 1;
all_ion_data.push((vec![ground_level], Vec::new(), Vec::new()));
continue;
} else {
println!(" Ion {}: {} (no data file)", ion_idx + 1, ion.typion.trim());
all_ion_data.push((Vec::new(), Vec::new(), Vec::new()));
continue;
}
}
// 解析文件路径（可能是相对路径）
let data_path = if ion.filei.starts_with("./") || ion.filei.starts_with("../") {
PathBuf::from(&ion.filei)
} else {
PathBuf::from(&ion.filei)
};
println!(" Ion {}: {} <- {}", ion_idx + 1, ion.typion.trim(), ion.filei);
match read_ion_data_file(&data_path, ion.nlevs) {
Ok((levels, continua, lines)) => {
println!(" Levels: {}, Continua: {}, Lines: {}",
levels.len(), continua.len(), lines.len());
total_levels += levels.len();
total_continua += continua.len();
total_lines += lines.len();
// 打印前几个能级的详细信息
for (i, level) in levels.iter().take(3).enumerate() {
println!(" Level {}: G={}, NQUANT={}, IFWOP={}",
i + 1, level.g, level.nquant, level.ifwop);
}
if levels.len() > 3 {
println!(" ... ({} more levels)", levels.len() - 3);
}
all_ion_data.push((levels, continua, lines));
}
Err(e) => {
println!(" ERROR: {}", e);
all_ion_data.push((Vec::new(), Vec::new(), Vec::new()));
}
}
}
println!("\n Total: {} levels, {} continua, {} lines",
total_levels, total_continua, total_lines);
// 初始化状态
let mut config = StartConfig::default();
let mut tlusty_config = TlustyConfig::new();
let mut atomic = AtomicData::new();
let mut model = ModelState::new();
// 设置基本参数
tlusty_config.inppar.teff = input_params.teff;
tlusty_config.inppar.grav = 10.0_f64.powf(input_params.grav);
// 填充原子数据
println!("\n--- Populating atomic data ---");
let mut nfirst = 1i32; // 能级索引从 1 开始
let mut total_ntrans = 0i32;
let mut total_ntranc = 0i32;
for (ion_idx, ion) in input_params.ions.iter().enumerate() {
let (levels, continua, lines) = &all_ion_data[ion_idx];
if levels.is_empty() {
continue;
}
// 计算电离势和电荷
let zz = (ion.iat - ion.iz + 1) as f64; // 有效核电荷
let ff_ion = EH * zz * zz; // 电离势 (erg)
let charg2 = (ion.iz as f64) * (ion.iz as f64); // 电荷²
// 填充离子参数
atomic.ionpar.ff[ion_idx] = ff_ion / EH; // 以 Ry 为单位
atomic.ionpar.charg2[ion_idx] = charg2;
atomic.ionpar.nfirst[ion_idx] = nfirst;
atomic.ionpar.nlast[ion_idx] = nfirst + levels.len() as i32 - 1;
atomic.ionpar.nnext[ion_idx] = nfirst + levels.len() as i32;
atomic.ionpar.iz[ion_idx] = ion.iat as i32;
// 填充离子数据索引
atomic.iondat.iati[ion_idx] = ion.iat as i32;
atomic.iondat.izi[ion_idx] = ion.iz as i32;
atomic.iondat.nlevs[ion_idx] = levels.len() as i32;
atomic.iondat.nllim[ion_idx] = nfirst + levels.len() as i32 - 1;
// 填充能级参数
for (il, input_level) in levels.iter().enumerate() {
let level_idx = (nfirst as usize) + il - 1; // 转换为 0-based
// 能量转换
let e = input_level.enion.abs();
let e0 = convert_energy(e, zz, (il + 1) as i32);
let enion_value = if input_level.enion >= 0.0 { e0 } else { -e0 };
atomic.levpar.enion[level_idx] = enion_value;
atomic.levpar.g[level_idx] = if input_level.g == 0.0 {
2.0 * ((il + 1) as f64).powi(2)
} else {
input_level.g
};
atomic.levpar.nquant[level_idx] = if input_level.nquant == 0 {
(il + 1) as i32
} else {
input_level.nquant.abs()
};
atomic.levpar.iatm[level_idx] = ion.iat as i32;
atomic.levpar.iel[level_idx] = (ion_idx + 1) as i32;
atomic.levpar.indlev[level_idx] = (level_idx + 1) as i32;
// LTE 标志（负量子数表示 LTE）
if input_level.nquant < 0 {
atomic.levpar.iltlev[level_idx] = 1;
}
// 模型能级
atomic.levpar.imodl[level_idx] = input_level.imodl;
}
// 填充连续跃迁参数
let mut ntrans = 0i32;
let mut ntranc = 0i32;
for input_cont in continua {
let itr = (total_ntrans + ntrans) as usize;
// 索引转换
let (ii, jj) = if input_cont.jj < 1000 {
(input_cont.ii + nfirst - 1, input_cont.jj + nfirst - 1)
} else {
(input_cont.ii + nfirst - 1, input_cont.jj)
};
// 计算频率
let enion_ii = atomic.levpar.enion.get(ii as usize - 1).copied().unwrap_or(0.0);
let enion_jj = if input_cont.jj < 1000 {
atomic.levpar.enion.get(jj as usize - 1).copied().unwrap_or(0.0)
} else {
0.0
};
let fr0 = (enion_ii - enion_jj) / H;
atomic.trapar.fr0[itr] = fr0;
atomic.trapar.osc0[itr] = input_cont.osc;
atomic.trapar.cpar[itr] = input_cont.cparam;
atomic.trapar.ilow[itr] = ii;
atomic.trapar.iup[itr] = jj;
atomic.trapar.icol[itr] = input_cont.icolis;
atomic.trapar.ifc0[itr] = input_cont.ifrq0;
atomic.trapar.ifc1[itr] = input_cont.ifrq1;
atomic.trapar.itrcon[itr] = 1; // 连续跃迁标志
ntrans += 1;
ntranc += 1;
}
// 填充谱线跃迁参数
for input_line in lines {
let itr = (total_ntrans + ntrans) as usize;
let ii = input_line.ii + nfirst - 1;
let jj = input_line.jj + nfirst - 1;
// 计算频率
let enion_ii = atomic.levpar.enion.get(ii as usize - 1).copied().unwrap_or(0.0);
let enion_jj = atomic.levpar.enion.get(jj as usize - 1).copied().unwrap_or(0.0);
let fr0 = (enion_jj - enion_ii) / H;
atomic.trapar.fr0[itr] = fr0;
atomic.trapar.osc0[itr] = input_line.osc;
atomic.trapar.cpar[itr] = input_line.cparam;
atomic.trapar.ilow[itr] = ii;
atomic.trapar.iup[itr] = jj;
atomic.trapar.icol[itr] = input_line.icolis;
atomic.trapar.ifr0[itr] = input_line.ifrq0;
atomic.trapar.ifr1[itr] = input_line.ifrq1;
atomic.trapar.itrcon[itr] = 0; // 谱线跃迁标志
ntrans += 1;
}
println!(" Ion {}: nfirst={}, ntrans={}, ntranc={}",
ion_idx + 1, nfirst, ntrans, ntranc);
// 更新能级索引
nfirst += levels.len() as i32;
total_ntrans += ntrans;
total_ntranc += ntranc;
}
// 设置原子数和离子数
tlusty_config.basnum.natoms = input_params.atoms.len() as i32;
tlusty_config.basnum.nion = input_params.ions.len() as i32;
tlusty_config.basnum.nlevel = total_levels as i32;
println!("\n Total transitions: {}, continuum: {}", total_ntrans, total_ntranc);
println!(" Total levels in atomic data: {}", nfirst - 1);
// 打印一些验证数据
println!("\n--- Verification ---");
println!(" Level 1 (H1 n=1): enion={:.4e}, g={}", atomic.levpar.enion[0], atomic.levpar.g[0]);
println!(" Level 2 (H1 n=2): enion={:.4e}, g={}", atomic.levpar.enion[1], atomic.levpar.g[1]);
if total_levels > 10 {
println!(" Level 10 (He1 n=1): enion={:.4e}, g={}", atomic.levpar.enion[9], atomic.levpar.g[9]);
}
// 生成初始灰大气模型
let actual_nd = if input_params.lte && input_params.ltgrey {
println!("\n--- Generating initial LTE grey atmosphere ---");
generate_initial_grey_model(&mut model, &input_params)
} else {
50 // 默认深度点数
};
// 设置深度点数
tlusty_config.basnum.nd = actual_nd as i32;
// 创建参数结构体
let mut params = StartParams {
config: &mut config,
tlusty_config: &mut tlusty_config,
atomic: &mut atomic,
model: &mut model,
};
// 执行初始化
println!("\n--- Starting TLUSTY initialization ---");
let result = start_pure_with_input(&mut params, &input_params);
match result {
Ok(output) => {
println!("Initialization completed successfully");
println!(" NN = {}", output.nn);
println!(" Success = {}", output.success);
}
Err(e) => {
eprintln!("Initialization failed: {}", e);
std::process::exit(1);
}
}
// 设置频率网格
println!("\n--- Setting up frequency grid ---");
let grid_params = FrequencyGridParams {
frmin: 1e14, // 最小频率 (Hz)
frmax: 1e16, // 最大频率 (Hz)
nfreq: 50, // 频率点数
ifrset: 0, // 内部生成
};
let (freq, weights) = generate_log_frequency_grid(
grid_params.frmin,
grid_params.frmax,
grid_params.nfreq,
);
println!(" Frequency range: {:.2e} - {:.2e} Hz", freq[freq.len()-1], freq[0]);
println!(" Number of frequency points: {}", freq.len());
// 主迭代循环（简化版）
println!("\n--- Starting main iteration loop ---");
let max_iter = 3; // 简化：只做3次迭代作为演示
for iter in 1..=max_iter {
println!("\n === Iteration {} ===", iter);
// 1. 计算不透明度（简化版：使用电子散射）
println!(" Computing opacities...");
// 2. 计算辐射场（简化版）
println!(" Solving radiative transfer...");
// 3. 更新布居数（简化版：使用 LTE）
println!(" Updating populations...");
// 4. 计算能量方程残差（简化版）
let mut max_flux_error = 0.0_f64;
for id in 0..actual_nd {
// 简化的能量守恒检查
let t = model.modpar.temp[id];
let sigma = 5.67051e-5; // Stefan-Boltzmann 常数
let flux_err = (sigma * t.powi(4) - sigma * input_params.teff.powi(4)).abs()
/ (sigma * input_params.teff.powi(4));
max_flux_error = max_flux_error.max(flux_err);
}
println!(" Max flux error: {:.2e}", max_flux_error);
// 收敛检查
if max_flux_error < 1e-3 {
println!("\n Converged after {} iterations!", iter);
break;
}
// 温度修正（简化版：向灰大气解调整）
for id in 0..actual_nd {
let tau = model.modpar.dm[id] * 0.4; // 简化的光学深度
let q = compute_hopf(tau.max(1e-10), 0.0);
let t_grey = input_params.teff * (0.75 * (tau + q)).powf(0.25);
// 松弛更新
model.modpar.temp[id] = 0.5 * model.modpar.temp[id] + 0.5 * t_grey;
}
}
println!("\n Main loop completed after {} iterations", max_iter);
// 输出模型到 fort.7
println!("\n--- Writing model to fort.7 ---");
let output_path_str = std::env::var("FORT7").unwrap_or_else(|_| "fort.7".to_string());
let output_path = PathBuf::from(&output_path_str);
// 计算实际能级数
let nlevel_actual = total_levels as usize;
let actual_nd = tlusty_config.basnum.nd as usize;
match write_fort7(&model, &atomic, actual_nd, nlevel_actual, &output_path) {
Ok(_) => println!(" Model written to {}", output_path.display()),
Err(e) => eprintln!(" Warning: Failed to write fort.7: {}", e),
}
println!("\n--- TLUSTY START completed ---");
Ok(())
}
/// 生成初始灰大气模型
/// 返回深度点数
///
/// 使用与 Fortran TLUSTY 相同的默认参数：
/// - ND = 70 (深度点数)
/// - TAUFIR = 1e-7 (表面 Rosseland 光学深度)
/// - TAULAS = 316 (底部 Rosseland 光学深度)
/// - ABROS0 = 初始 Rosseland 不透明度估计 (通过物理公式计算)
/// - DION0 = 1.0 (初始电离度估计，完全电离)
fn generate_initial_grey_model(model: &mut ModelState, input: &InputParams) -> usize {
// Fortran 默认值 (来自 nstpar.f)
let nd = 70; // ND = 70 (Fortran 默认)
let taufir = 1e-7; // TAUFIR = 1e-7
let taulas = 316.0; // TAULAS = 316.0
let dion0 = 1.0; // DION0 = 1.0 (完全电离)
// 计算 Hopf 函数 q(τ) 的简单近似
// T(τ) = Teff * (3/4 * (τ + q(τ)))^0.25
let teff = input.teff;
let t4 = teff.powi(4);
let grav = 10.0_f64.powf(input.grav); // log g -> g
// 生成 Rosseland 光学深度网格（对数等距）
let tau_min: f64 = taufir;
let tau_max: f64 = taulas;
let log_tau_min = tau_min.ln();
let log_tau_max = tau_max.ln();
// 常数
let dprad = 1.891204931e-15 * t4; // 辐射压力项
let prad0 = dprad / 1.732; // 表面辐射压力
// 电子密度计算配置
// 注意：dion0 = 1.0 表示完全电离，用于热星的初始估计
let eldens_config = EldensConfig {
ifmol: 0,
tmolim: 1e10,
ioptab: -1, // 简单模式
iath: 1,
iatref: 1,
ihm: 0,
ih2: 0,
ih2p: 0,
pfhyd: 2.0, // 氢配分函数 (不是电离度)
};
// 平均分子量（纯 H-He 混合）
let wmm = 1.0; // 简化：假设纯氢
// STATE 函数所需的原子数据数组
// 对于简单的 H-He 模型：
// - H: 丰度 = 1.0, 最高电离级 = 2 (H I, H II)
// - He: 丰度 = 0.1, 最高电离级 = 3 (He I, He II, He III)
let abndd: [f64; MATOM] = {
let mut arr = [0.0; MATOM];
arr[0] = 1.0; // H 丰度 (相对于 H = 1.0)
arr[1] = 0.1; // He 丰度
arr
};
let ioniz: [i32; MATOM] = {
let mut arr = [0; MATOM];
arr[0] = 2; // H: 2 个电离级
arr[1] = 3; // He: 3 个电离级
arr
};
let lgr: [bool; MATOM] = [false; MATOM]; // 没有显式能级
let lrm: [bool; MATOM] = {
let mut arr = [false; MATOM];
arr[0] = true; // H: 使用隐式能级
arr[1] = true; // He: 使用隐式能级
arr
};
// 预测-校正积分的压力历史
let mut plog1 = 0.0;
let mut plog2 = 0.0;
let mut plog3 = 0.0;
let mut plog4 = 0.0;
let mut dplog1 = 0.0;
let mut dplog2 = 0.0;
let mut dplog3 = 0.0;
let dlgm = (log_tau_max - log_tau_min) / (nd - 1) as f64;
// 计算初始不透明度估计 (使用 Teff 和典型大气参数)
let initial_opacity_params = LteOpacityParams {
t: teff,
ne: 1e12, // 典型电子密度估计
nh_total: 1e12,
np: 1e12, // 假设完全电离
nh_neutral: 0.0,
nhm: 0.0,
rho: 1e-12, // 典型密度估计
uh: 2.0,
uhe: 1.0,
uhep: 2.0,
xh: 0.70,
xhe: 0.28,
};
let mut abros = quick_lte_rosseland(&initial_opacity_params);
println!(" Initial opacity estimate: {:.4} cm²/g (computed from Teff={:.0}K)", abros, teff);
for id in 0..nd {
let frac = id as f64 / (nd - 1) as f64;
let log_tau = log_tau_min + frac * (log_tau_max - log_tau_min);
let tau = log_tau.exp();
// 使用精确 Hopf 函数
let q = compute_hopf(tau, 0.0);
// 温度: T = (0.75 * Teff^4 * (tau + q))^0.25
let temp = (0.75 * t4 * (tau + q)).powf(0.25);
if id == 0 {
eprintln!("DEBUG temp: tau={:.6e}, q={:.6e}, t4={:.6e}, temp={:.1}", tau, q, t4, temp);
}
model.modpar.temp[id] = temp;
// 使用 quick_lte_rosseland 估算当前深度点的不透明度
// 这是在压力计算之前，使用温度和典型大气参数
let estimated_dens = if id > 0 {
model.modpar.dens[id - 1] // 使用前一个深度点的密度作为估计
} else {
grav * taufir * wmm * HMASS / (abros * BOLK * teff) // 表面密度估计 (正确物理公式)
};
let estimated_ne = if id > 0 {
model.modpar.elec[id - 1]
} else {
estimated_dens / HMASS // 假设完全电离
};
// 使用温度估计电离度
// 对于热星 (Teff > 20000K)，表面温度约 0.75*Teff
// 在这个温度下，氢部分电离
let estimated_ion_frac = if temp > 15000.0 {
0.95 // 高温几乎完全电离
} else if temp > 10000.0 {
0.80 // 中等温度部分电离
} else if temp > 7000.0 {
0.30 // 较低温度
} else {
0.001 // 低温几乎不电离
};
let estimated_nh_total = estimated_dens / HMASS;
let estimated_np = estimated_nh_total * estimated_ion_frac;
let estimated_nh_neutral = estimated_nh_total * (1.0 - estimated_ion_frac);
let quick_opacity_params = LteOpacityParams {
t: temp,
ne: estimated_np, // 电子密度 = 质子密度
nh_total: estimated_nh_total,
np: estimated_np,
nh_neutral: estimated_nh_neutral,
nhm: 0.0,
rho: estimated_dens,
uh: 2.0,
uhe: 1.0,
uhep: 2.0,
xh: 0.70,
xhe: 0.28,
};
let mut current_abros = quick_lte_rosseland(&quick_opacity_params);
// 流体静力学平衡计算压力
// 预测步
let mut plog = if id == 0 {
(grav / current_abros * tau + prad0).ln()
} else if id <= 3 {
plog1 + dplog1
} else {
(3.0 * plog4 + 8.0 * dplog1 - 4.0 * dplog2 + 8.0 * dplog3) / 3.0
};
// 校正步迭代 - 与 Fortran LTEGR/ROSSOP 一致，在迭代中更新不透明度
let mut ptot = plog.exp();
let mut p = ptot - tau * dprad - prad0;
let mut an = p / (BOLK * temp);
let mut ane = estimated_ne;
let mut anp = estimated_np;
let mut ahtot = estimated_nh_total;
let mut nh_neutral = estimated_nh_neutral;
let mut dens = estimated_dens;
for j in 0..10 {
// 校正步计算 (与 Fortran LTEGR 一致)
let plog_new = if id == 0 {
(grav / current_abros * tau + prad0).ln()
} else if id <= 3 {
(plog + 2.0 * plog1 + dplog1 + dplog1) / 3.0
} else {
(126.0 * plog1 - 14.0 * plog3 + 9.0 * plog4
+ 42.0 * dplog1 + 108.0 * dplog2 - 54.0 * dplog3 + 24.0 * dplog3) / 121.0
};
let err = (plog_new - plog).abs();
plog = plog_new;
ptot = plog.exp();
// 计算气体压力和粒子数密度 (与 Fortran ROSSOP 一致)
p = ptot - tau * dprad - prad0;
an = p / (BOLK * temp);
// 计算电子密度 (与 Fortran ROSSOP 调用 ELDENS 一致)
// 创建 STATE 函数所需的参数
// STATE 需要: mode=1 (LTEGR 模式), 温度, 电子密度, 原子数据
let state_params = StateParams {
mode: 1, // MODE=1 用于 LTEGR (包含显式和非显式化学物种)
id: id + 1,
t: temp,
ane, // 使用当前电子密度估计
natoms: 2, // H 和 He
hpop: an * 0.9, // 氢数密度估计 (假设大部分是 H)
dens: estimated_dens,
wmm: wmm,
ytot: 1.1, // 总原子数/氢原子数 ≈ 1 + 0.1 (He)
abndd: &abndd,
ioniz: &ioniz,
irefa: 1, // 氢是参考原子
lgr: &lgr,
lrm: &lrm,
};
let eldens_params = EldensParams {
id: id + 1,
t: temp,
an,
ytot: 1.1, // H + He 丰度因子
qref: 0.0,
dqnr: 0.0,
wmy: wmm,
config: eldens_config.clone(),
state_params: Some(state_params),
molecule_data: None,
};
let eldens_output = eldens_pure(&eldens_params, 0);
ane = eldens_output.ane;
anp = eldens_output.anp;
ahtot = eldens_output.ahtot;
nh_neutral = (ahtot - anp).max(0.0);
// 密度计算 (与 Fortran ROSSOP 一致)
// Fortran: RHO = WMM * (AN - ANE)
// WMM 在 Fortran ELDENS 中计算: wmm = dens / (an - ane)
// 由于 WMM 以 g 为单位 (不是 amu), 我们需要乘以 HMASS
// 对于氢主导气体, wmm ≈ 1 amu, 所以 RHO ≈ (AN - ANE) * HMASS
dens = (an - ane) * HMASS; // 密度 (g/cm³)
if id == 0 && j == 0 {
eprintln!("DEBUG after ELDENS: ane={:.6e}, an={:.6e}, dens={:.6e}, ahtot={:.6e}", ane, an, dens, ahtot);
}
// 更新不透明度 (与 Fortran ROSSOP 调用 OPACF0 + MEANOP 一致)
// 使用频率积分计算完整 LTE Rosseland 不透明度
let lte_params_iter = LteOpacityParams {
t: temp,
ne: ane,
nh_total: ahtot,
np: anp,
nh_neutral,
nhm: 0.0,
rho: dens,
uh: 2.0,
uhe: 1.0,
uhep: 2.0,
xh: 0.70,
xhe: 0.28,
};
let grid_iter = generate_lte_frequency_grid(input.teff, 50);
let lte_output_iter = lte_meanopt(&lte_params_iter, &grid_iter);
current_abros = lte_output_iter.opros;
// 收敛检查 (在更新不透明度之后)
// 注意：对于表面点 (id=0)，需要至少迭代一次来更新不透明度
if err <= 1e-4 && j > 0 {
break;
}
// 防止无限循环的安全检查
if j == 9 {
if id == 0 {
println!(" Warning: Corrector iteration did not converge at depth {}", id + 1);
}
}
}
// 更新压力历史 (使用最终收敛的不透明度)
if id == 0 {
eprintln!("DEBUG before pressure update: ane={:.6e}, an={:.6e}", ane, an);
}
let dplog = grav / current_abros * tau / ptot * dlgm;
plog4 = plog3;
plog3 = plog2;
plog2 = plog1;
plog1 = plog;
dplog3 = dplog2;
dplog2 = dplog1;
dplog1 = dplog;
// 计算柱质量密度
let depth = (ptot - prad0) / grav;
model.modpar.dm[id] = depth;
// 计算密度
let dens = wmm * (an - ane) * HMASS;
model.modpar.dens[id] = dens;
model.modpar.elec[id] = ane;
model.modpar.totn[id] = an;
model.modpar.anto[id] = an - ane; // 总原子密度
// 辅助量
let t = temp;
let h = 6.62620e-27_f64;
model.modpar.sqt1[id] = t.sqrt();
model.modpar.hkt1[id] = h / (BOLK * t);
model.modpar.tk1[id] = 1.0 / t;
// 计算 Rosseland 不透明度
// 使用完整的 LTE 不透明度计算，包含电子散射、束缚-自由、自由-自由和 H-
let lte_params = LteOpacityParams {
t: temp,
ne: ane,
nh_total: ahtot,
np: anp,
nh_neutral,
nhm: 0.0, // H- 密度暂时为 0
rho: dens,
uh: 2.0, // 氢配分函数
uhe: 1.0, // 氦配分函数
uhep: 2.0, // He+ 配分函数
xh: 0.70, // 氢丰度
xhe: 0.28, // 氦丰度
};
// 使用频率积分计算完整不透明度
let grid = generate_lte_frequency_grid(input.teff, 100);
let lte_output = lte_meanopt(&lte_params, &grid);
abros = lte_output.opros;
// Debug output for first and last points
if id == 0 || id == nd - 1 {
let opes = 6.6524e-25 * ane / dens.max(1e-20); // 电子散射
let ion_frac = anp / ahtot.max(1e-30); // 电离度
let opbf = 4.3e-25 * (1.0 - ion_frac) * (temp/1e4).powf(-3.5); // 束缚-自由
println!(" Depth {}: T={:.0}K, ne={:.2e}, nH={:.2e}, rho={:.2e}", id + 1, temp, ane, nh_neutral, dens);
println!(" Quick estimate κ_R={:.4e}, Full LTE κ_R={:.4e}, κ_P={:.4e}",
current_abros, lte_output.opros, lte_output.oppla);
println!(" Components: κ_es={:.4e}, κ_bf={:.4e}, κ_ff={:.4e}, κ_H-={:.4e}",
lte_output.opes, lte_output.opbf, lte_output.opff, lte_output.ophm);
}
}
// 设置深度点数
model.modpar.dmtot = model.modpar.dm[nd - 1];
println!(" Generated {} depth points", nd);
println!(" Temperature range: {:.0} K (surface) to {:.0} K (bottom)",
model.modpar.temp[0], model.modpar.temp[nd - 1]);
println!(" Electron density range: {:.2e} to {:.2e} cm^-3",
model.modpar.elec[0], model.modpar.elec[nd - 1]);
nd
}
/// 带 InputParams 的 start 函数
fn start_pure_with_input(
params: &mut StartParams,
input: &InputParams,
) -> anyhow::Result<StartOutput> {
// 设置频率参数
params.tlusty_config.basnum.nfread = input.frequencies.nfread;
// 设置原子数
params.tlusty_config.basnum.natoms = input.atoms.len() as i32;
// 设置离子数
params.tlusty_config.basnum.nion = input.ions.len() as i32;
// 计算总能级数
let nlevel: i32 = input.ions.iter().map(|ion| ion.nlevs as i32).sum();
params.tlusty_config.basnum.nlevel = nlevel;
// 调用纯计算版本的 start
Ok(start_pure(params))
}
/// 打印输入参数摘要
fn print_input_summary(params: &InputParams) {
println!("\n================================");
println!(" M O D E L A T M O S P H E R E");
println!("================================\n");
println!(" TEFF = {:>12.1}", params.teff);
println!(" LOG G = {:>12.2}", params.grav);
println!(" LTE = {}", if params.lte { "T" } else { "F" });
println!(" LTGRAY = {}", if params.ltgrey { "T" } else { "F" });
println!("\n FREQUENCIES:");
println!(" NFREAD = {}", params.frequencies.nfread);
println!("\n ATOMS: {} elements configured", params.atoms.len());
println!("\n IONS:");
for ion in &params.ions {
println!(
" {:3} (Z={:2}, ion={}) - {} levels, file: {}",
ion.typion.trim(),
ion.iat,
ion.iz,
ion.nlevs,
if ion.filei.trim().is_empty() { "(none)" } else { &ion.filei }
);
}
}
/// 解析命令行参数
fn parse_args(args: &[String]) -> anyhow::Result<(Option<PathBuf>, Option<PathBuf>)> {
let mut input_path: Option<PathBuf> = None;
let mut output_path: Option<PathBuf> = None;
let mut i = 1;
while i < args.len() {
match args[i].as_str() {
"-i" | "--input" => {
i += 1;
if i < args.len() {
input_path = Some(PathBuf::from(&args[i]));
}
}
"-o" | "--output" => {
i += 1;
if i < args.len() {
output_path = Some(PathBuf::from(&args[i]));
}
}
"-h" | "--help" => {
print_usage();
std::process::exit(0);
}
_ => {
// 位置参数：第一个是输入文件
if input_path.is_none() {
input_path = Some(PathBuf::from(&args[i]));
}
}
}
i += 1;
}
Ok((input_path, output_path))
}
fn print_usage() {
println!("TLUSTY - Non-LTE Stellar Atmosphere Calculator");
println!();
println!("Usage:");
println!(" tlusty [OPTIONS] [INPUT_FILE]");
println!();
println!("Options:");
println!(" -i, --input <FILE> Input file (default: stdin)");
println!(" -o, --output <FILE> Output file (default: stdout)");
println!(" -h, --help Show this help message");
println!();
println!("Example:");
println!(" tlusty hhe35lt.5 > hhe35lt.6");
}
/// 写入 fort.7 格式模型文件
fn write_fort7(
model: &tlusty_rust::tlusty::state::model::ModelState,
atomic: &tlusty_rust::tlusty::state::atomic::AtomicData,
nd: usize,
nlevel_actual: usize,
path: &std::path::Path,
) -> anyhow::Result<()> {
use std::fs::File;
use std::io::{BufWriter, Write};
let file = File::create(path)?;
let mut writer = BufWriter::new(file);
// NUMPAR = 3 (T, NE, RHO) + nlevel_actual
let numpar = 3 + nlevel_actual as i32;
// 写入头部：ND NUMPAR
writeln!(writer, "{:4}{:5}", nd, numpar)?;
// 写入质量深度数组（每行 6 个值）
for i in 0..nd {
write!(writer, "{:13.6E}", model.modpar.dm[i])?;
if (i + 1) % 6 == 0 || i == nd - 1 {
writeln!(writer)?;
}
}
// 写入每个深度点的数据
for id in 0..nd {
// 温度
write!(writer, "{:15.7E}", model.modpar.temp[id])?;
// 电子密度
write!(writer, "{:15.7E}", model.modpar.elec[id])?;
// 质量密度
write!(writer, "{:15.7E}", model.modpar.dens[id])?;
// 能级占据数（LTE 模型使用 Boltzmann 分布）
for ilev in 0..nlevel_actual {
let enion = atomic.levpar.enion[ilev];
let g = atomic.levpar.g[ilev];
let hkt = model.modpar.hkt1[id];
// Boltzmann 分布
let pop = g * (-enion * hkt / H).exp() * model.modpar.elec[id];
write!(writer, "{:15.7E}", pop)?;
}
writeln!(writer)?;
}
writer.flush()?;
Ok(())
}
src/tlusty/io/iroset.rs
+87 -6
View File
@@ -135,6 +135,46 @@ pub struct IrosetOutput {
pub nftt: i32,
}
// ============================================================================
// Callback 接口
// ============================================================================
/// IROSET 子程序回调接口。
///
/// 用于在 IROSET 内部调用 LEVCD、INKUL、IJALI2 等子程序。
/// 这允许调用者提供具体实现，同时保持 IROSET 的流程与 Fortran 一致。
pub trait IrosetCallbacks {
/// 调用 LEVCD(ION, IOBS) - 设置超级能级
///
/// # 参数
/// * `ion` - 离子索引 (1-based)
/// * `iobs` - 观测标志 (0=标准, 1=使用观测能级, 2=使用所有能级)
fn call_levcd(&mut self, ion: usize, iobs: i32);
/// 调用 INKUL(ION, IOBS) - 读取谱线数据
///
/// # 参数
/// * `ion` - 离子索引 (1-based)
/// * `iobs` - 观测标志
fn call_inkul(&mut self, ion: usize, iobs: i32);
/// 调用 IJALI2() - 设置 ALI 频率索引
///
/// 在完全混合 CL/ALI 方案中，设置个别跃迁的 ALI 处理标志。
/// 对应 Fortran line 170: CALL IJALI2
fn call_ijali2(&mut self);
}
/// 空回调实现（默认不做任何操作）
#[derive(Debug, Clone, Default)]
pub struct NoOpCallbacks;
impl IrosetCallbacks for NoOpCallbacks {
fn call_levcd(&mut self, _ion: usize, _iobs: i32) {}
fn call_inkul(&mut self, _ion: usize, _iobs: i32) {}
fn call_ijali2(&mut self) {}
}
// ============================================================================
// 核心计算函数
// ============================================================================
@@ -216,14 +256,16 @@ fn setup_depth_interpolation(
/// - `lined`: 谱线数据 (可变，由 LEVCD/INKUL 填充)
/// - `colkur`: 碰撞强度数据 (可变)
/// - `wop`: 占据概率数组 (可变)
/// - `callbacks`: 子程序回调接口
///
/// # 返回
/// 计算结果
pub fn iroset_pure(
pub fn iroset_pure<C: IrosetCallbacks>(
params: &mut IrosetParams,
lined: &mut Lined,
_colkur: &mut ColKur,
wop: &mut [Vec<f64>],
callbacks: &mut C,
) -> IrosetOutput {
let nd = params.nd as usize;
let nion = params.nion as usize;
@@ -277,9 +319,16 @@ pub fn iroset_pure(
}
// 设置超级能级并读取谱线数据
// 注意: 实际实现中需要调用 LEVCD 和 INKUL
// 这里简化处理，假设数据已经填充
let _iobs = odfion.ikobs[ion];
// 对应 Fortran: CALL LEVCD(ION,IOBS)
let iobs = odfion.ikobs[ion];
// 调用 LEVCD 设置超级能级
// Fortran line 109: CALL LEVCD(ION,IOBS)
callbacks.call_levcd(ion + 1, iobs);
// 对应 Fortran: CALL INKUL(ION,IOBS)
// Fortran line 110: CALL INKUL(ION,IOBS)
callbacks.call_inkul(ion + 1, iobs);
// 输出进度信息 (对应 WRITE(6,610))
#[cfg(feature = "debug_output")]
@@ -319,6 +368,10 @@ pub fn iroset_pure(
splcom.nftt = nftt;
// 对应 Fortran line 170: CALL IJALI2
// 设置 ALI 频率索引
callbacks.call_ijali2();
IrosetOutput { nftmx, nftt }
}
@@ -339,7 +392,7 @@ pub fn iroset_pure(
///
/// # 返回
/// 计算结果
pub fn iroset<W6: Write, W10: Write, W41: Write>(
pub fn iroset<W6: Write, W10: Write, W41: Write, C: IrosetCallbacks>(
params: &mut IrosetParams,
lined: &mut Lined,
colkur: &mut ColKur,
@@ -347,6 +400,7 @@ pub fn iroset<W6: Write, W10: Write, W41: Write>(
writer6: &mut W6,
writer10: &mut W10,
writer41: &mut W41,
callbacks: &mut C,
) -> Result<IrosetOutput> {
let nd = params.nd as usize;
let nfreq = params.nfreq as usize;
@@ -403,6 +457,18 @@ pub fn iroset<W6: Write, W10: Write, W41: Write>(
continue;
}
// 设置超级能级并读取谱线数据
// 对应 Fortran: CALL LEVCD(ION,IOBS)
let iobs = odfion.ikobs[ion];
// 调用 LEVCD 设置超级能级
// Fortran line 109: CALL LEVCD(ION,IOBS)
callbacks.call_levcd(ion + 1, iobs);
// 对应 Fortran: CALL INKUL(ION,IOBS)
// Fortran line 110: CALL INKUL(ION,IOBS)
callbacks.call_inkul(ion + 1, iobs);
// 进度输出
writeln!(
writer6,
@@ -556,6 +622,14 @@ pub fn iroset<W6: Write, W10: Write, W41: Write>(
// 存储截面对数
for ij in ifrku..(ifrku + nft as usize) {
let kj = ij - ifrku + nftt as usize - nft as usize + 1;
// Fortran line 189-191: 检查是否超过 MCFE 限制
if kj > MCFE {
quit_func(
" Too many Fe cross-sect. to store",
kj as i32,
MCFE as i32,
);
}
for i in 0..jidn as usize {
let sxx = (sigt[i][ij] + 1e-40f64).ln();
if kj < splcom.sigfe[0].len() && i < splcom.sigfe[0][kj].len() {
@@ -600,6 +674,10 @@ pub fn iroset<W6: Write, W10: Write, W41: Write>(
splcom.nftt = nftt;
// 对应 Fortran line 170: CALL IJALI2
// 设置 ALI 频率索引
callbacks.call_ijali2();
Ok(IrosetOutput { nftmx, nftt })
}
@@ -698,8 +776,9 @@ mod tests {
let mut lined = Lined::new(100, 3);
let mut colkur = ColKur::default();
let mut wop = vec![vec![0.0; 3]; 10];
let mut callbacks = NoOpCallbacks;
let result = iroset_pure(&mut params, &mut lined, &mut colkur, &mut wop);
let result = iroset_pure(&mut params, &mut lined, &mut colkur, &mut wop, &mut callbacks);
// 由于所有能级在 LTE，应该跳过处理
assert!(result.nftmx >= 0);
@@ -740,6 +819,7 @@ mod tests {
let mut lined = Lined::new(100, 3);
let mut colkur = ColKur::default();
let mut wop = vec![vec![0.0; 3]; 10];
let mut callbacks = NoOpCallbacks;
let mut writer6 = Cursor::new(Vec::new());
let mut writer10 = Cursor::new(Vec::new());
@@ -753,6 +833,7 @@ mod tests {
&mut writer6,
&mut writer10,
&mut writer41,
&mut callbacks,
)
.unwrap();
src/tlusty/io/linset.rs
+6 -4
View File
@@ -581,10 +581,12 @@ fn compute_profiles_and_weights(
update_intmod(params, state, itr);
// 处理 INDEXP
if trapar.indexp[itr_idx] != 0 {
// 调用 IJALIS
// TODO: 实现 IJALIS 调用
}
// Fortran: CALL IJALIS(ITR,IFRQ0,IFRQ1)
// 注意: IJALIS 需要复杂的参数结构 (TraPar, LevPar, AtoPar, TraAli, FreAux, TraCor)
// 完整实现需要通过回调或外部处理
// if trapar.indexp[itr_idx] != 0 {
// callbacks.call_ijalis(itr, ij0, ij1);
// }
Ok(())
}
src/tlusty/io/mod.rs
+3 -2
View File
@@ -60,12 +60,13 @@ pub use initia::{
FrequencyGridParams, FrequencyGridOutput,
};
pub use inpmod::{inpmod, read_tlusty_model, InputModelData, InpmodParams, InpmodOutput};
pub use input::{InputParams, read_input_file};
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 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 ltegr::{ltegr, LtegrConfig, LtegrParams, LtegrOutput};
pub use ltegrd::{ltegrd_pure, LtegrdConfig, LtegrdParams, LtegrdOutput, LtegrdAtomicData};
pub use model::{ModelFile, ModelState, read_model, write_model};
pub use nstout::{nstout, NstoutParams, NstoutOutput};
@@ -85,7 +86,7 @@ pub use rayini::{rayini, rayini_pure, rayini_with_rayleigh, read_rayleigh_table,
pub use srtfrq::{srtfrq_pure, 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, StartConfig, StartParams, StartOutput};
pub use start::{start, start_pure, start_with_callbacks, StartConfig, StartParams, StartOutput, StartCallbacks, NoOpStartCallbacks};
/// 文件单元号常量（与 Fortran 保持一致）
pub mod units {
src/tlusty/io/reader.rs
+7
View File
@@ -88,6 +88,13 @@ impl<R: BufRead> FortranReader<R> {
Ok(())
}
/// 跳过当前行（读取并丢弃）
pub fn skip_line(&mut self) -> Result<()> {
self.read_line()?;
self.remaining.clear();
Ok(())
}
/// 读取并解析一个值
///
/// 从当前行读取下一个值，如果当前行为空则读取下一行。
src/tlusty/io/start.rs
+157 -33
View File
@@ -8,27 +8,83 @@
//! 1. 读取基本配置（idisk - 大气/盘模式）
//! 2. 调用 INITIA 进行完整初始化
//! 3. 可选调用 HEDIF（He 扩散）
//! 4. 调用 COMSET 设置 COMMON 块
//! 5. 调用 PRDINI 初始化 PRD（部分重分布）
//! 4. 保存 NN0 = NN
//! 5. 调用 COMSET 设置 COMMON 块
//! 6. 调用 PRDINI 初始化 PRD（部分重分布）
//!
//! # I/O 操作
//!
//! - fort.1: 读取 idisk 参数
//!
//! # Fortran 原始代码
//!
//! ```fortran
//! SUBROUTINE START
//! INCLUDE 'IMPLIC.FOR'
//! INCLUDE 'BASICS.FOR'
//! common/hediff/ hcmass,radstr
//!
//! read(1,*,end=10,err=10) idisk
//! 10 continue
//! call initia
//! if(hcmass.gt.0.) call hedif
//! nn0=nn
//! CALL COMSET
//! call prdini
//! return
//! end
//! ```
use super::FortranReader;
use crate::tlusty::math::{comset, ComsetParams};
use crate::tlusty::state::config::TlustyConfig;
use crate::tlusty::state::atomic::AtomicData;
use crate::tlusty::state::model::ModelState;
use crate::tlusty::state::constants::{MTRANS, MLEVEL, MATOM};
// ============================================================================
// 配置参数
// 回调接口 (用于 INITIA 和 PRDINI)
// ============================================================================
/// START 模块回调接口。
///
/// 用于在 START 内部调用 INITIA、HEDIF、PRDINI 等子程序。
/// 这允许调用者提供具体实现，同时保持 START 的流程与 Fortran 一致。
pub trait StartCallbacks {
/// 调用 INITIA - 完整初始化过程
fn call_initia(&mut self);
/// 调用 HEDIF - He 扩散计算
fn call_hedif(&mut self);
/// 调用 PRDINI - PRD 初始化
fn call_prdini(&mut self);
}
/// 空回调实现（默认不做任何操作）
#[derive(Debug, Clone, Default)]
pub struct NoOpStartCallbacks;
impl StartCallbacks for NoOpStartCallbacks {
fn call_initia(&mut self) {}
fn call_hedif(&mut self) {}
fn call_prdini(&mut self) {}
}
// ============================================================================
// 配置参数 (对应 Fortran common/hediff/)
// ============================================================================
/// START 配置参数。
///
/// 对应 Fortran:
/// ```fortran
/// common/hediff/ hcmass,radstr
/// ```
#[derive(Debug, Clone)]
pub struct StartConfig {
/// 盘模型标志 (0=大气, 1=盘)
/// 从 fort.1 读取
pub idisk: i32,
/// He 扩散质量 (HCMASS)
pub hcmass: f64,
@@ -77,23 +133,42 @@ pub struct StartOutput {
/// 执行 START 初始化过程。
///
/// # 参数
/// - `params`: 输入参数
/// - `reader`: 可选的输入读取器（用于读取 idisk）
/// 严格按照 Fortran start.f 的逻辑流程实现：
///
/// # 返回值
/// 初始化结果
/// ```fortran
/// read(1,*,end=10,err=10) idisk
/// 10 continue
/// call initia
/// if(hcmass.gt.0.) call hedif
/// nn0=nn
/// CALL COMSET
/// call prdini
/// return
/// ```
pub fn start<R: std::io::BufRead>(
params: &mut StartParams,
reader: Option<&mut FortranReader<R>>,
) -> StartOutput {
start_with_callbacks(params, reader, &mut NoOpStartCallbacks)
}
/// 执行 START 初始化过程（带回调）。
///
/// 使用回调接口调用 INITIA、HEDIF、PRDINI 等子程序。
pub fn start_with_callbacks<R: std::io::BufRead, C: StartCallbacks>(
params: &mut StartParams,
reader: Option<&mut FortranReader<R>>,
callbacks: &mut C,
) -> StartOutput {
let config = &mut params.config;
// -----------------------------------------------------------
// ========================================
// Step 1: 读取 idisk
// -----------------------------------------------------------
// 对应 Fortran: read(1,*,end=10,err=10) idisk
// 10 continue
// ========================================
if let Some(r) = reader {
// 尝试读取 idisk
// 对应 end=10,err=10，出错时跳过
if let Ok(idisk_val) = r.read_value::<i32>() {
config.idisk = idisk_val;
}
@@ -102,45 +177,48 @@ pub fn start<R: std::io::BufRead>(
// 更新 TLUSTY 配置中的 idisk
params.tlusty_config.basnum.idisk = config.idisk;
// -----------------------------------------------------------
// ========================================
// Step 2: 调用 INITIA
// -----------------------------------------------------------
// initia(params.tlusty_config, params.atomic, params.model);
// 简化实现：INITIA 尚未完全实现
// 对应 Fortran: call initia
// ========================================
callbacks.call_initia();
// -----------------------------------------------------------
// ========================================
// Step 3: 可选调用 HEDIF（He 扩散）
// -----------------------------------------------------------
// 对应 Fortran: if(hcmass.gt.0.) call hedif
// ========================================
if config.hcmass > 0.0 {
// 调用 HEDIF
// 需要完整的参数设置，这里简化处理
// let hedif_params = HedifParams {
// config: params.tlusty_config,
// atomic: params.atomic,
// model: params.model,
// };
// let _hedif_result = hedif(&mut hedif_params);
callbacks.call_hedif();
}
// -----------------------------------------------------------
// ========================================
// Step 4: 保存 NN0 = NN
// -----------------------------------------------------------
// 对应 Fortran: nn0=nn
// ========================================
let nn = params.tlusty_config.matkey.nn;
params.tlusty_config.matkey.nn0 = nn;
// -----------------------------------------------------------
// ========================================
// Step 5: 调用 COMSET
// -----------------------------------------------------------
let nd = params.model.modpar.dm.len();
// 对应 Fortran: CALL COMSET
// ========================================
let nd = if params.tlusty_config.basnum.nd > 0 {
params.tlusty_config.basnum.nd as usize
} else {
params.model.modpar.dm.len()
};
let comset_params = ComsetParams {
nd,
..Default::default()
};
let _comset_result = comset(&comset_params);
// -----------------------------------------------------------
// ========================================
// Step 6: 调用 PRDINI
// -----------------------------------------------------------
// prdini(params.tlusty_config, params.atomic, params.model);
// 对应 Fortran: call prdini
// ========================================
callbacks.call_prdini();
StartOutput {
success: true,
@@ -234,4 +312,50 @@ mod tests {
assert!(result.success);
assert_eq!(params.tlusty_config.basnum.idisk, 1);
}
#[test]
fn test_start_reads_idisk_from_input() {
let mut config = StartConfig::default();
let mut tlusty_config = TlustyConfig::default();
let mut atomic = AtomicData::default();
let mut model = ModelState::new();
let mut params = StartParams {
config: &mut config,
tlusty_config: &mut tlusty_config,
atomic: &mut atomic,
model: &mut model,
};
// 模拟 fort.1 输入: idisk = 1
let input_data = b"1\n";
let mut reader = FortranReader::new(&input_data[..]);
let result = start(&mut params, Some(&mut reader));
assert!(result.success);
assert_eq!(params.tlusty_config.basnum.idisk, 1);
}
#[test]
fn test_start_sets_nn0() {
let mut config = StartConfig::default();
let mut tlusty_config = TlustyConfig::default();
tlusty_config.matkey.nn = 100; // 设置测试值
let mut atomic = AtomicData::default();
let mut model = ModelState::new();
let mut params = StartParams {
config: &mut config,
tlusty_config: &mut tlusty_config,
atomic: &mut atomic,
model: &mut model,
};
let result = start_pure(&mut params);
assert!(result.success);
assert_eq!(result.nn, 100);
assert_eq!(params.tlusty_config.matkey.nn0, 100);
}
}
src/tlusty/main.rs
+600
View File
@@ -0,0 +1,600 @@
//! TLUSTY 主程序入口。
//!
//! 重构自 TLUSTY `tlusty.f` 主程序。
//!
//! # 算法概述
//!
//! TLUSTY 使用混合 Complete Linearization (CL) 和 Accelerated Lambda Iteration (ALI) 方法
//! 计算非LTE恒星大气模型。
//!
//! # Fortran 原始代码
//!
//! ```fortran
//! PROGRAM TLUSTY
//! INCLUDE 'IMPLIC.FOR'
//! INCLUDE 'BASICS.FOR'
//! INCLUDE 'ITERAT.FOR'
//! INCLUDE 'ALIPAR.FOR'
//!
//! OPEN(UNIT=91,STATUS='SCRATCH',FORM='UNFORMATTED')
//! OPEN(UNIT=92,STATUS='SCRATCH',FORM='UNFORMATTED')
//! OPEN(UNIT=93,STATUS='SCRATCH',FORM='UNFORMATTED')
//!
//! INIT=1
//! ITER=0
//! CALL START
//! LFIN=.FALSE.
//! IF(NITER.EQ.0) LFIN=.TRUE.
//!
//! 10 ITER=ITER+1
//! CALL RESOLV
//! INIT=0
//! IF(LFIN) GO TO 20
//!
//! IF(IACC.GT.0) CALL ACCEL2
//!
//! IF(IFRYB.EQ.0) THEN
//! IF(NN.GT.MSMX) THEN
//! CALL SOLVE
//! ELSE
//! CALL SOLVES
//! END IF
//! ELSE
//! CALL RYBSOL
//! END IF
//!
//! CALL TIMING(2,ITER)
//! GO TO 10
//! 20 CONTINUE
//! STOP
//! END
//! ```
use std::io;
use std::time::Instant;
use super::io::{
FortranReader, FortranWriter,
start, StartConfig, StartParams, StartOutput,
resolv, ResolvConfig, ResolvParams,
};
use super::math::solvers::{
accel2_pure, Accel2Config, Accel2Params,
solve_pure, SolveConfig, DepthMatrices,
solves_pure,
};
use super::math::io::{timing, TimingParams, TimingMode, reset_timer};
use super::state::config::TlustyConfig;
use super::state::atomic::AtomicData;
use super::state::model::ModelState;
use super::state::constants::{MDEPTH, MTOT};
// ============================================================================
// 常量 (从 Fortran 移植)
// ============================================================================
/// 最大简化矩阵维度 (与 Fortran MSMX 相同)
pub const MSMX: usize = 2000;
// ============================================================================
// 运行状态
// ============================================================================
/// TLUSTY 运行状态。
/// 对应 Fortran 中的 INIT, ITER, LFIN 等变量
#[derive(Debug, Clone)]
pub struct TlustyState {
/// 当前迭代次数 (ITER)
pub iter: i32,
/// 最大迭代次数 (NITER)
pub niter: i32,
/// 是否初始化阶段 (INIT: 1=第一次迭代前, 0=迭代中)
pub init: i32,
/// 是否最终迭代 (LFIN)
pub lfin: bool,
/// 系统维度 (NN)
pub nn: usize,
/// 深度点数 (ND)
pub nd: usize,
}
impl Default for TlustyState {
fn default() -> Self {
Self {
iter: 0,
niter: 100,
init: 1,
lfin: false,
nn: 0,
nd: 50,
}
}
}
// ============================================================================
// 运行结果
// ============================================================================
/// TLUSTY 运行结果。
#[derive(Debug, Clone)]
pub struct TlustyResult {
/// 最终迭代次数
pub total_iterations: usize,
/// 是否收敛
pub converged: bool,
/// 总运行时间（秒）
pub total_time_secs: f64,
}
// ============================================================================
// 解法类型
// ============================================================================
/// 解法类型。
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum SolverType {
/// 标准解法 (SOLVE) - NN > MSMX
Standard,
/// 简化解法 (SOLVES) - NN <= MSMX
Simple,
/// Ryan 解法 (RYBSOL) - IFRYB != 0
Ryan,
}
/// 选择解法。
///
/// 对应 Fortran:
/// ```fortran
/// IF(IFRYB.EQ.0) THEN
/// IF(NN.GT.MSMX) THEN
/// CALL SOLVE
/// ELSE
/// CALL SOLVES
/// END IF
/// ELSE
/// CALL RYBSOL
/// END IF
/// ```
pub fn select_solver(nn: usize, msmx: usize, ifryb: i32) -> SolverType {
if ifryb != 0 {
SolverType::Ryan
} else if nn > msmx {
SolverType::Standard
} else {
SolverType::Simple
}
}
/// 检查收敛性。
pub fn check_convergence(iter: usize, max_change: f64, tolerance: f64) -> bool {
iter > 0 && max_change < tolerance
}
// ============================================================================
// 工作数组
// ============================================================================
/// TLUSTY 工作数组。
/// 对应 Fortran 中的 PSY0, PSY1, PSY2, PSY3 等
pub struct TlustyWorkArrays {
/// PSY0 数组（当前解）[nn][nd]
pub psy0: Vec<Vec<f64>>,
/// PSY1 数组（前一次迭代）
pub psy1: Vec<Vec<f64>>,
/// PSY2 数组（前两次迭代）
pub psy2: Vec<Vec<f64>>,
/// PSY3 数组（前三次迭代）
pub psy3: Vec<Vec<f64>>,
/// 深度矩阵（用于 SOLVE）
pub depth_matrices: Vec<DepthMatrices>,
}
impl TlustyWorkArrays {
pub fn new(nn: usize, nd: usize) -> Self {
Self {
psy0: vec![vec![0.0; nd]; nn],
psy1: vec![vec![0.0; nd]; nn],
psy2: vec![vec![0.0; nd]; nn],
psy3: vec![vec![0.0; nd]; nn],
depth_matrices: vec![DepthMatrices::new(nn); nd],
}
}
}
// ============================================================================
// 临时文件（对应 Fortran UNIT=91,92,93）
// ============================================================================
/// 临时存储缓冲区（对应 Fortran UNIT=91,92,93 的 SCRATCH 文件）
pub struct ScratchFiles {
/// Unit 91 - 矩阵存储
pub unit91: Vec<u8>,
/// Unit 92 - ALI 数据
pub unit92: Vec<u8>,
/// Unit 93 - 其他中间数据
pub unit93: Vec<u8>,
}
impl Default for ScratchFiles {
fn default() -> Self {
Self {
unit91: Vec::new(),
unit92: Vec::new(),
unit93: Vec::new(),
}
}
}
// ============================================================================
// 主运行函数
// ============================================================================
/// 运行 TLUSTY 计算。
///
/// 严格按照 Fortran tlusty.f 的逻辑流程实现。
///
/// # 参数
/// - `config`: TLUSTY 配置
/// - `input_reader`: 输入读取器（fort.5）
/// - `output_writer`: 输出写入器（fort.6）
///
/// # 返回值
/// 计算结果
pub fn run_tlusty<R: io::BufRead, W: io::Write>(
config: &mut TlustyConfig,
input_reader: &mut FortranReader<R>,
output_writer: &mut FortranWriter<W>,
) -> TlustyResult {
let start_time = Instant::now();
// 重置计时器
reset_timer();
// ========================================
// 初始化（对应 Fortran COMMON 块）
// ========================================
let mut atomic = AtomicData::default();
let mut model = ModelState::new();
let mut state = TlustyState::default();
// 临时文件（对应 OPEN(UNIT=91,92,93)）
let mut _scratch = ScratchFiles::default();
// 从配置中获取 NITER
state.niter = if config.runkey.niter > 0 { config.runkey.niter } else { 100 };
// ========================================
// 初始化阶段
// 对应 Fortran:
// INIT=1
// ITER=0
// CALL START
// LFIN=.FALSE.
// IF(NITER.EQ.0) LFIN=.TRUE.
// ========================================
// INIT=1, ITER=0 已在 TlustyState::default() 中设置
// CALL START
let mut start_config = StartConfig {
idisk: config.basnum.idisk,
hcmass: 0.0,
radstr: 0.0,
};
let start_output: StartOutput = {
let mut start_params = StartParams {
config: &mut start_config,
tlusty_config: config,
atomic: &mut atomic,
model: &mut model,
};
start(&mut start_params, Some(input_reader))
};
state.nn = start_output.nn as usize;
state.nd = if config.basnum.nd > 0 { config.basnum.nd as usize } else { 50 };
// LFIN=.FALSE.
state.lfin = false;
// IF(NITER.EQ.0) LFIN=.TRUE.
if config.runkey.niter == 0 {
state.lfin = true;
}
// 初始化工作数组
let nn = state.nn.max(10);
let nd = state.nd;
let mut work_arrays = TlustyWorkArrays::new(nn, nd);
// ========================================
// 主迭代循环
// 对应 Fortran:
// 10 ITER=ITER+1
// CALL RESOLV
// INIT=0
// IF(LFIN) GO TO 20
// IF(IACC.GT.0) CALL ACCEL2
// IF(IFRYB.EQ.0) THEN
// IF(NN.GT.MSMX) THEN; CALL SOLVE; ELSE; CALL SOLVES; END IF
// ELSE; CALL RYBSOL; END IF
// CALL TIMING(2,ITER)
// GO TO 10
// 20 CONTINUE
// ========================================
// 10 ITER=ITER+1
while !state.lfin && state.iter < state.niter {
state.iter += 1;
// CALL RESOLV
let resolv_config = ResolvConfig {
iter: state.iter,
init: state.init,
lfin: state.lfin,
lte: config.inppar.lte,
ioptab: config.basnum.ioptab,
icompt: config.compti.icompt,
ifprec: 0,
hmix0: config.conkey.hmix0,
iprind: config.prints.iprind,
iacpp: config.acclp.iacpp,
ielcor: config.lambda.ielcor,
nitzer: config.runkey.nitzer,
iheso6: config.basnum.iheso6,
ihecor: 0,
izscal: config.basnum.izscal,
idisk: config.basnum.idisk,
ifryb: config.basnum.ifryb,
icoolp: config.prints.icoolp,
ipopac: config.prints.ipopac,
ichckp: config.prints.ichckp,
intens: config.basnum.intens,
lchc: config.inppar.lchc,
iconrs: 1,
iconre: config.conkey.iconre,
ipconf: config.iprkey.ipconf,
iacd: config.accel.iacd,
iacc: config.accel.iacc,
lres2: config.accel.lres2 != 0,
ifpopr: 0,
inzd: config.matkey.inzd,
nfreq: config.basnum.nfreq as usize,
nfreqe: config.basnum.nfreqe as usize,
nd,
nlevel: config.basnum.nlevel as usize,
ntrans: config.basnum.ntrans as usize,
teff: config.inppar.teff,
irder: config.basnum.irder,
};
{
let mut resolv_params = ResolvParams {
config: resolv_config,
tlusty_config: config,
atomic: &mut atomic,
model: &mut model,
};
resolv(&mut resolv_params, Some(output_writer));
}
// INIT=0
state.init = 0;
// IF(LFIN) GO TO 20
if state.lfin {
break;
}
// IF(IACC.GT.0) CALL ACCEL2
if config.accel.iacc > 0 {
let mut accel_config = Accel2Config {
iter: state.iter,
niter: state.niter,
iacc: config.accel.iacc,
iacc0: config.accel.iacc0,
iacd: config.accel.iacd,
lac2: config.accel.lac2 != 0,
lres2: config.accel.lres2 != 0,
lsng: (0..MTOT).map(|i| config.accel.lsng.get(i).copied().unwrap_or(0) != 0).collect(),
nd,
nn,
};
{
let mut accel_params = Accel2Params {
config: &mut accel_config,
model: &mut model,
psy0: &mut work_arrays.psy0,
psy1: &mut work_arrays.psy1,
psy2: &mut work_arrays.psy2,
psy3: &mut work_arrays.psy3,
};
let accel_output = accel2_pure(&mut accel_params);
// 更新状态
config.accel.lac2 = if accel_output.accelerated { 1 } else { 0 };
config.accel.lres2 = if accel_output.need_resolv { 1 } else { 0 };
}
}
// IF(IFRYB.EQ.0) THEN
// IF(NN.GT.MSMX) THEN
// CALL SOLVE
// ELSE
// CALL SOLVES
// END IF
// ELSE
// CALL RYBSOL
// END IF
let solver_type = select_solver(nn, MSMX, config.basnum.ifryb);
match solver_type {
SolverType::Standard => {
// CALL SOLVE
let solve_config = SolveConfig {
nn,
nfreqe: config.basnum.nfreqe as usize,
nd,
iter: state.iter,
niter: state.niter,
iconv: config.conkey.iconv,
nretc: config.basnum.nretc,
ifali: 5,
kant: config.accel.kant.clone(),
orelax: config.accel.orelax,
chmax: config.runkey.chmax,
chmaxt: config.chnad.chmaxt,
ispodf: config.basnum.ispodf,
inhe: config.matkey.inhe,
inre: config.matkey.inre,
inpc: config.matkey.inpc,
indl: config.conkey.indl,
inzd: config.matkey.inzd,
inse: config.matkey.inse,
inmp: config.matkey.inmp,
};
// 更新深度矩阵
for id in 0..nd.min(work_arrays.depth_matrices.len()) {
for i in 0..nn.min(MTOT) {
work_arrays.depth_matrices[id].psi0[i] = work_arrays.psy0[i][id];
}
}
let solve_output = solve_pure(
&solve_config,
&work_arrays.depth_matrices,
(0.0, 0.0, 0.0, 0.0),
);
// 更新状态
state.lfin = solve_output.lfin;
// 更新 PSY0
for id in 0..nd.min(solve_output.psy0.len()) {
for i in 0..nn.min(solve_output.psy0[id].len()) {
work_arrays.psy0[i][id] = solve_output.psy0[id][i];
}
}
}
SolverType::Simple => {
// CALL SOLVES
let solves_output = solves_pure(
nn, nd,
config.basnum.nfreqe as usize,
5, config.conkey.iconv, config.basnum.nretc,
state.iter, state.niter,
&config.accel.kant,
config.accel.orelax,
10.0, 3.0, 10.0, 10.0,
config.matkey.inre as usize,
config.matkey.inhe as usize,
config.matkey.inpc as usize,
config.conkey.indl as usize,
config.matkey.inse as usize,
config.matkey.inzd as usize,
config.matkey.inmp as usize,
config.chnad.chmaxt,
config.basnum.ispodf,
config.runkey.chmax,
);
state.lfin = solves_output.lfin;
}
SolverType::Ryan => {
// CALL RYBSOL
// 简化实现：直接检查收敛
if state.iter >= state.niter {
state.lfin = true;
}
}
}
// CALL TIMING(2,ITER)
let timing_params = TimingParams {
mode: TimingMode::Linearization,
iter: state.iter,
};
let timing_output = timing(&timing_params);
// 输出时间信息
let timing_msg = format!(
" {:4}{:4}{:11.2}{:11.2} {}\n",
timing_output.iter,
timing_output.mode,
timing_output.time,
timing_output.dt,
timing_output.route
);
let _ = output_writer.write_raw(&timing_msg);
// GO TO 10 (循环继续)
}
// 20 CONTINUE
// STOP (返回结果)
let total_time = start_time.elapsed().as_secs_f64();
TlustyResult {
total_iterations: state.iter as usize,
converged: state.lfin,
total_time_secs: total_time,
}
}
// ============================================================================
// 测试
// ============================================================================
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_tlusty_state_default() {
let state = TlustyState::default();
assert_eq!(state.iter, 0);
assert_eq!(state.init, 1);
assert!(!state.lfin);
}
#[test]
fn test_select_solver() {
assert_eq!(select_solver(100, 2000, 0), SolverType::Simple);
assert_eq!(select_solver(3000, 2000, 0), SolverType::Standard);
assert_eq!(select_solver(100, 2000, 1), SolverType::Ryan);
}
#[test]
fn test_check_convergence() {
assert!(check_convergence(5, 0.001, 0.01));
assert!(!check_convergence(0, 0.001, 0.01));
assert!(!check_convergence(5, 0.1, 0.01));
}
#[test]
fn test_work_arrays_creation() {
let work = TlustyWorkArrays::new(10, 50);
assert_eq!(work.psy0.len(), 10);
assert_eq!(work.psy0[0].len(), 50);
assert_eq!(work.depth_matrices.len(), 50);
}
#[test]
fn test_run_tlusty_basic() {
let mut config = TlustyConfig::default();
config.runkey.niter = 0; // 只做初始模型
let input_data = b"0\n";
let mut input_reader = FortranReader::new(&input_data[..]);
let mut output_buffer = Vec::new();
let mut output_writer = FortranWriter::new(&mut output_buffer);
let result = run_tlusty(&mut config, &mut input_reader, &mut output_writer);
assert!(result.converged);
assert_eq!(result.total_iterations, 0);
}
}
src/tlusty/math/continuum/lte_opacity.rs
+441
View File
@@ -0,0 +1,441 @@
//! LTE 不透明度的简化物理计算。
//!
//! 使用物理公式计算 Rosseland 和 Planck 平均不透明度，
//! 作为完整表插值方法的替代。
//!
//! # 不透明度来源
//!
//! 1. 电子散射 (Thomson 散射)
//! 2. 束缚-自由跃迁 (氢光致电离，Kramers 截面)
//! 3. 自由-自由跃迁 (氢轫致辐射)
//! 4. H- 不透明度 (负氢离子)
//!
//! # 参考
//!
//! - TLUSTY opacfl.f, meanopt.f
//! - Mihalas (1978) Stellar Atmospheres
use crate::tlusty::state::constants::{H, HK, BOLK, SIGE, EH, UN, TWO};
// ============================================================================
// 物理常数
// ============================================================================
/// 光速 (cm/s)
const CLIGHT: f64 = 2.99792458e10;
/// 氢电离频率 (Hz)
const FRH: f64 = 3.28805e15;
/// H- 自由-自由系数
const CFF1: f64 = 1.3727e-25;
const CFF2: f64 = 4.3748e-10;
const CFF3: f64 = 2.5993e-7;
/// 自由-自由基准截面
const SGFF0: f64 = 3.694e8;
/// H- 电离阈值频率 (Hz)
const FRHM: f64 = 1.82e15;
// ============================================================================
// 数据结构
// ============================================================================
/// LTE 不透明度输入参数
#[derive(Debug, Clone)]
pub struct LteOpacityParams {
/// 温度 (K)
pub t: f64,
/// 电子密度 (cm⁻³)
pub ne: f64,
/// 总氢密度 (中性 + 电离) (cm⁻³)
pub nh_total: f64,
/// 质子密度 (cm⁻³)
pub np: f64,
/// 中性氢密度 (cm⁻³)
pub nh_neutral: f64,
/// H- 密度 (cm⁻³)
pub nhm: f64,
/// 密度 (g/cm³)
pub rho: f64,
/// 氢配分函数
pub uh: f64,
/// 氦配分函数
pub uhe: f64,
/// He+ 配分函数
pub uhep: f64,
/// 氢丰度 (质量分数)
pub xh: f64,
/// 氦丰度 (质量分数)
pub xhe: f64,
}
impl Default for LteOpacityParams {
fn default() -> Self {
Self {
t: 10000.0,
ne: 1e12,
nh_total: 1e12,
np: 5e11,
nh_neutral: 5e11,
nhm: 0.0,
rho: 1e-12,
uh: 2.0,
uhe: 1.0,
uhep: 2.0,
xh: 0.70,
xhe: 0.28,
}
}
}
/// LTE 不透明度输出
#[derive(Debug, Clone)]
pub struct LteOpacityOutput {
/// Rosseland 平均不透明度 (cm²/g)
pub opros: f64,
/// Planck 平均不透明度 (cm²/g)
pub oppla: f64,
/// 电子散射不透明度 (cm²/g)
pub opes: f64,
/// 束缚-自由不透明度 (cm²/g)
pub opbf: f64,
/// 自由-自由不透明度 (cm²/g)
pub opff: f64,
/// H- 不透明度 (cm²/g)
pub ophm: f64,
}
/// LTE 频率网格
#[derive(Debug, Clone)]
pub struct LteFrequencyGrid {
/// 频率数组 (Hz)
pub freq: Vec<f64>,
/// 权重数组
pub weights: Vec<f64>,
/// Planck 函数
pub bnue: Vec<f64>,
}
// ============================================================================
// 核心计算函数
// ============================================================================
/// 生成用于 LTE 不透明度积分的频率网格。
pub fn generate_lte_frequency_grid(teff: f64, nfreq: usize) -> LteFrequencyGrid {
let frmin: f64 = 1e13;
let frmax: f64 = 3e16;
let log_frmin = frmin.ln();
let log_frmax = frmax.ln();
let dlog = (log_frmax - log_frmin) / (nfreq - 1) as f64;
let mut freq = Vec::with_capacity(nfreq);
let mut weights = Vec::with_capacity(nfreq);
let mut bnue = Vec::with_capacity(nfreq);
let c1 = 2.0 * H / (CLIGHT * CLIGHT);
for i in 0..nfreq {
let log_fr = log_frmin + i as f64 * dlog;
let fr = log_fr.exp();
freq.push(fr);
let w = if i == 0 || i == nfreq - 1 {
0.5 * dlog * fr
} else {
dlog * fr
};
weights.push(w);
let x = HK * fr / teff;
let ex = if x < 150.0 { x.exp() } else { 1e150 };
let bn = c1 * fr.powi(3) / (ex - 1.0);
bnue.push(bn);
}
LteFrequencyGrid { freq, weights, bnue }
}
/// 计算 LTE 模式的完整不透明度。
pub fn lte_meanopt(params: &LteOpacityParams, grid: &LteFrequencyGrid) -> LteOpacityOutput {
let t = params.t;
let ne = params.ne;
let nh = params.nh_neutral;
let np = params.np;
let nhm = params.nhm;
let rho = params.rho;
if rho <= 0.0 {
return LteOpacityOutput {
opros: 0.4,
oppla: 0.4,
opes: 0.0,
opbf: 0.0,
opff: 0.0,
ophm: 0.0,
};
}
let hkt = HK / t;
let sqrt_t = t.sqrt();
let sgff = SGFF0 / sqrt_t * ne;
let mut abr = 0.0;
let mut sumdb = 0.0;
let mut abp = 0.0;
let mut sumb = 0.0;
for (ij, &fr) in grid.freq.iter().enumerate() {
let w = grid.weights[ij];
let bnue = grid.bnue[ij];
let x = hkt * fr;
let x_clamped = x.min(150.0);
let ex = x_clamped.exp();
let e1 = 1.0 / (ex - 1.0);
let plan = bnue * e1 * w;
let dplan = plan * hkt * fr * ex * e1 * e1;
let (ab, sct) = compute_opacity_at_frequency(fr, t, ne, nh, np, nhm, hkt, sgff, params);
let total = ab + sct;
if total > 0.0 {
abr = abr + dplan / total;
}
sumdb = sumdb + dplan;
abp = abp + plan * ab;
sumb = sumb + plan;
}
let oprol = if abr > 0.0 { sumdb / abr } else { 0.0 };
let opplal = if sumb > 0.0 { abp / sumb } else { 0.0 };
let opros = oprol / rho;
let oppla = opplal / rho;
let opes = SIGE * ne / rho;
let (opbf, opff, ophm) = compute_mean_opacities_per_gram(params, sqrt_t);
LteOpacityOutput {
opros,
oppla,
opes,
opbf,
opff,
ophm,
}
}
/// 计算给定频率点的吸收和散射系数 (per cm³)。
fn compute_opacity_at_frequency(
fr: f64,
t: f64,
ne: f64,
nh: f64,
np: f64,
nhm: f64,
hkt: f64,
sgff: f64,
params: &LteOpacityParams,
) -> (f64, f64) {
let mut ab = 0.0;
let sct = SIGE * ne;
// 1. 氢束缚-自由
if fr >= FRH && nh > 0.0 {
let sigma_bf0 = 6.3e-18;
let sigma_bf = sigma_bf0 * (FRH / fr).powi(3);
let gaunt_bf = hydrogen_gaunt_bf(fr);
ab = ab + sigma_bf * gaunt_bf * nh;
}
// 2. 氢自由-自由
if np > 0.0 && ne > 0.0 {
let frinv = 1.0 / fr;
let fr3inv = frinv * frinv * frinv;
let sf1 = sgff * fr3inv;
let exp_factor = (-hkt * fr).exp();
let sf2 = 1.0 / (1.0 - exp_factor).max(1e-30);
let absoff = sf1 * sf2 * np;
ab = ab + absoff;
}
// 3. H- 自由-自由
if nhm > 0.0 && ne > 0.0 {
let frinv = 1.0 / fr;
let cfft = CFF2 - CFF3 / t;
let abhm_ff = (CFF1 + cfft * frinv) * nhm * ne * frinv;
ab = ab + abhm_ff;
}
// 4. H- 束缚-自由
if nhm > 0.0 && fr >= FRHM {
let sigma_hm = compute_hm_photodetachment_cross_section(fr);
ab = ab + sigma_hm * nhm;
}
// 5. He 束缚-自由
if params.xhe > 0.0 {
let he_abundance = params.xhe / 4.0 * params.nh_total / params.xh.max(0.1);
if fr >= 1.81e15 && he_abundance > 0.0 {
let sigma_he = 7.83e-18 * (1.81e15 / fr).powi(3);
let he_neutral = he_abundance * 0.9;
ab = ab + sigma_he * he_neutral;
}
}
(ab, sct)
}
/// 计算氢束缚-自由 Gaunt 因子。
fn hydrogen_gaunt_bf(fr: f64) -> f64 {
let u = fr / FRH;
if u < 1.0 {
0.0
} else if u < 2.0 {
0.9
} else if u < 10.0 {
0.85
} else {
0.8
}
}
/// 计算 H- 光致分离截面。
fn compute_hm_photodetachment_cross_section(fr: f64) -> f64 {
if fr < FRHM {
return 0.0;
}
let x = fr / FRHM - 1.0;
if x <= 0.0 {
return 0.0;
}
let sqrt_x = x.sqrt();
let sigma_0 = 4.0e-17;
let a = 1.0 + 0.5 * x - 0.1 * x * x;
sigma_0 * sqrt_x * a
}
/// 计算平均不透明度分量 (每克)。
fn compute_mean_opacities_per_gram(params: &LteOpacityParams, sqrt_t: f64) -> (f64, f64, f64) {
let rho = params.rho;
if rho <= 0.0 {
return (0.0, 0.0, 0.0);
}
let t = params.t;
let t4 = t / 1e4;
let t_factor = t4.powf(-3.5);
let ionization = if params.nh_total > 0.0 {
(params.np / params.nh_total).min(1.0)
} else {
1.0
};
let kappa_bf_h = 4.3e-25 * (1.0 - ionization) * t_factor * params.xh;
let kappa_bf_he = 1.0e-25 * (1.0 - ionization) * t_factor * params.xhe;
let opbf = kappa_bf_h + kappa_bf_he;
let kappa_ff = 1.0e-26 * ionization * (1.0 + ionization) * t_factor;
let opff = kappa_ff;
let ophm = if params.nhm > 0.0 && t < 10000.0 {
let t4_inv = 1e4 / t;
let sigma_hm = 4e-17;
sigma_hm * params.nhm / rho * t4_inv * t4_inv
} else {
0.0
};
(opbf, opff, ophm)
}
/// 快速计算 LTE Rosseland 平均不透明度（解析近似）。
pub fn quick_lte_rosseland(params: &LteOpacityParams) -> f64 {
let rho = params.rho;
if rho <= 0.0 {
return 0.4;
}
let t = params.t;
let ne = params.ne;
let np = params.np;
let nh_neutral = params.nh_neutral;
let kappa_es = SIGE * ne / rho;
let t4 = t / 1e4;
let t_factor = t4.powf(-3.5);
let nh_total = np + nh_neutral;
let ionization = if nh_total > 0.0 {
(np / nh_total).min(1.0)
} else {
1.0
};
let kramer_bf = 4.3e-25 * (1.0 - ionization) * t_factor;
let kramer_ff = 1.0e-26 * ionization * (1.0 + ionization) * t_factor;
let nh_factor = nh_total / rho.max(1e-30);
kappa_es + (kramer_bf + kramer_ff) * nh_factor
}
// ============================================================================
// 测试
// ============================================================================
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_lte_opacity_hot_star() {
let params = LteOpacityParams {
t: 30000.0,
ne: 1e14,
nh_total: 1e14,
np: 9e13,
nh_neutral: 1e13,
nhm: 0.0,
rho: 1e-10,
uh: 2.0,
uhe: 1.0,
uhep: 2.0,
xh: 0.70,
xhe: 0.28,
};
let grid = generate_lte_frequency_grid(35000.0, 100);
let result = lte_meanopt(&params, &grid);
assert!(result.opros > 0.0);
assert!(result.opes / result.opros > 0.3);
}
#[test]
fn test_quick_lte_rosseland() {
let params = LteOpacityParams {
t: 10000.0,
ne: 1e13,
nh_total: 1e15,
np: 5e12,
nh_neutral: 5e14,
nhm: 0.0,
rho: 1e-10,
..Default::default()
};
let kappar = quick_lte_rosseland(&params);
assert!(kappar > 0.0);
}
}
src/tlusty/math/continuum/mod.rs
+10
View File
@@ -1,5 +1,6 @@
//! continuum module
mod lte_opacity;
mod opacf0;
mod opacf1;
mod opacfa;
@@ -15,7 +16,15 @@ mod opaini;
mod opctab;
mod opdata;
mod opfrac;
mod opacity_table;
pub use lte_opacity::{
LteOpacityParams, LteOpacityOutput, LteFrequencyGrid,
lte_meanopt, generate_lte_frequency_grid, quick_lte_rosseland,
};
pub use lte_opacity::{
LteOpacityParams as LteOpacityParamsOld, LteOpacityOutput as LteOpacityOutputOld,
};
pub use opacf0::*;
pub use opacf1::*;
pub use opacfa::*;
@@ -31,3 +40,4 @@ pub use opaini::*;
pub use opctab::*;
pub use opdata::*;
pub use opfrac::*;
pub use opacity_table::{OwnedOpacityTable, create_simple_table};
src/tlusty/math/continuum/opacf0.rs
+276 -29
View File
@@ -15,7 +15,11 @@
//! 6. 初始化谱线不透明度
//! 7. 循环频率点计算总不透明度
use crate::tlusty::state::constants::{HK, H, UN, SIGE, NLMX, MFREQ, MFREQL, MLEVEL, MTRANS, MION, MMER};
use crate::tlusty::state::constants::{HK, H, UN, SIGE, NLMX, MFREQ, MFREQL, MLEVEL, MTRANS, MION, MMER, MDEPTH};
use crate::tlusty::state::{GffPar, DwnPar, ModPar, InpPar};
use crate::tlusty::math::atomic::{gfree0, gfree1};
use crate::tlusty::math::opacity::dwnfr0;
use crate::tlusty::math::opacity::dwnfr1;
// 物理常数 (来自 opacf0.f)
/// Rydberg 频率
@@ -78,6 +82,21 @@ impl Default for Opacf0Config {
}
}
/// OPACF0 上下文，包含所有必要的依赖结构体
///
/// 这个结构体用于传递 GFREE0、DWNFR0 等函数所需的状态
#[derive(Debug)]
pub struct Opacf0Context<'a> {
/// 自由-自由 Gaunt 因子参数
pub gffpar: &'a mut GffPar,
/// 溶解分数参数
pub dwnpar: &'a mut DwnPar,
/// 模型参数 (用于 DWNFR0)
pub modpar: &'a ModPar,
/// 输入参数 (用于 DWNFR1)
pub inppar: &'a InpPar,
}
/// OPACF0 模型状态参数
#[derive(Debug)]
pub struct Opacf0ModelState<'a> {
@@ -105,6 +124,8 @@ pub struct Opacf0ModelState<'a> {
pub sqt1: &'a mut [f64],
/// TEMP1 (nd) - 1/T
pub temp1: &'a mut [f64],
/// EMEL1 (nd) - 电子发射因子
pub emel1: &'a mut [f64],
/// ELEC1 (nd) - 1/ne
pub elec1: &'a mut [f64],
/// DENS1 (nd) - 1/n
@@ -284,6 +305,65 @@ pub struct Opacf0Output<'a> {
pub gmer: &'a [f64],
/// SGMG (mmer × nd) - Mermerges 截面
pub sgmg: &'a mut [f64],
// 束缚-自由截面数据
/// BFCS - 光电离截面表 (mcross × nfreqc)
pub bfcs: &'a [f32],
/// IJBF - 频率插值索引 (nfreq)
pub ijbf: &'a [i32],
/// AIJBF - 频率插值系数 (nfreq)
pub aijbf: &'a [f64],
}
// ============================================================================
// 截面计算辅助函数
// ============================================================================
/// 计算光电离截面 CROSS(IBFT, IJ)
///
/// 从预计算的截面表 BFCS 中通过线性插值获取截面值。
///
/// # 参数
/// * `ibft` - 束缚-自由跃迁索引 (0-indexed)
/// * `ij` - 频率索引 (0-indexed)
/// * `output` - 输出结构体，包含 BFCS, IJBF, AIJBF
///
/// # Fortran 对应
/// ```fortran
/// SG = CROSS(IBFT, IJ)
/// ```
#[inline]
fn cross(ibft: usize, ij: usize, output: &Opacf0Output) -> f64 {
let ij0 = output.ijbf[ij] as usize;
let a1 = output.aijbf[ij];
// BFCS 是 (mcross × nfreqc) 数组
let sig0 = output.bfcs[ibft * MFREQ + ij0] as f64;
let sig1 = output.bfcs[ibft * MFREQ + ij0 + 1] as f64;
a1 * sig0 + (UN - a1) * sig1
}
/// 计算含双电子复合的光电离截面 CROSSD(IBFT, IJ, ID)
///
/// 与 CROSS 类似，但考虑了双电子复合的深度相关修正。
/// 目前简化为调用 CROSS。
///
/// # 参数
/// * `ibft` - 束缚-自由跃迁索引 (0-indexed)
/// * `ij` - 频率索引 (0-indexed)
/// * `_id` - 深度索引 (0-indexed)，目前未使用
/// * `output` - 输出结构体
///
/// # Fortran 对应
/// ```fortran
/// SG = CROSSD(IBFT, IJ, ID)
/// ```
#[inline]
fn crossd(ibft: usize, ij: usize, _id: usize, output: &Opacf0Output) -> f64 {
// 简化版本：直接调用 cross
// 完整实现需要考虑双电子复合的深度相关修正
cross(ibft, ij, output)
}
/// 束缚-自由截面函数类型
@@ -292,6 +372,58 @@ pub type CrossFn = fn(ibft: usize, ij: usize) -> f64;
/// 双电子截面函数类型
pub type CrossDFn = fn(ibft: usize, ij: usize, id: usize) -> f64;
// ============================================================================
// 回调接口 (类似 OPACF1 的模式)
// ============================================================================
/// OPADD 回调结果
#[derive(Debug, Clone, Default)]
pub struct OpaddResult {
/// 吸收系数
pub abad: f64,
/// 发射系数
pub emad: f64,
/// 散射系数
pub scad: f64,
}
/// 子程序回调接口
///
/// 用于在 OPACF0 内部调用 WNSTOR、SABOLF、LINPRO、OPADD、OPACT1 等子程序。
/// 这允许调用者提供具体实现，同时保持 OPACF0 的流程与 Fortran 一致。
pub trait Opacf0Callbacks {
/// 调用 WNSTOR(ID) - 存储氢积分
fn call_wnstor(&mut self, id: usize);
/// 调用 SABOLF(ID) - Saha-Boltzmann 因子
fn call_sabolf(&mut self, id: usize);
/// 调用 LINPRO(ITR, ID, PRF) - 谱线轮廓
/// 返回 PRF 数组
fn call_linpro(&mut self, itr: usize, id: usize, prf: &mut [f64]);
/// 调用 OPADD(MODE, ICALL, IJ, ID) - 附加不透明度
/// 返回 (abad, emad, scad)
fn call_opadd(&mut self, mode: i32, icall: i32, ij: usize, id: usize) -> OpaddResult;
/// 调用 OPACT1(IJ) - 表格不透明度
fn call_opact1(&mut self, ij: usize);
}
/// 空回调实现（默认不做任何操作）
#[derive(Debug, Clone, Default)]
pub struct NoOpCallbacks;
impl Opacf0Callbacks for NoOpCallbacks {
fn call_wnstor(&mut self, _id: usize) {}
fn call_sabolf(&mut self, _id: usize) {}
fn call_linpro(&mut self, _itr: usize, _id: usize, _prf: &mut [f64]) {}
fn call_opadd(&mut self, _mode: i32, _icall: i32, _ij: usize, _id: usize) -> OpaddResult {
OpaddResult::default()
}
fn call_opact1(&mut self, _ij: usize) {}
}
// ============================================================================
// 主函数
// ============================================================================
@@ -309,7 +441,8 @@ pub type CrossDFn = fn(ibft: usize, ij: usize, id: usize) -> f64;
/// * `atomic` - 原子数据
/// * `freq_params` - 频率数据
/// * `output` - 输出数组
pub fn opacf0(
/// * `context` - 上下文，包含依赖结构体
pub fn opacf0<C: Opacf0Callbacks>(
id: usize,
nfrq: usize,
config: &Opacf0Config,
@@ -317,6 +450,8 @@ pub fn opacf0(
atomic: &mut Opacf0AtomicParams,
freq_params: &Opacf0FreqParams,
output: &mut Opacf0Output,
context: &mut Opacf0Context,
callbacks: &mut C,
) {
let id_idx = id - 1; // 转换为 0-indexed
let nd = model.nd;
@@ -333,8 +468,17 @@ pub fn opacf0(
model.sqt1[id_idx] = t.sqrt();
model.temp1[id_idx] = t1;
// 调用 GFREE0 初始化自由-自由 Gaunt 因子
// CALL GFREE0(ID) - 由外部调用或在此调用
// ========================================================================
// 1.1 调用 GFREE0 初始化自由-自由 Gaunt 因子
// 对应 Fortran: CALL GFREE0(ID)
// ========================================================================
gfree0(id_idx, model.temp, context.gffpar);
// ========================================================================
// 1.2 初始化电子发射因子
// 对应 Fortran line 38: EMEL1(ID)=UN
// ========================================================================
model.emel1[id_idx] = UN;
// ========================================================================
// 2. 初始化电子密度相关量 (类似 OPAINI)
@@ -354,10 +498,23 @@ pub fn opacf0(
model.elscat[id_idx] = ane * SIGE;
// 调用辅助函数
// CALL DWNFR0(ID) - 下沉修正初始化
// CALL WNSTOR(ID) - 氢积分存储
// CALL SABOLF(ID) - 束缚-自由 Sa Boltzmann 因子
// ========================================================================
// 2.1 调用 DWNFR0 初始化下沉修正
// 对应 Fortran: CALL DWNFR0(ID)
// ========================================================================
dwnfr0(id_idx, context.modpar, context.dwnpar);
// ========================================================================
// 2.2 调用 WNSTOR 存储氢积分
// 对应 Fortran: CALL WNSTOR(ID)
// ========================================================================
callbacks.call_wnstor(id);
// ========================================================================
// 2.3 调用 SABOLF 计算 Saha-Boltzmann 因子
// 对应 Fortran: CALL SABOLF(ID)
// ========================================================================
callbacks.call_sabolf(id);
// ========================================================================
// 3. 计算束缚-自由不透明度预备量
@@ -465,6 +622,10 @@ pub fn opacf0(
let laser = config.iter > config.itlas;
// 保存循环内最后有效的谱线索引用于后续计算
// (Fortran 中循环外的代码使用循环内最后的变量值)
let mut last_valid_itr: Option<usize> = None;
if nfrq > atomic.nfreqc {
// 初始化主谱线轮廓
for itr in 0..atomic.ntrans {
@@ -475,6 +636,8 @@ pub fn opacf0(
continue;
}
last_valid_itr = Some(itr);
let indxa = atomic.indexp[itr].abs();
let ijl0 = if config.ispodf >= 1 {
atomic.kfr0[itr] as usize
@@ -489,13 +652,81 @@ pub fn opacf0(
if indxa < 2 || indxa > 4 {
// 调用 LINPRO 计算谱线轮廓
// CALL LINPRO(ITR,ID,PRF)
// 这里需要外部提供 LINPRO 实现
// 对应 Fortran: CALL LINPRO(ITR,ID,PRF)
let mut prf = vec![0.0; MFREQL];
callbacks.call_linpro(itr + 1, id, &mut prf);
// 将 PRF 复制到 PRFLIN
for ij in ijl0..=ijl1 {
if ij - ijl0 < prf.len() && id_idx * MFREQL + ij < output.prflin.len() {
output.prflin[id_idx * MFREQL + ij] = prf[ij - ijl0] as f32;
}
}
}
}
// 计算谱线吸收/发射系数
// (这部分在原代码中有 bug - 循环外的代码使用了循环内的变量)
// 对应 Fortran lines 143-161
if let Some(itr) = last_valid_itr {
let ii = atomic.ilow[itr] as usize - 1;
let jj = atomic.iup[itr] as usize - 1;
let gg = atomic.g[ii] / atomic.g[jj];
let (pi, pj) = if atomic.ifwop[jj] >= 0 {
// PI = POPUL(II,ID)*WOP(JJ,ID)
// PJ = POPUL(JJ,ID)*WOP(II,ID)*GG
let pi_val = get_popul(atomic.nlevel, id_idx, ii, model.popul)
* get_wop(atomic.nlevel, id_idx, jj, atomic.wop);
let pj_val = get_popul(atomic.nlevel, id_idx, jj, model.popul)
* get_wop(atomic.nlevel, id_idx, ii, atomic.wop)
* gg;
(pi_val, pj_val)
} else {
// PI = POPUL(II,ID)
// PJ = POPUL(JJ,ID)*WOP(II,ID)*G(II)/GMER(IMRG(JJ),ID)
let imrg_jj = atomic.imrg[jj] as usize - 1;
let gmer_val = if imrg_jj < MMER {
output.gmer[imrg_jj * nd + id_idx]
} else {
1.0
};
let pi_val = get_popul(atomic.nlevel, id_idx, ii, model.popul);
let pj_val = get_popul(atomic.nlevel, id_idx, jj, model.popul)
* get_wop(atomic.nlevel, id_idx, ii, atomic.wop)
* atomic.g[ii]
/ gmer_val;
(pi_val, pj_val)
};
// ABTRA(ITR,ID) = PI
// EMTRA(ITR,ID) = PJ * EXP(FR0(ITR)*HKT1(ID))
output.abtra[itr * nd + id_idx] = pi;
let fr0_itr = atomic.fr0[itr];
output.emtra[itr * nd + id_idx] = pj * (fr0_itr * model.hkt1[id_idx]).exp();
// 激光抑制逻辑 (Fortran lines 153-161)
// IF(LASER) THEN ...
if laser {
let mut qtt = 0.0;
if (pi - pj).abs() > 1e-30 {
qtt = pj / (pi - pj) * ((fr0_itr * model.hkt1[id_idx]).exp() - UN);
}
// lfr = fr0(itr).lt.frtabm.and.iadop(iatm(ii)).gt.0
let lfr = fr0_itr < freq_params.frtabm
&& atomic.iatm[ii] > 0
&& {
let iatm_idx = atomic.iatm[ii] as usize - 1;
iatm_idx < atomic.iadop.len() && atomic.iadop[iatm_idx] > 0
};
// IF(QTT.LT.0. .OR. QTT.GT.QTLAS .or. lfr) THEN
if qtt < 0.0 || qtt > config.qtlas || lfr {
output.abtra[itr * nd + id_idx] = 0.0;
output.emtra[itr * nd + id_idx] = 0.0;
}
}
}
}
// ========================================================================
@@ -541,20 +772,29 @@ pub fn opacf0(
}
// 获取截面
let sg = if config.ifdiel == 0 {
// 对应 Fortran lines 191-195
let mut sg = if config.ifdiel == 0 {
// SG = CROSS(IBFT,IJ)
0.0 // 需要外部截面函数
cross(ibft, ij_idx, output)
} else {
// SG = CROSSD(IBFT,IJ,ID)
0.0 // 需要外部截面函数
crossd(ibft, ij_idx, id_idx, output)
};
// Mermerges 处理
// 对应 Fortran lines 196-201
if atomic.ifwop[ii] < 0 {
let imer = atomic.imrg[ii] as usize - 1;
// 调用 SGMER1
// CALL SGMER1(FRINV,FR3INV,IMER,ID,SGME1)
// output.sgmg[imer * nd + id_idx] = sgme1;
// 调用 SGMER1 计算 Mermerges 截面
// 对应 Fortran: CALL SGMER1(FRINV,FR3INV,IMER,ID,SGME1)
// ISU = INT(SQRT(FRCH(IMER)*FRINV)) + 1
let isu = ((output.frch[imer] * frinv).sqrt() as usize).min(NLMX - 1);
// SGME1 = SGMSUM(ISU,IMER,ID) * FR3INV
// SGMSUM 索引: (isu, imer, id) -> isu * MMER * MDEPTH + imer * MDEPTH + id
let sgme1 = output.sgmsum[isu * MMER * MDEPTH + imer * nd + id_idx] * fr3inv;
output.sgmg[imer * nd + id_idx] = sgme1;
// SG = SGME1 (替换原来的截面值)
sg = sgme1;
}
if sg <= 0.0 {
@@ -562,13 +802,16 @@ pub fn opacf0(
}
// Macfarlane 下沉修正
// 对应 Fortran lines 203-208
if atomic.mcdw[itr] > 0 {
let izz = atomic.iz[atomic.iel[ii] as usize - 1];
// 调用 DWNFR1
// CALL DWNFR1(FR,FR0(ITR),ID,IZZ,DW1)
// let dw1 = ...;
// output.dwf1[(atomic.mcdw[itr] - 1) as usize * nd + id_idx] = dw1;
// sg = sg * dw1;
let izz = atomic.iz[atomic.iel[ii] as usize - 1] as usize - 1;
let fr0_itr = atomic.fr0[itr];
// 调用 DWNFR1 计算下沉修正因子
// 对应 Fortran: CALL DWNFR1(FR,FR0(ITR),ID,IZZ,DW1)
let dw1 = dwnfr1(fr, fr0_itr, id_idx, izz, context.inppar, context.dwnpar);
output.dwf1[(atomic.mcdw[itr] - 1) as usize * nd + id_idx] = dw1;
// SG = SG * DW1
sg *= dw1;
}
let emis_bf = sg * output.emtra[itr * nd + id_idx];
@@ -605,14 +848,16 @@ pub fn opacf0(
}
2 => {
// 氢型精确 Gaunt
// 对应 Fortran lines 232-240
let sf1 = output.sff3[ion * nd + id_idx] * fr3inv;
let sf2 = if fr < atomic.ff[ion] {
let mut sf2 = if fr < atomic.ff[ion] {
UN / output.xkf[id_idx]
} else {
output.sff2[ion * nd + id_idx]
};
let x = C14 * atomic.charg2[ion] as f64 / fr;
// sf2 = sf2 - UN + GFREE1(ID,X)
sf2 = sf2 - UN + gfree1(id_idx, x, context.gffpar);
sf1 * sf2
}
3 => {
@@ -644,10 +889,11 @@ pub fn opacf0(
if config.iopadd != 0 {
// 调用 OPADD
// CALL OPADD(0,ICALL,IJ,ID)
// output.abso[ij_idx] += abad;
// output.emis[ij_idx] += emad;
// output.scat[ij_idx] += scad;
// 对应 Fortran: CALL OPADD(0,ICALL,IJ,ID)
let result = callbacks.call_opadd(0, icall, ij + 1, id);
output.abso[ij_idx] += result.abad;
output.emis[ij_idx] += result.emad;
output.scat[ij_idx] += result.scad;
}
// --------------------------------------------------------------------
@@ -777,7 +1023,8 @@ pub fn opacf0(
if config.ioptab > 0 {
// 调用 OPACT1
// CALL OPACT1(IJ)
// 对应 Fortran: call opact1(ij)
callbacks.call_opact1(ij + 1);
}
}
}
src/tlusty/math/continuum/opacf1.rs
+186 -31
View File
@@ -6,6 +6,8 @@
//! 这与 OPACF0 互补：OPACF0 计算单深度所有频率，OPACF1 计算单频率所有深度。
use crate::tlusty::state::constants::{UN, SIGE, NLMX, MFREQ, MLEVEL, MTRANS, MION, MMER, MDEPTH};
use crate::tlusty::state::config::InpPar;
use crate::tlusty::state::model::DwnPar;
// 物理常数
/// 光速 × 1e14
@@ -38,6 +40,61 @@ pub struct Opacf1Config {
pub ioplym: i32,
/// PRD 标志 (>0: 计算 PRD)
pub ifprd: i32,
/// 光致电离截面特殊深度点 (0: 无, >0: 深度索引)
pub iprcrs: i32,
/// 光致电离截面能级偏移 (nprcrs)
pub nprcrs: i32,
}
/// OPADD 回调结果
#[derive(Debug, Clone, Default)]
pub struct OpaddResult {
/// 吸收系数
pub abad: f64,
/// 发射系数
pub emad: f64,
/// 散射系数
pub scad: f64,
}
/// 子程序回调接口
///
/// 用于在 OPACF1 内部调用 QUASIM、GHYDOP、LYMLIN、OPACT1、OPADD、PRD 等子程序。
/// 这允许调用者提供具体实现，同时保持 OPACF1 的流程与 Fortran 一致。
pub trait Opacf1Callbacks {
/// 调用 QUASIM(IJ) - 准分子不透明度
fn call_quasim(&mut self, ij: usize);
/// 调用 GHYDOP(IJ) - 氢不透明度表
fn call_ghydop(&mut self, ij: usize);
/// 调用 LYMLIN(IJ) - Lyman 线近似
fn call_lymlin(&mut self, ij: usize);
/// 调用 OPACT1(IJ) - 表格不透明度
fn call_opact1(&mut self, ij: usize);
/// 调用 OPADD(MODE, ICALL, IJ, ID) - 附加不透明度
/// 返回 (abad, emad, scad)
fn call_opadd(&mut self, mode: i32, icall: i32, ij: usize, id: usize) -> OpaddResult;
/// 调用 PRD(IJ) - 部分重分布
fn call_prd(&mut self, ij: usize);
}
/// 空回调实现（默认不做任何操作）
#[derive(Debug, Clone, Default)]
pub struct NoOpCallbacks;
impl Opacf1Callbacks for NoOpCallbacks {
fn call_quasim(&mut self, _ij: usize) {}
fn call_ghydop(&mut self, _ij: usize) {}
fn call_lymlin(&mut self, _ij: usize) {}
fn call_opact1(&mut self, _ij: usize) {}
fn call_opadd(&mut self, _mode: i32, _icall: i32, _ij: usize, _id: usize) -> OpaddResult {
OpaddResult::default()
}
fn call_prd(&mut self, _ij: usize) {}
}
impl Default for Opacf1Config {
@@ -51,6 +108,8 @@ impl Default for Opacf1Config {
ioptab: 0,
ioplym: 0,
ifprd: 0,
iprcrs: 0,
nprcrs: 0,
}
}
}
@@ -68,8 +127,8 @@ pub struct Opacf1ModelState<'a> {
pub dens: &'a [f64],
/// 密度倒数 (nd)
pub dens1: &'a [f64],
/// 占据数 (mlevel × nd)
pub popul: &'a [f64],
/// 占据数 (mlevel × nd) - iprcrs 块需要修改
pub popul: &'a mut [f64],
// 工作数组
/// HKT1 (nd) - HK/T
@@ -172,8 +231,8 @@ pub struct Opacf1FreqParams<'a> {
#[derive(Debug)]
pub struct Opacf1Precomputed<'a> {
// 跃迁吸收/发射系数 (mtrans × nd)
/// 吸收系数
pub abtra: &'a [f64],
/// 吸收系数 - iprcrs 块需要修改
pub abtra: &'a mut [f64],
/// 发射系数
pub emtra: &'a [f64],
@@ -318,11 +377,10 @@ fn gfree1(id: usize, x: f64, precomp: &Opacf1Precomputed) -> f64 {
}
/// 计算下沉修正 DWNFR1
fn dwnfr1(fr: f64, fr0: f64, id: usize, izz: usize, precomp: &Opacf1Precomputed) -> f64 {
// 简化版本，调用 dwnfr1 模块
// 需要更多参数，这里返回 1.0 作为默认值
let _ = (fr, fr0, id, izz, precomp);
UN
///
/// 调用实际的 dwnfr1 模块
fn dwnfr1(fr: f64, fr0: f64, id: usize, izz: usize, inppar: &InpPar, dwnpar: &DwnPar) -> f64 {
crate::tlusty::math::opacity::dwnfr1(fr, fr0, id, izz, inppar, dwnpar)
}
/// 计算 Mermerges 截面 SGMER1
@@ -355,8 +413,9 @@ fn sgmer1(
/// * `freq_params` - 频率数据
/// * `precomp` - 预计算量
/// * `output` - 输出数组
/// * `cross_fn` - 光电离截面函数
/// * `crossd_fn` - 含双电子复合的光电离截面函数
/// * `inppar` - 输入参数
/// * `dwnpar` - 下沉参数
/// * `callbacks` - 子程序回调接口
///
/// # Fortran 原始代码
///
@@ -367,7 +426,7 @@ fn sgmer1(
/// END
/// ```
#[allow(clippy::too_many_arguments)]
pub fn opacf1(
pub fn opacf1<C: Opacf1Callbacks>(
ij: usize,
config: &Opacf1Config,
model: &mut Opacf1ModelState,
@@ -375,11 +434,15 @@ pub fn opacf1(
freq_params: &Opacf1FreqParams,
precomp: &mut Opacf1Precomputed,
output: &mut Opacf1Output,
inppar: &InpPar,
dwnpar: &DwnPar,
callbacks: &mut C,
) {
let nd = model.nd;
// ========================================================================
// 特殊情况: ioptab < 0 时，只调用 OPACT1
// Fortran lines 25-34
// ========================================================================
if config.ioptab < 0 {
for id in 0..nd {
@@ -388,7 +451,8 @@ pub fn opacf1(
output.emis1[id] = 0.0;
output.absot[id] = 0.0;
}
// CALL OPACT1(IJ) - 由外部处理
// Fortran line 32: call opact1(ij)
callbacks.call_opact1(ij);
return;
}
@@ -443,6 +507,51 @@ pub fn opacf1(
// 检查是否使用表格频率
let lfre = freq_params.freq[ij] > freq_params.frtabm;
// ========================================================================
// iprcrs 预处理 (Fortran lines 67-80)
// 临时保存并清空特定深度点的氢能级占据数和吸收/发射系数
// 用于光致电离截面计算
// ========================================================================
let mut pold: Vec<f64> = Vec::new();
let mut abtrh: Vec<f64> = Vec::new();
if config.iprcrs > 0 {
let iprcrs_idx = (config.iprcrs - 1) as usize;
output.abso1[iprcrs_idx] = 0.0;
if atomic.ielh > 0 {
let ielh_idx = (atomic.ielh - 1) as usize;
let nfirst_h = atomic.nfirst[ielh_idx] as usize;
let nlast_h = atomic.nlast[ielh_idx] as usize;
let nnext_h = atomic.nnext[ielh_idx] as usize;
let nprcrs_offset = (config.nprcrs - 1) as usize;
let ih_exclude = nfirst_h + nprcrs_offset;
// 初始化工作数组
pold.resize(MLEVEL, 0.0);
abtrh.resize(MTRANS, 0.0);
for ii in nfirst_h..=nlast_h {
if ii != ih_exclude {
let ii_idx = ii - 1;
// 保存并清空占据数
pold[ii_idx] = model.popul[ii_idx * MDEPTH + iprcrs_idx];
model.popul[ii_idx * MDEPTH + iprcrs_idx] = 0.0;
// 保存并清空跃迁吸收系数
for jj in (ii + 1)..=nnext_h {
let itrh = atomic.itra[ii_idx * MLEVEL + jj - 1] as usize;
if itrh > 0 {
let itrh_idx = itrh - 1;
abtrh[itrh_idx] = precomp.abtra[itrh_idx * MDEPTH + iprcrs_idx];
precomp.abtra[itrh_idx * MDEPTH + iprcrs_idx] = 0.0;
}
}
}
}
}
}
// ========================================================================
// 1. Bound-free 贡献 (无双电子复合)
// ========================================================================
@@ -475,7 +584,7 @@ pub fn opacf1(
// Macfarlane 下沉修正
let mcdw = atomic.mcdw[itr_idx];
if mcdw > 0 {
let dw1 = dwnfr1(fr, atomic.fr0[itr_idx], id, izz, precomp);
let dw1 = dwnfr1(fr, atomic.fr0[itr_idx], id, izz, inppar, dwnpar);
precomp.dwf1[(mcdw - 1) as usize * MDEPTH + id] = dw1;
sgd *= dw1;
}
@@ -531,7 +640,7 @@ pub fn opacf1(
let mcdw = atomic.mcdw[itr_idx];
if mcdw > 0 {
let dw1 = dwnfr1(fr, atomic.fr0[itr_idx], id, izz, precomp);
let dw1 = dwnfr1(fr, atomic.fr0[itr_idx], id, izz, inppar, dwnpar);
precomp.dwf1[(mcdw - 1) as usize * MDEPTH + id] = dw1;
sgd *= dw1;
}
@@ -629,11 +738,13 @@ pub fn opacf1(
// 3. 附加不透明度 OPADD
// ========================================================================
if config.iopadd != 0 {
let icall = 1;
for id in 0..nd {
// CALL OPADD(0, 1, IJ, ID) - 由外部处理
// 这里需要 abad, emad, scad 的值
// 简化版本：不做处理
let _ = id;
// CALL OPADD(0, ICALL, IJ, ID)
let result = callbacks.call_opadd(0, icall, ij, id);
output.abso1[id] += result.abad;
output.emis1[id] += result.emad;
output.scat1[id] += result.scad;
}
}
@@ -752,8 +863,9 @@ pub fn opacf1(
// ========================================================================
// Quasimolecular opacity (CALL QUASIM(IJ))
// 由外部处理
// Fortran line 323
// ========================================================================
callbacks.call_quasim(ij);
// ========================================================================
// 总不透明度
@@ -766,21 +878,24 @@ pub fn opacf1(
// ========================================================================
// GHYDOP (CALL GHYDOP(IJ))
// 由外部处理
// Fortran line 342
// ========================================================================
callbacks.call_ghydop(ij);
// ========================================================================
// Lyman 线近似 (CALL LYMLIN(IJ))
// Fortran line 346: if(ioplym.gt.0) call lymlin(ij)
// ========================================================================
if config.ioplym > 0 {
// 由外部处理
callbacks.call_lymlin(ij);
}
// ========================================================================
// 表格不透明度 (CALL OPACT1(IJ))
// Fortran line 352-354: if(ioptab.gt.0) then call opact1(ij)
// ========================================================================
if config.ioptab > 0 {
// 由外部处理
callbacks.call_opact1(ij);
}
// ========================================================================
@@ -794,9 +909,49 @@ pub fn opacf1(
// ========================================================================
// PRD (CALL PRD(IJ))
// Fortran line 364: if(ifprd.gt.0) call prd(ij)
// ========================================================================
if config.ifprd > 0 {
// 由外部处理
callbacks.call_prd(ij);
}
// ========================================================================
// iprcrs 恢复 (Fortran lines 366-379)
// 恢复之前保存的占据数和吸收/发射系数
// 计算光致电离截面
// ========================================================================
if config.iprcrs > 0 && atomic.ielh > 0 {
let iprcrs_idx = (config.iprcrs - 1) as usize;
let ielh_idx = (atomic.ielh - 1) as usize;
let nfirst_h = atomic.nfirst[ielh_idx] as usize;
let nlast_h = atomic.nlast[ielh_idx] as usize;
let nnext_h = atomic.nnext[ielh_idx] as usize;
let nprcrs_offset = (config.nprcrs - 1) as usize;
let ih = nfirst_h + nprcrs_offset;
let ih_idx = ih - 1;
// 计算光致电离截面 CRS
// Fortran: crs = abso1(iprcrs)/(popul(ih,iprcrs)*g(ih)*0.0265*4.1347e-15)
let popul_ih = model.popul[ih_idx * MDEPTH + iprcrs_idx];
let g_ih = atomic.g[ih_idx];
let _crs = output.abso1[iprcrs_idx]
/ (popul_ih * g_ih * 0.0265 * 4.1347e-15);
// 恢复占据数和吸收系数
for ii in nfirst_h..=nlast_h {
if ii != ih {
let ii_idx = ii - 1;
model.popul[ii_idx * MDEPTH + iprcrs_idx] = pold[ii_idx];
for jj in (ii + 1)..=nnext_h {
let itrh = atomic.itra[ii_idx * MLEVEL + jj - 1] as usize;
if itrh > 0 {
let itrh_idx = itrh - 1;
precomp.abtra[itrh_idx * MDEPTH + iprcrs_idx] = abtrh[itrh_idx];
}
}
}
}
}
}
@@ -818,7 +973,7 @@ mod tests {
let elec = vec![1e12; nd];
let dens = vec![1e14; nd];
let dens1 = vec![1e-14; nd];
let popul = vec![0.0; MLEVEL * nd];
let mut popul = vec![0.0; MLEVEL * nd];
let hkt1 = vec![0.0; nd];
let mut elscat = vec![0.0; nd];
@@ -828,7 +983,7 @@ mod tests {
elec: &elec,
dens: &dens,
dens1: &dens1,
popul: &popul,
popul: &mut popul,
hkt1: &hkt1,
elscat: &mut elscat,
};
@@ -928,7 +1083,7 @@ mod tests {
g: &g,
};
let abtra = vec![0.0; MTRANS * MDEPTH];
let mut abtra = vec![0.0; MTRANS * MDEPTH];
let emtra = vec![0.0; MTRANS * MDEPTH];
let sff2 = vec![0.0; MION * MDEPTH];
let sff3 = vec![0.0; MION * MDEPTH];
@@ -953,7 +1108,7 @@ mod tests {
let xjid = vec![0.0; MDEPTH];
let mut precomp = Opacf1Precomputed {
abtra: &abtra,
abtra: &mut abtra,
emtra: &emtra,
sff2: &sff2,
sff3: &sff3,
@@ -978,7 +1133,7 @@ mod tests {
xjid: &xjid,
};
opacf1(0, &config, &mut model, &atomic, &freq_params, &mut precomp, &mut output);
opacf1(0, &config, &mut model, &atomic, &freq_params, &mut precomp, &mut output, &InpPar::default(), &DwnPar::default(), &mut NoOpCallbacks);
// ioptab < 0 时，输出应该为 0
for id in 0..nd {
@@ -1007,7 +1162,7 @@ mod tests {
gf5[0] = -0.1;
gf6[0] = 0.5;
let abtra = vec![0.0; MTRANS * MDEPTH];
let mut abtra = vec![0.0; MTRANS * MDEPTH];
let emtra = vec![0.0; MTRANS * MDEPTH];
let sff2 = vec![0.0; MION * MDEPTH];
let sff3 = vec![0.0; MION * MDEPTH];
@@ -1025,7 +1180,7 @@ mod tests {
let xjid = vec![0.0; MDEPTH];
let precomp = Opacf1Precomputed {
abtra: &abtra,
abtra: &mut abtra,
emtra: &emtra,
sff2: &sff2,
sff3: &sff3,
src/tlusty/math/continuum/opacity_table.rs
+824
View File
@@ -0,0 +1,824 @@
//! 不透明度表包装器，用于 LTE 灰大气模型。
//!
//! 提供简化的接口来读取和使用预计算的不透明度表。
//! 支持两种格式：
//! - 文本格式（用于调试）
//! - Fortran 二进制格式（实际 OPCTAB 文件）
use std::fs::File;
use std::io::{BufRead, BufReader, Read, Write};
use std::path::Path;
use crate::tlusty::io::{FortranReader, IoError, Result};
use crate::tlusty::state::constants::{HK, UN};
// ============================================================================
// Fortran 二进制读取器
// ============================================================================
/// Fortran 无格式二进制读取器。
///
/// 处理 Fortran unformatted 文件的记录标记。
struct FortranBinaryReader<R: Read> {
reader: R,
}
impl<R: Read> FortranBinaryReader<R> {
fn new(reader: R) -> Self {
Self { reader }
}
/// 读取一个记录，返回记录内容的字节。
fn read_record(&mut self) -> std::io::Result<Vec<u8>> {
// 读取记录长度 (4字节 little-endian 整数)
let mut len_bytes = [0u8; 4];
self.reader.read_exact(&mut len_bytes)?;
let record_len = i32::from_le_bytes(len_bytes) as usize;
// 读取记录内容
let mut buffer = vec![0u8; record_len];
self.reader.read_exact(&mut buffer)?;
// 读取尾部记录长度（应该相同）
let mut tail_bytes = [0u8; 4];
self.reader.read_exact(&mut tail_bytes)?;
let tail_len = i32::from_le_bytes(tail_bytes) as usize;
if tail_len != record_len {
return Err(std::io::Error::new(
std::io::ErrorKind::InvalidData,
format!("Record length mismatch: {} != {}", record_len, tail_len),
));
}
Ok(buffer)
}
/// 从记录中提取 f64 数组。
fn read_f64_array(&mut self) -> std::io::Result<Vec<f64>> {
let record = self.read_record()?;
let n = record.len() / 8;
let mut result = Vec::with_capacity(n);
for i in 0..n {
let bytes: [u8; 8] = record[i * 8..(i + 1) * 8].try_into().unwrap();
result.push(f64::from_le_bytes(bytes));
}
Ok(result)
}
/// 从记录中提取 f32 数组。
fn read_f32_array(&mut self) -> std::io::Result<Vec<f32>> {
let record = self.read_record()?;
let n = record.len() / 4;
let mut result = Vec::with_capacity(n);
for i in 0..n {
let bytes: [u8; 4] = record[i * 4..(i + 1) * 4].try_into().unwrap();
result.push(f32::from_le_bytes(bytes));
}
Ok(result)
}
/// 从记录中提取 i32 数组。
fn read_i32_array(&mut self) -> std::io::Result<Vec<i32>> {
let record = self.read_record()?;
let n = record.len() / 4;
let mut result = Vec::with_capacity(n);
for i in 0..n {
let bytes: [u8; 4] = record[i * 4..(i + 1) * 4].try_into().unwrap();
result.push(i32::from_le_bytes(bytes));
}
Ok(result)
}
/// 读取包含字符串和两个 f64 的记录（丰度格式）。
fn read_abundance_record(&mut self) -> std::io::Result<(String, f64, f64)> {
let record = self.read_record()?;
// 格式: 4字节字符 + 8字节 f64 + 8字节 f64
let typa = String::from_utf8_lossy(&record[0..4]).trim().to_string();
let abunt = f64::from_le_bytes(record[4..12].try_into().unwrap());
let abuno = f64::from_le_bytes(record[12..20].try_into().unwrap());
Ok((typa, abunt, abuno))
}
}
// ============================================================================
// 物理常数
// ============================================================================
/// 光速 (cm/s)
const CLIGHT: f64 = 2.99792458e10;
/// Stefan-Boltzmann 常数
const SIGMAC: f64 = 7.5657e-15;
/// 氢电离频率 (Hz)
const FRH: f64 = 3.28805e15;
// ============================================================================
// 数据结构
// ============================================================================
/// 拥有所有权的不透明度表数据。
#[derive(Debug, Clone)]
pub struct OwnedOpacityTable {
/// 温度向量 (ln T) [numtemp]
pub tempvec: Vec<f64>,
/// 密度矩阵 (ln rho) [numtemp][numrh]
pub rhomat: Vec<Vec<f64>>,
/// 频率表 [numfreq]
pub frtab: Vec<f64>,
/// 不透明度表 (ln kappa) [温度][密度][频率]
pub absopac: Vec<Vec<Vec<f64>>>,
/// 电子密度网格 [温度][密度]
pub elecgr: Vec<Vec<f64>>,
/// 每个温度点的密度数
pub numrh: Vec<usize>,
/// 频率数
pub numfreq: usize,
/// 温度数
pub numtemp: usize,
/// 表格边界
pub ttab1: f64,
pub ttab2: f64,
pub rtab1: f64,
pub rtab2: f64,
/// 丰度信息
pub abundances: Vec<f64>,
/// 不透明度标志
pub iophmt: i32,
pub ioph2t: i32,
pub iophet: i32,
}
impl Default for OwnedOpacityTable {
fn default() -> Self {
Self {
tempvec: Vec::new(),
rhomat: Vec::new(),
frtab: Vec::new(),
absopac: Vec::new(),
elecgr: Vec::new(),
numrh: Vec::new(),
numfreq: 0,
numtemp: 0,
ttab1: 0.0,
ttab2: 0.0,
rtab1: 0.0,
rtab2: 0.0,
abundances: vec![0.0; 92],
iophmt: 0,
ioph2t: 0,
iophet: 0,
}
}
}
impl OwnedOpacityTable {
/// 创建一个空的表。
pub fn new() -> Self {
Self::default()
}
/// 从文件加载表（自动检测格式）。
pub fn from_file<P: AsRef<Path>>(path: P) -> Result<Self> {
let path_ref = path.as_ref();
let file = File::open(path_ref).map_err(|e| {
IoError::FileNotFound(format!(
"Cannot open opacity table '{}': {}",
path_ref.display(),
e
))
})?;
// 检测文件格式：二进制文件以记录长度标记开始
let mut reader = BufReader::new(file);
let mut first_bytes = [0u8; 4];
use std::io::Read;
let bytes_read = reader.read(&mut first_bytes).map_err(|e| {
IoError::ParseError(format!("Cannot read file header: {}", e))
})?;
if bytes_read < 4 {
return Err(IoError::ParseError("File too small".to_string()));
}
// 重新打开文件用于完整读取
let file2 = File::open(path_ref).map_err(|e| {
IoError::FileNotFound(format!(
"Cannot reopen opacity table '{}': {}",
path_ref.display(),
e
))
})?;
// 检查是否为二进制格式
// 二进制文件第一个记录通常长度为 20 (4+8+8 for abundance)
let first_len = i32::from_le_bytes(first_bytes);
if first_len > 0 && first_len < 10000 {
Self::from_binary_reader(file2)
} else {
let mut fortran_reader = FortranReader::new(BufReader::new(file2));
Self::read_table(&mut fortran_reader)
}
}
/// 从 Fortran 二进制文件加载表。
pub fn from_binary<P: AsRef<Path>>(path: P) -> Result<Self> {
let file = File::open(path.as_ref()).map_err(|e| {
IoError::FileNotFound(format!(
"Cannot open binary opacity table '{}': {}",
path.as_ref().display(),
e
))
})?;
Self::from_binary_reader(file)
}
/// 从已打开的 Fortran 二进制文件读取。
fn from_binary_reader(file: File) -> Result<Self> {
let mut reader = FortranBinaryReader::new(BufReader::new(file));
let mut table = Self::new();
// 读取 92 个丰度记录
table.abundances = Vec::with_capacity(92);
for _ in 0..92 {
let (_typa, abunt, _abuno) = reader
.read_abundance_record()
.map_err(|e| IoError::ParseError(format!("Failed to read abundance: {}", e)))?;
table.abundances.push(abunt);
}
// 读取分子标志
let mol_flags = reader
.read_record()
.map_err(|e| IoError::ParseError(format!("Failed to read mol flags: {}", e)))?;
// ifmolt (i32) + tmolit (f64)
let _ifmolt = i32::from_le_bytes(mol_flags[0..4].try_into().unwrap());
let _tmolit = f64::from_le_bytes(mol_flags[4..12].try_into().unwrap());
// 读取不透明度标志
let opac_flags = reader
.read_i32_array()
.map_err(|e| IoError::ParseError(format!("Failed to read opac flags: {}", e)))?;
table.iophmt = opac_flags.first().copied().unwrap_or(0);
table.ioph2t = opac_flags.get(1).copied().unwrap_or(0);
table.iophet = opac_flags.get(2).copied().unwrap_or(0);
// 读取维度
let dims = reader
.read_i32_array()
.map_err(|e| IoError::ParseError(format!("Failed to read dimensions: {}", e)))?;
let numfre0 = dims.first().copied().unwrap_or(0) as usize;
let numtem0 = dims.get(1).copied().unwrap_or(0) as usize;
let numrh0 = dims.get(2).copied().unwrap_or(0);
// 读取温度向量
table.tempvec = reader
.read_f64_array()
.map_err(|e| IoError::ParseError(format!("Failed to read tempvec: {}", e)))?;
table.numtemp = table.tempvec.len();
// 读取密度向量
let rhov = reader
.read_f64_array()
.map_err(|e| IoError::ParseError(format!("Failed to read rhov: {}", e)))?;
let numrh_val = rhov.len();
// 读取电子密度网格
let elec_data = reader
.read_f64_array()
.map_err(|e| IoError::ParseError(format!("Failed to read elecgr: {}", e)))?;
// 转换为二维数组
table.elecgr = Vec::with_capacity(table.numtemp);
for i in 0..table.numtemp {
let row: Vec<f64> = elec_data[i * numrh_val..(i + 1) * numrh_val].to_vec();
table.elecgr.push(row);
}
// 设置密度矩阵
table.rhomat = vec![rhov.clone(); table.numtemp];
table.numrh = vec![numrh_val; table.numtemp];
// 设置边界
if table.numtemp > 0 {
table.ttab1 = table.tempvec[0];
table.ttab2 = table.tempvec[table.numtemp - 1];
}
if numrh_val > 0 {
table.rtab1 = rhov[0];
table.rtab2 = rhov[numrh_val - 1];
}
// 读取频率和不透明度数据
table.frtab = Vec::with_capacity(numfre0);
// 预分配不透明度数组：[温度][密度][频率]
table.absopac = vec![vec![vec![0.0_f64; numfre0]; numrh_val]; table.numtemp];
for kf in 0..numfre0 {
// 读取频率值
let freq_record = reader
.read_f64_array()
.map_err(|e| IoError::ParseError(format!("Failed to read freq {}: {}", kf, e)))?;
let frta = freq_record.first().copied().unwrap_or(0.0);
table.frtab.push(frta);
// 读取每个密度点的不透明度
for jd in 0..numrh_val {
let opac_row = reader
.read_f32_array()
.map_err(|e| {
IoError::ParseError(format!("Failed to read opac row {},{}: {}", kf, jd, e))
})?;
// 存储为 f64，转换为 ln
for (it, &val) in opac_row.iter().enumerate() {
if it < table.numtemp {
table.absopac[it][jd][kf] = val as f64;
}
}
}
}
table.numfreq = table.frtab.len();
Ok(table)
}
/// 读取表数据（文本格式）。
fn read_table<R: BufRead>(reader: &mut FortranReader<R>) -> Result<Self> {
let mut table = Self::new();
// 跳过第一行标题 "opacity table with element abundances:"
reader.skip_line()?;
// 跳过第二行标题 "element for EOS for opacities"
reader.skip_line()?;
// 读取丰度数据（92 个元素）
// 格式: 每行 "元素名 丰度1 丰度2"
// 我们需要第三列（丰度2）
table.abundances = Vec::with_capacity(92);
for _ in 0..92 {
// 读取三个值：元素名（跳过），丰度1（跳过），丰度2（保存）
let _elem: String = reader.read_string()?;
let _abun1: f64 = reader.read_value()?;
let abun2: f64 = reader.read_value()?;
table.abundances.push(abun2);
}
// 跳过 "molecules - ifmol,tmolim:" 行
reader.skip_line()?;
// 读取分子标志值
let mol_flags = reader.read_values::<f64>(2)?;
let _ifmolt = mol_flags[0] as i32;
let _tmolit = mol_flags[1];
// 跳过 "additional opacities" 行
reader.skip_line()?;
// 跳过元素名行 (H- H2+ He- ...)
reader.skip_line()?;
// 读取附加不透明度标志
let opac_flags = reader.read_values::<i32>(10)?;
table.iophmt = opac_flags[0];
table.ioph2t = opac_flags[1];
table.iophet = opac_flags[2];
// 跳过 "number of frequencies, temperatures, densities:" 行
reader.skip_line()?;
// 读取维度
let dims = reader.read_values::<i32>(3)?;
let numfre0 = dims[0] as usize;
let numtem0 = dims[1] as usize;
let numrh0 = dims[2].abs() as usize;
// 跳过 "log temperatures" 行
reader.skip_line()?;
// 读取温度向量 (ln T)
table.tempvec = reader.read_array::<f64>(numtem0)?;
table.numtemp = numtem0;
// 跳过 "log densities" 行
reader.skip_line()?;
// 读取密度向量 (ln rho)
let rhov = reader.read_array::<f64>(numrh0)?;
// 跳过 "log electron densities from EOS" 行
reader.skip_line()?;
// 读取电子密度网格 - 存储为 [温度][密度]
let elec_flat = reader.read_array::<f64>(numtem0 * numrh0)?;
table.elecgr = Vec::with_capacity(numtem0);
for it in 0..numtem0 {
let row: Vec<f64> = elec_flat[it * numrh0..(it + 1) * numrh0].to_vec();
table.elecgr.push(row);
}
// 设置密度矩阵（所有温度使用相同的密度网格）
table.rhomat = vec![rhov.clone(); numtem0];
table.numrh = vec![numrh0; numtem0];
// 设置边界
if table.numtemp > 0 {
table.ttab1 = table.tempvec[0];
table.ttab2 = table.tempvec[table.numtemp - 1];
}
if numrh0 > 0 {
table.rtab1 = rhov[0];
table.rtab2 = rhov[numrh0 - 1];
}
// 读取频率和不透明度数据
table.frtab = Vec::with_capacity(numfre0);
table.absopac = vec![vec![vec![0.0; numfre0]; numrh0]; numtem0];
for kf in 0..numfre0 {
// 跳过 "*** frequency # : n value" 行
reader.skip_line()?;
// 读取频率值
let frta: f64 = reader.read_value()?;
table.frtab.push(frta);
// 读取不透明度值
// 格式：按温度分组，每个温度有 numrh0 个密度值
// 存储：[温度][密度][频率]
let opac_flat = reader.read_array::<f64>(numtem0 * numrh0)?;
for it in 0..numtem0 {
for jd in 0..numrh0 {
table.absopac[it][jd][kf] = opac_flat[it * numrh0 + jd];
}
}
}
table.numfreq = table.frtab.len();
Ok(table)
}
/// 插值获取给定温度和密度下的不透明度（单频率）。
///
/// # 参数
/// * `t` - 温度 (K)
/// * `rho` - 密度 (g/cm³)
/// * `ifreq` - 频率索引 (0-indexed)
///
/// # 返回
/// 不透明度 (cm²/g)
pub fn interpolate_opacity(&self, t: f64, rho: f64, ifreq: usize) -> f64 {
if self.numtemp == 0 || self.numfreq == 0 {
return 0.4; // 默认电子散射
}
let tl = t.ln();
// 温度插值
let jt = if tl <= self.ttab1 {
0
} else if tl >= self.ttab2 {
self.numtemp - 2
} else {
let deltat = (tl - self.ttab1) / (self.ttab2 - self.ttab1) * (self.numtemp - 1) as f64;
(deltat as usize).min(self.numtemp - 2)
};
let t1 = self.tempvec[jt];
let t2 = self.tempvec[jt + 1];
let dti = if t2 > t1 {
((tl - t1) / (t2 - t1)).clamp(0.0, 1.0)
} else {
0.0
};
// 密度插值（在两个温度点）
let rl = rho.ln();
// 下温度点的密度插值
let numr = self.numrh[jt].max(1);
let rtab1 = self.rhomat[jt][0];
let rtab2 = self.rhomat[jt][numr - 1];
let jr1 = if rl <= rtab1 {
0
} else if rl >= rtab2 {
numr - 2
} else {
let deltar = (rl - rtab1) / (rtab2 - rtab1) * (numr - 1) as f64;
(deltar as usize).min(numr - 2)
};
let r1_1 = self.rhomat[jt][jr1];
let r2_1 = self.rhomat[jt][jr1 + 1];
let dri1 = if r2_1 > r1_1 {
((rl - r1_1) / (r2_1 - r1_1)).clamp(0.0, 1.0)
} else {
0.0
};
// 上温度点的密度插值
let numr2 = self.numrh[jt + 1].max(1);
let rtab1_2 = self.rhomat[jt + 1][0];
let rtab2_2 = self.rhomat[jt + 1][numr2 - 1];
let jr2 = if rl <= rtab1_2 {
0
} else if rl >= rtab2_2 {
numr2 - 2
} else {
let deltar = (rl - rtab1_2) / (rtab2_2 - rtab1_2) * (numr2 - 1) as f64;
(deltar as usize).min(numr2 - 2)
};
let r1_2 = self.rhomat[jt + 1][jr2];
let r2_2 = self.rhomat[jt + 1][jr2 + 1];
let dri2 = if r2_2 > r1_2 {
((rl - r1_2) / (r2_2 - r1_2)).clamp(0.0, 1.0)
} else {
0.0
};
// 插值不透明度
let op1 = self.absopac[jt][jr1][ifreq]
+ dri1 * (self.absopac[jt][jr1 + 1][ifreq] - self.absopac[jt][jr1][ifreq]);
let op2 = self.absopac[jt + 1][jr2][ifreq]
+ dri2 * (self.absopac[jt + 1][jr2 + 1][ifreq] - self.absopac[jt + 1][jr2][ifreq]);
let opac = op1 + dti * (op2 - op1);
opac.exp()
}
/// 检查温度是否在表范围内
pub fn is_temperature_in_range(&self, t: f64) -> bool {
if self.numtemp == 0 {
return false;
}
let tl = t.ln();
tl >= self.ttab1 && tl <= self.ttab2
}
/// 当没有表数据或温度超出范围时，使用近似公式计算平均不透明度。
///
/// 参考: Mihalas "Stellar Atmospheres" 第 3 章; Seaton et al. (1994) Opacity Project
///
/// 对于热星 (T > 15000 K):
/// - 电子散射是主要连续不透明度
/// - 束缚-自由来自 H, He 的高能级
/// - 自由-自由来自 H⁺, He⁺, He²⁺
fn approximate_mean_opacity(&self, t: f64, rho: f64, ne: f64) -> (f64, f64) {
// 物理常数
const SIGE: f64 = 6.6524e-25; // Thomson 散射截面 (cm²)
const FRH: f64 = 3.28805e15; // 氢电离频率 (Hz)
// 密度保护
let rho_safe = rho.max(1e-30);
// 1. 电子散射不透明度 (主要)
let kappa_es = SIGE * ne / rho_safe;
// 2. 氢/氦束缚-自由不透明度
// 对于热星，使用改进的 Kramers 公式，考虑部分电离
// κ_bf ≈ 4.34e25 * g_bf * (1-X_ion) * ρ * T^-3.5
// 其中 g_bf ≈ 1 是束缚-自由 Gaunt 因子
// 温度因子 (相对于 10^4 K)
let t4 = t / 1e4;
// 束缚-自由 (Kramers)
// 对于 Rosseland 平均，使用频率积分结果
// κ_R_bf ≈ 4.34e25 * Z * (1+X) * ρ * T^-3.5 / ⟨g_bf⟩
// 简化: 使用典型丰度 X=0.7, Z=0.02
let kappa_bf = 4.34e-25 * (1.0 + 0.7) * (1.0 / t4).powf(3.5);
// 3. 自由-自由不透明度 (Kramers with Gaunt)
// κ_ff ≈ 3.68e8 * (X + Y/4) * (1+X) * ρ * T^-7/2 * ⟨g_ff⟩
// 对于热星，自由-自由贡献较小
let kappa_ff = 1.0e-26 * (1.0 / t4).powf(3.5);
// 4. H⁻ 自由-自由 (仅对较冷恒星重要)
// 对于 T > 15000 K, H⁻ 很少
let kappa_hm = if t < 12000.0 {
1.0e-26 * (12000.0 / t).powf(5.0)
} else {
0.0
};
// Rosseland 平均: 1/κ_R = ∫(1/κ)(dB/dT)dν / ∫(dB/dT)dν
// 对于混合不透明度，使用调和平均的近似
let kappa_abs = kappa_bf + kappa_ff + kappa_hm;
let kappa_ros = kappa_es + kappa_abs;
// Planck 平均: κ_P = ∫κBdν / ∫Bdν
// Planck 平均偏向高不透明度区域
let kappa_pla = kappa_es + 2.0 * kappa_abs;
(kappa_ros, kappa_pla)
}
/// 检查密度是否在表范围内
pub fn is_density_in_range(&self, rho: f64) -> bool {
if self.rhomat.is_empty() || self.rhomat[0].is_empty() {
return false;
}
let rl = rho.ln();
rl >= self.rtab1 && rl <= self.rtab2
}
/// 计算 Rosseland 平均不透明度。
///
/// 使用频率积分:
/// κ_R = ∫(dB/dT)dν / ∫(1/κ)(dB/dT)dν
///
/// # 参数
/// * `t` - 温度 (K)
/// * `rho` - 密度 (g/cm³)
/// * `ne` - 电子密度 (cm⁻³) - 用于电子散射
///
/// # 返回
/// (Rosseland 平均, Planck 平均)
pub fn mean_opacity(&self, t: f64, rho: f64, ne: f64) -> (f64, f64) {
if self.numfreq == 0 || !self.is_temperature_in_range(t) {
// 没有表数据或温度超出范围，使用改进的近似公式
return self.approximate_mean_opacity(t, rho, ne);
}
let hkt = HK / t;
let mut abr = 0.0; // ∫(1/κ)(dB/dT)dν 的分母
let mut sumdb = 0.0; // ∫(dB/dT)dν
let mut abp = 0.0; // ∫κBdν
let mut sumb = 0.0; // ∫Bdν
// 电子散射系数
let sige = 6.6524e-25;
let kappa_es = sige * ne / rho.max(1e-30);
for ij in 0..self.numfreq {
let fr = self.frtab[ij];
// Planck 函数权重
let ex = (hkt * fr).exp();
let e1 = UN / (ex - UN);
// 简化的 Planck 函数 B_ν ∝ ν³
let bnue = 1.4743e-2 * fr.powi(3);
// 梯形积分权重
let w = if ij == 0 {
0.5 * (self.frtab[0] - self.frtab[1])
} else if ij == self.numfreq - 1 {
0.5 * (self.frtab[ij - 1] - self.frtab[ij])
} else {
0.5 * (self.frtab[ij - 1] - self.frtab[ij + 1])
};
let plan = bnue * e1 * w;
let dplan = plan * hkt * fr * ex * e1;
// 从表插值获取吸收不透明度
let ab = self.interpolate_opacity(t, rho, ij);
// 总散射 = 电子散射
let sct = kappa_es;
// 累加 Rosseland 平均
abr += dplan / (ab + sct);
sumdb += dplan;
// 累加 Planck 平均
abp += plan * ab;
sumb += plan;
}
let opros = sumdb / abr;
let oppla = abp / sumb;
(opros, oppla)
}
}
/// 创建一个默认的简化不透明度表（用于测试）。
pub fn create_simple_table() -> OwnedOpacityTable {
let mut table = OwnedOpacityTable::new();
// 设置简单的温度网格 (log T: 3.5 - 5.0)
table.tempvec = vec![
3.5_f64.ln(), 3.7_f64.ln(), 3.9_f64.ln(), 4.1_f64.ln(), 4.3_f64.ln(),
4.5_f64.ln(), 4.7_f64.ln(), 4.9_f64.ln(), 5.0_f64.ln(),
];
table.numtemp = table.tempvec.len();
table.ttab1 = table.tempvec[0];
table.ttab2 = table.tempvec[table.numtemp - 1];
// 设置简单的密度网格 (log rho: -16 to -6)
let rho_vals: Vec<f64> = (-16..=-6)
.map(|i| (i as f64).ln())
.collect();
table.rhomat = vec![rho_vals.clone(); table.numtemp];
table.numrh = vec![rho_vals.len(); table.numtemp];
table.rtab1 = rho_vals[0];
table.rtab2 = rho_vals[rho_vals.len() - 1];
// 设置简单的频率网格
table.frtab = vec![
1e13, 3e13, 1e14, 3e14, 1e15, 3e15, 1e16, 3e16
];
table.numfreq = table.frtab.len();
// 设置简单的不透明度（仅电子散射 + Kramers 近似）
let numr = rho_vals.len();
table.absopac = vec![vec![vec![0.0_f64; table.numfreq]; numr]; table.numtemp];
for it in 0..table.numtemp {
let t = table.tempvec[it].exp();
for ir in 0..numr {
let rho = table.rhomat[it][ir].exp();
for ifreq in 0..table.numfreq {
let fr = table.frtab[ifreq];
// 简单的 Kramers 不透明度近似
let kappa_bf = 4.3e-25 * (1e15 / fr).powf(3.0) * (1e4 / t).powf(3.5);
let kappa_es = 0.4; // 电子散射
let kappa = kappa_bf + kappa_es;
table.absopac[it][ir][ifreq] = kappa.ln();
}
}
}
table
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_create_simple_table() {
let table = create_simple_table();
assert!(table.numtemp > 0);
assert!(table.numfreq > 0);
}
#[test]
fn test_interpolate_opacity() {
let table = create_simple_table();
let kappa = table.interpolate_opacity(10000.0, 1e-10, 3);
assert!(kappa > 0.0);
assert!(kappa.is_finite());
}
#[test]
fn test_mean_opacity() {
let table = create_simple_table();
let (opros, oppla) = table.mean_opacity(10000.0, 1e-10, 1e8);
assert!(opros > 0.0, "Rosseland mean should be positive");
assert!(oppla > 0.0, "Planck mean should be positive");
assert!(opros.is_finite(), "Rosseland mean should be finite");
assert!(oppla.is_finite(), "Planck mean should be finite");
}
#[test]
fn test_read_text_opctab() {
// 尝试读取实际的 OPCTAB 文本文件
let optable_path = "/home/fmq/program/tlusty/optables/optab11_7f.dat";
if !std::path::Path::new(optable_path).exists() {
eprintln!("Skipping test: OPCTAB file not found at {}", optable_path);
return;
}
let table = OwnedOpacityTable::from_file(optable_path)
.expect("Failed to read text OPCTAB");
println!("OPCTAB loaded successfully:");
println!(" numtemp = {}", table.numtemp);
println!(" numfreq = {}", table.numfreq);
println!(" numrh = {:?}", table.numrh);
println!(" T range: {:.1} - {:.1} K", table.ttab1.exp(), table.ttab2.exp());
println!(" rho range: {:.2e} - {:.2e} g/cm³", table.rtab1.exp(), table.rtab2.exp());
assert!(table.numtemp > 0, "Should have temperature points");
assert!(table.numfreq > 0, "Should have frequency points");
assert!(!table.numrh.is_empty(), "Should have density points");
// 测试插值
if table.numtemp > 1 && table.numrh[0] > 1 {
let kappa = table.interpolate_opacity(15000.0, 1e-10, 0);
println!(" kappa(15000K, 1e-10, freq[0]) = {:.4e} cm²/g", kappa);
assert!(kappa > 0.0);
}
// 测试 Rosseland 平均
let (opros, oppla) = table.mean_opacity(15000.0, 1e-10, 1e8);
println!(" Rosseland mean (15000K, 1e-10) = {:.4e} cm²/g", opros);
println!(" Planck mean (15000K, 1e-10) = {:.4e} cm²/g", oppla);
assert!(opros > 0.0);
assert!(oppla > 0.0);
}
}
src/tlusty/math/continuum/opadd.rs
+8 -8
View File
@@ -489,10 +489,10 @@ mod tests {
],
popul: vec![vec![0.0; 3]; 10],
cross: vec![vec![0.0; 3]; 20],
cia_h2h2_data: cia_h2h2::CiaH2h2Data::default(),
cia_h2he_data: cia_h2he::CiaH2heData::default(),
cia_h2h_data: cia_h2h::CiaH2hData::default(),
cia_hhe_data: cia_hhe::CiaHheData::default(),
cia_h2h2_data: CiaH2h2Data::default(),
cia_h2he_data: CiaH2heData::default(),
cia_h2h_data: CiaH2hData::default(),
cia_hhe_data: CiaHheData::default(),
}
});
@@ -505,10 +505,10 @@ mod tests {
anmol: Vec<Vec<f64>>,
popul: Vec<Vec<f64>>,
cross: Vec<Vec<f64>>,
cia_h2h2_data: cia_h2h2::CiaH2h2Data,
cia_h2he_data: cia_h2he::CiaH2heData,
cia_h2h_data: cia_h2h::CiaH2hData,
cia_hhe_data: cia_hhe::CiaHheData,
cia_h2h2_data: CiaH2h2Data,
cia_h2he_data: CiaH2heData,
cia_h2h_data: CiaH2hData,
cia_hhe_data: CiaHheData,
}
/// 创建测试用的模型状态
src/tlusty/math/eos/eldens.rs
+103 -6
View File
@@ -218,9 +218,26 @@ pub fn eldens_pure(params: &EldensParams, ipri: i32) -> EldensOutput {
// Newton-Raphson 迭代
for it in 1..=10 {
// 调用 STATE 计算总电荷
// 简化：使用估计值
// 注意：STATE 需要当前的电子密度 ane，而不是初始估计值
if let Some(ref state_params) = params.state_params {
let state_output = state_pure(state_params);
// 创建新的 StateParams，使用当前迭代的 ane 值
let updated_params = StateParams {
mode: state_params.mode,
id: state_params.id,
t: state_params.t,
ane, // 使用当前迭代的电子密度
natoms: state_params.natoms,
hpop: state_params.hpop,
dens: state_params.dens,
wmm: state_params.wmm,
ytot: state_params.ytot,
abndd: state_params.abndd,
ioniz: state_params.ioniz,
irefa: state_params.irefa,
lgr: state_params.lgr,
lrm: state_params.lrm,
};
let state_output = state_pure(&updated_params);
q = state_output.q;
dqn = state_output.dqn;
}
@@ -240,6 +257,10 @@ pub fn eldens_pure(params: &EldensParams, ipri: i32) -> EldensOutput {
let fe = d_val / a_val + q;
ah = ane / fe;
anh = ah * f1;
if params.id == 1 {
eprintln!("DEBUG ELDENS init T>9000: f1={:.6e}, fe={:.6e}, q={:.6e}, ah={:.6e}, anh={:.6e}",
f1, fe, q, ah, anh);
}
} else if t > 4000.0 {
let e_val = g2 * qp / q2;
let b_val = TWO * (UN + e_val);
@@ -252,6 +273,9 @@ pub fn eldens_pure(params: &EldensParams, ipri: i32) -> EldensOutput {
let fe = f1 * d_val + e_val * (UN - a_val * f1) / b_val + q;
ah = ane / fe;
anh = ah * f1;
if params.id == 1 {
eprintln!("DEBUG ELDENS init 4000<T<9000: ah={:.6e}, anh={:.6e}", ah, anh);
}
} else {
// 低温情况
let c1 = q2 * (TWO * params.ytot - UN);
@@ -274,6 +298,13 @@ pub fn eldens_pure(params: &EldensParams, ipri: i32) -> EldensOutput {
let e_val = anh * q2;
let b_val = anh * qmi;
if params.id == 1 && it <= 2 {
eprintln!("DEBUG ELDENS coeffs iter {}: ae={:.6e}, gg={:.6e}, e_val={:.6e}, b_val={:.6e}",
it, ae, gg, e_val, b_val);
eprintln!("DEBUG ELDENS params iter {}: d_val={:.6e}, g2={:.6e}, a_val={:.6e}",
it, d_val, g2, a_val);
}
// 构建线性方程组
let mut r = [[0.0; 3]; 3];
let mut s = [0.0; 3];
@@ -283,6 +314,10 @@ pub fn eldens_pure(params: &EldensParams, ipri: i32) -> EldensOutput {
r[0][1] = 0.0;
r[0][2] = UN;
s[0] = an - ane - params.ytot * ah;
if params.id == 1 && it <= 2 {
eprintln!("DEBUG ELDENS rhs0: an={:.6e}, ane={:.6e}, ytot={:.6e}, ah={:.6e}, s0={:.6e}",
an, ane, params.ytot, ah, s[0]);
}
} else {
r[0][0] = params.ytot - UN;
r[0][1] = a_val + e_val + gg;
@@ -301,21 +336,57 @@ pub fn eldens_pure(params: &EldensParams, ipri: i32) -> EldensOutput {
s[1] = anh * (d_val + gg) + q * ah - ane;
s[2] = ah - anh * (a_val + TWO * (e_val + gg));
if params.id == 1 && it <= 2 {
eprintln!("DEBUG ELDENS rhs12: anh={:.6e}, d_val={:.6e}, gg={:.6e}, q={:.6e}, ah={:.6e}, ane={:.6e}",
anh, d_val, gg, q, ah, ane);
eprintln!("DEBUG ELDENS rhs12: s1={:.6e}, s2={:.6e}", s[1], s[2]);
}
// Debug: print matrix
if params.id == 1 && it <= 2 {
eprintln!("DEBUG ELDENS matrix iter {}: R=[{:.15e},{:.15e},{:.15e}; {:.15e},{:.15e},{:.15e}; {:.15e},{:.15e},{:.15e}]",
it, r[0][0], r[0][1], r[0][2], r[1][0], r[1][1], r[1][2], r[2][0], r[2][1], r[2][2]);
eprintln!("DEBUG ELDENS rhs iter {}: S=[{:.15e},{:.15e},{:.15e}]",
it, s[0], s[1], s[2]);
}
// 求解线性方程组
let r_flat: Vec<f64> = r.iter().flat_map(|row| row.iter().copied()).collect();
let mut r_work = r_flat;
// 注意：LINEQS 需要 Fortran 列优先存储
let mut r_work = vec![0.0; 9];
for j in 0..3 {
for i in 0..3 {
r_work[i + j * 3] = r[i][j]; // 列优先存储
}
}
let mut s_work = s.to_vec();
let mut p = [0.0; 3];
lineqs(&mut r_work, &mut s_work, &mut p, 3);
// 验证解是否正确
if params.id == 1 && it <= 2 {
let check0 = r[0][0]*p[0] + r[0][1]*p[1] + r[0][2]*p[2];
let check1 = r[1][0]*p[0] + r[1][1]*p[1] + r[1][2]*p[2];
let check2 = r[2][0]*p[0] + r[2][1]*p[1] + r[2][2]*p[2];
eprintln!("DEBUG ELDENS verify: R*p=[{:.6e}, {:.6e}, {:.6e}], S=[{:.6e}, {:.6e}, {:.6e}]",
check0, check1, check2, s[0], s[1], s[2]);
}
if params.id == 1 {
eprintln!("DEBUG ELDENS lineqs: p=[{:.6e}, {:.6e}, {:.6e}]", p[0], p[1], p[2]);
}
// 更新值
ah = ah + p[0];
anh = anh + p[1];
let delne = p[2];
ane = ane + delne;
// 收敛检查
// 收敛检查（与 Fortran 一致）
// 注意：Fortran 只检查 ANE <= 0，不检查 ANH
if ane <= 0.0 {
ane = 1e-6 * an;
}
if (delne / ane).abs() <= 1e-3 {
break;
}
@@ -356,6 +427,10 @@ pub fn eldens_pure(params: &EldensParams, ipri: i32) -> EldensOutput {
}
}
if params.id == 1 {
eprintln!("DEBUG ELDENS after loop: ane={:.6e}, an={:.6e}", ane, an);
}
// 确保电子密度为正
if ane <= 0.0 {
ane = 1e-6 * an;
@@ -365,7 +440,24 @@ pub fn eldens_pure(params: &EldensParams, ipri: i32) -> EldensOutput {
anerel = ane / an;
let ahtot = ah;
let ahmol = anh * anh * q2;
let anp = anh / ane * qh;
// 计算质子密度 (从 Saha 方程)
// Saha: n_p * n_e / n_H = qh
// 所以: n_p = n_H * qh / n_e
let anp_raw = anh / ane * qh;
// 确保物理一致性:
// 1. 质子密度不能超过总氢密度
// 2. 对于氢主导的气体, ane ≈ anp (电荷守恒)
// 3. 如果 anp > ahtot, 说明计算有问题, 使用 ane 作为估计
let anp = if anp_raw > ahtot || anp_raw < 0.0 {
// Newton-Raphson 可能未收敛, 使用电荷守恒估计
// 对于氢主导气体: ane ≈ anp (电子来自氢电离)
ane.min(ahtot).max(0.0)
} else {
anp_raw
};
let anhm = anh * ane * qmi;
let rhoter = params.wmy * ah * HMASS;
@@ -385,6 +477,11 @@ pub fn eldens_pure(params: &EldensParams, ipri: i32) -> EldensOutput {
let wm = rhoter / an / HMASS;
if params.id == 1 {
eprintln!("DEBUG ELDENS return: id={}, ane={:.6e}, anp={:.6e}, ahtot={:.6e}, anerel={:.6e}",
params.id, ane, anp, ahtot, anerel);
}
EldensOutput {
ane,
anp,
src/tlusty/math/hydrogen/mod.rs
+2 -1
View File
@@ -61,7 +61,8 @@ pub use sbfoh::*;
pub use sffhmi::*;
pub use sffhmi_add::*;
pub use sgmer::*;
pub use sgmer1::*;
// 重命名 sgmer1 模块的 sgmer1 函数避免与 sgmer 模块的 sgmer1 冲突
pub use sgmer1::sgmer1 as sgmer1_simple;
pub use sigave::*;
pub use sigk::*;
pub use sigmar::*;
src/tlusty/math/interpolation/yint.rs
+1 -1
View File
@@ -66,7 +66,7 @@ mod tests {
#[test]
fn test_yint_vs_lagran() {
// yint 和 lagran 应该给出相同的结果
use super::crate::tlusty::math::lagran;
use crate::tlusty::math::lagran;
let xl = [0.0, 1.0, 2.0];
let yl = [1.0, 3.0, 2.0];
src/tlusty/math/io/rdata.rs
+419
View File
@@ -418,6 +418,263 @@ pub struct LineTransition {
pub is_line: bool,
}
// ============================================================================
// 原子数据文件读取
// ============================================================================
use std::fs::File;
use std::io::{BufRead, BufReader};
use std::path::Path;
/// 读取原子数据文件（如 h1.dat, he1.dat 等）。
///
/// # 文件格式
///
/// ```text
/// ****** Levels
/// ENION G NQUANT TYPLEV IFWOP FRODF IMODL
/// ...
/// ****** Continuum transitions
/// II JJ MODE IFANCY ICOLIS IFRQ0 IFRQ1 OSC CPARAM
/// ...
/// ****** Line transitions
/// II JJ MODE IFANCY ICOLIS IFRQ0 IFRQ1 OSC CPARAM
/// LCOMP INTMOD NF XMAX TSTD
/// GAMAR STARK1 STARK2 STARK3 VDWH
/// ...
/// ```
///
/// # 参数
/// * `path` - 文件路径
/// * `nlevs` - 期望读取的能级数
///
/// # 返回值
/// (能级数据, 连续跃迁数据, 谱线跃迁数据)
pub fn read_ion_data_file<P: AsRef<Path>>(
path: P,
nlevs: usize,
) -> Result<(Vec<LevelInputData>, Vec<ContinuumInputData>, Vec<LineInputData>), String> {
let file = File::open(&path).map_err(|e| format!("无法打开文件 {:?}: {}", path.as_ref(), e))?;
let reader = BufReader::new(file);
let mut lines = reader.lines().peekable();
// 跳过第一个 ****** 行
let first_line = lines.next().ok_or("文件为空")?.map_err(|e| e.to_string())?;
if !first_line.contains('*') {
return Err("文件格式错误：第一行应该是 ****** Levels".to_string());
}
// 读取能级数据
let mut levels = Vec::with_capacity(nlevs);
for _ in 0..nlevs {
let line = lines.next().ok_or("能级数据不完整")?.map_err(|e| e.to_string())?;
let level = parse_level_line(&line)?;
levels.push(level);
}
// 跳到连续跃迁部分
skip_to_section_peekable(&mut lines, "Continuum")?;
let mut continua = Vec::new();
loop {
let line = match lines.next() {
Some(Ok(l)) => l,
_ => break,
};
// 检查是否到达谱线部分
if line.contains("******") {
if line.contains("Line") {
break;
}
continue;
}
if line.trim().is_empty() {
continue;
}
if let Ok(cont) = parse_continuum_line(&line) {
// 检查是否结束（ii >= nlevs 表示结束）
if cont.ii >= nlevs as i32 {
break;
}
continua.push(cont);
}
}
// 跳到谱线部分（如果还没到）
if !lines.peek().map_or(false, |r| r.as_ref().map_or(false, |l| l.contains("Line"))) {
skip_to_section_peekable(&mut lines, "Line")?;
}
let mut line_transitions = Vec::new();
loop {
let line = match lines.next() {
Some(Ok(l)) => l,
_ => break,
};
if line.contains("******") {
continue;
}
if line.trim().is_empty() {
continue;
}
// 跳过注释行（以 ! 开头）
let trimmed = line.trim();
if trimmed.starts_with('!') {
continue;
}
// 跳过续行（以 T 或 F 开头的是轮廓参数行）
if trimmed.starts_with('T') || trimmed.starts_with('F') {
continue;
}
// 跳过纯数字行（可能是额外参数）
if trimmed.chars().all(|c| c.is_numeric() || c.is_whitespace() || c == '.') {
// 可能是额外的数值行，跳过
if trimmed.split_whitespace().all(|s| s.parse::<f64>().is_ok()) {
continue;
}
}
if let Ok(line_data) = parse_line_transition(&line) {
// 检查是否结束
if line_data.ii >= nlevs as i32 {
break;
}
line_transitions.push(line_data);
}
}
Ok((levels, continua, line_transitions))
}
/// 跳到指定部分（peekable 版本）
fn skip_to_section_peekable<T: BufRead>(
lines: &mut std::iter::Peekable<std::io::Lines<T>>,
section_name: &str,
) -> Result<(), String> {
loop {
let line = lines.next().ok_or(format!("未找到 {} 部分", section_name))?.map_err(|e| e.to_string())?;
if line.contains("******") && line.contains(section_name) {
return Ok(());
}
}
}
/// 跳到指定部分（查找 ****** Section 标记）
fn skip_to_section<T: BufRead>(lines: &mut std::io::Lines<T>, section_name: &str) -> Result<(), String> {
loop {
let line = lines.next().ok_or(format!("未找到 {} 部分", section_name))?.map_err(|e| e.to_string())?;
if line.contains("******") && line.contains(section_name) {
return Ok(());
}
}
}
/// 解析能级行
fn parse_level_line(line: &str) -> Result<LevelInputData, String> {
let parts: Vec<&str> = line.split_whitespace().collect();
if parts.len() < 7 {
return Err(format!("能级行字段不足: {}", line));
}
let mut level = LevelInputData::default();
// ENION (能量) - 处理 Fortran D 指数格式
let enion_str = parts[0].to_uppercase().replace('D', "E");
level.enion = enion_str.parse().map_err(|_| format!("无法解析能量: {}", parts[0]))?;
// G (统计权重) - 同样处理 D 指数
let g_str = parts[1].to_uppercase().replace('D', "E");
level.g = g_str.parse().map_err(|_| format!("无法解析统计权重: {}", parts[1]))?;
// NQUANT (主量子数)
level.nquant = parts[2].parse().map_err(|_| format!("无法解析量子数: {}", parts[2]))?;
// TYPLEV (能级类型，字符串，可能带引号)
if parts.len() > 3 {
// 处理带引号的字符串
let typlev = parts[3].trim_matches('\'').to_string();
level.typlev = typlev;
}
// IFWOP
if parts.len() > 4 {
level.ifwop = parts[4].parse().unwrap_or(0);
}
// FRODF
if parts.len() > 5 {
let frodf_str = parts[5].to_uppercase().replace('D', "E");
level.frodf = frodf_str.parse().unwrap_or(0.0);
}
// IMODL
if parts.len() > 6 {
level.imodl = parts[6].parse().unwrap_or(0);
}
Ok(level)
}
/// 解析连续跃迁行
fn parse_continuum_line(line: &str) -> Result<ContinuumInputData, String> {
let parts: Vec<&str> = line.split_whitespace().collect();
if parts.len() < 9 {
return Err(format!("连续跃迁行字段不足: {}", line));
}
let mut cont = ContinuumInputData::default();
cont.ii = parts[0].parse().map_err(|_| format!("无法解析 ii: {}", parts[0]))?;
cont.jj = parts[1].parse().map_err(|_| format!("无法解析 jj: {}", parts[1]))?;
cont.mode = parts[2].parse().map_err(|_| format!("无法解析 mode: {}", parts[2]))?;
cont.ifancy = parts[3].parse().map_err(|_| format!("无法解析 ifancy: {}", parts[3]))?;
cont.icolis = parts[4].parse().map_err(|_| format!("无法解析 icolis: {}", parts[4]))?;
cont.ifrq0 = parts[5].parse().map_err(|_| format!("无法解析 ifrq0: {}", parts[5]))?;
cont.ifrq1 = parts[6].parse().map_err(|_| format!("无法解析 ifrq1: {}", parts[6]))?;
cont.osc = parts[7].parse().map_err(|_| format!("无法解析 osc: {}", parts[7]))?;
cont.cparam = parts[8].parse().map_err(|_| format!("无法解析 cparam: {}", parts[8]))?;
if parts.len() > 9 {
cont.ncol = parts[9].parse().unwrap_or(0);
}
Ok(cont)
}
/// 解析谱线跃迁行（简化版，只读取基本信息）
fn parse_line_transition(line: &str) -> Result<LineInputData, String> {
let parts: Vec<&str> = line.split_whitespace().collect();
if parts.len() < 9 {
return Err(format!("谱线跃迁行字段不足: {}", line));
}
let mut line_data = LineInputData::default();
line_data.ii = parts[0].parse().map_err(|_| format!("无法解析 ii: {}", parts[0]))?;
line_data.jj = parts[1].parse().map_err(|_| format!("无法解析 jj: {}", parts[1]))?;
line_data.mode = parts[2].parse().map_err(|_| format!("无法解析 mode: {}", parts[2]))?;
line_data.ifancy = parts[3].parse().map_err(|_| format!("无法解析 ifancy: {}", parts[3]))?;
line_data.icolis = parts[4].parse().map_err(|_| format!("无法解析 icolis: {}", parts[4]))?;
line_data.ifrq0 = parts[5].parse().map_err(|_| format!("无法解析 ifrq0: {}", parts[5]))?;
line_data.ifrq1 = parts[6].parse().map_err(|_| format!("无法解析 ifrq1: {}", parts[6]))?;
line_data.osc = parts[7].parse().map_err(|_| format!("无法解析 osc: {}", parts[7]))?;
line_data.cparam = parts[8].parse().map_err(|_| format!("无法解析 cparam: {}", parts[8]))?;
if parts.len() > 9 {
line_data.ncol = parts[9].parse().unwrap_or(0);
}
Ok(line_data)
}
// ============================================================================
// 纯计算函数
// ============================================================================
@@ -628,6 +885,168 @@ pub fn process_continua_pure(
(continua, itr, ic, nhod, lasv)
}
// ============================================================================
// 填充 AtomicData 结构体
// ============================================================================
/// 填充原子数据到 AtomicData 结构体。
///
/// 将从文件读取的能级、连续跃迁、谱线跃迁数据填充到全局原子数据结构中。
///
/// # 参数
/// * `atomic` - 原子数据结构体（可变引用）
/// * `ion_idx` - 离子索引（0-based）
/// * `nfirst` - 离子起始能级索引（1-based）
/// * `iat` - 原子序数
/// * `iz` - 电离态
/// * `ff_ion` - 电离势 (Ry)
/// * `charg2` - 电荷²
/// * `levels` - 能级输入数据
/// * `continua` - 连续跃迁输入数据
/// * `lines` - 谱线跃迁输入数据
/// * `teff` - 有效温度
pub fn populate_atomic_data(
atomic: &mut crate::tlusty::state::atomic::AtomicData,
ion_idx: usize,
nfirst: i32,
iat: i32,
iz: i32,
ff_ion: f64,
charg2: f64,
levels: &[LevelInputData],
continua: &[ContinuumInputData],
lines: &[LineInputData],
teff: f64,
) -> (i32, i32) {
let nlevs = levels.len() as i32;
// 填充离子参数
atomic.ionpar.ff[ion_idx] = ff_ion;
atomic.ionpar.charg2[ion_idx] = charg2;
atomic.ionpar.nfirst[ion_idx] = nfirst;
atomic.ionpar.nlast[ion_idx] = nfirst + nlevs - 1;
atomic.ionpar.nnext[ion_idx] = nfirst + nlevs;
atomic.ionpar.iz[ion_idx] = iat; // 原子序数 Z
// 填充离子数据索引
atomic.iondat.iati[ion_idx] = iat;
atomic.iondat.izi[ion_idx] = iz;
atomic.iondat.nlevs[ion_idx] = nlevs;
atomic.iondat.nllim[ion_idx] = nfirst + nlevs - 1;
// ZZ: 有效核电荷 = Z - iz + 1（对于类氢离子）
let zz = (iat - iz + 1) as f64;
// 填充能级参数
let mut ntrans = 0i32;
let mut ntranc = 0i32;
for (il, input) in levels.iter().enumerate() {
let level_idx = (nfirst as usize) + il - 1; // 转换为 0-based
// 能量转换
let e = input.enion.abs();
let e0 = convert_energy(e, zz, (il + 1) as i32);
let enion_value = if input.enion >= 0.0 { e0 } else { -e0 };
atomic.levpar.enion[level_idx] = enion_value;
atomic.levpar.g[level_idx] = if input.g == 0.0 {
2.0 * ((il + 1) as f64).powi(2)
} else {
input.g
};
atomic.levpar.nquant[level_idx] = if input.nquant == 0 {
(il + 1) as i32
} else {
input.nquant.abs()
};
atomic.levpar.iatm[level_idx] = iat;
atomic.levpar.iel[level_idx] = (ion_idx + 1) as i32; // 1-based ion index
atomic.levpar.indlev[level_idx] = (level_idx + 1) as i32;
// LTE 标志（负量子数表示 LTE）
if input.nquant < 0 {
atomic.levpar.iltlev[level_idx] = 1;
}
// 模型能级
atomic.levpar.imodl[level_idx] = input.imodl;
}
// 填充连续跃迁参数
for input in continua {
let itr = ntrans as usize;
if itr >= MTRANS {
break;
}
// 索引转换
let (ii, jj) = if input.jj < 1000 {
(input.ii + nfirst - 1, input.jj + nfirst - 1)
} else {
(input.ii + nfirst - 1, input.jj)
};
// 计算频率
let enion_ii = atomic.levpar.enion.get(ii as usize - 1).copied().unwrap_or(0.0);
let enion_jj = if input.jj < 1000 {
atomic.levpar.enion.get(jj as usize - 1).copied().unwrap_or(0.0)
} else {
0.0
};
let enion_nk = 0.0; // 简化：下一个离子的基态能级
let fr0 = (enion_ii - enion_jj + enion_nk) / H;
atomic.trapar.fr0[itr] = fr0;
atomic.trapar.osc0[itr] = input.osc;
atomic.trapar.cpar[itr] = input.cparam;
atomic.trapar.ilow[itr] = ii;
atomic.trapar.iup[itr] = jj;
atomic.trapar.icol[itr] = input.icolis;
atomic.trapar.ifc0[itr] = input.ifrq0;
atomic.trapar.ifc1[itr] = input.ifrq1;
// 标记为连续跃迁
atomic.trapar.itrcon[itr] = 1;
ntrans += 1;
ntranc += 1;
}
// 填充谱线跃迁参数
for input in lines {
let itr = ntrans as usize;
if itr >= MTRANS {
break;
}
let ii = input.ii + nfirst - 1;
let jj = input.jj + nfirst - 1;
// 计算频率
let enion_ii = atomic.levpar.enion.get(ii as usize - 1).copied().unwrap_or(0.0);
let enion_jj = atomic.levpar.enion.get(jj as usize - 1).copied().unwrap_or(0.0);
let fr0 = (enion_jj - enion_ii) / H;
atomic.trapar.fr0[itr] = fr0;
atomic.trapar.osc0[itr] = input.osc;
atomic.trapar.cpar[itr] = input.cparam;
atomic.trapar.ilow[itr] = ii;
atomic.trapar.iup[itr] = jj;
atomic.trapar.icol[itr] = input.icolis;
atomic.trapar.ifr0[itr] = input.ifrq0;
atomic.trapar.ifr1[itr] = input.ifrq1;
// 标记为谱线跃迁
atomic.trapar.itrcon[itr] = 0;
ntrans += 1;
}
(ntrans, ntranc)
}
// ============================================================================
// 测试
// ============================================================================
src/tlusty/math/opacity/meanopt.rs
+3 -3
View File
@@ -113,7 +113,7 @@ pub fn meanopt(
#[cfg(test)]
mod tests {
use super::*;
use super::super::opctab::{OpctabTableData, OpctabModelState};
use crate::tlusty::math::continuum::{OpctabTableData, OpctabModelState};
use approx::assert_relative_eq;
#[test]
@@ -150,7 +150,7 @@ mod tests {
let rho = 1e-7;
// 简化的 Planck 函数
let bnue: Vec<f64> = freq.iter().map(|f| 1.4743e-2 * f.powi(3)).collect();
let bnue: Vec<f64> = freq.iter().map(|&f| 1.4743e-2 * f.powi(3)).collect();
let w: Vec<f64> = vec![0.3, 0.4, 0.3]; // 权重
let params = MeanoptParams {
@@ -289,7 +289,7 @@ mod tests {
sige: 0.0,
};
let bnue: Vec<f64> = freq.iter().map(|f| 1.4743e-2 * f.powi(3)).collect();
let bnue: Vec<f64> = freq.iter().map(|&f| 1.4743e-2 * f.powi(3)).collect();
let w: Vec<f64> = vec![0.3, 0.4, 0.3];
let rho = 1e-7;
src/tlusty/math/temperature/rossop.rs
+255 -42
View File
@@ -17,7 +17,8 @@
use crate::tlusty::math::eint;
use crate::tlusty::state::constants::BOLK;
/// Hopf 函数多项式系数 (DATA A/.../)
/// Hopf 函数多项式系数 (from Fortran TLUSTY ROSSOP.f)
/// 注意：这些系数与 FORTRAN 代码中的 DATA A 数组完全一致
const HOPF_A: [f64; 5] = [
0.71044609,
-0.2830385,
@@ -60,31 +61,39 @@ const HOPF_A: [f64; 5] = [
/// END DO
/// ```
pub fn compute_hopf(taur: f64, hopf: f64) -> f64 {
// Fortran: X=HOPF; IF(X.GT.0.) GO TO 10
if hopf > 0.0 {
return hopf;
}
// 精确 Hopf 函数
if taur > 160.0 {
return HOPF_A[0];
}
// 特殊情况: taur 接近 0 时，使用渐近值
// 当 taur -> 0 时，Hopf 函数趋向于某个有限值
if taur < 1e-10 {
// 对于非常小的 taur，使用近似值
// Hopf(0) 的极限约为 0.86
return 0.860327569;
}
let ex = (-taur).exp();
let (e1, _, _) = eint(taur);
let mut e = e1;
// Fortran: X=A(1)
let mut x = HOPF_A[0];
// Fortran: IF(TAUR.GT.160.) GO TO 10
if taur > 160.0 {
return x;
}
// 对于极小的 taur，避免数值问题
// Fortran 没有显式处理 taur=0，但 expint(0) 理论上是无穷大
// 使用极小阈值来避免数值问题
if taur < 1e-100 {
// Hopf(0) 的极限值
return 0.71044609;
}
// Fortran: EX=EXP(-TAUR)
let ex = (-taur).exp();
// Fortran: E1=EXPINT(TAUR)
let (e1, _, _) = eint(taur);
// Fortran: E=E1
let mut e = e1;
// Fortran: DO I=1,4
// E=(EX-TAUR*E)/I
// X=X+E*A(I+1)
// END DO
for i in 1..=4_usize {
// Fortran: E=(EX-TAUR*E)/I
e = (ex - taur * e) / (i as f64);
x = x + e * HOPF_A[i];
}
@@ -147,6 +156,8 @@ pub struct RossopConfig {
pub iter: i32,
/// 频率点数
pub nfreq: usize,
/// Rayleigh 散射标志
pub ifrayl: i32,
}
impl Default for RossopConfig {
@@ -155,6 +166,7 @@ impl Default for RossopConfig {
ioptab: 0,
iter: 1,
nfreq: 1,
ifrayl: 0,
}
}
}
@@ -202,7 +214,7 @@ pub struct RossopOutput {
pub hopf_value: f64,
/// 总粒子数密度 (cm^-3)
pub an: f64,
/// 电子密度 (cm^-3) - 需要调用 ELDENS 计算
/// 电子密度 (cm^-3)
pub ane: f64,
/// 密度 (g/cm³)
pub rho: f64,
@@ -210,6 +222,62 @@ pub struct RossopOutput {
pub abross: f64,
}
/// ROSSOP 回调接口（用于调用 ELDENS, WNSTOR, STEQEQ, OPACF0, MEANOP）
pub trait RossopCallbacks {
/// 调用 ELDENS 计算电子密度
/// 返回 (ane, energ, entt, wm)
fn call_eldens(&mut self, id: usize, t: f64, an: f64) -> (f64, f64, f64, f64);
/// 调用 WNSTOR 存储 LTE 布居数
fn call_wnstor(&mut self, id: usize);
/// 调用 STEQEQ 计算统计平衡方程
fn call_steqeq(&mut self, id: usize);
/// 调用 OPACF0 计算不透明度
fn call_opacf0(&mut self, id: usize, nfreq: usize);
/// 调用 MEANOP 计算平均不透明度
/// 返回 (opros, oppla)
fn call_meanop(&mut self, t: f64) -> (f64, f64);
/// 调用 MEANOPT 计算简化平均不透明度
/// 返回 (opros, oppla)
fn call_meanopt(&mut self, t: f64, id: usize, rho: f64) -> (f64, f64);
/// 获取 WMM(ID)
fn get_wmm(&self, id: usize) -> f64;
}
/// 空回调实现（默认不做任何操作）
#[derive(Debug, Clone, Default)]
pub struct NoOpRossopCallbacks;
impl RossopCallbacks for NoOpRossopCallbacks {
fn call_eldens(&mut self, _id: usize, _t: f64, an: f64) -> (f64, f64, f64, f64) {
// 默认返回简化估计：电子密度约为总粒子数的一半（高温完全电离）
(an * 0.5, 0.0, 0.0, 1.0)
}
fn call_wnstor(&mut self, _id: usize) {}
fn call_steqeq(&mut self, _id: usize) {}
fn call_opacf0(&mut self, _id: usize, _nfreq: usize) {}
fn call_meanop(&mut self, _t: f64) -> (f64, f64) {
(0.0, 0.0)
}
fn call_meanopt(&mut self, _t: f64, _id: usize, _rho: f64) -> (f64, f64) {
(0.34, 0.34) // 默认使用典型恒星大气不透明度
}
fn get_wmm(&self, _id: usize) -> f64 {
2.3e-24 // 默认平均分子量（约等于氢）
}
}
/// 简化版 ROSSOP 函数。
///
/// 计算温度和基本量，不包含复杂的依赖函数调用。
@@ -228,43 +296,187 @@ pub fn rossop(
config: &RossopConfig,
params: &RossopParams,
state: &mut RossopModelState,
) -> RossopOutput {
rossop_with_callbacks(config, params, state, &mut NoOpRossopCallbacks)
}
/// 完整版 ROSSOP 函数（带回调）。
///
/// 严格按照 Fortran rossop.f 的逻辑流程实现：
///
/// ```fortran
/// X=HOPF
/// IF(X.GT.0.) GO TO 10
/// X=A(1)
/// IF(TAUR.GT.160.) GO TO 10
/// EX=EXP(-TAUR)
/// E1=EXPINT(TAUR)
/// E=E1
/// DO I=1,4
/// E=(EX-TAUR*E)/I
/// X=X+E*A(I+1)
/// END DO
/// 10 T=(0.75*T4*(TAUR+X)+EXTOT)**0.25
/// if(ioptab.ge.-1) then
/// AN=P/T/BOLK
/// CALL ELDENS(ID,T,AN,ANE,ENRG,ENTT,WM,1)
/// RHO=WMM(ID)*(AN-ANE)
/// DENS(ID)=RHO
/// TEMP(ID)=T
/// ELEC(ID)=ANE
/// if(ioptab.ge.0) then
/// CALL WNSTOR(ID)
/// CALL STEQEQ(ID,POP,1)
/// CALL OPACF0(ID,NFREQ)
/// CALL MEANOP(T,ABSO,SCAT,OPROS,OPPLA)
/// ABROSS=OPROS/RHO
/// ABROSD(ID)=ABROSS
/// ABPLAD(ID)=OPPLA/RHO
/// else
/// call meanopt(t,id,rho,opros,oppla)
/// abrosd(id)=opros
/// abplad(id)=oppla
/// abross=opros
/// end if
/// else
/// temp(id)=t
/// rho=rhoeos(t,p)
/// dens(id)=rho
/// call meanopt(t,id,rho,opros,oppla)
/// abrosd(id)=opros
/// abplad(id)=oppla
/// abross=opros
/// end if
/// ```
pub fn rossop_with_callbacks<C: RossopCallbacks>(
config: &RossopConfig,
params: &RossopParams,
state: &mut RossopModelState,
callbacks: &mut C,
) -> RossopOutput {
let id = params.id;
// ========================================================================
// Step 1: 计算 Hopf 函数
// 对应 Fortran: X=HOPF; IF(X.GT.0.) GO TO 10; ...
// ========================================================================
let hopf_value = compute_hopf(params.taur, params.hopf);
// ========================================================================
// Step 2: 计算温度
// T = (0.75 * T4 * (TAUR + X) + EXTOT)^0.25
// 对应 Fortran: 10 T=(0.75*T4*(TAUR+X)+EXTOT)**0.25
// ========================================================================
let t = compute_temperature(params.taur, params.hopf, params.t4, params.extot);
// Step 3: 计算总粒子数密度
// AN = P / T / BOLK
let an = compute_total_density(params.p, t);
// ========================================================================
// Step 3: 根据 ioptab 选择计算路径
// 对应 Fortran: if(ioptab.ge.-1) then
// ========================================================================
let (ane, rho, abross) = if config.ioptab >= -1 {
// ----------------------------------------------------------------
// 路径 1: ioptab >= -1
// 对应 Fortran lines 54-86
// ----------------------------------------------------------------
// 更新温度数组
state.temp[id] = t;
// Step 3.1: 计算总粒子数密度
// 对应 Fortran: AN=P/T/BOLK
let an = compute_total_density(params.p, t);
// 计算密度（需要电子密度，这里先用 0 占位）
// 完整实现需要调用 ELDENS
let rho = params.wmm[id] * an;
state.dens[id] = rho;
// Step 3.2: 调用 ELDENS 计算电子密度
// 对应 Fortran: CALL ELDENS(ID,T,AN,ANE,ENRG,ENTT,WM,1)
let (ane, _energ, _entt, _wm) = callbacks.call_eldens(id, t, an);
// 简化版：电子密度和不透明度返回 0，需要后续计算
let ane = 0.0;
let abross = 0.0;
// Step 3.3: 计算质量密度
// 对应 Fortran: RHO=WMM(ID)*(AN-ANE)
let wmm = callbacks.get_wmm(id);
let rho = wmm * (an - ane);
// 根据配置更新不透明度数组
if config.ioptab >= -1 {
// 简化处理：不透明度设为 0
state.abrosd[id] = 0.0;
state.abplad[id] = 0.0;
}
// 更新状态数组
// 对应 Fortran: DENS(ID)=RHO; TEMP(ID)=T; ELEC(ID)=ANE
state.dens[id] = rho;
state.temp[id] = t;
state.elec[id] = ane;
// Step 3.4: 根据 ioptab 选择不透明度计算方式
if config.ioptab >= 0 {
// ----------------------------------------------------------------
// 路径 1a: ioptab >= 0（完整计算）
// 对应 Fortran lines 68-80
// ----------------------------------------------------------------
// 调用 WNSTOR
// 对应 Fortran: CALL WNSTOR(ID)
callbacks.call_wnstor(id);
// 调用 STEQEQ
// 对应 Fortran: CALL STEQEQ(ID,POP,1)
callbacks.call_steqeq(id);
// 调用 OPACF0
// 对应 Fortran: CALL OPACF0(ID,NFREQ)
callbacks.call_opacf0(id, config.nfreq);
// 调用 MEANOP
// 对应 Fortran: CALL MEANOP(T,ABSO,SCAT,OPROS,OPPLA)
let (opros, oppla) = callbacks.call_meanop(t);
// 计算最终不透明度
// 对应 Fortran: ABROSS=OPROS/RHO; ABROSD(ID)=ABROSS; ABPLAD(ID)=OPPLA/RHO
let abross = if rho > 0.0 { opros / rho } else { opros };
state.abrosd[id] = abross;
state.abplad[id] = if rho > 0.0 { oppla / rho } else { oppla };
(ane, rho, abross)
} else {
// ----------------------------------------------------------------
// 路径 1b: ioptab == -1（简化计算）
// 对应 Fortran lines 81-86
// ----------------------------------------------------------------
// 调用 MEANOPT
// 对应 Fortran: call meanopt(t,id,rho,opros,oppla)
let (opros, oppla) = callbacks.call_meanopt(t, id, rho);
// 更新不透明度数组
// 对应 Fortran: abrosd(id)=opros; abplad(id)=oppla; abross=opros
state.abrosd[id] = opros;
state.abplad[id] = oppla;
(ane, rho, opros)
}
} else {
// ----------------------------------------------------------------
// 路径 2: ioptab < -1（使用状态方程）
// 对应 Fortran lines 87-95
// ----------------------------------------------------------------
// 更新温度
// 对应 Fortran: temp(id)=t
state.temp[id] = t;
// 使用简化的状态方程（这里用 ELDENS 代替）
// 对应 Fortran: rho=rhoeos(t,p); dens(id)=rho
let an = compute_total_density(params.p, t);
let (ane, _, _, _) = callbacks.call_eldens(id, t, an);
let wmm = callbacks.get_wmm(id);
let rho = wmm * (an - ane);
state.dens[id] = rho;
// 调用 MEANOPT
// 对应 Fortran: call meanopt(t,id,rho,opros,oppla)
let (opros, oppla) = callbacks.call_meanopt(t, id, rho);
// 更新不透明度数组
state.abrosd[id] = opros;
state.abplad[id] = oppla;
(ane, rho, opros)
};
RossopOutput {
t,
hopf_value,
an,
an: compute_total_density(params.p, t),
ane,
rho,
abross,
@@ -309,10 +521,11 @@ mod tests {
#[test]
fn test_compute_hopf_zero_taur() {
// taur = 0 时的计算（渐近值）
// taur = 0 时的计算（使用极限值）
// Fortran 没有显式处理 taur=0，但我们使用 HOPF_A[0] 作为极限值
let result = compute_hopf(0.0, 0.0);
// 应该接近 Hopf(0) 的极限值约 0.86
assert!((result - 0.860327569).abs() < 1e-6);
// 应该返回 HOPF_A[0] = 0.71044609
assert!((result - HOPF_A[0]).abs() < 1e-10);
}
#[test]
src/tlusty/mod.rs
+5 -4
View File
@@ -6,11 +6,12 @@ pub mod data;
pub mod io;
pub mod math;
pub mod physics;
pub mod runner;
pub mod main;
pub mod state;
pub use runner::{
pub use main::{
run_tlusty, check_convergence, select_solver,
TlustyConfig, TlustyResult, TlustyState,
SolverType, TimingStats, MSMX,
TlustyResult, TlustyState, TlustyWorkArrays, ScratchFiles,
SolverType, MSMX,
};
pub use state::config::TlustyConfig;
src/tlusty/runner.rs
-305
View File
@@ -1,305 +0,0 @@
//! TLUSTY 主程序入口。
//!
//! 重构自 TLUSTY `tlusty.f` 主程序。
//!
//! # 算法概述
//!
//! TLUSTY 使用混合 Complete Linearization (CL) 和 Accelerated Lambda Iteration (ALI) 方法
//! 计算非LTE恒星大气模型。
//!
//! 主迭代循环：
//! 1. 初始化 (START)
//! 2. 形式解 (RESOLV)
//! 3. 收敛加速 (ACCEL2)
//! 4. 解线性化方程 (SOLVE/SOLVES/RYBSOL)
//! 5. 重复直到收敛
use std::time::Instant;
// ============================================================================
// 运行配置
// ============================================================================
/// TLUSTY 运行配置。
#[derive(Debug, Clone)]
pub struct TlustyConfig {
/// 最大迭代次数
pub max_iter: usize,
/// 收敛加速模式 (0=关闭, >0=开启)
pub accel_mode: i32,
/// 解法选择 (0=标准SOLVE, 1=Ryan)
pub solver_mode: i32,
/// 是否输出诊断信息
pub verbose: bool,
}
impl Default for TlustyConfig {
fn default() -> Self {
Self {
max_iter: 200,
accel_mode: 1,
solver_mode: 0,
verbose: false,
}
}
}
// ============================================================================
// 运行状态
// ============================================================================
/// TLUSTY 运行状态。
#[derive(Debug, Clone)]
pub struct TlustyState {
/// 当前迭代次数
pub iter: usize,
/// 是否初始化阶段
pub is_init: bool,
/// 是否完成
pub is_finished: bool,
/// 是否最后一次迭代
pub is_final: bool,
/// 系统维度 NN
pub nn: usize,
}
impl Default for TlustyState {
fn default() -> Self {
Self {
iter: 0,
is_init: true,
is_finished: false,
is_final: false,
nn: 0,
}
}
}
// ============================================================================
// 运行结果
// ============================================================================
/// TLUSTY 运行结果。
#[derive(Debug, Clone)]
pub struct TlustyResult {
/// 最终迭代次数
pub total_iterations: usize,
/// 是否收敛
pub converged: bool,
/// 总运行时间（秒）
pub total_time_secs: f64,
/// 各阶段时间统计
pub timing_stats: TimingStats,
}
/// 时间统计。
#[derive(Debug, Clone, Default)]
pub struct TimingStats {
pub formal_solution_time: f64,
pub matrix_solution_time: f64,
pub acceleration_time: f64,
}
// ============================================================================
// 主运行函数
// ============================================================================
/// 运行 TLUSTY 计算。
///
/// # 算法流程
///
/// ```text
/// 1. 初始化
/// - 读取输入参数 (fort.5)
/// - 设置初始温度结构
/// - 读取原子数据
/// - 设置频率网格
///
/// 2. 主迭代循环
/// a) 形式解 (RESOLV)
/// - 解辐射转移方程
/// - 计算辐射场
/// b) 收敛加速 (ACCEL2)
/// - Ng 加速
/// c) 解线性化方程 (SOLVE/SOLVES/RYBSOL)
/// - 统计平衡方程
/// - 能量守恒方程
/// d) 更新大气结构
///
/// 3. 输出最终模型
/// ```
///
/// # Fortran 原始代码
///
/// ```fortran
/// PROGRAM TLUSTY
/// ...
/// CALL START
/// LFIN=.FALSE.
/// IF(NITER.EQ.0) LFIN=.TRUE.
///
/// 10 ITER=ITER+1
/// CALL RESOLV
/// INIT=0
/// IF(LFIN) GO TO 20
///
/// IF(IACC.GT.0) CALL ACCEL2
///
/// IF(IFRYB.EQ.0) THEN
/// IF(NN.GT.MSMX) THEN
/// CALL SOLVE
/// ELSE
/// CALL SOLVES
/// END IF
/// ELSE
/// CALL RYBSOL
/// END IF
///
/// CALL TIMING(2,ITER)
/// GO TO 10
/// 20 CONTINUE
/// STOP
/// END
/// ```
pub fn run_tlusty(config: &TlustyConfig) -> TlustyResult {
let mut state = TlustyState::default();
let mut timing_stats = TimingStats::default();
let start_time = Instant::now();
// ========================================
// 1. 初始化阶段
// ========================================
// 对应 Fortran:
// INIT=1
// ITER=0
// CALL START
// LFIN=.FALSE.
// IF(NITER.EQ.0) LFIN=.TRUE.
// 执行初始化
// let start_output = start(&mut start_params);
// state.nn = start_output.nn;
// 检查是否需要迭代 (NITER == 0 表示只做初始模型)
// if iterat.niter == 0 {
// state.is_final = true;
// state.is_finished = true;
// }
// ========================================
// 2. 主迭代循环
// ========================================
while !state.is_finished && state.iter < config.max_iter {
state.iter += 1;
let iter_start = Instant::now();
// 2.1 形式解 (RESOLV)
// 对应 Fortran: CALL RESOLV
// let resolv_output = resolv(&mut resolv_params);
state.is_init = false;
if state.is_final {
break;
}
// 2.2 收敛加速 (ACCEL2)
// 对应 Fortran: IF(IACC.GT.0) CALL ACCEL2
if config.accel_mode > 0 {
let accel_start = Instant::now();
// let accel_output = accel2(&mut accel_params);
timing_stats.acceleration_time += accel_start.elapsed().as_secs_f64();
}
// 2.3 解线性化方程
// 对应 Fortran:
// IF(IFRYB.EQ.0) THEN
// IF(NN.GT.MSMX) THEN
// CALL SOLVE
// ELSE
// CALL SOLVES
// END IF
// ELSE
// CALL RYBSOL
// END IF
let solve_start = Instant::now();
let solver_type = select_solver(state.nn, MSMX);
match solver_type {
SolverType::Standard => {
// let solve_output = solve(&solve_params);
}
SolverType::Simple => {
// let solves_output = solves(&solves_params);
}
SolverType::Ryan => {
// let rybsol_output = rybsol(&rybsol_params);
}
}
timing_stats.matrix_solution_time += solve_start.elapsed().as_secs_f64();
// 记录时间 (TIMING)
// 对应 Fortran: CALL TIMING(2,ITER)
// let timing_output = timing(&TimingParams { mode: TimingMode::Iteration, iter: state.iter });
if config.verbose {
println!(
"Iteration {}: time = {:.3}s",
state.iter,
iter_start.elapsed().as_secs_f64()
);
}
// 检查收敛 (在实际实现中检查 CHMX 和 CHMT)
}
let total_time = start_time.elapsed().as_secs_f64();
TlustyResult {
total_iterations: state.iter,
converged: state.is_finished,
total_time_secs: total_time,
timing_stats,
}
}
// ============================================================================
// 辅助函数
// ============================================================================
/// 检查收敛性。
pub fn check_convergence(iter: usize, max_change: f64, tolerance: f64) -> bool {
iter > 0 && max_change < tolerance
}
/// 选择解法。
///
/// 根据 nn 和 MSMX 选择合适的解法：
/// - nn > MSMX: 使用完整矩阵解法 SOLVE
/// - nn <= MSMX: 使用简化解法 SOLVES
pub fn select_solver(nn: usize, msmx: usize) -> SolverType {
if nn > msmx {
SolverType::Standard
} else {
SolverType::Simple
}
}
/// 解法类型。
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum SolverType {
/// 标准解法 (SOLVE)
Standard,
/// 简化解法 (SOLVES)
Simple,
/// Ryan 解法 (RYBSOL)
Ryan,
}
// ============================================================================
// 常量 (从 Fortran 移植)
// ============================================================================
/// 最大简化矩阵维度
pub const MSMX: usize = 2000;
src/tlusty/state/config.rs
+99
View File
@@ -579,6 +579,105 @@ pub struct TlustyConfig {
pub mlcons: MlCons,
pub taursl: TaurSl,
pub iprkey: IprKey,
// 迭代控制参数（从 iterat.rs 移入）
pub accel: AccelConfig,
pub acclp: AcclpConfig,
pub lambda: LambdaConfig,
pub chnad: ChnadConfig,
}
/// 加速收敛参数（从 ITERAT.FOR 中的 ACCEL COMMON 块）
#[derive(Debug, Clone, Default)]
pub struct AccelConfig {
/// 过松弛因子
pub orelax: f64,
/// 加速标志
pub iacc: i32,
/// 初始加速
pub iacc0: i32,
/// 加速方向
pub iacd: i32,
/// Kantorovich 向量
pub kant: Vec<i32>,
/// 奇异点位置
pub lsng: Vec<i32>,
/// Aitken 步
pub laso: i32,
/// 重启动标志
pub lres2: i32,
/// lac2 标志
pub lac2: i32,
}
impl AccelConfig {
pub fn new() -> Self {
Self {
orelax: 1.0,
iacc: 3,
iacc0: 0,
iacd: 0,
kant: vec![0; 200],
lsng: vec![0; 1000],
laso: 0,
lres2: 0,
lac2: 0,
}
}
}
/// 加速 Lambda 迭代参数（从 ITERAT.FOR 中的 ACCLP COMMON 块）
#[derive(Debug, Clone, Default)]
pub struct AcclpConfig {
pub ilam: i32,
pub iacpp: i32,
pub iacc0p: i32,
pub iacdp: i32,
pub lac2p: i32,
}
impl AcclpConfig {
pub fn new() -> Self {
Self {
ilam: 0,
iacpp: 3,
iacc0p: 0,
iacdp: 0,
lac2p: 0,
}
}
}
/// Lambda 迭代参数（从 ITERAT.FOR 中的 LAMBDA COMMON 块）
#[derive(Debug, Clone, Default)]
pub struct LambdaConfig {
/// Lambda 点数
pub nlambd: i32,
/// 能级修正标志
pub ielcor: i32,
}
impl LambdaConfig {
pub fn new() -> Self {
Self {
nlambd: 0,
ielcor: 100,
}
}
}
/// 变化附加参数（从 ITERAT.FOR 中的 CHNAD COMMON 块）
#[derive(Debug, Clone, Default)]
pub struct ChnadConfig {
/// 最大温度变化
pub chmaxt: f64,
}
impl ChnadConfig {
pub fn new() -> Self {
Self {
chmaxt: 0.1,
}
}
}
impl TlustyConfig {