diff --git a/src/tlusty/math/convection/conref.rs b/src/tlusty/math/convection/conref.rs index 8fc8d7d..7b27421 100644 --- a/src/tlusty/math/convection/conref.rs +++ b/src/tlusty/math/convection/conref.rs @@ -232,10 +232,8 @@ pub fn conref_pure(params: &mut ConrefParams) -> ConrefOutput { let mut flxtt = vec![0.0f64; nd]; let mut delta0 = vec![0.0f64; nd]; - // 保存初始 delta - for id in 0..nd { - delta0[id] = params.delta[id]; - } + // Fortran: delta0(mdepth) 未初始化（-fno-automatic 下全为 0.0） + // Rust: vec![0.0; nd] 已初始化为 0.0，无需额外赋值 // 第一遍：计算对流梯度 let mut icbeg: usize = 0; @@ -362,7 +360,7 @@ pub fn conref_pure(params: &mut ConrefParams) -> ConrefOutput { } // Fortran: icbegp=icbeg0; icendp=icend — 保存状态用于下次调用 let icbegp_save = icbeg0; - let icendp_save = icend; + let mut icendp_save = icend; // 如果没有对流区，直接返回 if icbeg0 == 0 || icend == 0 { @@ -524,6 +522,7 @@ pub fn conref_pure(params: &mut ConrefParams) -> ConrefOutput { // Fortran: icbeg0=icbeg; icend=nd; icendp=nd let imucon_icbeg0 = icbeg; let imucon_icend = nd; + icendp_save = nd; // Fortran: icendp=nd (overwrites earlier icendp=icend) // FORMAT 674: modification with imucon eprintln!("\n modification with imucon: icbeg0,icend{:4}{:4}{:4}", config.imucon, imucon_icbeg0, imucon_icend); @@ -559,51 +558,71 @@ pub fn conref_pure(params: &mut ConrefParams) -> ConrefOutput { cubcon.b = bcnv; cubcon.grdadb = grdadb; - let alp = params.flrd[idx].min(flxtt[idx]) / t0.powi(4) / dlt; + let mut alp = params.flrd[idx].min(flxtt[idx]) / t0.powi(4) / dlt; let fcnv = flxtt[idx] - params.flrd[idx]; if fcnv <= 0.0 || fc0 <= 0.0 { // 辐射主导 — Fortran label 200-220 - if flxtt[idx] < params.flrd[idx] { - let mut alp_val = flxtt[idx] / params.flrd[idx] * t0.powi(4) * delta0[idx]; - let mut t_iter = t; - let mut t0_i = t0; - let mut dlt_i = dlt; + // Fortran: 同一作用域的 t, t0, dlt 在 Newton 循环中被修改 + let mut t_loc = t; + let mut t0_loc = t0; + let mut dlt_loc = dlt; + let mut last_dele = 0.0f64; - for iter_idx in 0..20 { + if flxtt[idx] < params.flrd[idx] { + // Fortran line 272: alp=flxtt(id)/flrd(id)*t0**4*delta0(id) + alp = flxtt[idx] / params.flrd[idx] * t0_loc.powi(4) * delta0[idx]; + let mut t_iter = t_loc; + let mut t0_i = t0_loc; + let mut dlt_i = dlt_loc; + + // Fortran label 210: Newton iteration (itrnrc=1..21) + for iter_idx in 0..21 { let t1 = UN / t_iter; let dltp = t1 / (p / pm).ln(); - // Fortran: dele=alp-t0**4*dlt (使用 t0 和 dlt 而非 t_iter) - let dele = alp_val - t0_i.powi(4) * dlt_i; + let dele = alp - t0_i.powi(4) * dlt_i; + last_dele = dele; let delep = -t0_i.powi(4) * dlt_i * (2.0 * t1 + dltp / dlt_i); let dt = -dele / delep * t1; + // Fortran: t=t*(un+dt) BEFORE convergence check + t_iter = t_iter * (UN + dt); - // FORMAT 645: Newton iteration progress eprintln!("{:4}{:4}{:11.3e}{:8.2}", id, iter_idx + 1, dt, t_iter); - if dt.abs() < 1e-6 { + if dt.abs() < 1e-6 || iter_idx >= 20 { + // Fortran: exit with updated t but old t0, dlt break; } - t_iter = t_iter * (UN + dt); - // Fortran: 更新 t0 和 dlt 用于下一次迭代 t0_i = (t_iter * tm).sqrt(); dlt_i = (t_iter / tm).ln() / (p / pm).ln(); } - // Fortran: alp=flrd(id)/t0**4/dlt - alp_val = params.flrd[idx] / t0_i.powi(4) / dlt_i; - params.delta[idx] = dlt_i; - params.temp[idx] = t_iter; - modified = true; + // Fortran line 289: alp=flrd(id)/t0**4/dlt (old t0, dlt) + alp = params.flrd[idx] / t0_i.powi(4) / dlt_i; + // Fortran label 230: temp(id)=t (new), delta(id)=dlt (old) + t_loc = t_iter; // new t + t0_loc = t0_i; // old t0 + dlt_loc = dlt_i; // old dlt } - // Fortran label 230: flr=alp*t0**4*dlt - let dlt_val = params.delta[idx]; - let t0_diag = (params.temp[idx] * tm).sqrt(); - let flr_diag = alp * t0_diag.powi(4) * dlt_val; - if dlt_val < cubcon.grdadb { - eprintln!("{:4}{:4}{:8.1}{:12.4e}{:12.4e}{:12.4e}{:12.4e}", - id, 0, params.temp[idx], dlt_val, cubcon.grdadb, - params.flrd[idx] / flxtt[idx], flr_diag / flxtt[idx]); + // Fortran label 230: delta(id)=dlt, temp(id)=t, flr=alp*t0**4*dlt + params.delta[idx] = dlt_loc; + params.temp[idx] = t_loc; + modified = true; + { + // Fortran label 230: flr=alp*t0**4*dlt (old t0, old dlt, new alp) + let flr = alp * t0_loc.powi(4) * dlt_loc; + if dlt_loc >= cubcon.grdadb { + // Fortran: flc=fc0*dele (dele from Newton iteration) + let flc = fc0 * last_dele; + eprintln!("{:4}{:4}{:8.1}{:12.4e}{:12.4e}{:12.4e}{:12.4e}{:12.4e}{:12.4e}{:12.4e}", + id, 0, params.temp[idx], dlt_loc, cubcon.grdadb, + params.flrd[idx] / flxtt[idx], flr / flxtt[idx], flc / flxtt[idx], + 0.0_f64, (flr + 0.0_f64) / flxtt[idx]); + } else { + eprintln!("{:4}{:4}{:8.1}{:12.4e}{:12.4e}{:12.4e}{:12.4e}", + id, 0, params.temp[idx], dlt_loc, cubcon.grdadb, + params.flrd[idx] / flxtt[idx], flr / flxtt[idx]); + } } } else { // 对流主导 — Fortran label 100 (Newton 迭代) @@ -616,6 +635,7 @@ pub fn conref_pure(params: &mut ConrefParams) -> ConrefOutput { let mut t0_i = (t_iter * tm).sqrt(); let mut dlt_i = (t_iter / tm).ln() / (p / pm).ln(); let mut fc0_i = fc0; + let mut flxcn0_conv = 0.0f64; let mut itrnrc = 0; // Fortran: label 100 Newton 迭代循环 @@ -639,7 +659,7 @@ pub fn conref_pure(params: &mut ConrefParams) -> ConrefOutput { t0_i = (t_iter * tm).sqrt(); dlt_i = (t_iter / tm).ln() / (p / pm).ln(); // Fortran: call convc1(id,t0,p0,pg0,prad0,ab0,dlt,flxcn0,fc0) - let (_, fc0_new, _, _) = compute_convc1_simple( + let (flxcn0_conv, fc0_new, _, _) = compute_convc1_simple( idx + 1, t0_i, p0, pg0, pr0, ab0, dlt_i, config, flxtt[idx], cubcon.gravd, ); fc0_i = fc0_new; @@ -654,13 +674,13 @@ pub fn conref_pure(params: &mut ConrefParams) -> ConrefOutput { modified = true; // Fortran label 230: flr=alp*t0**4*dlt let flr = alp * t0_i.powi(4) * dlt_i; - // Fortran: flc=fc0*dele (使用最后迭代值) - let flc = fc0_i * (flxtt[idx] - alp * t0_i.powi(4) * dlt_i) / fc0_i; + // Fortran: flc=fc0*dele + let flc = flxtt[idx] - alp * t0_i.powi(4) * dlt_i; // FORMAT 646: convection-dominant depth result eprintln!("{:4}{:4}{:8.1}{:12.4e}{:12.4e}{:12.4e}{:12.4e}{:12.4e}{:12.4e}{:12.4e}", id, itrnrc, params.temp[idx], dlt_i, cubcon.grdadb, params.flrd[idx] / flxtt[idx], flr / flxtt[idx], flc / flxtt[idx], - 0.0_f64, (flr + 0.0_f64) / flxtt[idx]); + flxcn0_conv / flxtt[idx], (flr + flxcn0_conv) / flxtt[idx]); } } } @@ -695,12 +715,13 @@ pub fn conref_pure(params: &mut ConrefParams) -> ConrefOutput { /// 简化的对流计算 (内部使用)。 /// -/// 返回 (flxcnv, vcon, grdadb, bcnv) +/// 使用理想气体近似替代 TRMDER 热力学导数。 +/// 返回 (flxcnv, vconv, grdadb, bcnv) fn compute_convection_simple( _id: usize, t0: f64, pt0: f64, - _pg0: f64, + pg0: f64, _pr0: f64, ab0: f64, dlt: f64, @@ -713,7 +734,13 @@ fn compute_convection_simple( return (0.0, 0.0, 0.0, 0.0); } - // 绝热梯度近似 (单原子理想气体 = 0.4) + // 理想气体近似热力学量 (替代 Fortran TRMDER 输出) + const KB: f64 = 1.380649e-16; // Boltzmann constant (erg/K) + const MH: f64 = 1.6726e-24; // proton mass (g) + const HEATCP_IDEAL: f64 = 2.07e8; // 5/2 * KB/MH (erg/K/g) + const DLRDLT_IDEAL: f64 = 1.0; // d(ln rho)/d(ln T) for ideal gas + + // 绝热梯度 (单原子理想气体, gamma=5/3, ∇_ad = 1 - 1/gamma = 0.4) let grdadb = 0.4; // 检查对流不稳定性 @@ -722,15 +749,18 @@ fn compute_convection_simple( return (0.0, 0.0, grdadb, 0.0); } - // 计算引力 - let grav = if config.idisk == 1 { gravd } else { config.grav }; - if grav == 0.0 { + // 计算引力 (Fortran: if(idisk.eq.0) use GRAV, else use GRAVD) + let grav_val = if config.idisk == 1 { gravd } else { config.grav }; + if grav_val == 0.0 { return (0.0, 0.0, grdadb, 0.0); } - // 粗略估计密度 - let rho = if t0 > 0.0 { - pt0 / (t0 * 1.38e-16 * grav) + // 密度估计: rho = P_gas * m_H / (k_B * T) (理想气体状态方程) + // Fortran TRMDER 从完整 EOS 计算密度 + let rho = if t0 > 0.0 && pg0 > 0.0 { + pg0 * MH / (KB * t0) + } else if t0 > 0.0 { + pt0 * MH / (KB * t0) } else { 1e-7 }; @@ -738,29 +768,29 @@ fn compute_convection_simple( // 混合长度 let hmix = if config.hmix0 == 0.0 { 1.0 } else { config.hmix0 }; - // 压力标高 - let hscale = pt0 / rho / grav; + // 压力标高: HSCALE = PTOT / (RHO * GRAV) + let hscale = pt0 / rho / grav_val; - // 对流速度 - let vco = hmix * (config.aconml * pt0 / rho).abs().sqrt(); + // 对流速度: VCO = HMIX * sqrt(|ACONML * PTOT / RHO * DLRDLT|) + let vco = hmix * (config.aconml * pt0 / rho * DLRDLT_IDEAL).abs().sqrt(); - // 对流系数 - let flco = config.bconml * rho * t0 * hmix / 12.5664; + // 对流系数: FLCO = BCONML * RHO * HEATCP * T * HMIX / (4*PI) + let flco = config.bconml * rho * HEATCP_IDEAL * t0 * hmix / 12.5664; - // 光学厚度 + // 光学厚度: TAUE = HMIX * ABROS * RHO * HSCALE let taue = hmix * ab0 * rho * hscale; - // 辐射耗散因子 + // 辐射耗散因子: FAC = TAUE / (1 + 0.5*TAUE^2) let fac = taue / (UN + HALF * taue * taue); - // 参数 B (参考 Mihalas) — 这是 COMMON/CUBCON/ 中的 BCNV - let b = 5.67e-5 * t0.powi(3) / (rho * vco) * fac * config.cconml * HALF; + // 参数 B: B = sigma * T^3 / (RHO * HEATCP * VCO) * FAC * CCONML * 0.5 + let b = 5.67e-5 * t0.powi(3) / (rho * HEATCP_IDEAL * vco) * fac * config.cconml * HALF; // 参数 D let d = b * b / 2.0; let disc = d / 2.0 + ddel; - // 计算有效 DLT + // 计算有效 DLT = D + DDEL - B * sqrt(D/2 + DDEL) let dlt_eff = if disc >= 0.0 { let val = d + ddel - b * disc.sqrt(); if val < 0.0 { 0.0 } else { val } @@ -777,12 +807,13 @@ fn compute_convection_simple( /// 简化的 CONVC1 计算 (返回 fc0)。 /// +/// 与 CONVEC 类似，但在检查不稳定性之前计算 FC0 = FLCO * VCO。 /// 返回 (flxcnv, fc0, grdadb, bcnv) fn compute_convc1_simple( _id: usize, t0: f64, pt0: f64, - _pg0: f64, + pg0: f64, _pr0: f64, ab0: f64, dlt: f64, @@ -795,18 +826,26 @@ fn compute_convc1_simple( return (0.0, 0.0, 0.0, 0.0); } - // 绝热梯度近似 + // 理想气体近似热力学量 + const KB: f64 = 1.380649e-16; + const MH: f64 = 1.6726e-24; + const HEATCP_IDEAL: f64 = 2.07e8; + const DLRDLT_IDEAL: f64 = 1.0; + + // 绝热梯度 let grdadb = 0.4; // 计算引力 - let grav = if config.idisk == 1 { gravd } else { config.grav }; - if grav == 0.0 { + let grav_val = if config.idisk == 1 { gravd } else { config.grav }; + if grav_val == 0.0 { return (0.0, 0.0, grdadb, 0.0); } - // 估计密度 - let rho = if t0 > 0.0 { - pt0 / (t0 * 1.38e-16 * grav) + // 密度估计 (理想气体状态方程) + let rho = if t0 > 0.0 && pg0 > 0.0 { + pg0 * MH / (KB * t0) + } else if t0 > 0.0 { + pt0 * MH / (KB * t0) } else { 1e-7 }; @@ -815,15 +854,15 @@ fn compute_convc1_simple( let hmix = if config.hmix0 == 0.0 { 1.0 } else { config.hmix0 }; // 压力标高 - let hscale = pt0 / rho / grav; + let hscale = pt0 / rho / grav_val; // 对流速度 - let vco = hmix * (config.aconml * pt0 / rho).abs().sqrt(); + let vco = hmix * (config.aconml * pt0 / rho * DLRDLT_IDEAL).abs().sqrt(); // 对流系数 - let flco = config.bconml * rho * t0 * hmix / 12.5664; + let flco = config.bconml * rho * HEATCP_IDEAL * t0 * hmix / 12.5664; - // FC0 = FLCO * VCO + // CONVC1 特有: FC0 = FLCO * VCO (在检查不稳定性之前计算) let fc0 = flco * vco; // 检查对流不稳定性 @@ -839,7 +878,7 @@ fn compute_convc1_simple( let fac = taue / (UN + HALF * taue * taue); // 参数 B - let b = 5.67e-5 * t0.powi(3) / (rho * vco) * fac * config.cconml * HALF; + let b = 5.67e-5 * t0.powi(3) / (rho * HEATCP_IDEAL * vco) * fac * config.cconml * HALF; // 参数 D let d = b * b / 2.0; diff --git a/src/tlusty/math/hydrogen/colis.rs b/src/tlusty/math/hydrogen/colis.rs index cedc0dc..fe79901 100644 --- a/src/tlusty/math/hydrogen/colis.rs +++ b/src/tlusty/math/hydrogen/colis.rs @@ -11,11 +11,12 @@ use crate::tlusty::math::cion; use crate::tlusty::math::{colh, ColhAtomicData, ColhOutput, ColhParams}; +use crate::tlusty::math::{colhe, ColheParams}; use crate::tlusty::math::cspec; use crate::tlusty::math::irc; use crate::tlusty::math::ylintp; -use crate::tlusty::state::constants::{EH, HK, TWO, UN}; -// f2r_depends: COLH, COLHE, CSPEC, IRC +use crate::tlusty::state::constants::{EH, HK, MLEVEL, TWO, UN}; +// f2r_depends: COLH, COLHE, CSPEC, IRC, CION, YLININTP // ============================================================================ // 常量 @@ -187,6 +188,8 @@ pub struct ColisParams<'a> { pub colh_params: Option, /// COLH 原子数据 pub colh_atomic: Option>, + /// COLHE 参数（如果需要调用 COLHE） + pub colhe_params: Option, } /// COLIS 输出结果 @@ -232,11 +235,27 @@ pub fn colis(params: &ColisParams) -> ColisOutput { }; // 调用 COLH 和 COLHE（如果有氢和氦） - if params.iath != 0 || params.iathe != 0 { - // 注意：COLH 和 COLHE 的结果需要乘以电子密度 - // 这里我们假设调用者已经处理了这些 + // Fortran: IF(IATH.NE.0) CALL COLH(ID,T,COL) + // IF(IATHE.NE.0) CALL COLHE(T,COL) + if params.iath != 0 { + if let (Some(colh_p), Some(colh_a)) = (¶ms.colh_params, ¶ms.colh_atomic) { + let mut colh_out = ColhOutput { col: &mut col[..] }; + colh(colh_p, colh_a, &mut colh_out); + } + } + if params.iathe != 0 { + if let Some(colhe_p) = ¶ms.colhe_params { + let colhe_result = colhe(colhe_p); + for (it, val) in colhe_result.col_rates.iter().enumerate() { + if it < col.len() { + col[it] += *val; + } + } + } + } - // 计算 CLOC 数组 + // Fortran: IF(IATH.NE.0.OR.IATHE.NE.0) THEN — 计算 CLOC 数组 + if params.iath != 0 || params.iathe != 0 { for it in 1..=params.ntrans { let it_idx = it - 1; if col[it_idx] != 0.0 { @@ -406,8 +425,11 @@ pub fn colis(params: &ColisParams) -> ColisOutput { let cs = creger(t32, u0_new, cc1, gg) * params.ane; col[it_idx] += cs; - // 注意：这里需要 WOP 数组 - // cloc[it_idx] += cs * params.ane * params.sbf[i - 1] * params.wop[i - 1 + (id - 1) * ?] * corr; + // Fortran: CLOC(IT)=CLOC(IT)+CS*ANE*SBF(I)*WOP(I,ID)*CORR + // WOP(MLEVEL,MDEPTH) column-major: linear index (ID-1)*MLEVEL + (I-1) + let wop_idx = (params.id - 1) * MLEVEL + (i - 1); + let wop_val = params.wop.get(wop_idx).copied().unwrap_or(0.0); + cloc[it_idx] += cs * params.ane * params.sbf[i - 1] * wop_val * corr; } } } else { @@ -439,7 +461,15 @@ pub fn colis(params: &ColisParams) -> ColisOutput { } else if ic == 4 { cupsx(srt, u0, c2 / params.g[i - 1]) * params.ane } else if ic == 9 { - params.omecol[it_idx] * srt0 * (-u0 * tt0).exp() * params.ane + // Fortran: CS=OMECOL(I,J)*SRT0*EXP(-U0*TT0) * ANE + // OMECOL(MLEVEL,MLEVEL) — column-major: (I-1) + (J-1)*MLEVEL + let omecol_idx = (i - 1) + (j - 1) * MLEVEL; + let omecol_val = if omecol_idx < params.omecol.len() { + params.omecol[omecol_idx] + } else { + 0.0 + }; + omecol_val * srt0 * (-u0 * tt0).exp() * params.ane } else if ic < 0 { cspec(i as i32, j as i32, ic, c1, c2, u0, t) * params.ane } else { @@ -457,10 +487,10 @@ pub fn colis(params: &ColisParams) -> ColisOutput { } /// 从 ITRA 数组获取跃迁索引 +/// Fortran: ITRA(I,J) — column-major array ITRA(MLEVEL,MLEVEL) +/// Linear index: (I-1) + (J-1)*MLEVEL fn get_itra(itra: &[i32], i: usize, j: usize, _nki: usize) -> i32 { - // 简化：假设 itra 是二维数组展平的 - // 实际索引方式取决于 Fortran 原始代码 - let idx = (i - 1) * 1000 + (j - 1); // 简化索引 + let idx = (i - 1) + (j - 1) * MLEVEL; if idx < itra.len() { itra[idx] } else { @@ -469,10 +499,10 @@ fn get_itra(itra: &[i32], i: usize, j: usize, _nki: usize) -> i32 { } /// 获取 OSH 振子强度 +/// Fortran: OSH(IQ,JQ) — column-major array OSH(20,20) +/// Linear index: (IQ-1) + (JQ-1)*20 fn get_osh(osh: &[f64], i: usize, j: usize) -> f64 { - // OSH(I, J) - 从能级 i 到 j 的振子强度 - // 假设是 (nlevel x 20) 的数组 - let idx = (i - 1) * 20 + (j - 1); + let idx = (i - 1) + (j - 1) * 20; if idx < osh.len() { osh[idx] } else { @@ -867,6 +897,7 @@ mod tests { ctemp_tab: &[[[0.0; MCFIT]; MXTCOL]], colh_params: None, colh_atomic: None, + colhe_params: None, }; let result = colis(¶ms); diff --git a/src/tlusty/math/hydrogen/sigave.rs b/src/tlusty/math/hydrogen/sigave.rs index fe666dc..e8ea0dc 100644 --- a/src/tlusty/math/hydrogen/sigave.rs +++ b/src/tlusty/math/hydrogen/sigave.rs @@ -88,23 +88,27 @@ fn read_element_data( itra: &[Vec], iup: &[i32], ilow: &[i32], - ibf: &[i32], + _ibf: &[i32], enion: &[f64], ierr: i32, izrr: i32, + ie: usize, ) -> Result)>> { let mut transitions = Vec::new(); - let mut itr = 0; for i in nl1..=nl2 { - itr += 1; + // Fortran: ITR = I in the inner loop (ITR starts at NL1-1, increments to NL1, ..., NL2) + let itr = i; // ITR = I always holds + if indexp.get(itr - 1).copied().unwrap_or(0) == 0 { continue; } + // Fortran: ITRA(IUP(ITR), ILOW(ITR)) — column-major + // Rust row-major: swap indices → itra[ILOW-1][IUP-1] let ic = itra - .get(iup.get(itr - 1).copied().unwrap_or(0) as usize - 1) - .and_then(|row| row.get(ilow.get(itr - 1).copied().unwrap_or(0) as usize - 1).copied()) + .get(ilow.get(itr - 1).copied().unwrap_or(0) as usize - 1) + .and_then(|row| row.get(iup.get(itr - 1).copied().unwrap_or(0) as usize - 1).copied()) .unwrap_or(0) as usize; // 读取跃迁数据 @@ -116,9 +120,22 @@ fn read_element_data( // Fe 特殊处理：能量校准 if ierr == 26 { ecmr = XIFE.get(izrr as usize - 1).copied().unwrap_or(0.0) - ecmr; + // Fortran: DE=ABS((ENION(ILOW(ITR))-HCCM*ECMR)/ENION(ILOW(ITR))) + let ilow_itr = ilow.get(itr - 1).copied().unwrap_or(0) as usize; + if ilow_itr > 0 { + let enion_val = enion.get(ilow_itr - 1).copied().unwrap_or(0.0); + if enion_val.abs() > 1e-30 { + let de = ((enion_val - HCCM * ecmr) / enion_val).abs(); + if de > 2e-2 { + eprintln!(" SIGAVE: energy discrepancy at ie={}, itr={}, i={}: DE={:.3e}", ie, itr, i, de); + } + } + } } - // 读取截面数据点（按频率递减顺序） + // 读取截面数据点 + // Fortran: DO IJ=1,NFIS; JI=NFIS-IJ+1; READ FRINSG(JI),CRIN(JI) + // Fortran reverses the data — read into descending frequency order let mut points = Vec::with_capacity(nfis); for _ in 0..nfis { let fr: f64 = reader.read_value()?; @@ -128,6 +145,8 @@ fn read_element_data( cross_section: cr, }); } + // Fortran reversal: file has ascending freq, Fortran stores descending + points.reverse(); transitions.push((ic, points)); } @@ -265,10 +284,12 @@ fn sigave_impl(params: &SigaveParams, mut reader: FortranReader) // 获取跃迁索引 let iup_val = params.iup.get(itr - 1).copied().unwrap_or(0) as usize; let ilow_val = params.ilow.get(itr - 1).copied().unwrap_or(0) as usize; + // Fortran: ITRA(IUP(ITR), ILOW(ITR)) — column-major + // Rust row-major: swap → itra[ILOW-1][IUP-1] let ic = params .itra - .get(iup_val - 1) - .and_then(|row| row.get(ilow_val - 1).copied()) + .get(ilow_val - 1) + .and_then(|row| row.get(iup_val - 1).copied()) .unwrap_or(0) as usize; if ic == 0 { @@ -318,6 +339,7 @@ fn sigave_impl(params: &SigaveParams, mut reader: FortranReader) params.enion, ierr, izrr, + ie, ) { Ok(t) => t, Err(e) => { @@ -343,6 +365,10 @@ fn sigave_impl(params: &SigaveParams, mut reader: FortranReader) ); transitions_read += transitions.len(); + + // Fortran: after inner DO 100 loop, ITR has been incremented (NL2-NL1+1) times + // Advance outer itr past all processed transitions + itr += nl2 - nl1 + 1; } SigaveOutput { diff --git a/src/tlusty/math/radiative/coolrt.rs b/src/tlusty/math/radiative/coolrt.rs index e14cb5d..ea9bf14 100644 --- a/src/tlusty/math/radiative/coolrt.rs +++ b/src/tlusty/math/radiative/coolrt.rs @@ -52,8 +52,8 @@ pub struct CoolrtParams<'a> { // 深度数据 /// 密度 (nd) pub dens: &'a [f64], - /// dm1 导数 (nd) - pub dedm1: &'a [f64], + /// dm1 导数（标量，来自 MODELQ COMMON/MODPAR/） + pub dedm1: f64, /// 深度间隔 (nd-1) pub deldmz: &'a [f64], @@ -192,6 +192,11 @@ pub fn coolrt_pure(params: &CoolrtParams) -> CoolrtOutput { // WRITE(87,686) ij,freq(ij) // FORMAT 686: i5,1pe15.7 eprintln!("{:5}{:15.7e}", ij + 1, params.freq[ij]); + // Fortran lines 53-57: taud 计算（结果未输出，仅内部使用） + let _taud = compute_taud( + nd, params.abso1[0], params.dedm1, + params.deldmz, params.absot, + ); } } @@ -293,6 +298,8 @@ pub fn find_fe2_ion( numat: &[i32], iz: &[i32], ) -> Option { + // Fortran iterates ALL ions and takes the LAST match (IOFE2 is overwritten) + let mut iofe2 = None; for ion in 0..nion { let nn1 = nfirst[ion]; if nn1 == 0 || nn1 > iatm.len() { @@ -300,10 +307,10 @@ pub fn find_fe2_ion( } let iat2 = iatm[nn1 - 1]; // Fortran 1-indexed -> Rust 0-indexed if iat2 > 0 && iat2 <= numat.len() && numat[iat2 - 1] == 26 && iz[ion] == 2 { - return Some(ion); + iofe2 = Some(ion); } } - None + iofe2 } // ============================================================================ @@ -347,7 +354,7 @@ mod tests { let w = vec![0.25, 0.25, 0.25, 0.25]; let freq = vec![1.0e14, 2.0e14, 3.0e14, 4.0e14]; let dens = vec![1.0e-6, 1.0e-5, 1.0e-4]; - let dedm1 = vec![1.0, 2.0, 3.0]; + let dedm1 = 1.0; let deldmz = vec![0.1, 0.2]; let abso1 = vec![0.1, 0.2, 0.3]; let emis1 = vec![0.01, 0.02, 0.03]; @@ -383,7 +390,7 @@ mod tests { w: Box::leak(w.into_boxed_slice()), freq: Box::leak(freq.into_boxed_slice()), dens: Box::leak(dens.into_boxed_slice()), - dedm1: Box::leak(dedm1.into_boxed_slice()), + dedm1, deldmz: Box::leak(deldmz.into_boxed_slice()), abso1: Box::leak(abso1.into_boxed_slice()), emis1: Box::leak(emis1.into_boxed_slice()),