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

重构4

Browse Source
This commit is contained in:
fmq
2026-03-22 00:30:20 +08:00
parent 497e62e13c
commit 59404207c3
27 changed files with 4035 additions and 4352 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/math/alifr1.rs
+1483
View File
File diff suppressed because it is too large Load Diff
src/math/cspec.rs
+199
View File
@@ -0,0 +1,199 @@
//! 碰撞强度计算。
//!
//! 重构自 TLUSTY `cspec.f`
//!
//! 使用 Van Regemorter 公式计算非标准的碰撞率，
//! 遵循 Mihalas (1978, Stellar Atmospheres, 2nd edition) 的建议。
//!
//! ## 碰撞类型 (IC)
//!
//! - IC = -1: 中性粒子（Neutrals）
//! - IC = -2: 离子（Ions）
//! - IC = -11: 特殊处理
//! - IC = -12: He I 禁戒跃迁
use crate::state::constants::UN;
/// He I 禁戒跃迁的系数 (CHE1FB)
/// 原始 Fortran DATA 语句:
/// CHE1FB(1,1-4) = 9.63675, -2.22941, -17.30103
/// CHE1FB(2,1-4) = 10.85578, -2.40931, -27.00903
/// CHE1FB(3,1-4) = 8.38043, -2.04791, -7.36621
const CHE1FB: [[f64; 4]; 3] = [
[9.63675, -2.22941, -17.30103, 0.0], // IFORB = 1
[10.85578, -2.40931, -27.00903, 0.0], // IFORB = 2
[8.38043, -2.04791, -7.36621, 0.0], // IFORB = 3
];
/// 指数积分近似系数 (Abramowitz & Stegun)
const EXPIA1: f64 = -0.57721566;
const EXPIA2: f64 = 0.99999193;
const EXPIA3: f64 = -0.24991055;
const EXPIA4: f64 = 0.05519968;
const EXPIA5: f64 = -0.00976004;
const EXPIA6: f64 = 0.00107857;
const EXPIB1: f64 = 0.2677737343;
const EXPIB2: f64 = 8.6347608925;
const EXPIB3: f64 = 18.059016973;
const EXPIB4: f64 = 8.5733287401;
const EXPIC1: f64 = 3.9584969228;
const EXPIC2: f64 = 21.0996530827;
const EXPIC3: f64 = 25.6329561486;
const EXPIC4: f64 = 9.5733223454;
/// 计算碰撞强度。
///
/// # 参数
///
/// * `i` - 下能级索引 (1-indexed)
/// * `j` - 上能级索引 (1-indexed)
/// * `ic` - 碰撞类型:
/// - -1: 中性粒子
/// - -2: 离子
/// - -11: 特殊处理
/// - -12: He I 禁戒跃迁
/// * `os` - 振子强度
/// * `cp` - 碰撞参数（用于离子）
/// * `u0` - 约化能量 (E/kT)
/// * `t` - 温度 (K)
///
/// # 返回值
///
/// 碰撞强度 CS
#[allow(clippy::too_many_arguments)]
pub fn cspec(i: i32, j: i32, ic: i32, os: f64, cp: f64, u0: f64, t: f64) -> f64 {
let mut cs = 0.0;
if ic > -10 {
// 计算指数积分 E1(U0)
let expiu0 = if u0 <= UN {
// U0 <= 1: 使用级数展开
(-u0.ln()) + EXPIA1
+ u0 * (EXPIA2 + u0 * (EXPIA3 + u0 * (EXPIA4 + u0 * (EXPIA5 + u0 * EXPIA6))))
} else {
// U0 > 1: 使用连分式近似
let num = EXPIB1 + u0 * (EXPIB2 + u0 * (EXPIB3 + u0 * EXPIB4));
let den = EXPIC1 + u0 * (EXPIC2 + u0 * (EXPIC3 + u0 * EXPIC4));
(-u0).exp() * num / den / u0
};
let gg = if ic == -1 {
// 中性粒子 (Auer & Mihalas 1973)
if u0 <= 14.0 {
0.276 * u0.exp() * expiu0
} else {
0.066 * (1.0 + 1.5 / u0) / u0.sqrt()
}
} else if ic == -2 {
// 离子 (Mihalas 1972)
let gg0 = 0.276 * u0.exp() * expiu0;
if gg0 > cp { gg0 } else { cp }
} else {
0.0
};
let t32 = t.powf(-1.5);
cs += 19.7363 * t32 * (-u0).exp() / u0 * gg * os;
return cs;
}
if ic == -11 {
// 特殊处理
let xr = -1.68_f64;
cs += 2.16 * u0.powf(xr) / t / t.sqrt() * (-u0).exp() * os;
return cs;
}
if ic == -12 {
// He I 禁戒跃迁 (from Klaus Werner)
// 确定跃迁类型
let iforb = match (i, j) {
(2, 3) => 1,
(2, 5) => 2,
(3, 4) => 3,
(4, 5) => 4,
_ => {
panic!("Inconsistent ICOL - CSPEC: i={}, j={}, iforb=0", i, j);
}
};
let xt = t.log10();
let gam = if iforb <= 3 {
CHE1FB[0][iforb as usize - 1]
+ CHE1FB[1][iforb as usize - 1] * xt
+ CHE1FB[2][iforb as usize - 1] / xt / xt
} else {
0.0
};
let gam = (2.30258509299405_f64 * gam).exp();
cs += 5.465e-11 * t.sqrt() * (-u0).exp() * gam;
}
cs
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
#[test]
fn test_cspec_neutral_low_u0() {
// 中性粒子，U0 <= 1
let cs = cspec(1, 2, -1, 0.5, 0.0, 0.5, 10000.0);
assert!(cs > 0.0);
}
#[test]
fn test_cspec_neutral_high_u0() {
// 中性粒子，U0 > 14
let cs = cspec(1, 2, -1, 0.5, 0.0, 15.0, 10000.0);
assert!(cs > 0.0);
}
#[test]
fn test_cspec_ion() {
// 离子
let cs = cspec(1, 2, -2, 0.5, 0.1, 2.0, 10000.0);
assert!(cs > 0.0);
}
#[test]
fn test_cspec_special() {
// IC = -11
let cs = cspec(1, 2, -11, 0.5, 0.0, 2.0, 10000.0);
assert!(cs > 0.0);
}
#[test]
fn test_cspec_hei_forbidden() {
// He I 禁戒跃迁 (2->3)
let cs = cspec(2, 3, -12, 0.0, 0.0, 2.0, 10000.0);
assert!(cs > 0.0);
}
#[test]
fn test_expiu0_low() {
// 测试指数积分 U0 <= 1
let u0: f64 = 0.5;
let expiu0 = (-u0.ln()) + EXPIA1
+ u0 * (EXPIA2 + u0 * (EXPIA3 + u0 * (EXPIA4 + u0 * (EXPIA5 + u0 * EXPIA6))));
// E1(0.5) ≈ 0.5598
assert_relative_eq!(expiu0, 0.5598, epsilon = 1e-3);
}
#[test]
fn test_expiu0_high() {
// 测试指数积分 U0 > 1
let u0: f64 = 2.0;
let num = EXPIB1 + u0 * (EXPIB2 + u0 * (EXPIB3 + u0 * EXPIB4));
let den = EXPIC1 + u0 * (EXPIC2 + u0 * (EXPIC3 + u0 * EXPIC4));
let expiu0 = (-u0).exp() * num / den / u0;
// E1(2.0) ≈ 0.0477 (Abramowitz & Stegun 近似)
assert_relative_eq!(expiu0, 0.0477, epsilon = 1e-3);
}
}
src/math/levset.rs
+427
View File
@@ -0,0 +1,427 @@
//! 能级参数设置。
//!
//! 重构自 TLUSTY `levset.f`
//!
//! 设置能级参数 IIEXP 和 IIFOR，控制能级的处理方式。
//!
//! ## 能级模式 (IMODL)
//!
//! - 0: 显式处理，完全 NLTE
//! - 1, 3: 跳过
//! - 4, 5: LTE
//! - -1, -3: 跳过，禁戒跃迁
//! - 6: 跳过
//! - -5, -6: 跳过，禁戒跃迁
//! - < -100: 分组处理
//! - < -200: 分组处理（不同方式）
use crate::state::constants::{MLEVEL, MDEPTH, MLVEXP};
/// LEVSET 的输入参数
pub struct LevsetParams {
/// 是否启用 LTE 模式
pub lte: bool,
/// 能级处理模式 (0: 自动, 其他: 由 IMODL 决定)
pub iflev: i32,
/// 选项表标志
pub ioptab: i32,
}
/// LEVSET 的模型状态
pub struct LevsetModelState<'a> {
/// 深度点数
pub nd: usize,
/// 原子数
pub natom: usize,
/// 能级数
pub nlevel: usize,
/// 能级模式 (NLEVEL)
pub imodl: &'a mut [i32],
/// 能级所属原子 (NLEVEL)
pub iatm: &'a [i32],
/// 能级所属元素/离子 (NLEVEL)
pub iel: &'a [i32],
/// 原子固定标志 (NATOM)
pub iifix: &'a [i32],
/// 原子起始能级 (NATOM)
pub n0a: &'a [i32],
/// 原子结束能级 (NATOM)
pub nka: &'a [i32],
/// 离子第一个能级 (NLEVEL)
pub nfirst: &'a [i32],
/// 离子最后一个能级 (NLEVEL)
pub nnext: &'a [i32],
/// LTE 能级标志 (NLEVEL)
pub iltlev: &'a [i32],
}
/// LEVSET 的输出状态
pub struct LevsetOutputState<'a> {
/// 显式能级索引 (NLEVEL) - 正值表示显式处理
pub iiexp: &'a mut [i32],
/// 禁戒能级索引 (NLEVEL)
pub iifor: &'a mut [i32],
/// 显式能级对应的实际能级 (MLVEXP)
pub indlev: &'a mut [i32],
/// LTE 参考能级 (NLEVEL × MDEPTH)
pub iltref: &'a mut [i32],
/// B 因子 (NLEVEL × MDEPTH)
pub bfac: &'a mut [f64],
/// 积分零标志 (MLVEXP × MDEPTH)
pub igzero: &'a mut [i32],
/// 种群零标志 (NLEVEL × MDEPTH)
pub ipzero: &'a mut [i32],
/// 显式能级数
pub nlvexp: &'a mut i32,
/// 禁戒能级数
pub nlvfor: &'a mut i32,
}
/// 设置能级参数。
///
/// # 参数
///
/// * `params` - 输入参数
/// * `model` - 模型状态
/// * `output` - 输出状态
pub fn levset(
params: &LevsetParams,
model: &mut LevsetModelState,
output: &mut LevsetOutputState,
) {
// 如果 IOPTAB < 0，直接返回
if params.ioptab < 0 {
return;
}
let nd = model.nd;
let natom = model.natom;
let nlevel = model.nlevel;
// 初始化 B 因子
for i in 0..nlevel {
for id in 0..nd {
output.bfac[i * nd + id] = 1.0;
}
}
if params.iflev == 0 {
// 情况 1: 由 IMODL 决定处理方式
process_by_imodl(params, model, output, natom, nlevel, nd);
} else {
// 情况 2: 自动处理 - 所有 ILK=0 的能级在更新的 LTE 模式下
process_automatic(params, model, output, nlevel, nd);
}
// 初始化 IGZERO 和 IPZERO
let nlvexp = *output.nlvexp as usize;
for ii in 0..nlvexp.min(MLVEXP) {
output.indlev[ii] = 0;
for id in 0..nd {
output.igzero[ii * nd + id] = 0;
}
}
for i in 0..nlevel {
for id in 0..nd {
output.ipzero[i * nd + id] = 0;
}
if model.imodl[i].abs() <= 6 {
if output.iiexp[i] > 0 {
let ii_idx = (output.iiexp[i] - 1) as usize;
if ii_idx < MLVEXP {
output.indlev[ii_idx] = (i + 1) as i32;
}
}
}
}
}
/// 按能级模式处理
#[allow(clippy::too_many_arguments)]
fn process_by_imodl(
params: &LevsetParams,
model: &mut LevsetModelState,
output: &mut LevsetOutputState,
natom: usize,
nlevel: usize,
nd: usize,
) {
// 初始化
for i in 0..nlevel {
output.iiexp[i] = 0;
output.iifor[i] = 0;
}
let mut iie: i32 = 0;
let mut iif: i32 = 0;
let mut igrp: i32 = 0;
for iat in 0..natom {
igrp = 0;
if model.iifix[iat] == 1 {
continue;
}
let n0a = model.n0a[iat] as usize;
let nka = model.nka[iat] as usize;
for i in n0a..=nka {
let imodl = model.imodl[i];
let mut inew = 1;
if imodl == 0 {
// 显式处理
iie += 1;
iif += 1;
output.iiexp[i] = iie;
output.iifor[i] = iif;
output.indlev[(iie - 1) as usize] = (i + 1) as i32;
} else if imodl > 0 {
// 正模式
iif += 1;
output.iifor[i] = iif;
if model.iltlev[i] >= 1 {
iie += 1;
output.iiexp[i] = iie;
}
let nfirst = model.nfirst[model.iel[i] as usize] as usize;
let nnext = model.nnext[model.iel[i] as usize] as usize;
if i + 1 == nfirst || i + 1 == nnext {
iie += 1;
output.iiexp[i] = iie;
}
} else if imodl < -100 {
// 分组处理 (IMODL < -100)
if i > n0a {
if model.imodl[i] == model.imodl[i - 1] {
inew = 0;
}
}
output.iiexp[i] = -iie;
if inew == 1 {
iie += 1;
output.iiexp[i] = -iie;
let nfirst = model.nfirst[model.iel[i] as usize] as usize;
let mut im = nfirst;
let mut lml = true;
while im < i && lml {
if model.imodl[i] == model.imodl[im] {
output.iiexp[i] = output.iiexp[im];
iie -= 1;
lml = false;
}
im += 1;
}
}
igrp = 1;
iif += 1;
output.iifor[i] = iif;
} else if imodl < -200 {
// 分组处理 (IMODL < -200)
if i > n0a {
if model.imodl[i] == model.imodl[i - 1] {
inew = 0;
}
}
if inew == 1 {
iie += 1;
}
if inew == 1 {
iif += 1;
}
output.iiexp[i] = -iie;
output.iifor[i] = -iif;
igrp = 1;
}
}
// 处理分组
if igrp == 1 {
for i in n0a..=nka {
if output.iiexp[i] > 0 {
output.iiexp[i] = -output.iiexp[i];
}
if model.imodl[i] == 0 {
model.imodl[i] = 7;
}
}
}
}
*output.nlvexp = iie.abs();
if *output.nlvexp > MLVEXP as i32 {
panic!("nlvexp.gt.mlvexp: {} > {}", *output.nlvexp, MLVEXP);
}
*output.nlvfor = iif.abs();
// 清理特定模式
for i in 0..nlevel {
let imodl = model.imodl[i];
if imodl == 1 || imodl == 3 {
output.iiexp[i] = 0;
} else if imodl == 4 || imodl == 5 {
output.iiexp[i] = 0;
} else if imodl == -1 || imodl == -3 {
output.iiexp[i] = 0;
output.iifor[i] = 0;
} else if imodl == 6 {
output.iiexp[i] = 0;
} else if imodl == -5 || imodl == -6 {
output.iiexp[i] = 0;
output.iifor[i] = 0;
} else if imodl < -100 {
model.imodl[i] = 7;
} else if imodl < -200 {
model.imodl[i] = -7;
}
// 设置 ILTREF
let nnext = model.nnext[model.iel[i] as usize];
for id in 0..nd {
output.iltref[i * nd + id] = nnext;
}
}
// 如果有分组，将所有能级设为模式 7
if igrp == 1 {
for iat in 0..natom {
if model.iifix[iat] == 1 {
continue;
}
let n0a = model.n0a[iat] as usize;
let nka = model.nka[iat] as usize;
for i in n0a..=nka {
model.imodl[i] = 7;
}
}
}
}
/// 自动处理模式
fn process_automatic(
params: &LevsetParams,
model: &mut LevsetModelState,
output: &mut LevsetOutputState,
nlevel: usize,
nd: usize,
) {
let mut iif: i32 = 0;
for i in 0..nlevel {
if model.iifix[model.iatm[i] as usize] == 1 {
continue;
}
model.imodl[i] = 5;
let nfirst = model.nfirst[model.iel[i] as usize] as usize;
let nnext = model.nnext[model.iel[i] as usize] as usize;
if i + 1 == nfirst || i + 1 == nnext {
iif += 1;
output.iifor[i] = iif;
}
}
*output.nlvexp = iif;
if *output.nlvexp > MLVEXP as i32 {
panic!("nlvexp.gt.mlvexp: {} > {}", *output.nlvexp, MLVEXP);
}
*output.nlvfor = iif;
// 清理非第一/最后能级的 IIFOR
for i in 0..nlevel {
let nfirst = model.nfirst[model.iel[i] as usize] as usize;
let nnext = model.nnext[model.iel[i] as usize] as usize;
if i + 1 != nfirst && i + 1 != nnext {
output.iifor[i] = 0;
}
}
// 设置 IIEXP 和 INDLEV
for i in 0..nlevel {
output.iiexp[i] = output.iifor[i];
if output.iiexp[i] > 0 {
output.indlev[(output.iiexp[i] - 1) as usize] = (i + 1) as i32;
}
let nnext = model.nnext[model.iel[i] as usize];
for id in 0..nd {
output.iltref[i * nd + id] = nnext;
}
}
// 非 LTE 模式下，所有能级都是显式的
if !params.lte {
for i in 0..nlevel {
output.iifor[i] = (i + 1) as i32;
}
*output.nlvfor = nlevel as i32;
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_levset_ioptab_negative() {
// IOPTAB < 0 时应直接返回
let params = LevsetParams {
lte: false,
iflev: 0,
ioptab: -1,
};
// 创建简单的模型状态
let mut imodl = vec![0i32; 10];
let iatm = vec![1i32; 10];
let iel = vec![1i32; 10];
let iifix = vec![0i32; 2];
let n0a = vec![0i32; 2];
let nka = vec![9i32; 2];
let nfirst = vec![0i32; 10];
let nnext = vec![9i32; 10];
let iltlev = vec![0i32; 10];
let mut model = LevsetModelState {
nd: 5,
natom: 2,
nlevel: 10,
imodl: &mut imodl,
iatm: &iatm,
iel: &iel,
iifix: &iifix,
n0a: &n0a,
nka: &nka,
nfirst: &nfirst,
nnext: &nnext,
iltlev: &iltlev,
};
let mut iiexp = vec![0i32; 10];
let mut iifor = vec![0i32; 10];
let mut indlev = vec![0i32; MLVEXP];
let mut iltref = vec![0i32; 10 * 5];
let mut bfac = vec![1.0; 10 * 5];
let mut igzero = vec![0i32; MLVEXP * 5];
let mut ipzero = vec![0i32; 10 * 5];
let mut nlvexp = 0i32;
let mut nlvfor = 0i32;
let mut output = LevsetOutputState {
iiexp: &mut iiexp,
iifor: &mut iifor,
indlev: &mut indlev,
iltref: &mut iltref,
bfac: &mut bfac,
igzero: &mut igzero,
ipzero: &mut ipzero,
nlvexp: &mut nlvexp,
nlvfor: &mut nlvfor,
};
levset(&params, &mut model, &mut output);
// IOPTAB < 0 时，nlvexp 应该保持为 0
assert_eq!(*output.nlvexp, 0);
}
}
src/math/mod.rs
+14
View File
@@ -1,5 +1,6 @@
//! 数学工具函数，重构自 TLUSTY Fortran。
mod alifr1;
mod alifr3;
mod alifr6;
mod alifrk;
@@ -18,6 +19,7 @@ mod collhe;
mod compt0;
mod comset;
mod cross;
mod cspec;
mod ctdata;
mod cubic;
mod dielrc;
@@ -54,6 +56,7 @@ mod irc;
mod interpolate;
mod laguer;
mod levsol;
mod levset;
mod levgrp;
mod lineqs;
mod linspl;
@@ -62,6 +65,9 @@ mod matinv;
mod meanop;
mod minv3;
mod odfhst;
mod odffr;
mod opadd0;
mod opctab;
mod pfcno;
mod pffe;
mod prdini;
@@ -76,6 +82,7 @@ mod quit;
mod reflev;
mod raph;
mod ratmal;
mod ratmat;
mod rayleigh;
mod rybmat;
mod rayset;
@@ -120,6 +127,7 @@ mod xk2dop;
mod ylintp;
mod zmrho;
pub use alifr1::{alifr1, Alifr1Params, Alifr1ModelState, Alifr1RadState};
pub use alifr3::{alifr3, Alifr3Params};
pub use alifr6::{alifr6, Alifr6Params, Alifr6State};
pub use alifrk::{alifrk, AlifrkParams, AlifrkState};
@@ -137,6 +145,7 @@ pub use ckoest::ckoest;
pub use collhe::collhe;
pub use comset::{comset, ComsetParams, ComsetResult};
pub use cross::{cross, crossd};
pub use cspec::cspec;
pub use ctdata::{hction, hctrecom, CTION, CTRECOMB};
pub use cubic::{cubic, CubicCon};
pub use dielrc::dielrc;
@@ -173,6 +182,7 @@ pub use irc::irc;
pub use interpolate::{lagran, yint};
pub use laguer::laguer;
pub use levsol::levsol;
pub use levset::{levset, LevsetParams, LevsetModelState, LevsetOutputState};
pub use levgrp::{levgrp, LevgrpParams, LevgrpResult};
pub use lineqs::{lineqs, lineqs_nr};
pub use linspl::{linspl, LinsplParams};
@@ -180,7 +190,10 @@ pub use locate::locate;
pub use matinv::matinv;
pub use meanop::meanop;
pub use minv3::minv3;
pub use opadd0::{opadd0, Opadd0Params, Opadd0FreqData, Opadd0OutputState};
pub use opctab::{opctab, OpctabParams, OpctabTableData, OpctabModelState, OpctabOutput};
pub use odfhst::odfhst;
pub use odffr::{odffr, OdffrParams, OdffrAtomicData, OdffrModelData, OdffrOutputState};
pub use pfcno::pfcno;
pub use pffe::pffe;
pub use prdini::prdini;
@@ -195,6 +208,7 @@ pub use reflev::reflev;
pub use quit::{quit, quit_error};
pub use raph::raph;
pub use ratmal::ratmal;
pub use ratmat::{ratmat, RatmatParams, RatmatOutput};
pub use rayleigh::{
rayleigh, rayleigh_h2_cross_section, rayleigh_h_cross_section, rayleigh_he_cross_section,
RayleighParams, RayleighResult,
src/math/odffr.rs
+307
View File
@@ -0,0 +1,307 @@
//! ODF 频率设置。
//!
//! 重构自 TLUSTY `odffr.f`
//!
//! 为 ODF（Opacity Distribution Function）设置内部频率，
//! 用于处理系限附近的重叠谱线。
//!
//! ## 算法说明
//!
//! 谱线向连续跃迁 (IL - IU) 的边缘收敛。
//! - IL: 下能级索引
//! - IU: 上能级索引（通常是下一离子的基态或 ODF 形式中的平均能级）
use crate::state::constants::{HALF, MFRO};
/// 常量 (从 Fortran PARAMETER 语句)
/// FRH = 3.28805D15
const FRH: f64 = 3.28805e15;
/// CDOP = 2.84511D-7
const CDOP: f64 = 2.84511e-7;
/// CDOM = 14.
const CDOM: f64 = 14.0;
/// SIX = 6.
const SIX: f64 = 6.0;
/// SEPT = 7.
const SEPT: f64 = 7.0;
/// ODFFR 输入参数
pub struct OdffrParams {
/// 下能级索引 (1-indexed)
pub il: usize,
/// 上能级索引 (1-indexed)
pub iu: usize,
/// 有效温度 (K)
pub teff: f64,
/// 最大主量子数
pub nlmx: usize,
}
/// ODFFR 原子数据
pub struct OdffrAtomicData<'a> {
/// 元素索引 (NLEVEL)
pub iel: &'a [i32],
/// 电荷数 (NATOM)
pub iz: &'a [f64],
/// 电离能 (NLEVEL)
pub enion: &'a [f64],
/// 主量子数 (NLEVEL)
pub nquant: &'a [i32],
}
/// ODFFR 模型数据
pub struct OdffrModelData<'a> {
/// 跃迁索引 (NLEVEL × NLEVEL)
pub itra: &'a [i32],
/// ODF 索引
pub jndodf: &'a [i32],
}
/// ODFFR 输出状态
pub struct OdffrOutputState<'a> {
/// ODF 频率点数 (MODF)
pub nfrodf: &'a mut [i32],
/// ODF 频率 (MFRO × MODF)
pub fros: &'a mut [f64],
/// ODF 权重 (MFRO × MODF)
pub wnus: &'a mut [f64],
}
/// 普朗克常数 (从 TLUSTY 常量)
const H: f64 = 6.626176e-27;
/// UN = 1.0 (里德伯常数的倍数)
const UN: f64 = 1.0;
/// 为 ODF 设置内部频率。
///
/// # 参数
///
/// * `params` - 输入参数
/// * `atomic` - 原子数据
/// * `model` - 模型数据
/// * `output` - 输出状态
///
/// # Panics
///
/// 当频率点数超过 MFRO 时 panic。
#[allow(clippy::too_many_arguments)]
pub fn odffr(
params: &OdffrParams,
atomic: &OdffrAtomicData,
model: &OdffrModelData,
output: &mut OdffrOutputState,
) {
let il = params.il;
let iu = params.iu;
let il_idx = il - 1;
let iu_idx = iu - 1;
// 计算电荷平方
let iel_idx = atomic.iel[il_idx] as usize;
let ch = atomic.iz[iel_idx - 1] * atomic.iz[iel_idx - 1];
let frion = ch * FRH;
// 电离频率
let fre = atomic.enion[il_idx] / H;
// 临时频率数组
let mut ffro = vec![0.0_f64; MFRO];
let mut nf: usize = 1;
let nq1 = atomic.nquant[iu_idx] as usize;
let xl2 = UN / (atomic.nquant[il_idx] as f64 * atomic.nquant[il_idx] as f64);
let xu1 = UN / ((atomic.nquant[iu_idx] - 1) as f64 * (atomic.nquant[iu_idx] - 1) as f64);
let xu2 = UN / (atomic.nquant[iu_idx] as f64 * atomic.nquant[iu_idx] as f64);
let frc = frion * (xl2 - xu2);
ffro[nf - 1] = HALF * (frc + frion * (xl2 - xu1));
let kt = model.itra[il_idx * atomic.iel.len() + iu_idx] as usize;
let kl = model.jndodf[kt - 1] as usize;
let dopo = CDOP * params.teff.sqrt() * frc;
let dopm = CDOM * dopo;
let mut fr1 = ffro[0];
// 遍历主量子数
for i in nq1..=params.nlmx {
let ii = (i * i) as f64;
let fr2 = fre - frion / ii;
let df = fr2 - fr1;
if df > dopm {
// 大间隔：添加密集点
for j in 1..=7 {
nf += 1;
if nf > MFRO {
nf -= 1;
goto_end(&ffro, nf, kl, output);
return;
}
ffro[nf - 1] = fr1 + j as f64 * dopo;
}
let df_inner = fr2 - SEPT * dopo - ffro[nf - 1];
let ni = (df_inner / (SIX * dopo)) as usize;
let ddf = df_inner / ((ni + 1) as f64);
for j in 1..=ni {
nf += 1;
if nf > MFRO - 3 {
nf -= 1;
goto_end(&ffro, nf, kl, output);
return;
}
ffro[nf - 1] = fr1 + SEPT * dopo + j as f64 * ddf;
}
for j in (0..=7).rev() {
nf += 1;
if nf > MFRO {
nf -= 1;
goto_end(&ffro, nf, kl, output);
return;
}
ffro[nf - 1] = fr2 - j as f64 * dopo;
}
fr1 = fr2;
} else {
// 小间隔：均匀分布
let ni = (df / dopo) as usize;
let ddf = df / ((ni + 1) as f64);
for j in 1..=ni {
nf += 1;
if nf > MFRO - 3 {
nf -= 1;
goto_end(&ffro, nf, kl, output);
return;
}
ffro[nf - 1] = fr1 + j as f64 * ddf;
}
nf += 1;
if nf > MFRO {
nf -= 1;
goto_end(&ffro, nf, kl, output);
return;
}
ffro[nf - 1] = fr2;
fr1 = fr2;
}
}
goto_end(&ffro, nf, kl, output);
}
/// 完成频率设置并计算权重
fn goto_end(ffro: &[f64], nf: usize, kl: usize, output: &mut OdffrOutputState) {
let mut nf = nf;
nf += 1;
if nf > MFRO {
panic!(
"too many points for hydrogen ODF - nf.gt.mfro: {} > {}",
nf, MFRO
);
}
// 添加最后一个点（略低于电离限）
let kl_idx = kl - 1;
let mut fros_local = vec![0.0; nf];
// 逆序复制
let fre = ffro[0] / 0.999999999_f64; // 从第一个点估算 fre
fros_local[nf - 1] = fre * 0.999999999;
for i in 0..nf - 1 {
fros_local[i] = ffro[nf - 2 - i];
}
output.nfrodf[kl_idx] = nf as i32;
if nf > MFRO {
panic!(
"too many points for hydrogen ODF - nf.gt.mfro: {} > {}",
nf, MFRO
);
}
// 复制到输出
for i in 0..nf {
output.fros[i * output.nfrodf.len() + kl_idx] = fros_local[i];
}
// 计算权重
// 第一个点
output.wnus[kl_idx] = HALF * (fros_local[0] - fros_local[1]);
// 最后一个点
output.wnus[(nf - 1) * output.nfrodf.len() + kl_idx] =
HALF * (fros_local[nf - 2] - fros_local[nf - 1]);
// 中间点
for i in 2..nf {
output.wnus[(i - 1) * output.nfrodf.len() + kl_idx] =
HALF * (fros_local[i - 2] - fros_local[i]);
}
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
#[test]
fn test_odffr_constants() {
// 验证常量
assert_relative_eq!(FRH, 3.28805e15, epsilon = 1e10);
assert_relative_eq!(CDOP, 2.84511e-7, epsilon = 1e-12);
assert_relative_eq!(CDOM, 14.0, epsilon = 1e-10);
}
#[test]
fn test_odffr_basic() {
// 基本测试
let params = OdffrParams {
il: 2,
iu: 10,
teff: 10000.0,
nlmx: 10,
};
let iel = vec![1, 1, 1, 1, 1, 1, 1, 1, 1, 1];
let iz = vec![1.0];
let enion = vec![13.6, 13.5, 13.4, 13.3, 13.2, 13.1, 13.0, 12.9, 12.8, 12.7];
let nquant = vec![1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
let atomic = OdffrAtomicData {
iel: &iel,
iz: &iz,
enion: &enion,
nquant: &nquant,
};
let itra = vec![1i32; 100]; // 所有跃迁索引设为 1
let jndodf = vec![1i32]; // ODF 索引
let model = OdffrModelData {
itra: &itra,
jndodf: &jndodf,
};
let mut nfrodf = vec![0i32; 10];
let mut fros = vec![0.0; MFRO * 10];
let mut wnus = vec![0.0; MFRO * 10];
let mut output = OdffrOutputState {
nfrodf: &mut nfrodf,
fros: &mut fros,
wnus: &mut wnus,
};
odffr(&params, &atomic, &model, &mut output);
// 验证输出
assert!(output.nfrodf[0] > 0);
}
}
src/math/opadd0.rs
+271
View File
@@ -0,0 +1,271 @@
//! 附加不透明度源截面设置。
//!
//! 重构自 TLUSTY `opadd0.f`
//!
//! 设置各种附加不透明度源的截面，包括：
//! - H I Rayleigh 散射
//! - He I Rayleigh 散射
//! - H2 Rayleigh 散射
//! - H- 束缚-自由和自由-自由
//! - H2+ 束缚-自由和自由-自由
//! - He- 自由-自由
use crate::math::sbfhmi::sbfhmi;
/// 常量 (从 Fortran PARAMETER 语句)
/// H I Rayleigh 散射阈值频率
const FRRAY: f64 = 2.463e15;
/// He I Rayleigh 散射阈值频率
const FRAYHE: f64 = 5.150e15;
/// H2 Rayleigh 散射阈值频率
const FRAYH2: f64 = 2.922e15;
/// 光速 (Angstrom/s)
const CLS: f64 = 2.997925e18;
/// H I Rayleigh 散射系数
const CR0: f64 = 5.799e-13;
const CR1: f64 = 1.422e-6;
const CR2: f64 = 2.784;
/// OPADD0 输入参数
pub struct Opadd0Params {
/// 频率索引 (1-indexed)
pub ij: usize,
/// 连续谱数
pub ncon: usize,
/// 最大截面数
pub mcross: usize,
/// H I Rayleigh 散射标志
pub irsct: i32,
/// He I Rayleigh 散射标志
pub irsche: i32,
/// H2 Rayleigh 散射标志
pub irsch2: i32,
/// 分子标志
pub ifmol: i32,
/// H- 不透明度标志
pub iophmi: i32,
/// H2+ 不透明度标志
pub ioph2p: i32,
/// He- 不透明度标志
pub iophem: i32,
/// ODF 标志
pub ispodf: i32,
}
/// OPADD0 频率数据
pub struct Opadd0FreqData<'a> {
/// 频率数组
pub freq: &'a [f64],
/// ODF 频率映射 (可选)
pub ifreqb: Option<&'a [i32]>,
}
/// OPADD0 输出状态
pub struct Opadd0OutputState<'a> {
/// 截面数组 (MCROSS × MFREQ)
pub bfcs: &'a mut [f32],
/// 频率维度
pub mfreq: usize,
}
/// 设置附加不透明度源的截面。
///
/// # 参数
///
/// * `params` - 输入参数
/// * `freq_data` - 频率数据
/// * `output` - 输出状态
///
/// # Panics
///
/// 当 IT > MCROSS 时 panic。
pub fn opadd0(
params: &Opadd0Params,
freq_data: &Opadd0FreqData,
output: &mut Opadd0OutputState,
) {
let ij = params.ij;
let ij_idx = ij - 1;
// 获取频率
let fr = if params.ispodf >= 1 {
let ifreqb = freq_data.ifreqb.expect("ifreqb required for ODF mode");
freq_data.freq[ifreqb[ij_idx] as usize - 1]
} else {
freq_data.freq[ij_idx]
};
let mut it = params.ncon;
// H I Rayleigh 散射
if params.irsct != 0 {
it += 1;
if it > params.mcross {
panic!("it.gt.mcross in opadd: {} > {}", it, params.mcross);
}
let frm = fr.min(FRRAY);
let x = (CLS / frm).powi(2);
let cs = (CR0 + (CR1 + CR2 / x) / x) / x / x / x;
output.bfcs[(it - 1) * output.mfreq + ij_idx] = cs as f32;
}
// He I Rayleigh 散射
if params.irsche != 0 {
it += 1;
if it > params.mcross {
panic!("it.gt.mcross in opadd: {} > {}", it, params.mcross);
}
let x = (CLS / fr.min(FRAYHE)).powi(2);
let cs = 5.484e-14 / x / x * (1.0 + (2.44e5 + 5.94e10 / (x - 2.90e5)) / x).powi(2);
output.bfcs[(it - 1) * output.mfreq + ij_idx] = cs as f32;
}
// H2 Rayleigh 散射
if params.irsch2 != 0 && params.ifmol > 0 {
it += 1;
if it > params.mcross {
panic!("it.gt.mcross in opadd: {} > {}", it, params.mcross);
}
let x = (CLS / fr.min(FRAYH2)).powi(2);
let x2 = 1.0 / x / x;
let cs = (8.14e-13 + 1.28e-6 / x + 1.61 * x2) * x2;
output.bfcs[(it - 1) * output.mfreq + ij_idx] = cs as f32;
}
// H- 束缚-自由和自由-自由
if params.iophmi > 0 {
it += 1;
if it > params.mcross {
panic!("it.gt.mcross in opadd: {} > {}", it, params.mcross);
}
let cs = sbfhmi(fr);
output.bfcs[(it - 1) * output.mfreq + ij_idx] = cs as f32;
}
// H2+ 束缚-自由和自由-自由
if params.ioph2p > 0 {
it += 1;
if it + 1 > params.mcross {
panic!("it.gt.mcross in opadd: {} > {}", it, params.mcross);
}
let x = fr * 1e-15;
// H2+ bound-free
let cs_bf = (-7.342e-3 + (-2.409 + (1.028 + (-4.23e-1 + (1.224e-1 - 1.351e-2 * x) * x) * x) * x) * x)
* 1.602e-12
/ 1.3806e-16; // BOLK
output.bfcs[(it - 1) * output.mfreq + ij_idx] = cs_bf as f32;
it += 1;
let x = fr.ln();
// H2+ free-free
let cs_ff = -3.0233e3 + (3.7797e2 + (-1.82496e1 + (3.9207e-1 - 3.1672e-3 * x) * x) * x) * x;
output.bfcs[(it - 1) * output.mfreq + ij_idx] = cs_ff as f32;
}
// He- 自由-自由
if params.iophem > 0 {
it += 1;
if it + 2 > params.mcross {
panic!("it.gt.mcross in opadd: {} > {}", it, params.mcross);
}
let a = 3.397e-46 + (-5.216e-31 + 7.039e-15 / fr) / fr;
let b = -4.116e-42 + (1.067e-26 + 8.135e-11 / fr) / fr;
let c = 5.081e-37 + (-8.724e-23 - 5.659e-8 / fr) / fr;
output.bfcs[(it - 1) * output.mfreq + ij_idx] = a as f32;
output.bfcs[it * output.mfreq + ij_idx] = b as f32;
output.bfcs[(it + 1) * output.mfreq + ij_idx] = c as f32;
}
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
#[test]
fn test_opadd0_constants() {
assert_relative_eq!(FRRAY, 2.463e15, epsilon = 1e10);
assert_relative_eq!(FRAYHE, 5.150e15, epsilon = 1e10);
assert_relative_eq!(FRAYH2, 2.922e15, epsilon = 1e10);
assert_relative_eq!(CLS, 2.997925e18, epsilon = 1e13);
}
#[test]
fn test_opadd0_hi_rayleigh() {
let params = Opadd0Params {
ij: 1,
ncon: 0,
mcross: 10,
irsct: 1,
irsche: 0,
irsch2: 0,
ifmol: 0,
iophmi: 0,
ioph2p: 0,
iophem: 0,
ispodf: 0,
};
let freq = vec![1e15];
let freq_data = Opadd0FreqData {
freq: &freq,
ifreqb: None,
};
let mut bfcs = vec![0.0f32; 10 * 1];
let mut output = Opadd0OutputState {
bfcs: &mut bfcs,
mfreq: 1,
};
opadd0(&params, &freq_data, &mut output);
// 验证 H I Rayleigh 截面被设置
assert!(output.bfcs[0] > 0.0);
}
#[test]
fn test_opadd0_all_sources() {
let params = Opadd0Params {
ij: 1,
ncon: 0,
mcross: 20,
irsct: 1,
irsche: 1,
irsch2: 1,
ifmol: 1,
iophmi: 1,
ioph2p: 1,
iophem: 1,
ispodf: 0,
};
let freq = vec![1e15];
let freq_data = Opadd0FreqData {
freq: &freq,
ifreqb: None,
};
let mut bfcs = vec![0.0f32; 20 * 1];
let mut output = Opadd0OutputState {
bfcs: &mut bfcs,
mfreq: 1,
};
opadd0(&params, &freq_data, &mut output);
// 验证所有截面都被设置
// H I Rayleigh (1)
assert!(output.bfcs[0] > 0.0);
// He I Rayleigh (2)
assert!(output.bfcs[1] > 0.0);
// H2 Rayleigh (3)
assert!(output.bfcs[2] > 0.0);
// H- (4)
assert!(output.bfcs[3] > 0.0);
// H2+ bf (5)
// H2+ ff (6)
// He- (7,8,9)
}
}
src/math/opctab.rs
+418
View File
@@ -0,0 +1,418 @@
//! 不透明度表插值与散射计算。
//!
//! 重构自 TLUSTY `opctab.f`
//!
//! 通过对预计算的不透明度表进行二维线性插值（温度-密度）来计算给定温度和密度下的吸收不透明度。
//! 同时计算散射不透明度，包括：
//! - Rayleigh 散射（简单公式或完整 rayleigh 函数）
//! - 电子散射
//!
//! 这是一个简化版本，所有插值都是线性的，fortran中也是线性插值。
use crate::math::rayleigh::{rayleigh, RayleighParams};
use crate::state::model::{EosPar, RaySct};
/// 参考频率 (FRRAY0)
const FRRAY0: f64 = 5.0872638e14;
/// OPCTAB 输入参数
pub struct OpctabParams<'a> {
/// 频率 (Hz)
pub fr: f64,
/// 频率索引 (1-indexed)
pub ij: usize,
/// 深度索引 (1-indexed)
pub id: usize,
/// 温度 (K)
pub t: f64,
/// 密度 (g cm^-3)
pub rho: f64,
/// 克拉马标志 (0: 返回每克, 1: 返回每体积)
pub igram: i32,
/// 迭代次数
pub iter: i32,
/// Rayleigh 散射标志 (<0: 使用简单公式, >0: 调用 rayleigh, =0: 关闭)
pub ifrayl: i32,
/// 选项表标志 (<0: 添加电子散射)
pub ioptab: i32,
/// Rayleigh 参数 (当 ifrayl > 0 时需要)
pub rayleigh_params: Option<&'a RayleighParams<'a>>,
}
/// OPCTAB 表数据
pub struct OpctabTableData<'a> {
/// 温度数
pub numtemp: usize,
/// 深度点数
pub nd: usize,
/// 频率数
pub nfreq: usize,
/// 最大密度点数 (用于数组维度)
pub max_numrh: usize,
/// 温度向量 (numtemp)
pub tempvec: &'a [f64],
/// 温度下限 (ln T)
pub ttab1: f64,
/// 温度上限 (ln T)
pub ttab2: f64,
/// 每个温度的密度数 (numtemp)
pub numrh: &'a [i32],
/// 密度矩阵 (numtemp × max_numrh) - 存储 ln(rho)
pub rhomat: &'a [f64],
/// 吸收不透明度表 (numtemp × max_numrh × nfreq) - 存储 ln(opacity)
pub absopac: &'a [f64],
/// Rayleigh 散射系数 (nd)
pub raysc: &'a [f64],
/// 频率数组 (nfreq)
pub freq: &'a [f64],
/// 电子散射系数
pub sige: f64,
}
/// OPCTAB 模型状态
pub struct OpctabModelState<'a> {
/// 电子密度 (nd)
pub elec: &'a [f64],
/// 密度 (nd)
pub dens: &'a [f64],
/// Rayleigh 散射截面 (当 ifrayl > 0 时需要)
pub raysct: Option<&'a mut RaySct>,
/// EOS 粒子数密度 (当 ifrayl > 0 时需要)
pub eospar: Option<&'a EosPar>,
}
/// OPCTAB 输出
pub struct OpctabOutput {
/// 吸收不透明度
pub ab: f64,
/// 散射不透明度
pub sc: f64,
/// 总散射
pub sct: f64,
}
/// 通过插值计算不透明度。
///
/// # 参数
///
/// * `params` - 输入参数
/// * `table` - 表数据
/// * `model` - 模型状态
///
/// # 返回值
///
/// 不透明度输出结构
///
/// # Panics
///
/// 当 `ifrayl > 0` 但 `rayleigh_params`、`raysct` 或 `eospar` 为 `None` 时 panic。
pub fn opctab(
params: &OpctabParams,
table: &OpctabTableData,
model: &mut OpctabModelState,
) -> OpctabOutput {
let jf = params.ij;
let id = params.id;
let id_idx = id - 1;
// 检查是否直接使用预计算的值 (numtemp == nd)
// Fortran: if(numtemp.eq.nd) then
let opac = if table.numtemp == table.nd {
// 直接使用当前深度点的值
// Fortran: opac=absopac(id,1,jf)
table.absopac[id_idx * table.max_numrh * table.nfreq + jf - 1]
} else {
compute_interpolated_opacity(params, table, jf)
};
let ab = opac.exp();
// 散射计算
let mut sct = 0.0;
// 1. Rayleigh 散射
if params.ifrayl < 0 {
// 简单公式: sct = raysc(id) * (freq(jf) / frray0)^4
sct = table.raysc[id_idx] * (table.freq[jf - 1] / FRRAY0).powi(4);
} else if params.ifrayl > 0 {
// 调用完整的 rayleigh 函数
let rayleigh_params = params
.rayleigh_params
.expect("rayleigh_params required when ifrayl > 0");
let raysct = model
.raysct
.as_mut()
.expect("raysct required when ifrayl > 0");
let eospar = model
.eospar
.as_ref()
.expect("eospar required when ifrayl > 0");
let scr = rayleigh(1, params.ij, params.id, rayleigh_params, raysct, eospar);
sct = scr / model.dens[id_idx];
}
// 电子散射
if params.ioptab < 0 {
sct = sct + table.sige * model.elec[id_idx] / params.rho;
}
// 如果迭代次数 <= 0，直接返回
if params.iter <= 0 {
return OpctabOutput { ab, sc: sct, sct };
}
// 根据 IGRAM 标志调整
let (ab_out, sc_out, sct_out) = if params.igram == 0 {
(ab * params.rho, sct * params.rho, sct * params.rho)
} else {
(ab, sct, sct)
};
OpctabOutput {
ab: ab_out,
sc: sc_out,
sct: sct_out,
}
}
/// 计算插值不透明度（内部函数）
fn compute_interpolated_opacity(params: &OpctabParams, table: &OpctabTableData, jf: usize) -> f64 {
let tl = params.t.ln();
let deltat = (tl - table.ttab1) / (table.ttab2 - table.ttab1) * (table.numtemp - 1) as f64;
let mut jt = (1.0 + deltat.floor()) as isize;
if jt < 1 {
jt = 1;
}
if jt > table.numtemp as isize - 1 {
jt = table.numtemp as isize - 1;
}
let jt = jt as usize;
let jt_idx = jt - 1;
let t1i = table.tempvec[jt_idx];
let t2i = table.tempvec[jt_idx + 1];
let mut dti = (tl - t1i) / (t2i - t1i);
if deltat < 0.0 {
dti = 0.0;
}
// 检查是否需要密度插值
// Fortran: if(numrh(1).ne.1) then
if table.numrh[0] != 1 {
compute_2d_interpolation(params, table, jf, jt_idx, dti)
} else {
// 单密度情况: jr = 1 (Fortran), 所以 jr_idx = 0
// Fortran: opac=absopac(jt,jr,jf)+(absopac(ju,jr,jf)-absopac(jt,jr,jf))*dti
let jr_idx = 0; // jr = 1 in Fortran
let idx_jt = jt_idx * table.max_numrh * table.nfreq + jr_idx * table.nfreq + (jf - 1);
let idx_ju = (jt_idx + 1) * table.max_numrh * table.nfreq + jr_idx * table.nfreq + (jf - 1);
table.absopac[idx_jt] + (table.absopac[idx_ju] - table.absopac[idx_jt]) * dti
}
}
/// 二维插值（温度和密度）
fn compute_2d_interpolation(
params: &OpctabParams,
table: &OpctabTableData,
jf: usize,
jt_idx: usize,
dti: f64,
) -> f64 {
let rl = params.rho.ln();
// 低温下的密度插值
// Fortran: numrho=numrh(jt), rtab1=rhomat(jt,1), rtab2=rhomat(jt,numrho)
let numrho = table.numrh[jt_idx] as usize;
let rtab1 = table.rhomat[jt_idx * table.max_numrh];
let rtab2 = table.rhomat[jt_idx * table.max_numrh + numrho - 1];
let deltar = (rl - rtab1) / (rtab2 - rtab1) * (numrho - 1) as f64;
let mut jr = (1.0 + deltar.floor()) as isize;
if jr < 1 {
jr = 1;
}
if jr > numrho as isize - 1 {
jr = numrho as isize - 1;
}
let jr = jr as usize;
let jr_idx = jr - 1;
// Fortran: r1i=rhomat(jt,jr), r2i=rhomat(jt,jr+1)
let r1i = table.rhomat[jt_idx * table.max_numrh + jr_idx];
let r2i = table.rhomat[jt_idx * table.max_numrh + jr_idx + 1];
let mut dri = (rl - r1i) / (r2i - r1i);
if deltar < 0.0 {
dri = 0.0;
}
// Fortran: opr1=absopac(jt,jr,jf)+dri*(absopac(jt,jr+1,jf)-absopac(jt,jr,jf))
let idx_jt_jr = jt_idx * table.max_numrh * table.nfreq + jr_idx * table.nfreq + (jf - 1);
let idx_jt_jrp1 = jt_idx * table.max_numrh * table.nfreq + (jr_idx + 1) * table.nfreq + (jf - 1);
let opr1 =
table.absopac[idx_jt_jr] + dri * (table.absopac[idx_jt_jrp1] - table.absopac[idx_jt_jr]);
// 高温下的密度插值
// Fortran: ju=jt+1, numrho=numrh(ju)
let ju_idx = jt_idx + 1;
let numrho_h = table.numrh[ju_idx] as usize;
let rtab1_h = table.rhomat[ju_idx * table.max_numrh];
let rtab2_h = table.rhomat[ju_idx * table.max_numrh + numrho_h - 1];
let deltar_h = (rl - rtab1_h) / (rtab2_h - rtab1_h) * (numrho_h - 1) as f64;
let mut jr_h = (1.0 + deltar_h.floor()) as isize;
if jr_h < 1 {
jr_h = 1;
}
if jr_h > numrho_h as isize - 1 {
jr_h = numrho_h as isize - 1;
}
let jr_h = jr_h as usize;
let jr_h_idx = jr_h - 1;
let r1i_h = table.rhomat[ju_idx * table.max_numrh + jr_h_idx];
let r2i_h = table.rhomat[ju_idx * table.max_numrh + jr_h_idx + 1];
let mut dri_h = (rl - r1i_h) / (r2i_h - r1i_h);
if deltar_h < 0.0 {
dri_h = 0.0;
}
// Fortran: opr2=absopac(ju,jr,jf)+dri*(absopac(ju,jr+1,jf)-absopac(ju,jr,jf))
let idx_ju_jr = ju_idx * table.max_numrh * table.nfreq + jr_h_idx * table.nfreq + (jf - 1);
let idx_ju_jrp1 = ju_idx * table.max_numrh * table.nfreq + (jr_h_idx + 1) * table.nfreq + (jf - 1);
let opr2 = table.absopac[idx_ju_jr]
+ dri_h * (table.absopac[idx_ju_jrp1] - table.absopac[idx_ju_jr]);
// Fortran: opac=opr1+(opr2-opr1)*dti
opr1 + (opr2 - opr1) * dti
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
#[test]
fn test_opctab_constants() {
assert_relative_eq!(FRRAY0, 5.0872638e14, epsilon = 1e9);
}
#[test]
fn test_opctab_direct_value() {
let params = OpctabParams {
fr: 1e15,
ij: 1,
id: 1,
t: 10000.0,
rho: 1e-7,
igram: 0,
iter: 1,
ifrayl: -1,
ioptab: 0,
rayleigh_params: None,
};
let numtemp = 1;
let nd = 1;
let nfreq = 1;
let max_numrh = 1;
let tempvec = vec![9.2103];
let numrh = vec![1];
let rhomat = vec![-16.1181];
let absopac = vec![0.0];
let raysc = vec![1e-20];
let freq = vec![5.0872638e14];
let table = OpctabTableData {
numtemp,
nd,
nfreq,
max_numrh,
tempvec: &tempvec,
ttab1: 8.0,
ttab2: 11.0,
numrh: &numrh,
rhomat: &rhomat,
absopac: &absopac,
raysc: &raysc,
freq: &freq,
sige: 6.65e-25,
};
let elec = vec![1e-10];
let dens = vec![1e-7];
let mut model = OpctabModelState {
elec: &elec,
dens: &dens,
raysct: None,
eospar: None,
};
let result = opctab(&params, &table, &mut model);
// 当 numtemp == nd 时，直接使用 absopac(id, 1, jf) = 0.0
// ab = exp(0.0) * rho = 1.0 * 1e-7 = 1e-7
assert_relative_eq!(result.ab, 1e-7, epsilon = 1e-10);
}
#[test]
fn test_opctab_rayleigh_scattering() {
let params = OpctabParams {
fr: 1e15,
ij: 1,
id: 1,
t: 10000.0,
rho: 1e-7,
igram: 0,
iter: 1,
ifrayl: -1,
ioptab: 0,
rayleigh_params: None,
};
let numtemp = 1;
let nd = 1;
let nfreq = 1;
let max_numrh = 1;
let tempvec = vec![9.2103];
let numrh = vec![1];
let rhomat = vec![-16.1181];
let absopac = vec![0.0];
let raysc = vec![1e-20];
let freq = vec![5.0872638e14];
let table = OpctabTableData {
numtemp,
nd,
nfreq,
max_numrh,
tempvec: &tempvec,
ttab1: 8.0,
ttab2: 11.0,
numrh: &numrh,
rhomat: &rhomat,
absopac: &absopac,
raysc: &raysc,
freq: &freq,
sige: 6.65e-25,
};
let elec = vec![1e-10];
let dens = vec![1e-7];
let mut model = OpctabModelState {
elec: &elec,
dens: &dens,
raysct: None,
eospar: None,
};
let result = opctab(&params, &table, &mut model);
// Rayleigh 散射 = raysc * (freq/frray0)^4 * rho
// = 1e-20 * 1.0^4 * 1e-7 = 1e-27
assert_relative_eq!(result.sct, 1e-27, epsilon = 1e-30);
}
}
src/math/ratmat.rs
+369
View File
@@ -0,0 +1,369 @@
//! 速率矩阵计算 - RATMAT。
//!
//! 重构自 TLUSTY `ratmat.f`
//!
//! 计算统计平衡方程的速率矩阵 A 和右端向量 B。
//!
//! 方程形式为：A * n = B，其中 n 是能级占据数向量。
//!
//! 参考：Mihalas, 1978, pp.138-139
use crate::math::reflev::reflev;
use crate::state::atomic::AtomicData;
use crate::state::config::TlustyConfig;
use crate::state::constants::{HK, MLEVEL, UN};
use crate::state::iterat::IterControl;
use crate::state::model::ModelState;
/// RATMAT 输入参数
pub struct RatmatParams<'a> {
/// 深度索引 (1-indexed)
pub id: usize,
/// 能级索引数组 (可能为负)
pub iical: &'a mut [i32],
/// 模式标志 (<0: 特殊模式, 0: 正常, >0: 调用 REFLEV)
pub imode: i32,
}
/// RATMAT 输出
pub struct RatmatOutput {
/// 速率矩阵 A (MLEVEL × MLEVEL)
pub a: Vec<Vec<f64>>,
/// 右端向量 B (MLEVEL)
pub b: Vec<f64>,
}
/// 计算速率矩阵。
///
/// # 参数
///
/// * `params` - 输入参数
/// * `config` - 配置
/// * `atomic` - 原子数据
/// * `model` - 模型状态
/// * `iterat` - 迭代控制
///
/// # 返回值
///
/// 速率矩阵 A 和右端向量 B
pub fn ratmat(
params: &mut RatmatParams,
config: &mut TlustyConfig,
atomic: &mut AtomicData,
model: &mut ModelState,
iterat: &IterControl,
) -> RatmatOutput {
let id = params.id;
let id_idx = id - 1;
// 检查选项表标志
if config.basnum.ioptab < 0 {
return RatmatOutput {
a: vec![vec![0.0; MLEVEL]; MLEVEL],
b: vec![0.0; MLEVEL],
};
}
let nlevel = config.basnum.nlevel as usize;
let ntrans = config.basnum.ntrans as usize;
let natom = config.basnum.natom as usize;
let lte = config.inppar.lte;
let ipslte = config.inppar.ipslte;
// 如果使用 Slater 迭代，清零辐射速率
if ipslte != 0 {
for itr in 0..ntrans {
model.rrrates.rru[itr][id_idx] = 0.0;
model.rrrates.rrd[itr][id_idx] = 0.0;
model.rrrates.drdt[itr][id_idx] = 0.0;
}
}
// 初始化输出
let mut a = vec![vec![0.0; MLEVEL]; MLEVEL];
let mut b = vec![0.0; MLEVEL];
// 工作数组
let mut aij = vec![0.0; ntrans.max(1)];
let mut aji = vec![0.0; ntrans.max(1)];
let mut sbw = vec![0.0; MLEVEL];
let mut llte = vec![false; MLEVEL];
// 获取物理量
let t = model.modpar.temp[id_idx];
let ane = model.modpar.elec[id_idx];
let hkt = HK / t;
let tk = hkt / HK;
// 初始化
for i in 0..nlevel {
b[i] = 0.0;
sbw[i] = ane * model.levpop.sbf[i] * model.wmcomp.wop[i][id_idx];
llte[i] = false;
// 设置 LTE 标志
let iel_i = atomic.levpar.iel[i] as usize;
let ilt = atomic.ionpar.iltion[iel_i];
let llt = ilt == 1 && params.imode == 0;
llte[i] = llt || lte || atomic.ionpar.iltlev[i] >= 1 || ilt >= 2;
llte[i] = llte[i] || id >= config.inppar.idlte as usize;
for j in 0..nlevel {
a[j][i] = 0.0;
}
}
// 确定参考能级
if params.imode != 0 {
reflev(id, params.imode.abs(), config, atomic, model, iterat);
}
// 第一部分：LTE 情况下的简化表达式
for iat in 0..natom {
if atomic.atopar.iifix[iat] != 1 || params.imode < 0 {
let nrefi = atomic.atopar.nrefs[iat][id_idx] as usize;
let n0 = atomic.atopar.n0a[iat] as usize;
let nk = atomic.atopar.nka[iat] as usize;
for i in n0..=nk {
let ii = params.iical[i].abs() as usize;
if i != nrefi && ii != 0 && llte[i] {
a[ii][ii] = UN;
let n = params.iical[model.levref.ilterf[i][id_idx] as usize].abs() as usize;
a[ii][n] = a[ii][n] - model.levref.sblpsi[i][id_idx];
}
}
}
}
// 第二部分：NLTE 速率方程
if !lte {
// 计算跃迁速率
for itr in 0..ntrans {
let i = atomic.trapar.ilow[itr] as usize;
if atomic.atopar.iifix[atomic.trapar.iatm[i] as usize] == 1 {
continue;
}
let j = atomic.trapar.iup[itr] as usize;
let nke = atomic.ionpar.nnext[atomic.levpar.iel[i] as usize] as usize;
// 向上总速率
aij[itr] = (model.rrrates.colrat[itr][id_idx] + model.rrrates.rru[itr][id_idx])
* model.wmcomp.wop[j][id_idx];
// 向下总速率
if atomic.trapar.line[itr] {
// 束缚-束缚跃迁
aji[itr] = (model.rrrates.coltar[itr][id_idx]
+ model.rrrates.rrd[itr][id_idx]
* atomic.levpar.g[i] / atomic.levpar.g[j]
* (hkt * atomic.trapar.fr0[itr]).exp())
* model.wmcomp.wop[i][id_idx];
} else {
// 束缚-自由跃迁
let corr = if nke != j {
atomic.levpar.g[nke] / atomic.levpar.g[j]
* ((atomic.levpar.enion[nke] - atomic.levpar.enion[j]) * tk).exp()
} else {
UN
};
aji[itr] = model.rrrates.coltar[itr][id_idx] * model.wmcomp.wop[i][id_idx]
+ model.rrrates.rrd[itr][id_idx] * sbw[i] * corr;
}
// 处理负索引
if params.iical[i] < 0 && params.iical[j] < 0 {
aij[itr] *= model.levref.sbpsi[i][id_idx];
aji[itr] *= model.levref.sbpsi[j][id_idx];
}
}
// 填充速率矩阵
for itr in 0..ntrans {
let i = atomic.trapar.ilow[itr] as usize;
if atomic.atopar.iifix[atomic.trapar.iatm[i] as usize] == 1 {
continue;
}
let nrefi = atomic.atopar.nrefs[atomic.trapar.iatm[i] as usize][id_idx] as usize;
let j = atomic.trapar.iup[itr] as usize;
let ii = params.iical[i].abs() as usize;
let jj = params.iical[j].abs() as usize;
if model.popzr0.ipzero[i][id_idx] != 0 || model.popzr0.ipzero[j][id_idx] != 0 {
continue;
}
// 下能级贡献
if i != nrefi && ii > 0 && !llte[i] {
a[ii][ii] += aij[itr];
if jj > 0 {
a[ii][jj] -= aji[itr];
} else {
let jjj = params.iical[model.levref.iltref[j][id_idx] as usize] as usize;
a[ii][jjj] -= aji[itr] * model.levref.sbpsi[j][id_idx];
}
}
// 上能级贡献
if j != nrefi && jj > 0 && !llte[j] {
a[jj][jj] += aji[itr];
if ii > 0 {
a[jj][ii] -= aij[itr];
} else {
let iii = params.iical[model.levref.iltref[i][id_idx] as usize] as usize;
a[jj][iii] -= aij[itr] * model.levref.sbpsi[i][id_idx];
}
}
}
}
// 重置"小"占据数的速率矩阵元素
for i in 0..nlevel {
let ii = params.iical[i];
if ii > 0 {
let ii_idx = ii as usize;
if model.popzr0.ipzero[i][id_idx] > 0 {
for j in 0..nlevel {
a[ii_idx][j] = 0.0;
}
a[ii_idx][ii_idx] = 1.0;
}
} else if ii < 0 {
let ii_idx = (-ii) as usize;
if model.popzr0.igzero[ii_idx - 1][id_idx] > 0 {
for j in 0..nlevel {
a[ii_idx - 1][j] = 0.0;
}
a[ii_idx - 1][ii_idx - 1] = 1.0;
}
}
}
// 第三部分：丰度定义方程
for iat in 0..natom {
if atomic.atopar.iifix[iat] == 1 && params.imode >= 0 {
continue;
}
let nrefii = params.iical[atomic.atopar.nrefs[iat][id_idx] as usize].abs() as usize;
let n0 = atomic.atopar.n0a[iat] as usize;
let nk = atomic.atopar.nka[iat] as usize;
for i in n0..=nk {
let il = atomic.levpar.ilk[i];
let ii = params.iical[i];
if ii > 0 {
let ii_idx = ii as usize;
a[nrefii][ii_idx] += UN;
if il != 0 {
a[nrefii][ii_idx] += ane * model.levpop.usum[il as usize];
}
} else if ii < 0 {
let ii_idx = (-ii) as usize;
a[nrefii][ii_idx - 1] += model.levref.sbpsi[i][id_idx];
if il != 0 {
a[nrefii][ii_idx - 1] +=
ane * model.levpop.usum[il as usize] * model.levref.sbpsi[i][id_idx];
}
} else {
let iii = params.iical[model.levref.iltref[i][id_idx] as usize] as usize;
if iii > 0 && iii <= nlevel {
a[nrefii][iii - 1] += model.levref.sbpsi[i][id_idx];
}
}
}
b[nrefii] += model.modpar.dens[id_idx]
/ model.modpar.wmm[id_idx]
/ model.modpar.ytot[id_idx]
* atomic.atopar.abund[iat][id_idx];
}
RatmatOutput { a, b }
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_ratmat_ioptab_negative() {
let mut config = TlustyConfig::default();
let mut atomic = AtomicData::default();
let mut model = ModelState::default();
let iterat = IterControl::default();
config.basnum.ioptab = -1;
config.basnum.nlevel = 10;
config.basnum.ntrans = 5;
config.basnum.natom = 1;
let mut iical = vec![0; MLEVEL];
let mut params = RatmatParams {
id: 1,
iical: &mut iical,
imode: 0,
};
let result = ratmat(&mut params, &mut config, &mut atomic, &mut model, &iterat);
// 当 ioptab < 0 时，应返回零矩阵
assert_eq!(result.a.len(), MLEVEL);
assert_eq!(result.b.len(), MLEVEL);
for i in 0..MLEVEL {
assert_relative_eq!(result.b[i], 0.0);
}
}
#[test]
fn test_ratmat_basic_lte() {
let mut config = TlustyConfig::default();
let mut atomic = AtomicData::default();
let mut model = ModelState::default();
let iterat = IterControl::default();
config.basnum.ioptab = 0;
config.basnum.nlevel = 5;
config.basnum.ntrans = 3;
config.basnum.natom = 1;
config.inppar.lte = true;
atomic.atopar.n0a[0] = 0;
atomic.atopar.nka[0] = 4;
atomic.atopar.nrefs[0][0] = 2;
atomic.atopar.iifix[0] = 0;
atomic.levpar.iel[0] = 0;
atomic.levpar.iel[1] = 0;
atomic.levpar.iel[2] = 1;
atomic.levpar.iel[3] = 1;
atomic.levpar.iel[4] = 1;
atomic.levpar.ilk[0] = 0;
atomic.levpar.ilk[1] = 0;
atomic.levpar.ilk[2] = 1;
atomic.levpar.ilk[3] = 0;
atomic.levpar.ilk[4] = 0;
atomic.ionpar.iltion[0] = 1;
atomic.ionpar.iltion[1] = 1;
model.modpar.temp[0] = 10000.0;
model.modpar.elec[0] = 1e12;
model.modpar.dens[0] = 1e14;
model.modpar.wmm[0] = 1.0;
model.modpar.ytot[0] = 1.0;
atomic.atopar.abund[0][0] = 1.0;
let mut iical = vec![1, 2, 3, 4, 5];
let mut params = RatmatParams {
id: 1,
iical: &mut iical,
imode: 0,
};
let result = ratmat(&mut params, &mut config, &mut atomic, &mut model, &iterat);
// 验证矩阵维度
assert_eq!(result.a.len(), MLEVEL);
assert_eq!(result.b.len(), MLEVEL);
}
}