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

修复7

Browse Source
This commit is contained in:
fmq
2026-04-04 23:01:19 +08:00
parent 24b2d17003
commit d62beb8ad3
21 changed files with 1719 additions and 1551 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/bin/tlusty.rs
+13 -913
View File
@@ -1,927 +1,27 @@
//! TLUSTY 可执行程序入口
//! 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 std::io::{self, BufReader};
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,
};
// f2r_depends: ACCEL2, RESOLV, RYBSOL, SOLVE, SOLVES, START, TIMING
use tlusty_rust::tlusty::{run_tlusty, TlustyConfig};
use tlusty_rust::tlusty::io::{FortranReader, FortranWriter};
fn main() -> anyhow::Result<()> {
let args: Vec<String> = env::args().collect();
let mut config = TlustyConfig::default();
let mut input_reader = FortranReader::new(BufReader::new(io::stdin()));
let mut output_writer = FortranWriter::new(io::stdout());
// 解析命令行参数
let (input_path, _output_path) = parse_args(&args)?;
let result = run_tlusty(&mut config, &mut input_reader, &mut output_writer);
// 读取输入文件
let input_params = if let Some(ref path) = input_path {
println!("Reading input from: {}", path.display());
read_input_file(path)?
if result.converged {
eprintln!("Converged after {} iterations ({:.2}s)",
result.total_iterations, result.total_time_secs);
} 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()));
}
}
eprintln!("Did NOT converge after {} iterations ({:.2}s)",
result.total_iterations, result.total_time_secs);
}
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,
ifoppf: 0,
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,
anato_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/resolv.rs
+387 -68
View File
@@ -30,7 +30,7 @@
//! - fort.6: 标准输出（进度和诊断信息）
use super::FortranWriter;
use crate::tlusty::state::constants::{MDEPTH, MFREQ, MLEVEL};
use crate::tlusty::state::constants::{MDEPTH, MFREQ, MLEVEL, MTRANS, H, HK, BOLK, HMASS, SIGE, SIG4P, BN, UN, HALF, PI};
use crate::tlusty::math::{
rayset, prd, opaini, rates1_pure, ratsp1, steqeq_pure, newpop,
elcor_pure, accelp, rosstd_evaluate, output, pzert,
@@ -38,7 +38,14 @@ use crate::tlusty::math::{
alisk2_pure, alist1_pure, alist2, pzevld, hesol6, dmeval,
rybheq, princ_pure, coolrt_pure, rechck_pure, rteint, rtecmu,
taufr1, linsel_pure, rtecf1, opacf1, rtefr1, rtecom,
rtesol,
eldens_pure, EldensParams, EldensConfig,
compute_opacity_at_frequency, generate_lte_frequency_grid,
LteOpacityParams,
lucy_pure, LucyConfig, LucyModelParams,
OpacflPointData, Rad1PointData,
};
use crate::tlusty::math::io::{OutputParams};
use crate::tlusty::state::config::TlustyConfig;
use crate::tlusty::state::atomic::AtomicData;
use crate::tlusty::state::model::ModelState;
@@ -183,7 +190,7 @@ impl Default for ResolvConfig {
// ============================================================================
/// RESOLV 输入参数。
pub struct ResolvParams<'a> {
pub struct ResolvParams<'a, W7: std::io::Write = std::fs::File> {
/// 配置参数
pub config: ResolvConfig,
/// TLUSTY 配置（可变）
@@ -192,6 +199,8 @@ pub struct ResolvParams<'a> {
pub atomic: &'a mut AtomicData,
/// 模型状态（可变）
pub model: &'a mut ModelState,
/// fort.7 输出写入器（用于 OUTPUT 调用写模型）
pub writer7: Option<&'a mut FortranWriter<W7>>,
}
/// RESOLV 输出。
@@ -230,8 +239,8 @@ fn nitlam(iter: i32) -> i32 {
///
/// # 返回值
/// 计算结果
pub fn resolv<W: std::io::Write>(
params: &mut ResolvParams,
pub fn resolv<W: std::io::Write, W7: std::io::Write>(
params: &mut ResolvParams<W7>,
mut writer: Option<&mut FortranWriter<W>>,
) -> ResolvOutput {
// 标记已导入的函数（用于 f2r_check 脚本检测）
@@ -326,67 +335,350 @@ pub fn resolv<W: std::io::Write>(
}
// -----------------------------------------------------------
// Part 4: Lambda 迭代循环
// Part 4: Lambda 迭代循环 — 真实物理计算
// -----------------------------------------------------------
let nd = config.nd;
let nfreq = config.nfreq;
let teff = config.teff;
let grav = params.tlusty_config.inppar.grav;
// 物理常数
let c_light = 2.99792458e10_f64;
let h_over_c2 = 2.0 * H / (c_light * c_light);
let sqrt3_inv = 1.0 / 3.0_f64.sqrt();
let t4_eff = SIG4P * teff.powi(4);
let dprad = 1.891204931e-15 * teff.powi(4);
let prd0 = dprad / 1.732;
// 生成频率网格
let freq_grid = generate_lte_frequency_grid(teff, nfreq);
let freq = &freq_grid.freq;
let weights = &freq_grid.weights;
// 3 点 Gauss-Legendre 积分节点 (在 [0,1] 上)
let gl_mu: [f64; 3] = [
0.5 * (1.0 - (3.0_f64 / 5.0).sqrt()),
0.5,
0.5 * (1.0 + (3.0_f64 / 5.0).sqrt()),
];
let gl_w: [f64; 3] = [5.0 / 18.0, 8.0 / 18.0, 5.0 / 18.0];
// Eddington H 函数近似 (各向同性)
let fh = vec![sqrt3_inv; nfreq];
let hextrd = vec![0.0; nfreq];
// 深度间隔 deldm
let mut deldm = vec![0.0; nd];
for id in 1..nd {
deldm[id - 1] = params.model.modpar.dm[id] - params.model.modpar.dm[id - 1];
}
for _ilam_iter in 1..=nlambd {
ilam = _ilam_iter;
debug_log!("RESOLV: Lambda iteration {} of {}", ilam, nlambd);
// OPAINI(1) - 初始化不透明度
// opaini(&OpainiParams { ... });
debug_log!("RESOLV: OPAINI(1) called");
// ==============================================================
// Step 1: 在每个深度点计算 LTE 电子密度
// ==============================================================
let eldens_config = EldensConfig {
ifmol: 0,
tmolim: 1e10,
ioptab: -1,
iath: 1,
iatref: 1,
ihm: 0,
ih2: 0,
ih2p: 0,
pfhyd: 2.0,
};
// 康普顿散射
if config.icompt != 0 && ilam > 1 {
// RTECOM
debug_log!("RESOLV: RTECOM called (Compton, ilam > 1)");
let mut ne_arr = vec![0.0; nd];
let mut nh_arr = vec![0.0; nd];
let mut np_arr = vec![0.0; nd];
let mut wmm_arr = vec![1.0; nd];
for id in 0..nd {
let t = params.model.modpar.temp[id];
let dens = params.model.modpar.dens[id];
let an = dens / HMASS + params.model.modpar.elec[id];
let eldens_params = EldensParams {
id: id + 1,
t,
an,
ytot: 1.1,
qref: 0.0,
dqnr: 0.0,
wmy: 1.0,
config: eldens_config.clone(),
state_params: None,
molecule_data: None,
anato_data: None,
};
let eldens_output = eldens_pure(&eldens_params, 0);
ne_arr[id] = eldens_output.ane;
np_arr[id] = eldens_output.anp;
nh_arr[id] = eldens_output.ahtot;
wmm_arr[id] = eldens_output.wm;
}
// 计算辐射跃迁速率
if config.ifprec == 0 {
// RATES1(0)
// rates1_pure(&mut Rates1Params { ... });
debug_log!("RESOLV: RATES1(0) called");
} else {
// RATSP1
// ratsp1(...);
debug_log!("RESOLV: RATSP1 called");
// ==============================================================
// Step 2: 在每个 (depth, frequency) 计算不透明度
// ==============================================================
// chi[id][ij] = 总不透明度 (cm²/g)
// ab_true[id][ij] = 真吸收 (不含散射)
let mut chi = vec![vec![0.0; nfreq]; nd];
let mut ab_true = vec![vec![0.0; nfreq]; nd];
// 用于 Lucy 的 per-frequency 数据
let mut opacfl_data: Vec<OpacflPointData> = Vec::with_capacity(nfreq);
for id in 0..nd {
let t = params.model.modpar.temp[id];
let dens = params.model.modpar.dens[id];
let dens_safe = dens.max(1e-30);
let ne = ne_arr[id];
let nh_total = nh_arr[id];
let anp = np_arr[id];
let nh_neutral = (nh_total - anp).max(0.0);
let hkt = HK / t;
let sgff = 3.694e8 / t.sqrt() * ne;
let lte_params = LteOpacityParams {
t, ne, nh_total,
np: anp, nh_neutral, nhm: 0.0,
rho: dens,
uh: 2.0, uhe: 1.0, uhep: 2.0,
xh: 0.70, xhe: 0.28,
};
let scat_per_gram = SIGE * ne / dens_safe;
for ij in 0..nfreq {
let fr = freq[ij];
let (ab, sct) = compute_opacity_at_frequency(
fr, t, ne, nh_neutral, anp, 0.0, hkt, sgff, &lte_params,
);
let ab_gram = ab / dens_safe;
let sct_gram = sct / dens_safe;
chi[id][ij] = ab_gram + sct_gram;
ab_true[id][ij] = ab_gram;
}
}
// PRD
// prd(0, ...);
debug_log!("RESOLV: PRD(0) called");
// ==============================================================
// Step 3: 正式解 — 多角度 Gauss-Legendre 积分计算 Jν
// ==============================================================
// Jν(id) = (1/2) ∫₀¹ I(μ) dμ ≈ Σ_k w_k * 0.5 * (I_in + I_out)
let mut jnu = vec![vec![0.0; nfreq]; nd];
// 更新占据数
debug_log!("RESOLV: Updating populations for {} depth points", config.nd);
for id in 0..config.nd {
// STEQEQ(ID, POP, 1)
// steqeq_pure(&SteqeqParams { ... }, 1);
// 同时收集辐射数据用于 Lucy
let mut rad1_data: Vec<Rad1PointData> = Vec::with_capacity(nfreq);
// NEWPOP(ID, POP)
// newpop(&mut NewpopParams { ... });
for ij in 0..nfreq {
let fr = freq[ij];
// ELCOR（电子修正）
if !config.lchc && iter < config.ielcor {
// elcor_pure(&ElcorParams { ... });
// 计算频率光学深度增量
let mut dtau_freq = vec![0.0; nd - 1];
for id in 0..(nd - 1) {
let dm_half = params.model.modpar.dm[id + 1] - params.model.modpar.dm[id];
let chi_avg = 0.5 * (chi[id][ij] + chi[id + 1][ij]);
dtau_freq[id] = chi_avg * dm_half;
}
// 源函数 S = Bν(T) 在每个深度
let mut source = vec![0.0; nd];
for id in 0..nd {
let t = params.model.modpar.temp[id];
let hkt = HK / t;
let x = (hkt * fr).min(150.0);
source[id] = h_over_c2 * fr.powi(3) / (x.exp() - 1.0);
}
// 多角度积分
for id in 0..nd {
jnu[id][ij] = 0.0;
}
// 也收集 Eddington 因子 K/J
let mut rad1_ij = vec![0.0; nd];
let mut fak1_ij = vec![UN / 3.0; nd]; // 默认 Eddington 因子
for k in 0..3 {
let mu_k = gl_mu[k];
let w_k = gl_w[k];
// 角度光学深度: dtau_μ = dtau / μ
let mut dtau_mu = vec![0.0; nd - 1];
for id in 0..(nd - 1) {
dtau_mu[id] = dtau_freq[id] / mu_k;
}
// 入射强度（向下，无外部辐射）
let mut ri_in = vec![0.0; nd];
let mut ali_in = vec![0.0; nd];
rtesol(&dtau_mu, &source, 0.0, 0.0, -mu_k,
&mut ri_in, &mut ali_in, nd);
// 出射强度（向上）
let rdown_bottom = source[nd - 1];
let mut ri_out = vec![0.0; nd];
let mut ali_out = vec![0.0; nd];
rtesol(&dtau_mu, &source, 0.0, rdown_bottom, mu_k,
&mut ri_out, &mut ali_out, nd);
// 累加 Jν += w_k * (I_in + I_out) / 2
for id in 0..nd {
jnu[id][ij] += w_k * 0.5 * (ri_in[id] + ri_out[id]);
}
}
// 安全夹紧
for id in 0..nd {
if jnu[id][ij] < 0.0 || jnu[id][ij].is_nan() {
jnu[id][ij] = source[id];
}
rad1_ij[id] = jnu[id][ij];
}
// 近似 Eddington 因子: f = K/J ≈ 1/3 (各向同性)
// 在深处修正: f ≈ J/(4πB) 或更精确地由形式解计算
// 对于 LTE 初始迭代，使用 1/3 即可
rad1_data.push(Rad1PointData {
rad1: rad1_ij,
fak1: fak1_ij,
});
// 构建不透明度数据用于 Lucy
let mut abso1_ij = vec![0.0; nd];
let mut abso1l_ij = vec![0.0; nd];
let mut emis1l_ij = vec![0.0; nd];
let mut scat1_ij = vec![0.0; nd];
for id in 0..nd {
let dens = params.model.modpar.dens[id];
let dens_safe = dens.max(1e-30);
let ne = ne_arr[id];
abso1_ij[id] = chi[id][ij] * dens_safe; // cm⁻¹
abso1l_ij[id] = abso1_ij[id]; // LTE: 全部为真吸收 + 散射
emis1l_ij[id] = ab_true[id][ij] * dens_safe * source[id]; // 发射 = 真吸收 × Bν
scat1_ij[id] = SIGE * ne; // 电子散射 (cm⁻¹)
}
opacfl_data.push(OpacflPointData {
abso1: abso1_ij,
abso1l: abso1l_ij,
emis1l: emis1l_ij,
scat1: scat1_ij,
});
}
// ==============================================================
// Step 4: 计算辐射压力梯度 pradt
// ==============================================================
// dP_rad/dm = (4π/c) Σ_ν κ_ν^true × (J_ν - B_ν) × w_ν
// Lucy 用 pradt 来修正流体静力学平衡
let mut pradt = vec![0.0; nd];
for id in 0..nd {
let t = params.model.modpar.temp[id];
let hkt = HK / t;
let mut sum_prad = 0.0;
for ij in 0..nfreq {
let fr = freq[ij];
let x = (hkt * fr).min(150.0);
let bnu = h_over_c2 * fr.powi(3) / (x.exp() - 1.0);
let jnu_val = rad1_data[ij].rad1[id];
// 真吸收系数 (不含散射) per gram
let ab_true_gram = opacfl_data[ij].abso1l[id]
/ params.model.modpar.dens[id].max(1e-30);
sum_prad += ab_true_gram * (jnu_val - bnu) * weights[ij];
}
pradt[id] = 4.0 * PI / c_light * sum_prad;
}
// 辅助数组
let dens1: Vec<f64> = params.model.modpar.dens.iter()
.map(|d| 1.0 / d.max(1e-30))
.collect();
let vturb = vec![0.0; nd];
let mut pgs = vec![0.0; nd];
let lucy_config = LucyConfig {
itlucy: 1, // 单次 Lambda 迭代只做一步 Lucy
iaclt: 10,
iacldt: 1,
ihecor: 1,
lte: config.lte,
lchc: config.lchc,
ielcor: config.ielcor,
iter,
ntrl: 0,
iluctr: vec![0; MTRANS],
};
let lucy_model = LucyModelParams {
nd,
nfreq,
ntrans: config.ntrans,
teff,
grav,
prd0,
dm: &params.model.modpar.dm[..nd],
temp: &params.model.modpar.temp[..nd],
elec: &params.model.modpar.elec[..nd],
dens: &params.model.modpar.dens[..nd],
dens1: &dens1,
wmm: &wmm_arr,
deldm: &deldm,
vturb: &vturb,
pradt: &pradt,
pgs: &mut pgs,
freq: &freq[..nfreq],
w: &weights[..nfreq],
fh: &fh,
hextrd: &hextrd,
lac2t: false,
};
let lucy_output = lucy_pure(&lucy_config, &lucy_model, &opacfl_data, &rad1_data);
// ==============================================================
// Step 5: 更新模型状态
// ==============================================================
for id in 0..nd {
let t_new = lucy_output.temp[id].max(3000.0).min(200000.0);
params.model.modpar.temp[id] = t_new;
params.model.modpar.elec[id] = lucy_output.elec[id];
// Density from Lucy may be wrong at surface (nearly fully ionized gas gives rho~0)
// Keep original density; will be updated properly when SOLVE is implemented
if lucy_output.dens[id] > 1e-25 {
params.model.modpar.dens[id] = lucy_output.dens[id];
} params.model.modpar.sqt1[id] = t_new.sqrt();
params.model.modpar.hkt1[id] = HK / t_new;
params.model.modpar.tk1[id] = 1.0 / t_new;
}
// 存储辐射场到 TotRad (供后续 SOLVE 使用)
for ij in 0..nfreq {
for id in 0..nd {
params.model.totrad.rad[ij][id] = rad1_data[ij].rad1[id];
params.model.totrad.fak[ij][id] = rad1_data[ij].fak1[id];
}
}
// 诊断输出
if config.iprind == 2 {
// output(writer, &OutputParams { ... });
debug_log!("RESOLV: OUTPUT called (iprind=2)");
if let Some(w) = writer.as_mut() {
let _ = w.write_raw(&format!(
" Lambda iter {}: max dH/H = {:.4e}, T_surf = {:.0}, T_bottom = {:.0}\n",
ilam, lucy_output.dhhmx1,
params.model.modpar.temp[0], params.model.modpar.temp[nd - 1],
));
}
}
debug_log!("RESOLV: Lambda iter {} done, dhhmx={:.4e}", ilam, lucy_output.dhhmx1);
// 加速收敛
if config.iacpp > 0 {
// accelp(&mut AccelpParams { ... });
debug_log!("RESOLV: ACCELP called (iacpp={})", config.iacpp);
}
// Lucy 迭代
// lucy_pure(&LucyParams { ... });
debug_log!("RESOLV: LUCY called");
}
// -----------------------------------------------------------
@@ -401,9 +693,14 @@ pub fn resolv<W: std::io::Write>(
debug_log!("RESOLV: ROSSTD(0) called");
}
// 输出模型
// output(writer, &OutputParams { ... });
debug_log!("RESOLV: OUTPUT called");
// 输出模型 — 对应 Fortran CALL OUTPUT (line 84)
if let Some(w7) = params.writer7.as_mut() {
let output_params = OutputParams {
config: params.tlusty_config,
model: params.model,
};
let _ = output(w7, &output_params, None, None);
}
// -----------------------------------------------------------
// Part 6: 压力评估
@@ -539,10 +836,15 @@ pub fn resolv<W: std::io::Write>(
}
// -----------------------------------------------------------
// Part 13: 输出压缩模型到 fort.7
// Part 13: 输出压缩模型到 fort.7 — 对应 Fortran CALL OUTPUT (line 159)
// -----------------------------------------------------------
// output(writer, &OutputParams { ... });
debug_log!("RESOLV: OUTPUT called (model to fort.7)");
if let Some(w7) = params.writer7.as_mut() {
let output_params = OutputParams {
config: params.tlusty_config,
model: params.model,
};
let _ = output(w7, &output_params, None, None);
}
// -----------------------------------------------------------
// Part 14: 最终输出
@@ -576,8 +878,8 @@ pub fn resolv<W: std::io::Write>(
}
/// 最终输出处理。
fn final_output<W: std::io::Write>(
params: &mut ResolvParams,
fn final_output<W: std::io::Write, W7: std::io::Write>(
params: &mut ResolvParams<W7>,
_writer: Option<&mut FortranWriter<W>>,
) -> ResolvOutput {
let config = &params.config;
@@ -636,7 +938,7 @@ fn final_output<W: std::io::Write>(
/// 纯计算版本的 RESOLV（无 I/O 操作）。
///
/// 用于测试和嵌入式使用。
pub fn resolv_pure(params: &mut ResolvParams) -> ResolvOutput {
pub fn resolv_pure<W7: std::io::Write>(params: &mut ResolvParams<W7>) -> ResolvOutput {
resolv(params, None::<&mut FortranWriter<std::io::Empty>>)
}
@@ -667,31 +969,29 @@ mod tests {
#[test]
fn test_resolv_pure_basic() {
// 创建默认配置
// 创建默认配置 - 使用小网格避免 OOM
let config = ResolvConfig {
iter: 1,
init: 1,
lfin: false,
nd: 10,
nfreq: 100,
nd: 5,
nfreq: 10,
..Default::default()
};
// 创建最小化的状态
let mut tlusty_config = TlustyConfig::default();
tlusty_config.inppar.teff = 10000.0;
tlusty_config.inppar.grav = 1e4;
let mut atomic = AtomicData::default();
let mut model = ModelState::new();
// 初始化模型温度
for i in 0..10 {
model.modpar.temp[i] = 10000.0 - i as f64 * 500.0;
}
let mut model = create_minimal_model(5);
let mut params = ResolvParams {
config,
tlusty_config: &mut tlusty_config,
atomic: &mut atomic,
model: &mut model,
writer7: None,
};
// 执行 RESOLV
@@ -701,6 +1001,19 @@ mod tests {
assert_eq!(result.iter, 1);
}
/// 辅助函数：创建具有最小有效物理数据的 ModelState
fn create_minimal_model(nd: usize) -> ModelState {
let mut model = ModelState::new();
for i in 0..nd {
let t = 10000.0 - i as f64 * 500.0;
model.modpar.temp[i] = t;
model.modpar.elec[i] = 1e12;
model.modpar.dens[i] = 1e-12;
model.modpar.dm[i] = (i + 1) as f64 * 1e-6;
}
model
}
#[test]
fn test_resolv_final_iteration() {
// 测试最终迭代
@@ -708,20 +1021,22 @@ mod tests {
iter: 5,
init: 0,
lfin: true,
nd: 10,
nfreq: 100,
nd: 5,
nfreq: 10,
..Default::default()
};
let mut tlusty_config = TlustyConfig::default();
tlusty_config.inppar.teff = 10000.0;
let mut atomic = AtomicData::default();
let mut model = ModelState::new();
let mut model = create_minimal_model(5);
let mut params = ResolvParams {
config,
tlusty_config: &mut tlusty_config,
atomic: &mut atomic,
model: &mut model,
writer7: None,
};
let result = resolv_pure(&mut params);
@@ -738,20 +1053,22 @@ mod tests {
init: 1,
lfin: false,
lte: true,
nd: 10,
nfreq: 100,
nd: 5,
nfreq: 10,
..Default::default()
};
let mut tlusty_config = TlustyConfig::default();
tlusty_config.inppar.teff = 10000.0;
let mut atomic = AtomicData::default();
let mut model = ModelState::new();
let mut model = create_minimal_model(5);
let mut params = ResolvParams {
config,
tlusty_config: &mut tlusty_config,
atomic: &mut atomic,
model: &mut model,
writer7: None,
};
let result = resolv_pure(&mut params);
@@ -769,20 +1086,22 @@ mod tests {
iconre: 5,
iconrs: 1,
ipconf: 1,
nd: 10,
nfreq: 100,
nd: 5,
nfreq: 10,
..Default::default()
};
let mut tlusty_config = TlustyConfig::default();
tlusty_config.inppar.teff = 10000.0;
let mut atomic = AtomicData::default();
let mut model = ModelState::new();
let mut model = create_minimal_model(5);
let mut params = ResolvParams {
config,
tlusty_config: &mut tlusty_config,
atomic: &mut atomic,
model: &mut model,
writer7: None,
};
let result = resolv_pure(&mut params);
src/tlusty/main.rs
+10
View File
@@ -52,6 +52,7 @@
use std::io;
use std::time::Instant;
use std::fs::File;
use super::io::{
FortranReader, FortranWriter,
@@ -260,6 +261,10 @@ pub fn run_tlusty<R: io::BufRead, W: io::Write>(
// 临时文件（对应 OPEN(UNIT=91,92,93)）
let mut _scratch = ScratchFiles::default();
// fort.7 输出文件（Fortran OUTPUT 写入 UNIT=7）
let fort7_file = File::create("fort.7").expect("Cannot create fort.7");
let mut fort7_writer = FortranWriter::new(fort7_file);
// 从配置中获取 NITER
state.niter = if config.runkey.niter > 0 { config.runkey.niter } else { 100 };
@@ -307,6 +312,10 @@ pub fn run_tlusty<R: io::BufRead, W: io::Write>(
let nd = state.nd;
let mut work_arrays = TlustyWorkArrays::new(nn, nd);
// 打开 fort.7 输出文件（对应 Fortran OPEN(UNIT=7,...) 在 OUTPUT 中使用）
let fort7_file = File::create("fort.7").expect("Cannot create fort.7");
let mut fort7_writer = FortranWriter::new(fort7_file);
// ========================================
// 主迭代循环
// 对应 Fortran:
@@ -374,6 +383,7 @@ pub fn run_tlusty<R: io::BufRead, W: io::Write>(
tlusty_config: config,
atomic: &mut atomic,
model: &mut model,
writer7: Some(&mut fort7_writer),
};
resolv(&mut resolv_params, Some(output_writer));
}
src/tlusty/math/continuum/lte_opacity.rs
+1 -1
View File
@@ -231,7 +231,7 @@ pub fn lte_meanopt(params: &LteOpacityParams, grid: &LteFrequencyGrid) -> LteOpa
}
/// 计算给定频率点的吸收和散射系数 (per cm³)。
fn compute_opacity_at_frequency(
pub fn compute_opacity_at_frequency(
fr: f64,
t: f64,
ne: f64,
src/tlusty/math/continuum/mod.rs
+1
View File
@@ -31,6 +31,7 @@ pub use crate::tlusty::math::opacity::{cia_h2h, cia_h2h2, cia_h2he, cia_hhe};
pub use lte_opacity::{
LteOpacityParams, LteOpacityOutput, LteFrequencyGrid,
lte_meanopt, generate_lte_frequency_grid, quick_lte_rosseland,
compute_opacity_at_frequency,
};
pub use lte_opacity::{
LteOpacityParams as LteOpacityParamsOld, LteOpacityOutput as LteOpacityOutputOld,
src/tlusty/math/opacity/corrwm.rs
+12 -12
View File
@@ -470,14 +470,14 @@ mod tests {
#[test]
fn test_corrwm_basic() {
let (mut basnum, trapar, mut frqall, mut freaux, phoexp) = create_test_state();
let (mut basnum, trapar, mut frqall, mut freaux, mut phoexp) = create_test_state();
let mut params = CorrwmParams {
basnum: &mut basnum,
trapar: &trapar,
frqall: &mut frqall,
freaux: &mut freaux,
phoexp: &phoexp,
phoexp: &mut phoexp,
};
corrwm(&mut params);
@@ -496,14 +496,14 @@ mod tests {
#[test]
fn test_corrwm_lskip_radiation_pressure() {
let (mut basnum, trapar, mut frqall, mut freaux, phoexp) = create_test_state();
let (mut basnum, trapar, mut frqall, mut freaux, mut phoexp) = create_test_state();
let mut params = CorrwmParams {
basnum: &mut basnum,
trapar: &trapar,
frqall: &mut frqall,
freaux: &mut freaux,
phoexp: &phoexp,
phoexp: &mut phoexp,
};
corrwm(&mut params);
@@ -534,14 +534,14 @@ mod tests {
#[test]
fn test_corrwm_w0e_bnue_wc() {
let (mut basnum, trapar, mut frqall, mut freaux, phoexp) = create_test_state();
let (mut basnum, trapar, mut frqall, mut freaux, mut phoexp) = create_test_state();
let mut params = CorrwmParams {
basnum: &mut basnum,
trapar: &trapar,
frqall: &mut frqall,
freaux: &mut freaux,
phoexp: &phoexp,
phoexp: &mut phoexp,
};
corrwm(&mut params);
@@ -565,7 +565,7 @@ mod tests {
#[test]
fn test_corrwm_rybicki_mode() {
let (mut basnum, trapar, mut frqall, mut freaux, phoexp) = create_test_state();
let (mut basnum, trapar, mut frqall, mut freaux, mut phoexp) = create_test_state();
// 启用 Rybicki 模式
basnum.ifryb = 1;
@@ -575,7 +575,7 @@ mod tests {
trapar: &trapar,
frqall: &mut frqall,
freaux: &mut freaux,
phoexp: &phoexp,
phoexp: &mut phoexp,
};
corrwm(&mut params);
@@ -591,7 +591,7 @@ mod tests {
#[test]
fn test_corrwm_skip_all_radiation_pressure() {
let (mut basnum, trapar, mut frqall, mut freaux, phoexp) = create_test_state();
let (mut basnum, trapar, mut frqall, mut freaux, mut phoexp) = create_test_state();
// 跳过所有辐射压力
basnum.ifprad = 0;
@@ -601,7 +601,7 @@ mod tests {
trapar: &trapar,
frqall: &mut frqall,
freaux: &mut freaux,
phoexp: &phoexp,
phoexp: &mut phoexp,
};
corrwm(&mut params);
@@ -619,14 +619,14 @@ mod tests {
#[test]
fn test_corrwm_io() {
let (mut basnum, trapar, mut frqall, mut freaux, phoexp) = create_test_state();
let (mut basnum, trapar, mut frqall, mut freaux, mut phoexp) = create_test_state();
let mut params = CorrwmParams {
basnum: &mut basnum,
trapar: &trapar,
frqall: &mut frqall,
freaux: &mut freaux,
phoexp: &phoexp,
phoexp: &mut phoexp,
};
let mut writer = FortranWriter::to_memory();
src/tlusty/math/opacity/quasim.rs
+4 -4
View File
@@ -167,8 +167,8 @@ mod tests {
#[test]
fn test_quasim_disabled() {
let (model, atomic, basnum, freq) = create_test_data();
let result = quasim(0, &model, &atomic, &basnum, &freq);
let (mut model, atomic, basnum, freq) = create_test_data();
let result = quasim(0, &mut model, &atomic, &basnum, &freq);
// 当禁用时，所有值应为 0
for &val in &result.sgd {
assert_relative_eq!(val, 0.0, epsilon = 1e-30);
@@ -180,7 +180,7 @@ mod tests {
let (mut model, atomic, basnum, mut freq) = create_test_data();
model.quasun.iquasi = 1;
freq[0] = 1e15; // 波长 < 911 Å
let result = quasim(0, &model, &atomic, &basnum, &freq);
let result = quasim(0, &mut model, &atomic, &basnum, &freq);
for &val in &result.sgd {
assert_relative_eq!(val, 0.0, epsilon = 1e-30);
}
@@ -207,7 +207,7 @@ mod tests {
}
}
let result = quasim(0, &model, &atomic, &basnum, &freq);
let result = quasim(0, &mut model, &atomic, &basnum, &freq);
// 应该返回非零值
assert!(result.sgd.iter().any(|&v| v > 0.0));
}
src/tlusty/math/solvers/accel2.rs
+2 -1
View File
@@ -119,7 +119,8 @@ pub struct Accel2Output {
/// 当 `need_resolv` 为 true 时，调用者需要执行 RESOLV 重新计算
pub fn accel2_pure(params: &mut Accel2Params) -> Accel2Output {
// 标记 RESOLV 依赖（Fortran 在加速后调用 RESOLV）
let _ = resolv_pure;
// 标记依赖（确保编译器知道 accel2 依赖 resolv_pure）
let _ = "accel2 depends on resolv_pure";
let config = &mut params.config;
let iter = config.iter;
src/tlusty/math/temperature/lucy.rs
+4 -2
View File
@@ -393,7 +393,8 @@ pub fn lucy_pure(
let mut xx1 = 0.0;
// 表面点
let tp3 = model.temp[0].powi(3);
// Fortran: TP3 = TEMP1(ID)**3, where TEMP1 = 1/T
let tp3 = model.temp[0].powi(-3);
xx = state.eddf[0] / state.eddh * state.delh[0];
xx1 = xx;
state.dt1[0] = state.heat[0] / 16.0 / SIG4P * tp3 / state.absp[0];
@@ -402,7 +403,8 @@ pub fn lucy_pure(
// 内部点
for id in 1..nd {
let tp3 = model.temp[id].powi(3);
// Fortran: TP3 = TEMP1(ID)**3, where TEMP1 = 1/T
let tp3 = model.temp[id].powi(-3);
xx += model.deldm[id - 1]
* (state.absh[id - 1] * model.dens1[id - 1] * state.delh[id - 1]
+ state.absh[id] * model.dens1[id] * state.delh[id]);
src/tlusty/math/utils/change.rs
+1
View File
@@ -312,6 +312,7 @@ pub fn saha_factor(s: f64, t: f64, g1: f64, g2: f64, enion: f64) -> f64 {
#[cfg(test)]
mod tests {
use super::*;
use crate::tlusty::state::constants::MLEVEL;
#[test]
fn test_saha_factor() {
src/tlusty/math/utils/comset.rs
+5 -5
View File
@@ -394,7 +394,7 @@ mod tests {
params.kij[i] = 5 - i;
}
let result = comset(&params);
let result = comset(&params, None);
// 当 ICOMPT = 0 时，只计算 SIGEC
// 检查 SIGEC 有限
@@ -406,7 +406,7 @@ mod tests {
#[test]
fn test_comset_basic() {
let params = create_test_params();
let result = comset(&params);
let result = comset(&params, None);
// 检查 IJORIG 映射 (只检查前 nfreq 个元素)
for i in 0..params.nfreq {
@@ -441,7 +441,7 @@ mod tests {
let mut params = create_test_params();
params.ichcoo = 1; // 高阶模式
let result = comset(&params);
let result = comset(&params, None);
// 检查结果有限
for i in 0..params.nfreq {
@@ -455,7 +455,7 @@ mod tests {
let mut params = create_test_params();
params.knish = 1; // 完整 Klein-Nishina
let result = comset(&params);
let result = comset(&params, None);
// 检查 SIGEC 有限且为正
for i in 0..params.nfreq {
@@ -477,7 +477,7 @@ mod tests {
// 高频 (xf > 1000)
params.freq[2] = 1e24; // xf ≈ 8e3
let result = comset(&params);
let result = comset(&params, None);
// 所有 SIGEC 应为正且有限
for i in 0..3 {