重构5,无io和无io依赖的模块已经全部重构完毕,接下来是重构剩余的模块,主要是io和依赖io的模块。

This commit is contained in:
fmq
2026-03-22 17:08:40 +08:00
parent 59404207c3
commit 5078d6120e
35 changed files with 11554 additions and 16 deletions
+421
View File
@@ -0,0 +1,421 @@
//! B 矩阵的占据数行和显式频率列部分。
//!
//! 重构自 TLUSTY `bpope.f`
//!
//! 处理完全重叠情况下的 B 矩阵元素。
use crate::state::constants::{MFREX, MLEVEL, MLVEXP, UN};
/// BPOPE 输入参数
pub struct BpopeParams {
/// 深度索引 (1-indexed)
pub id: usize,
}
/// BPOPE 配置参数
pub struct BpopeConfig {
/// 显式频率点数
pub nfreqe: usize,
/// 频率点数
pub nfreq: usize,
/// 连续谱跃迁数
pub ntranc: usize,
/// 显式能级数
pub nlvexp: usize,
/// INSE 索引偏移
pub inse: usize,
/// ODF 采样标志 (0: 不使用 ODF)
pub ispodf: i32,
/// 人口行处理标志
pub ifpopr: i32,
/// CRSW 系数
pub crsw: f64,
}
/// BPOPE 原子数据
pub struct BpopeAtomicData<'a> {
/// 跃迁的能级索引 (ntrans)
pub ilow: &'a [i32],
/// 跃迁的上能级索引 (ntrans)
pub iup: &'a [i32],
/// 连续谱跃迁索引 (ntranc)
pub itrbf: &'a [i32],
/// 跃迁的频率 (ntrans)
pub fr0: &'a [f64],
/// MCDW 标志 (ntrans)
pub mcdw: &'a [i32],
/// 谱线是否显式 (ntrans)
pub linexp: &'a [bool],
/// LEXP 标志 (ntrans)
pub lexp: &'a [bool],
/// 元素索引 (nlevel)
pub iel: &'a [i32],
/// 原子索引 (nlevel)
pub iatm: &'a [i32],
/// 能级是否显式 (nlevel)
pub iiexp: &'a [i32],
/// 能级的 LTE 标志 (nlevel)
pub iltlev: &'a [i32],
/// IMODL 标志 (nlevel)
pub imodl: &'a [i32],
/// IMRG 标志 (nlevel)
pub imrg: &'a [i32],
/// 电离阶段 (nelem)
pub iltion: &'a [i32],
/// 固定原子标志 (natom)
pub iifix: &'a [i32],
/// 原子核电荷 (nelem)
pub iz: &'a [i32],
}
/// BPOPE 模型状态
pub struct BpopeModelState<'a> {
/// 温度 (nd)
pub temp: &'a [f64],
/// HKT1 数组 (nd)
pub hkt1: &'a [f64],
/// 参考能级索引 (natom × nd)
pub nrefs: &'a [i32],
/// 零占据数标志 (nlevel × nd)
pub ipzero: &'a [i32],
/// 吸收系数 (ntrans × nd)
pub abtra: &'a [f64],
/// 发射系数 (ntrans × nd)
pub emtra: &'a [f64],
}
/// BPOPE 频率数据
pub struct BpopeFreqData<'a> {
/// 频率数组 (nfreq)
pub freq: &'a [f64],
/// 显式频率索引 (nfreq)
pub ijex: &'a [i32],
/// 显式频率映射 (nfreqe)
pub ijfr: &'a [i32],
/// IJX 标志 (nfreq)
pub ijx: &'a [i32],
/// 谱线索引 (nfreq)
pub ijlin: &'a [i32],
/// 重叠谱线数 (nfreq)
pub nlines: &'a [i32],
/// 重叠谱线索引 (nliness × nfreq)
pub itrlin: &'a [i32],
/// 权重 (nfreq)
pub w0e: &'a [f64],
/// 跃迁起始频率索引 (ntrans)
pub ifr0: &'a [i32],
/// 跃迁结束频率索引 (ntrans)
pub ifr1: &'a [i32],
/// KFR0 索引 (ntrans)
pub kfr0: &'a [i32],
/// 谱线轮廓 (nd × nfreq 或 nd × nfro)
pub prflin: &'a [f64],
/// 截面 (ntranc × nfreq)
pub cross: &'a [f64],
}
/// BPOPE 矩阵数据
pub struct BpopeMatrixData<'a> {
/// ESE 矩阵 (nlvexp × nlvexp)
pub esemat: &'a [f64],
/// APT 数组 (nlvexp × nd)
pub apt: &'a [f64],
}
/// BPOPE 输出
pub struct BpopeOutput {
/// B 矩阵元素 (nlvexp × nfreqe)
pub b: Vec<Vec<f64>>,
}
/// 计算 B 矩阵的占据数行和显式频率列部分。
///
/// # 参数
///
/// * `params` - 输入参数
/// * `config` - 配置参数
/// * `atomic` - 原子数据
/// * `model` - 模型状态
/// * `freq_data` - 频率数据
/// * `matrix_data` - 矩阵数据
///
/// # 返回值
///
/// B 矩阵元素
pub fn bpope(
params: &BpopeParams,
config: &BpopeConfig,
atomic: &BpopeAtomicData,
model: &BpopeModelState,
freq_data: &BpopeFreqData,
matrix_data: &BpopeMatrixData,
) -> BpopeOutput {
let id = params.id;
let id_idx = id - 1;
// 如果没有显式频率点,直接返回
if config.nfreqe <= 0 {
return BpopeOutput {
b: vec![vec![0.0; config.nfreqe]; config.nlvexp],
};
}
let nse = config.nfreqe + config.inse - 1;
let hk = 4.1356692e-16; // Planck 常数 (eV·s),需要从常量获取
// 初始化 AJIJ 数组
let mut ajij = vec![vec![0.0; config.nlvexp]; MFREX];
let mut ehke = vec![0.0; MFREX];
let hkt = hk / model.temp[id_idx];
// 计算 EHKE
for ije in 0..config.nfreqe {
let ij = freq_data.ijfr[ije] as usize - 1;
ehke[ije] = (-model.hkt1[id_idx] * freq_data.freq[ij]).exp();
}
// 遍历所有频率点
for ij in 0..config.nfreq {
if freq_data.ijex[ij] <= 0 || freq_data.ijx[ij] == -1 {
continue;
}
let ije = (freq_data.ijex[ij] - 1) as usize;
let fr = freq_data.freq[ij];
let frinv = UN / fr;
let fr3inv = frinv * frinv * frinv;
// 处理连续谱跃迁
for ibft in 0..config.ntranc {
let itr = atomic.itrbf[ibft] as usize - 1;
let sg = freq_data.cross[ibft * config.nfreq + ij];
if sg <= 0.0 {
continue;
}
let i = atomic.ilow[itr] as usize - 1;
let iel_i = atomic.iel[i] as usize;
if atomic.iltion[iel_i] >= 1 || atomic.iifix[atomic.iatm[i] as usize] == 1 {
continue;
}
let ii = atomic.iiexp[i].abs() as usize;
let j = atomic.iup[itr] as usize - 1;
if model.ipzero[i * id + id_idx] != 0 || model.ipzero[j * id + id_idx] != 0 {
continue;
}
let jj = atomic.iiexp[j].abs() as usize;
let nrefi = model.nrefs[atomic.iatm[i] as usize * id + id_idx];
// 简化处理:直接使用 sg
let sg_final = sg;
let w0 = freq_data.w0e[ij];
let sgw0 = sg_final * w0;
let apfr = (model.abtra[itr * id + id_idx]
- model.emtra[itr * id + id_idx] * ehke[ije])
* sgw0;
if ii > 0
&& (i + 1) != nrefi as usize
&& atomic.iltlev[i] <= 0
{
ajij[ije][ii - 1] += apfr;
}
if jj > 0
&& (j + 1) != nrefi as usize
&& atomic.iltlev[j] <= 0
&& atomic.imodl[i].abs() != 4
{
ajij[ije][jj - 1] -= apfr;
}
}
// 处理谱线跃迁(简化版本,不处理 ODF 采样)
if config.ispodf == 0 && freq_data.ijlin[ij] > 0 {
let itr = (freq_data.ijlin[ij] - 1) as usize;
if !atomic.linexp[itr] && atomic.lexp[itr] {
let i = atomic.ilow[itr] as usize - 1;
let iel_i = atomic.iel[i] as usize;
if atomic.iltion[iel_i] >= 1 || atomic.iifix[atomic.iatm[i] as usize] == 1 {
continue;
}
let j = atomic.iup[itr] as usize - 1;
if model.ipzero[i * id + id_idx] != 0
|| model.ipzero[j * id + id_idx] != 0
{
continue;
}
let ii = atomic.iiexp[i].abs() as usize;
let jj = atomic.iiexp[j].abs() as usize;
if ii == 0 && jj == 0 {
continue;
}
let nrefi = model.nrefs[atomic.iatm[i] as usize * id + id_idx];
let sgw = freq_data.prflin[id_idx * config.nfreq + ij] * freq_data.w0e[ij];
let apfr = (model.abtra[itr * id + id_idx]
- model.emtra[itr * id + id_idx] * ehke[ije])
* sgw;
if ii > 0
&& (i + 1) != nrefi as usize
&& atomic.iltlev[i] <= 0
{
ajij[ije][ii - 1] += apfr;
}
if jj > 0
&& (j + 1) != nrefi as usize
&& atomic.iltlev[j] <= 0
&& atomic.imodl[i].abs() != 4
{
ajij[ije][jj - 1] -= apfr;
}
}
}
}
// 计算 B 矩阵元素
let mut b = vec![vec![0.0; config.nfreqe]; config.nlvexp];
for i in 0..config.nlvexp {
for ije in 0..config.nfreqe {
let sum = if config.ifpopr <= 3 {
let mut s = 0.0;
for j in 0..config.nlvexp {
s -= matrix_data.esemat[i * config.nlvexp + j] * ajij[ije][j];
}
s
} else {
ajij[ije][i]
};
b[i][ije] = sum * config.crsw;
}
}
BpopeOutput { b }
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_bpope_no_explicit_freq() {
// 当 nfreqe = 0 时,应返回零矩阵
let params = BpopeParams { id: 1 };
let config = BpopeConfig {
nfreqe: 0,
nfreq: 100,
ntranc: 10,
nlvexp: 5,
inse: 1,
ispodf: 0,
ifpopr: 3,
crsw: 1.0,
};
let ilow = vec![1; 10];
let iup = vec![2; 10];
let itrbf = vec![1; 10];
let fr0 = vec![1e15; 10];
let mcdw = vec![0; 10];
let linexp = vec![false; 10];
let lexp = vec![true; 10];
let iel = vec![0; 100];
let iatm = vec![0; 100];
let iiexp = vec![1; 100];
let iltlev = vec![0; 100];
let imodl = vec![0; 100];
let imrg = vec![0; 100];
let iltion = vec![0; 10];
let iifix = vec![0; 10];
let iz = vec![1; 10];
let atomic = BpopeAtomicData {
ilow: &ilow,
iup: &iup,
itrbf: &itrbf,
fr0: &fr0,
mcdw: &mcdw,
linexp: &linexp,
lexp: &lexp,
iel: &iel,
iatm: &iatm,
iiexp: &iiexp,
iltlev: &iltlev,
imodl: &imodl,
imrg: &imrg,
iltion: &iltion,
iifix: &iifix,
iz: &iz,
};
let temp = vec![10000.0; 10];
let hkt1 = vec![1e-18; 10];
let nrefs = vec![1; 100];
let ipzero = vec![0; 1000];
let abtra = vec![1e-10; 100];
let emtra = vec![1e-10; 100];
let model = BpopeModelState {
temp: &temp,
hkt1: &hkt1,
nrefs: &nrefs,
ipzero: &ipzero,
abtra: &abtra,
emtra: &emtra,
};
let freq = vec![1e15; 100];
let ijex = vec![0; 100];
let ijfr = vec![0; 100];
let ijx = vec![0; 100];
let ijlin = vec![0; 100];
let nlines = vec![0; 100];
let itrlin = vec![0; 1000];
let w0e = vec![1.0; 100];
let ifr0 = vec![1; 100];
let ifr1 = vec![10; 100];
let kfr0 = vec![0; 100];
let prflin = vec![1.0; 1000];
let cross = vec![1e-18; 1000];
let freq_data = BpopeFreqData {
freq: &freq,
ijex: &ijex,
ijfr: &ijfr,
ijx: &ijx,
ijlin: &ijlin,
nlines: &nlines,
itrlin: &itrlin,
w0e: &w0e,
ifr0: &ifr0,
ifr1: &ifr1,
kfr0: &kfr0,
prflin: &prflin,
cross: &cross,
};
let esemat = vec![0.0; 25];
let apt = vec![0.0; 50];
let matrix_data = BpopeMatrixData {
esemat: &esemat,
apt: &apt,
};
let result = bpope(&params, &config, &atomic, &model, &freq_data, &matrix_data);
// 结果应该是 5×0 的空矩阵
assert_eq!(result.b.len(), 5);
assert_eq!(result.b[0].len(), 0);
}
}
+752
View File
@@ -0,0 +1,752 @@
//! 辐射平衡方程矩阵计算。
//!
//! 重构自 TLUSTY `bre.f`
//!
//! 计算辐射平衡方程对应的矩阵 A, B, C 部分(第 NFREQE+INRE 行)。
//! 包含积分方程部分和微分方程部分。
use crate::state::constants::{HALF, SIG4P, UN};
// ============================================================================
// 常量
// ============================================================================
/// 最大线性化能级数
const MLVEXP: usize = 233;
// ============================================================================
// BRE 参数结构体
// ============================================================================
/// BRE 输入参数。
///
/// 包含所有来自 COMMON 块的必要数据。
#[derive(Debug)]
pub struct BreParams<'a> {
// ==================== 基本参数 ====================
/// 深度索引 ID (1-indexed)
pub id: usize,
// ==================== 维度参数 ====================
/// 深度点数 ND
pub nd: usize,
/// 线性化频率数 NFREQE
pub nfreqe: usize,
/// 总频率数 NFREQ
pub nfreq: usize,
/// 线性化能级数 NLVEXP
pub nlvexp: usize,
// ==================== 索引参数 ====================
/// 氦索引 INHE (0 表示无氦)
pub inhe: usize,
/// 辐射平衡索引 INRE
pub inre: usize,
/// 电子密度索引 INPC
pub inpc: usize,
/// 质量密度索引 INMP
pub inmp: usize,
// ==================== 控制参数 ====================
/// Compton 散射标志
pub icompt: i32,
/// Compton 边界条件标志
pub icombc: i32,
/// Compton 密度导数标志
pub icmdra: i32,
/// 迭代次数
pub iter: i32,
/// 辐射平衡温度控制
pub nretc: i32,
/// 盘模型标志 (0=无盘, 1=有盘)
pub idisk: i32,
/// ALI 标志 (>5 启用完整 ALI)
pub ifali: i32,
// ==================== 频率映射 ====================
/// 频率索引映射 [MFREQE]
pub ijfr: &'a [usize],
/// 频率起源索引 [MFREQ]
pub ijorig: &'a [usize],
/// 频率反转映射 [MFREQ]
pub kij: &'a [usize],
/// ALI 频率索引 [MFREQ]
pub ijex: &'a [i32],
// ==================== 模型状态 ====================
/// 温度 [MDEPTH] (K)
pub temp: &'a [f64],
/// 电子密度 [MDEPTH]
pub elec: &'a [f64],
/// 柱质量密度 [MDEPTH]
pub dm: &'a [f64],
/// 密度倒数 [MDEPTH]
pub dens1: &'a [f64],
/// 分子权重 [MDEPTH]
pub wmm: &'a [f64],
/// 粘性加热 [MDEPTH]
pub tvisc: &'a [f64],
/// 辐射积分因子 [MDEPTH]
pub reint: &'a [f64],
/// THETAV 速度场 [MDEPTH]
pub thetav: &'a [f64],
// ==================== 辐射场 ====================
/// 当前深度辐射 [MFREQE]
pub rad0: &'a [f64],
/// 前深度辐射 [MFREQE]
pub radm: &'a [f64],
/// FK 当前 [MFREQE]
pub fk0: &'a [f64],
/// FK 前深度 [MFREQE]
pub fkm: &'a [f64],
// ==================== 吸收/发射系数 ====================
/// 吸收系数当前 [MFREQE]
pub abso0: &'a [f64],
/// 吸收系数前深度 [MFREQE]
pub absom: &'a [f64],
/// 发射系数当前 [MFREQE]
pub emis0: &'a [f64],
/// 发射系数前深度 [MFREQE]
pub emism: &'a [f64],
/// 散射系数当前 [MFREQE]
pub scat0: &'a [f64],
// ==================== 导数 ====================
/// 吸收系数 T 导数当前 [MFREQE]
pub dabt0: &'a [f64],
/// 吸收系数 T 导数前深度 [MFREQE]
pub dabtm: &'a [f64],
/// 吸收系数 N 导数当前 [MFREQE]
pub dabn0: &'a [f64],
/// 吸收系数 N 导数前深度 [MFREQE]
pub dabnm: &'a [f64],
/// 吸收系数 M 导数当前 [MFREQE]
pub dabm0: &'a [f64],
/// 发射系数 T 导数当前 [MFREQE]
pub demt0: &'a [f64],
/// 发射系数 T 导数前深度 [MFREQE]
pub demtm: &'a [f64],
/// 发射系数 N 导数当前 [MFREQE]
pub demn0: &'a [f64],
/// 发射系数 N 导数前深度 [MFREQE]
pub demnm: &'a [f64],
/// 发射系数 M 导数当前 [MFREQE]
pub demm0: &'a [f64],
// ==================== 能级导数 ====================
/// 能级导数当前 [MLVEXP × MFREQE]
pub drch0: &'a [Vec<f64>],
/// 能级导数前深度 [MLVEXP × MFREQE]
pub drchm: &'a [Vec<f64>],
/// 能级发射导数当前 [MLVEXP × MFREQE]
pub dret0: &'a [Vec<f64>],
// ==================== 深度权重 ====================
/// 深度权重当前 [MFREQE]
pub wdep0: &'a [f64],
// ==================== 冷却率 ====================
/// ALI 冷却率 [MDEPTH]
pub fcool: &'a [f64],
// ==================== ALI 辐射等效项 ====================
/// REIT [MDEPTH]
pub reit: &'a [f64],
/// REIN [MDEPTH]
pub rein: &'a [f64],
/// REIM [MDEPTH]
pub reim: &'a [f64],
/// REIX [MDEPTH]
pub reix: &'a [f64],
/// REIP [MLVEXP × MDEPTH]
pub reip: &'a [Vec<f64>],
// ==================== ALI A 矩阵项 ====================
/// AREIT [MDEPTH]
pub areit: &'a [f64],
/// AREIN [MDEPTH]
pub arein: &'a [f64],
/// AREIM [MDEPTH]
pub areim: &'a [f64],
/// AREIP [MLVEXP × MDEPTH]
pub areip: &'a [Vec<f64>],
// ==================== ALI C 矩阵项 ====================
/// CREIT [MDEPTH]
pub creit: &'a [f64],
/// CREIN [MDEPTH]
pub crein: &'a [f64],
/// CREIM [MDEPTH]
pub creim: &'a [f64],
/// CREIX [MDEPTH]
pub creix: &'a [f64],
/// CREIP [MLVEXP × MDEPTH]
pub creip: &'a [Vec<f64>],
// ==================== RED 微分方程项 ====================
/// REDIF [MDEPTH]
pub redif: &'a [f64],
/// REDT [MDEPTH]
pub redt: &'a [f64],
/// REDTP [MDEPTH]
pub redtp: &'a [f64],
/// REDX [MDEPTH]
pub redx: &'a [f64],
/// REDXP [MDEPTH]
pub redxp: &'a [f64],
/// REDN [MDEPTH]
pub redn: &'a [f64],
/// REDP [MLVEXP × MDEPTH]
pub redp: &'a [Vec<f64>],
// ==================== RED 前深度项 ====================
/// REDTM [MDEPTH]
pub redtm: &'a [f64],
/// REDXM [MDEPTH]
pub redxm: &'a [f64],
/// REDNM [MDEPTH]
pub rednm: &'a [f64],
/// REDPM [MLVEXP × MDEPTH]
pub redpm: &'a [Vec<f64>],
// ==================== RED 后深度项 ====================
/// REDNP [MDEPTH]
pub rednp: &'a [f64],
/// REDPP [MLVEXP × MDEPTH]
pub redpp: &'a [Vec<f64>],
// ==================== 盘模型项 ====================
/// DTVIST [MDEPTH]
pub dtvist: &'a [f64],
/// DTVISR [MDEPTH]
pub dtvisr: &'a [f64],
/// DTVISN [MDEPTH]
pub dtvisn: &'a [f64],
// ==================== 边界条件 ====================
/// FH [MFREQ]
pub fh: &'a [f64],
/// HEXTRD [MFREQ]
pub hextrd: &'a [f64],
// ==================== 有效温度 ====================
pub teff: f64,
/// 氢原子质量
pub hmass: f64,
// ==================== 电子散射截面 ====================
/// SIGEC [MFREQ]
pub sigec: &'a [f64],
/// 汤姆逊散射截面
pub sige: f64,
/// CMD - Compton 密度导数
pub cmd: f64,
}
// ============================================================================
// BRE 可变状态结构体
// ============================================================================
/// BRE 可变状态(矩阵和向量)。
#[derive(Debug)]
pub struct BreState<'a> {
/// A 矩阵 [MTOT × MTOT]
pub a: &'a mut [Vec<f64>],
/// B 矩阵 [MTOT × MTOT]
pub b: &'a mut [Vec<f64>],
/// C 矩阵 [MTOT × MTOT]
pub c: &'a mut [Vec<f64>],
/// 左向量 VECL [MTOT]
pub vecl: &'a mut [f64],
/// REX 辅助数组 [MLEVEL]
pub rex: &'a mut [f64],
}
// ============================================================================
// Compton 辅助计算
// ============================================================================
/// 计算 Compton 散射辅助量(简化版)。
fn compt0_bre(
_ij: usize,
_id: usize,
ab: f64,
nfreq: usize,
kij: &[usize],
elec: &[f64],
sige: f64,
ij_idx: usize,
id_idx: usize,
) -> (f64, f64, f64, f64, f64, f64) {
// IJI = NFREQ - KIJ(IJ) + 1
let iji = nfreq - kij[ij_idx] + 1;
if iji == 1 {
return (0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
}
// 简化计算 - 完整实现需要调用 compt0 函数
let ss0 = elec[id_idx] * sige / ab;
// 返回 (CMA, CMB, CMC, CME, CMS, CMD)
(0.0, 0.0, 0.0, 0.0, ss0, 0.0)
}
// ============================================================================
// BRE 主函数
// ============================================================================
/// 计算辐射平衡方程的矩阵 A, B, C 部分。
///
/// # 参数
///
/// * `params` - 输入参数
/// * `state` - 可变状态(矩阵 A, B, C, VECL
///
/// # 说明
///
/// 此函数修改矩阵 A, B, C 的第 (NFREQE+INRE) 行,
/// 对应辐射平衡方程的积分部分和微分部分。
pub fn bre(params: &BreParams, state: &mut BreState) {
let id = params.id;
let id_idx = id - 1; // 0-indexed
// 计算矩阵列索引
let nhe = params.nfreqe + params.inhe;
let nre = params.nfreqe + params.inre;
let npc = params.nfreqe + params.inpc;
let nmp = params.nfreqe + params.inmp;
let nse = params.nfreqe + params.inse() - 1;
// IJ1 = 1 或 2 (Compton 散射时从 2 开始)
let mut ij1 = 1;
if params.icompt > 0 && params.icombc > 0 && !params.ijex.is_empty() && params.ijex[0] > 0 {
ij1 = 2;
}
// 检查温度控制
let ittc = (params.nretc.abs() / 100) as i32;
if params.iter > ittc {
let mod_val = (params.nretc.abs() % 100) as usize;
if id <= mod_val {
state.b[nre - 1][nre - 1] = 1.0;
if params.nretc < 0 {
state.c[nre - 1][nre - 1] = -1.0;
if id_idx + 1 < params.temp.len() {
state.vecl[nre - 1] = params.temp[id_idx + 1] - params.temp[id_idx];
}
}
return;
}
}
// RHS 向量初始化(ALI 冷却率)
state.vecl[nre - 1] = params.fcool[id_idx];
if params.idisk == 1 {
state.vecl[nre - 1] -= params.reint[id_idx] * params.tvisc[id_idx];
}
if params.reint[id_idx] <= 0.0 {
// 跳转到微分方程部分
bre_differential(params, state, id, nre, nhe, npc, nmp, nse);
return;
}
// ==================== 积分方程部分 ====================
let mut brepc = 0.0;
let mut bremp = 0.0;
// 初始化 REX
for i in 0..params.nlvexp.min(state.rex.len()) {
state.rex[i] = 0.0;
}
if params.nfreqe > 0 {
for ij in ij1..=params.nfreqe {
let ij_idx = ij - 1;
let ijt = params.ijfr[ij_idx];
// 累积积分项
let sigec_val = if ijt > 0 && ijt <= params.sigec.len() {
params.sigec[ijt - 1]
} else {
0.0
};
brepc += ((params.dabn0[ij_idx] - sigec_val) * params.rad0[ij_idx]
- params.demn0[ij_idx])
* params.wdep0[ij_idx];
bremp += (params.dabm0[ij_idx] * params.rad0[ij_idx] - params.demm0[ij_idx])
* params.wdep0[ij_idx];
for i in 0..params.nlvexp.min(state.rex.len()) {
state.rex[i] += (params.drch0[i][ij_idx] * params.rad0[ij_idx]
- params.dret0[i][ij_idx])
* params.wdep0[ij_idx];
}
// B 矩阵对角项
state.b[nre - 1][nre - 1] += (params.dabt0[ij_idx] * params.rad0[ij_idx]
- params.demt0[ij_idx])
* params.wdep0[ij_idx]
* params.reint[id_idx];
// 加热项
let heat = params.abso0[ij_idx] - params.scat0[ij_idx];
state.b[nre - 1][ij - 1] = params.wdep0[ij_idx] * heat * params.reint[id_idx];
// RHS 向量
state.vecl[nre - 1] -= (heat * params.rad0[ij_idx] - params.emis0[ij_idx])
* params.wdep0[ij_idx]
* params.reint[id_idx];
// Compton 散射项
if params.icompt > 5 {
let (_cma, cmb, _cmc, cme, cms, _cmd) = compt0_bre(
ijt,
id,
params.abso0[ij_idx],
params.nfreq,
params.kij,
params.elec,
params.sige,
ij_idx,
id_idx,
);
state.vecl[nre - 1] +=
params.abso0[ij_idx] * cms * params.wdep0[ij_idx] * params.reint[id_idx];
if params.icompt > 6 {
if params.icmdra > 0 {
state.b[nre - 1][ij - 1] -=
params.abso0[ij_idx] * (cmb + cme) * params.wdep0[ij_idx]
* params.reint[id_idx];
} else {
state.b[nre - 1][ij - 1] -=
params.abso0[ij_idx] * (cmb + cme) * params.reint[id_idx];
}
// 注:完整的 Compton 处理需要更多代码
}
}
}
}
// ALI 修正项
state.b[nre - 1][nre - 1] += params.reit[id_idx] * params.reint[id_idx];
if params.inpc > 0 {
state.b[nre - 1][npc - 1] += (brepc + params.rein[id_idx]) * params.reint[id_idx];
}
if params.inmp > 0 {
state.b[nre - 1][nmp - 1] += (bremp + params.reim[id_idx]) * params.reint[id_idx];
}
if params.inhe > 0 {
state.b[nre - 1][nhe - 1] = params.reix[id_idx] * params.reint[id_idx];
}
// IFALI > 5 时的 A 和 C 矩阵项
if params.ifali > 5 {
state.a[nre - 1][nre - 1] = params.areit[id_idx] * params.reint[id_idx];
if params.inpc > 0 {
state.a[nre - 1][npc - 1] = params.arein[id_idx] * params.reint[id_idx];
}
if params.inmp > 0 {
state.a[nre - 1][nmp - 1] = params.areim[id_idx] * params.reint[id_idx];
}
state.c[nre - 1][nre - 1] = params.creit[id_idx] * params.reint[id_idx];
if params.inpc > 0 {
state.c[nre - 1][npc - 1] = params.crein[id_idx] * params.reint[id_idx];
}
if params.inmp > 0 {
state.c[nre - 1][nmp - 1] = params.creim[id_idx] * params.reint[id_idx];
}
if params.inhe > 0 {
state.c[nre - 1][nhe - 1] = params.creix[id_idx] * params.reint[id_idx];
}
}
// 盘模型项
if params.idisk == 1 {
state.b[nre - 1][nre - 1] += params.dtvist[id_idx] * params.reint[id_idx];
if params.inpc > 0 {
state.b[nre - 1][npc - 1] -= params.dtvisr[id_idx] * params.reint[id_idx];
}
if params.inhe > 0 {
state.b[nre - 1][nhe - 1] =
(params.dtvisr[id_idx] + params.dtvisn[id_idx]) * params.reint[id_idx];
}
if params.inmp > 0 {
state.b[nre - 1][nmp - 1] = params.dtvisr[id_idx] * params.hmass
/ params.wmm[id_idx]
* params.reint[id_idx];
}
}
// 能级相关项
for ii in 0..params.nlvexp {
if nse + ii < state.b[nre - 1].len() {
state.b[nre - 1][nse + ii] +=
(state.rex[ii] + params.reip[ii][id_idx]) * params.reint[id_idx];
}
}
if params.ifali > 5 && id > 1 {
for ii in 0..params.nlvexp {
if nse + ii < state.a[nre - 1].len() {
state.a[nre - 1][nse + ii] += params.areip[ii][id_idx] * params.reint[id_idx];
}
}
}
if params.ifali > 5 && id < params.nd {
for ii in 0..params.nlvexp {
if nse + ii < state.c[nre - 1].len() {
state.c[nre - 1][nse + ii] += params.creip[ii][id_idx] * params.reint[id_idx];
}
}
}
// ==================== 微分方程部分 ====================
bre_differential(params, state, id, nre, nhe, npc, nmp, nse);
}
/// 计算微分方程部分。
fn bre_differential(
params: &BreParams,
state: &mut BreState,
id: usize,
nre: usize,
nhe: usize,
npc: usize,
nmp: usize,
nse: usize,
) {
let id_idx = id - 1;
if params.redif[id_idx] == 0.0 {
return;
}
// TEFF^4 项
let mut teffd = params.teff.powi(4);
if params.idisk == 1 {
teffd *= UN - params.thetav[id_idx];
}
state.vecl[nre - 1] += SIG4P * teffd * params.redif[id_idx];
if id == 1 {
// 上边界条件
bre_upper_boundary(params, state, nre, nhe, npc, nse);
return;
}
// ==================== 内部深度点 ====================
let ddm = (params.dm[id_idx] - params.dm[id_idx - 1]) * HALF;
let mut aren = 0.0;
let mut bren = 0.0;
let mut arepc = 0.0;
let mut brepc = 0.0;
// GP, GN 几何因子
let (gp, gn) = if params.inmp > 0 { (UN, 0.0) } else { (0.0, UN) };
// 初始化辅助数组
let mut rexa = vec![0.0; params.nlvexp];
let mut rexb = vec![0.0; params.nlvexp];
if params.nfreqe > 0 {
for ij in 1..=params.nfreqe {
let ij_idx = ij - 1;
let omeg0 = params.abso0[ij_idx] * params.dens1[id_idx];
let omegm = params.absom[ij_idx] * params.dens1[id_idx - 1];
let dtaum = (omeg0 + omegm) * ddm;
if dtaum.abs() < 1e-30 {
continue;
}
let frd = params.fk0[ij_idx] * params.rad0[ij_idx]
- params.fkm[ij_idx] * params.radm[ij_idx];
let gamr = frd / dtaum;
let a1 = gamr / (omeg0 + omegm);
let a3r = a1 * params.dens1[id_idx - 1] * params.wdep0[ij_idx];
let b3r = a1 * params.dens1[id_idx] * params.wdep0[ij_idx];
// A 矩阵项
state.a[nre - 1][ij - 1] =
-params.wdep0[ij_idx] * params.fkm[ij_idx] / dtaum * params.redif[id_idx];
let rtr = omegm * params.wmm[id_idx - 1] * a3r;
aren += rtr * gn;
arepc -= a3r * params.dabnm[ij_idx] + rtr * gn;
if params.inmp != 0 {
state.a[nre - 1][nmp - 1] += rtr * gp * params.redif[id_idx];
}
state.a[nre - 1][nre - 1] -= a3r * params.dabtm[ij_idx] * params.redif[id_idx];
// B 矩阵项
state.b[nre - 1][ij - 1] +=
params.wdep0[ij_idx] * params.fk0[ij_idx] / dtaum * params.redif[id_idx];
let rtr = omeg0 * params.wmm[id_idx] * b3r;
bren += rtr * gn;
brepc -= b3r * params.dabn0[ij_idx] - rtr * gn;
if params.inmp != 0 {
state.b[nre - 1][nmp - 1] +=
(rtr + params.redx[id_idx]) * gp * params.redif[id_idx];
}
// 温度列
state.b[nre - 1][nre - 1] -= b3r * params.dabt0[ij_idx] * params.redif[id_idx];
// 能级导数
for i in 0..params.nlvexp {
rexa[i] -= a3r * params.drchm[i][ij_idx];
rexb[i] -= b3r * params.drch0[i][ij_idx];
}
// RHS
state.vecl[nre - 1] -= params.wdep0[ij_idx] * gamr * params.redif[id_idx];
}
}
// N 列(氦)
if params.inhe != 0 {
state.a[nre - 1][nhe - 1] = (aren + params.redxm[id_idx]) * params.redif[id_idx];
state.b[nre - 1][nhe - 1] += (bren + params.redx[id_idx]) * params.redif[id_idx];
}
// 温度列
state.a[nre - 1][nre - 1] += params.redtm[id_idx] * params.redif[id_idx];
state.b[nre - 1][nre - 1] += params.redt[id_idx] * params.redif[id_idx];
state.c[nre - 1][nre - 1] += params.redtp[id_idx] * params.redif[id_idx];
// 电子密度列
if params.inpc != 0 {
state.a[nre - 1][npc - 1] +=
(arepc + params.rednm[id_idx] - params.redxm[id_idx]) * params.redif[id_idx];
state.b[nre - 1][npc - 1] +=
(brepc + params.redn[id_idx] - params.redx[id_idx]) * params.redif[id_idx];
state.c[nre - 1][npc - 1] += params.rednp[id_idx] * params.redif[id_idx];
}
// 能级列
for ii in 0..params.nlvexp {
if nse + ii < state.a[nre - 1].len() {
state.a[nre - 1][nse + ii] +=
(rexa[ii] + params.redpm[ii][id_idx]) * params.redif[id_idx];
}
if nse + ii < state.b[nre - 1].len() {
state.b[nre - 1][nse + ii] +=
(rexb[ii] + params.redp[ii][id_idx]) * params.redif[id_idx];
}
if nse + ii < state.c[nre - 1].len() {
state.c[nre - 1][nse + ii] += params.redpp[ii][id_idx] * params.redif[id_idx];
}
}
}
/// 上边界条件(ID = 1)。
fn bre_upper_boundary(
params: &BreParams,
state: &mut BreState,
nre: usize,
nhe: usize,
npc: usize,
nse: usize,
) {
let id_idx = 0; // ID = 1
if params.nfreqe > 0 {
for ij in 1..=params.nfreqe {
let ij_idx = ij - 1;
let ijt = params.ijfr[ij_idx];
let fh_val = if ijt > 0 && ijt <= params.fh.len() {
params.fh[ijt - 1]
} else {
0.0
};
let hextrd_val = if ijt > 0 && ijt <= params.hextrd.len() {
params.hextrd[ijt - 1]
} else {
0.0
};
let wf = params.wdep0[ij_idx] * fh_val * params.redif[id_idx];
state.b[nre - 1][ij - 1] += wf;
state.vecl[nre - 1] -=
wf * params.rad0[ij_idx] + params.wdep0[ij_idx] * hextrd_val * params.redif[id_idx];
}
}
// 温度列
state.b[nre - 1][nre - 1] += params.redt[id_idx] * params.redif[id_idx];
state.c[nre - 1][nre - 1] += params.redtp[id_idx] * params.redif[id_idx];
// N 和电子密度列
if params.inhe != 0 {
state.b[nre - 1][nhe - 1] += params.redx[id_idx] * params.redif[id_idx];
}
if params.inpc != 0 {
state.b[nre - 1][npc - 1] += params.redn[id_idx] * params.redif[id_idx];
}
if params.inhe != 0 {
state.c[nre - 1][nhe - 1] += params.redxp[id_idx] * params.redif[id_idx];
}
if params.inpc != 0 {
state.c[nre - 1][npc - 1] += params.rednp[id_idx] * params.redif[id_idx];
}
// 能级列
for ii in 0..params.nlvexp {
if nse + ii < state.b[nre - 1].len() {
state.b[nre - 1][nse + ii] += params.redp[ii][id_idx] * params.redif[id_idx];
}
if nse + ii < state.c[nre - 1].len() {
state.c[nre - 1][nse + ii] += params.redpp[ii][id_idx] * params.redif[id_idx];
}
}
}
// ============================================================================
// 辅助方法
// ============================================================================
impl<'a> BreParams<'a> {
/// 计算 INSE(谱线起始索引)
fn inse(&self) -> usize {
1
}
}
// ============================================================================
// 测试
// ============================================================================
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_bre_compile() {
// 编译测试 - 验证类型签名正确
assert!(true);
}
#[test]
fn test_compt0_bre_iji_1() {
// 当 IJI = 1 时,所有输出应为 0
let kij = vec![100]; // NFREQ - 100 + 1 = 1
let elec = vec![1e12];
let result = compt0_bre(1, 1, 1e-8, 100, &kij, &elec, 6.6516e-25, 0, 0);
assert!((result.0).abs() < 1e-15);
assert!((result.1).abs() < 1e-15);
assert!((result.2).abs() < 1e-15);
}
}
+656
View File
@@ -0,0 +1,656 @@
//! 辐射平衡方程矩阵计算(几何深度版本)。
//!
//! 重构自 TLUSTY `brez.f`
//!
//! 与 BRE 类似,但使用几何深度 ZD 而非柱质量密度 DM 进行差分。
//! 用于球对称或柱对称几何配置。
use crate::state::constants::{HALF, SIG4P, UN};
/// 最大线性化能级数
const MLVEXP: usize = 233;
// ============================================================================
// BREZ 参数结构体
// ============================================================================
/// BREZ 输入参数。
#[derive(Debug)]
pub struct BrezParams<'a> {
// ==================== 基本参数 ====================
/// 深度索引 ID (1-indexed)
pub id: usize,
// ==================== 维度参数 ====================
pub nd: usize,
pub nfreqe: usize,
pub nfreq: usize,
pub nlvexp: usize,
// ==================== 索引参数 ====================
pub inhe: usize,
pub inre: usize,
pub inpc: usize,
pub inmp: usize,
// ==================== 控制参数 ====================
pub icompt: i32,
pub icombc: i32,
pub icmdra: i32,
pub iter: i32,
pub nretc: i32,
pub idisk: i32,
pub ifali: i32,
// ==================== 频率映射 ====================
pub ijfr: &'a [usize],
pub ijorig: &'a [usize],
pub kij: &'a [usize],
pub ijex: &'a [i32],
// ==================== 模型状态 ====================
pub temp: &'a [f64],
pub elec: &'a [f64],
pub zd: &'a [f64], // 几何深度(BREZ 特有)
pub dens1: &'a [f64],
pub wmm: &'a [f64],
pub tvisc: &'a [f64],
pub reint: &'a [f64],
pub thetav: &'a [f64],
// ==================== 辐射场 ====================
pub rad0: &'a [f64],
pub radm: &'a [f64], // 前深度辐射(BREZ 使用)
pub fk0: &'a [f64],
pub fkm: &'a [f64],
// ==================== 吸收/发射系数 ====================
pub abso0: &'a [f64],
pub absom: &'a [f64], // 前深度吸收(BREZ 使用)
pub emis0: &'a [f64],
pub emism: &'a [f64],
pub scat0: &'a [f64],
// ==================== 导数 ====================
pub dabt0: &'a [f64],
pub dabtm: &'a [f64],
pub dabn0: &'a [f64],
pub dabnm: &'a [f64],
pub dabm0: &'a [f64],
pub demt0: &'a [f64],
pub demtm: &'a [f64],
pub demn0: &'a [f64],
pub demnm: &'a [f64],
pub demm0: &'a [f64],
// ==================== 能级导数 ====================
pub drch0: &'a [Vec<f64>],
pub drchm: &'a [Vec<f64>],
pub dret0: &'a [Vec<f64>],
// ==================== 深度权重 ====================
pub wdep0: &'a [f64],
// ==================== 冷却率 ====================
pub fcool: &'a [f64],
// ==================== ALI 辐射等效项 ====================
pub reit: &'a [f64],
pub rein: &'a [f64],
pub reim: &'a [f64],
pub reix: &'a [f64],
pub reip: &'a [Vec<f64>],
// ==================== ALI A 矩阵项 ====================
pub areit: &'a [f64],
pub arein: &'a [f64],
pub areip: &'a [Vec<f64>],
// ==================== ALI C 矩阵项 ====================
pub creit: &'a [f64],
pub crein: &'a [f64],
pub creim: &'a [f64],
pub creix: &'a [f64],
pub creip: &'a [Vec<f64>],
// ==================== RED 微分方程项 ====================
pub redif: &'a [f64],
pub redt: &'a [f64],
pub redtp: &'a [f64],
pub redn: &'a [f64],
pub rednm: &'a [f64],
pub redp: &'a [Vec<f64>],
pub redpm: &'a [Vec<f64>],
pub rednp: &'a [f64],
pub redpp: &'a [Vec<f64>],
pub redtm: &'a [f64],
// ==================== 盘模型项 ====================
pub dtvist: &'a [f64],
pub dtvisr: &'a [f64],
pub dtvisn: &'a [f64],
// ==================== 边界条件 ====================
pub fh: &'a [f64],
// ==================== 有效温度 ====================
pub teff: f64,
pub hmass: f64,
// ==================== 电子散射截面 ====================
pub sigec: &'a [f64],
pub sige: f64,
pub cmd: f64,
}
/// BREZ 可变状态(矩阵和向量)。
#[derive(Debug)]
pub struct BrezState<'a> {
pub a: &'a mut [Vec<f64>],
pub b: &'a mut [Vec<f64>],
pub c: &'a mut [Vec<f64>],
pub vecl: &'a mut [f64],
pub rex: &'a mut [f64],
}
// ============================================================================
// Compton 辅助计算
// ============================================================================
/// 简化版 Compton 计算。
fn compt0_brez(
ij: usize,
id: usize,
ab: f64,
nfreq: usize,
kij: &[usize],
elec: &[f64],
sige: f64,
) -> (f64, f64, f64, f64, f64, f64) {
let iji = nfreq - kij[ij - 1] + 1;
if iji == 1 {
return (0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
}
let ss0 = elec[id - 1] * sige / ab;
(0.0, 0.0, 0.0, 0.0, ss0, 0.0)
}
// ============================================================================
// BREZ 主函数
// ============================================================================
/// 计算辐射平衡方程的矩阵 A, B, C 部分(几何深度版本)。
pub fn brez(params: &BrezParams, state: &mut BrezState) {
let id = params.id;
let id_idx = id - 1;
// 计算矩阵列索引
let nre = params.nfreqe + params.inre;
let nhe = params.nfreqe + params.inhe;
let npc = params.nfreqe + params.inpc;
let nmp = params.nfreqe + params.inmp;
let nse = params.nfreqe + params.inse() - 1;
// IJ1 = 1 或 2
let mut ij1 = 1;
if params.icompt > 0 && params.icombc > 0 && !params.ijex.is_empty() && params.ijex[0] > 0 {
ij1 = 2;
}
// 温度控制检查
let ittc = (params.nretc.abs() / 100) as i32;
if params.iter > ittc {
let mod_val = (params.nretc.abs() % 100) as usize;
if id <= mod_val {
state.b[nre - 1][nre - 1] = 1.0;
if params.nretc < 0 {
state.c[nre - 1][nre - 1] = -1.0;
if id_idx + 1 < params.temp.len() {
state.vecl[nre - 1] = params.temp[id_idx + 1] - params.temp[id_idx];
}
}
return;
}
}
// RHS 向量初始化
state.vecl[nre - 1] = params.fcool[id_idx] - params.reint[id_idx] * params.tvisc[id_idx];
if params.reint[id_idx] <= 0.0 {
brez_differential(params, state, id, nre, nhe, npc, nse);
return;
}
// ==================== 积分方程部分 ====================
let mut brepc = 0.0;
let mut bremp = 0.0;
for i in 0..params.nlvexp.min(state.rex.len()) {
state.rex[i] = 0.0;
}
if params.nfreqe > 0 {
for ij in ij1..=params.nfreqe {
let ij_idx = ij - 1;
let ijt = params.ijfr[ij_idx];
let sigec_val = if ijt > 0 && ijt <= params.sigec.len() {
params.sigec[ijt - 1]
} else {
0.0
};
brepc += ((params.dabn0[ij_idx] - sigec_val) * params.rad0[ij_idx]
- params.demn0[ij_idx])
* params.wdep0[ij_idx];
bremp += (params.dabm0[ij_idx] * params.rad0[ij_idx] - params.demm0[ij_idx])
* params.wdep0[ij_idx];
for i in 0..params.nlvexp.min(state.rex.len()) {
state.rex[i] += (params.drch0[i][ij_idx] * params.rad0[ij_idx]
- params.dret0[i][ij_idx])
* params.wdep0[ij_idx];
}
// B 矩阵对角项
state.b[nre - 1][nre - 1] += (params.dabt0[ij_idx] * params.rad0[ij_idx]
- params.demt0[ij_idx])
* params.wdep0[ij_idx]
* params.reint[id_idx];
// 加热项(HEAT = ABSO0 - SCAT0
let heat = params.abso0[ij_idx] - params.scat0[ij_idx];
state.b[nre - 1][ij - 1] = params.wdep0[ij_idx] * heat * params.reint[id_idx];
// RHS 向量
state.vecl[nre - 1] -= (heat * params.rad0[ij_idx] - params.emis0[ij_idx])
* params.wdep0[ij_idx]
* params.reint[id_idx];
// Compton 散射项
if params.icompt > 5 {
let (_cma, cmb, _cmc, cme, cms, cmd) = compt0_brez(
ijt,
id,
params.abso0[ij_idx],
params.nfreq,
params.kij,
params.elec,
params.sige,
);
state.vecl[nre - 1] +=
params.abso0[ij_idx] * cms * params.wdep0[ij_idx] * params.reint[id_idx];
if params.icompt > 6 {
if params.icmdra > 0 {
state.b[nre - 1][ij - 1] -=
params.abso0[ij_idx] * (cmb + cme) * params.wdep0[ij_idx]
* params.reint[id_idx];
} else {
state.b[nre - 1][ij - 1] -=
params.abso0[ij_idx] * (cmb + cme) * params.reint[id_idx];
}
// 邻近频率项
let iji = params.nfreq - params.kij[ijt - 1] + 1;
if iji > 1 {
let ijm = params.ijex[params.ijorig[iji - 2]] as usize;
if ijm > 0 {
if params.icmdra > 0 {
state.b[nre - 1][ijm - 1] -=
params.abso0[ij_idx] * _cma * params.wdep0[ij_idx]
* params.reint[id_idx];
} else {
state.b[nre - 1][ijm - 1] -=
params.abso0[ij_idx] * _cma * params.reint[id_idx];
}
}
}
if iji < params.nfreq {
let ijp = params.ijex[params.ijorig[iji]] as usize;
if ijp > 0 {
if params.icmdra > 0 {
state.b[nre - 1][ijp - 1] -=
params.abso0[ij_idx] * _cmc * params.wdep0[ij_idx]
* params.reint[id_idx];
} else {
state.b[nre - 1][ijp - 1] -=
params.abso0[ij_idx] * _cmc * params.reint[id_idx];
}
}
}
state.b[nre - 1][nre - 1] -=
cmd * params.abso0[ij_idx] * params.wdep0[ij_idx] * params.reint[id_idx];
state.b[nre - 1][npc - 1] -= cms * params.abso0[ij_idx] / params.elec[id_idx]
* params.wdep0[ij_idx]
* params.reint[id_idx];
}
}
}
}
// ALI 修正项
state.b[nre - 1][nre - 1] += params.reit[id_idx] * params.reint[id_idx];
if params.inpc > 0 {
state.b[nre - 1][npc - 1] += (brepc + params.rein[id_idx]) * params.reint[id_idx];
}
if params.inmp > 0 {
state.b[nre - 1][nmp - 1] += (bremp + params.reim[id_idx]) * params.reint[id_idx];
}
if params.inhe > 0 {
state.b[nre - 1][nhe - 1] = params.reix[id_idx] * params.reint[id_idx];
}
// A 和 C 矩阵项(BREZ 总是设置,不像 BRE 需要 IFALI > 5
state.a[nre - 1][nre - 1] = params.areit[id_idx] * params.reint[id_idx];
if params.inpc > 0 {
state.a[nre - 1][npc - 1] = params.arein[id_idx] * params.reint[id_idx];
}
state.c[nre - 1][nre - 1] = params.creit[id_idx] * params.reint[id_idx];
if params.inpc > 0 {
state.c[nre - 1][npc - 1] = params.crein[id_idx] * params.reint[id_idx];
}
if params.inmp > 0 {
state.c[nre - 1][nmp - 1] = params.creim[id_idx] * params.reint[id_idx];
}
if params.inhe > 0 {
state.c[nre - 1][nhe - 1] = params.creix[id_idx] * params.reint[id_idx];
}
// 盘模型项
state.b[nre - 1][nre - 1] += params.dtvist[id_idx] * params.reint[id_idx];
if params.inpc > 0 {
state.b[nre - 1][npc - 1] -= params.dtvisr[id_idx] * params.reint[id_idx];
}
if params.inhe > 0 {
state.b[nre - 1][nhe - 1] =
(params.dtvisr[id_idx] + params.dtvisn[id_idx]) * params.reint[id_idx];
}
if params.inmp > 0 {
state.b[nre - 1][nmp - 1] =
params.dtvisr[id_idx] * params.hmass / params.wmm[id_idx] * params.reint[id_idx];
}
// 能级相关项
for ii in 0..params.nlvexp {
if nse + ii < state.b[nre - 1].len() {
state.b[nre - 1][nse + ii] +=
(state.rex[ii] + params.reip[ii][id_idx]) * params.reint[id_idx];
}
}
if params.ifali > 5 && id > 1 {
for ii in 0..params.nlvexp {
if nse + ii < state.a[nre - 1].len() {
state.a[nre - 1][nse + ii] += params.areip[ii][id_idx] * params.reint[id_idx];
}
}
}
if params.ifali > 5 && id < params.nd {
for ii in 0..params.nlvexp {
if nse + ii < state.c[nre - 1].len() {
state.c[nre - 1][nse + ii] += params.creip[ii][id_idx] * params.reint[id_idx];
}
}
}
// ==================== 微分方程部分 ====================
brez_differential(params, state, id, nre, nhe, npc, nse);
}
/// 计算微分方程部分(使用 ZD 几何深度)。
fn brez_differential(
params: &BrezParams,
state: &mut BrezState,
id: usize,
nre: usize,
_nhe: usize,
npc: usize,
nse: usize,
) {
let id_idx = id - 1;
if params.redif[id_idx] == 0.0 {
return;
}
// TEFF^4 项(盘模型总是有 THETAV
let teffd = params.teff.powi(4) * (UN - params.thetav[id_idx]);
state.vecl[nre - 1] += SIG4P * teffd * params.redif[id_idx];
if id == 1 {
brez_upper_boundary(params, state, nre, npc, nse);
return;
}
// ==================== 内部深度点 ====================
// 使用 ZD 差分(BREZ 特有)
let ddm = (params.zd[id_idx - 1] - params.zd[id_idx]) * HALF;
let mut arepc = 0.0;
let mut brepc = 0.0;
let mut rexa = vec![0.0; params.nlvexp];
let mut rexb = vec![0.0; params.nlvexp];
if params.nfreqe > 0 {
for ij in 1..=params.nfreqe {
let ij_idx = ij - 1;
let omeg0 = params.abso0[ij_idx];
let omegm = params.absom[ij_idx];
let dtaum = (omeg0 + omegm) * ddm;
if dtaum.abs() < 1e-30 {
continue;
}
let frd = params.fk0[ij_idx] * params.rad0[ij_idx]
- params.fkm[ij_idx] * params.radm[ij_idx];
let gamr = frd / dtaum;
let a1 = gamr / (omeg0 + omegm) * params.wdep0[ij_idx];
// A 矩阵项
state.a[nre - 1][ij - 1] =
-params.wdep0[ij_idx] * params.fkm[ij_idx] / dtaum * params.redif[id_idx];
arepc -= a1 * params.dabnm[ij_idx];
state.a[nre - 1][nre - 1] -= a1 * params.dabtm[ij_idx] * params.redif[id_idx];
// B 矩阵项
state.b[nre - 1][ij - 1] +=
params.wdep0[ij_idx] * params.fk0[ij_idx] / dtaum * params.redif[id_idx];
brepc -= a1 * params.dabn0[ij_idx];
// 温度列
state.b[nre - 1][nre - 1] -= a1 * params.dabt0[ij_idx] * params.redif[id_idx];
// 能级导数
for i in 0..params.nlvexp {
rexa[i] -= a1 * params.drchm[i][ij_idx];
rexb[i] -= a1 * params.drch0[i][ij_idx];
}
// RHS
state.vecl[nre - 1] -= params.wdep0[ij_idx] * gamr * params.redif[id_idx];
}
}
// 温度列
state.a[nre - 1][nre - 1] += params.redtm[id_idx] * params.redif[id_idx];
state.b[nre - 1][nre - 1] += params.redt[id_idx] * params.redif[id_idx];
state.c[nre - 1][nre - 1] += params.redtp[id_idx] * params.redif[id_idx];
// 电子密度列
if params.inpc != 0 {
state.a[nre - 1][npc - 1] +=
(arepc + params.rednm[id_idx]) * params.redif[id_idx];
state.b[nre - 1][npc - 1] +=
(brepc + params.redn[id_idx]) * params.redif[id_idx];
state.c[nre - 1][npc - 1] += params.rednp[id_idx] * params.redif[id_idx];
}
// 能级列
for ii in 0..params.nlvexp {
if nse + ii < state.a[nre - 1].len() {
state.a[nre - 1][nse + ii] +=
(rexa[ii] + params.redpm[ii][id_idx]) * params.redif[id_idx];
}
if nse + ii < state.b[nre - 1].len() {
state.b[nre - 1][nse + ii] +=
(rexb[ii] + params.redp[ii][id_idx]) * params.redif[id_idx];
}
if nse + ii < state.c[nre - 1].len() {
state.c[nre - 1][nse + ii] += params.redpp[ii][id_idx] * params.redif[id_idx];
}
}
}
/// 上边界条件(ID = 1)。
fn brez_upper_boundary(
params: &BrezParams,
state: &mut BrezState,
nre: usize,
npc: usize,
nse: usize,
) {
let id_idx = 0;
if params.nfreqe > 0 {
for ij in 1..=params.nfreqe {
let ij_idx = ij - 1;
let ijt = params.ijfr[ij_idx];
let fh_val = if ijt > 0 && ijt <= params.fh.len() {
params.fh[ijt - 1]
} else {
0.0
};
let wf = params.wdep0[ij_idx] * fh_val * params.redif[id_idx];
state.b[nre - 1][ij - 1] += wf;
state.vecl[nre - 1] -= wf * params.rad0[ij_idx];
}
}
// 温度列
state.b[nre - 1][nre - 1] += params.redt[id_idx] * params.redif[id_idx];
// 电子密度列
if params.inpc != 0 {
state.b[nre - 1][npc - 1] += params.redn[id_idx] * params.redif[id_idx];
}
// 能级列
for ii in 0..params.nlvexp {
if nse + ii < state.b[nre - 1].len() {
state.b[nre - 1][nse + ii] += params.redp[ii][id_idx] * params.redif[id_idx];
}
}
}
// ============================================================================
// 辅助方法
// ============================================================================
impl<'a> BrezParams<'a> {
fn inse(&self) -> usize {
1
}
}
// ============================================================================
// 测试
// ============================================================================
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_brez_compile() {
// 编译测试 - 验证类型签名正确
assert!(true);
}
#[test]
fn test_brez_inse() {
// INSE 应该返回 1
assert_eq!(BrezParams::inse(&BrezParams {
id: 1, nd: 2, nfreqe: 1, nfreq: 1, nlvexp: 1,
inhe: 0, inre: 1, inpc: 1, inmp: 0,
icompt: 0, icombc: 0, icmdra: 0, iter: 1, nretc: 0,
idisk: 0, ifali: 0,
ijfr: &[], ijorig: &[], kij: &[], ijex: &[],
temp: &[], elec: &[], zd: &[], dens1: &[], wmm: &[],
tvisc: &[], reint: &[], thetav: &[],
rad0: &[], radm: &[], fk0: &[], fkm: &[],
abso0: &[], absom: &[], emis0: &[], emism: &[], scat0: &[],
dabt0: &[], dabtm: &[], dabn0: &[], dabnm: &[], dabm0: &[],
demt0: &[], demtm: &[], demn0: &[], demnm: &[], demm0: &[],
drch0: &[], drchm: &[], dret0: &[], wdep0: &[], fcool: &[],
reit: &[], rein: &[], reim: &[], reix: &[], reip: &[],
areit: &[], arein: &[], areip: &[],
creit: &[], crein: &[], creim: &[], creix: &[], creip: &[],
redif: &[], redt: &[], redtp: &[], redn: &[], rednm: &[],
redp: &[], redpm: &[], rednp: &[], redpp: &[], redtm: &[],
dtvist: &[], dtvisr: &[], dtvisn: &[], fh: &[],
teff: 10000.0, hmass: 1.67333e-24,
sigec: &[], sige: 6.6516e-25, cmd: 0.0,
}), 1);
}
#[test]
fn test_compt0_brez_iji_1() {
// 当 IJI = 1 时,所有输出应为 0
let kij = vec![100]; // NFREQ - 100 + 1 = 1
let elec = vec![1e12];
let result = compt0_brez(1, 1, 1e-8, 100, &kij, &elec, 6.6516e-25);
assert!((result.0).abs() < 1e-15);
assert!((result.1).abs() < 1e-15);
assert!((result.2).abs() < 1e-15);
}
#[test]
fn test_brez_constants() {
// 验证常量值
assert!((HALF - 0.5).abs() < 1e-15);
assert!((UN - 1.0).abs() < 1e-15);
assert!(SIG4P > 0.0);
}
#[test]
fn test_brez_matrix_indices() {
// 测试矩阵索引计算
let nfreqe = 5;
let inre = 1;
let inpc = 2;
let inhe = 3;
let inmp = 4;
// NRE = NFREQE + INRE
let nre = nfreqe + inre;
assert_eq!(nre, 6);
// NPC = NFREQE + INPC
let npc = nfreqe + inpc;
assert_eq!(npc, 7);
// NHE = NFREQE + INHE
let nhe = nfreqe + inhe;
assert_eq!(nhe, 8);
// NMP = NFREQE + INMP
let nmp = nfreqe + inmp;
assert_eq!(nmp, 9);
}
#[test]
fn test_brez_differential_ddm() {
// 测试 DDM 计算(几何深度差分)
let zd = vec![1.0, 0.5, 0.0];
let id: usize = 2;
let ddm = (zd[id - 2] - zd[id - 1]) * HALF;
assert!((ddm - 0.25).abs() < 1e-15);
}
}
+929
View File
@@ -0,0 +1,929 @@
//! 辐射转移方程矩阵计算。
//!
//! 重构自 TLUSTY `brte.f`
//!
//! 计算线性化辐射转移方程的矩阵 A, B, C 部分(前 NFREQE 行)。
//! 处理三种深度情况:上边界、内部点、下边界。
use crate::state::constants::{HALF, UN};
// ============================================================================
// 常量
// ============================================================================
/// Compton 常量
const XCON: f64 = 8.0935e-21;
const YCON: f64 = 1.68638e-10;
/// 1/6
const SIXTH: f64 = 1.0 / 6.0;
/// 1/3
const THIRD: f64 = 1.0 / 3.0;
// ============================================================================
// BRTE 参数结构体
// ============================================================================
/// BRTE 输入参数。
#[derive(Debug)]
pub struct BrteParams<'a> {
// ==================== 基本参数 ====================
pub id: usize,
pub nd: usize,
pub nfreqe: usize,
pub nfreq: usize,
pub nlvexp: usize,
// ==================== 索引参数 ====================
pub inhe: usize,
pub inre: usize,
pub inpc: usize,
pub inse: usize,
pub inmp: usize,
// ==================== 控制参数 ====================
pub icompt: i32,
pub icombc: i32,
pub ichcoo: i32,
pub isplin: i32,
pub idisk: i32,
pub ifz0: i32,
pub ibc: i32,
pub iwinbl: i32,
pub inre_idx: i32,
pub inpc_idx: i32,
pub ndre: usize,
pub radzer: f64,
// ==================== 频率映射 ====================
pub ijfr: &'a [usize],
pub ijorig: &'a [usize],
pub kij: &'a [usize],
pub ijex: &'a [i32],
pub kijt: &'a [usize],
// ==================== 模型状态 ====================
pub temp: &'a [f64],
pub tempbd: f64,
pub dm: &'a [f64],
pub dens: &'a [f64],
pub wmm: &'a [f64],
pub elec: &'a [f64],
// ==================== 频率相关 ====================
pub freq: &'a [f64],
pub dlnfr: &'a [f64],
pub delj: &'a [Vec<f64>],
// ==================== 辐射场 ====================
pub rad0: &'a [f64],
pub radm: &'a [f64],
pub radp: &'a [f64],
pub radex: &'a [Vec<f64>],
// ==================== FK 系数 ====================
pub fk0: &'a [f64],
pub fkm: &'a [f64],
pub fkp: &'a [f64],
// ==================== 吸收/发射系数 ====================
pub abso0: &'a [f64],
pub absom: &'a [f64],
pub absop: &'a [f64],
pub emis0: &'a [f64],
pub emism: &'a [f64],
pub emisp: &'a [f64],
pub scat0: &'a [f64],
pub scatm: &'a [f64],
pub scatp: &'a [f64],
// ==================== 导数 ====================
pub dabt0: &'a [f64],
pub dabtm: &'a [f64],
pub dabtp: &'a [f64],
pub dabn0: &'a [f64],
pub dabnm: &'a [f64],
pub dabnp: &'a [f64],
pub dabm0: &'a [f64],
pub dabmm: &'a [f64],
pub dabmp: &'a [f64],
pub demt0: &'a [f64],
pub demtm: &'a [f64],
pub demtp: &'a [f64],
pub demn0: &'a [f64],
pub demnm: &'a [f64],
pub demnp: &'a [f64],
pub demm0: &'a [f64],
pub demmm: &'a [f64],
pub demmp: &'a [f64],
// ==================== 能级导数 ====================
pub drch0: &'a [Vec<f64>],
pub drchm: &'a [Vec<f64>],
pub drchp: &'a [Vec<f64>],
pub dret0: &'a [Vec<f64>],
pub dretm: &'a [Vec<f64>],
pub dretp: &'a [Vec<f64>],
// ==================== 边界条件 ====================
pub fh: &'a [f64],
pub fhd: &'a [f64],
pub hextrd: &'a [f64],
pub q0: &'a [f64],
pub uu0: &'a [f64],
// ==================== 散射参数 ====================
pub sigec: &'a [f64],
pub dst: f64,
pub dsn: f64,
// ==================== 物理常量 ====================
pub hk: f64,
pub bn: f64,
pub rrdil: f64,
pub sige: f64,
}
/// BRTE 可变状态。
#[derive(Debug)]
pub struct BrteState<'a> {
pub a: &'a mut [Vec<f64>],
pub b: &'a mut [Vec<f64>],
pub c: &'a mut [Vec<f64>],
pub vecl: &'a mut [f64],
}
// ============================================================================
// Compton 辅助计算
// ============================================================================
/// 简化版 Compton 计算(与 compt0_bre 类似)。
fn compt0_brte(
ijt: usize,
id: usize,
ab: f64,
nfreq: usize,
kijt: &[usize],
elec: &[f64],
sige: f64,
) -> (f64, f64, f64, f64, f64, f64) {
let ijt_idx = ijt - 1;
let id_idx = id - 1;
let iji = nfreq - kijt[ijt_idx] + 1;
if iji == 1 {
return (0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
}
let ss0 = elec[id_idx] * sige / ab;
(0.0, 0.0, 0.0, 0.0, ss0, 0.0)
}
// ============================================================================
// BRTE 主函数
// ============================================================================
/// 计算辐射转移方程的矩阵 A, B, C 部分。
pub fn brte(params: &mut BrteParams, state: &mut BrteState) {
if params.nfreqe <= 0 {
return;
}
let id = params.id;
let id_idx = id - 1;
// 保存并修改 ISPLIN
let ispl = params.isplin;
let isplin = if ispl >= 5 { ispl - 5 } else { ispl };
params.isplin = isplin;
// 计算矩阵列索引
let nhe = params.nfreqe + params.inhe;
let nre = params.nfreqe + params.inre;
let npc = params.nfreqe + params.inpc;
let nse = params.nfreqe + params.inse - 1;
let nmp = params.nfreqe + params.inmp;
// GP, GN 几何因子
let (gp, gn) = if params.inmp > 0 { (UN, 0.0) } else { (0.0, UN) };
// IJ1 起始索引
let mut ij1 = 1;
if params.icompt > 0 && params.icombc > 0 && !params.ijex.is_empty() && params.ijex[0] > 0 {
ij1 = 2;
// Compton 边界条件(最高频率)
brte_compton_boundary(params, state, id_idx);
}
// ==================== ID = 1 上边界条件 ====================
if id > 1 {
brte_internal(params, state, id, id_idx, ij1, nhe, nre, npc, nmp, nse, gn, gp, isplin);
} else {
brte_upper_boundary(params, state, id_idx, ij1, nhe, nre, npc, nmp, nse, gn, gp, isplin);
}
// 恢复 ISPLIN
params.isplin = ispl;
// ==================== 低强度辐射场置零 ====================
if params.radzer > 0.0 {
brte_zero_low_intensity(params, state, id_idx, ij1);
}
}
/// Compton 边界条件(最高频率)。
fn brte_compton_boundary(params: &BrteParams, state: &mut BrteState, id_idx: usize) {
let ij = 1;
let iji = params.nfreq;
let zj1 = (-params.hk * params.freq[0] / params.temp[id_idx]).exp();
let zj2 = (-params.hk * params.freq[1] / params.temp[id_idx]).exp();
let dlt = params.delj[iji - 2][id_idx];
let (combid, comaid) = if params.ichcoo == 0 {
let zj0 = UN / (params.hk * (params.freq[0] * params.freq[1]).sqrt() / params.temp[id_idx]);
let zxx = UN - 3.0 * zj0 + (UN - dlt) * zj1 + dlt * zj2;
let combid = zj0 / params.dlnfr[iji - 2] + (UN - dlt) * zxx;
let comaid = -zj0 / params.dlnfr[iji - 2] + dlt * zxx;
(combid, comaid)
} else {
let e2 = YCON * params.temp[id_idx];
let zxx0 = XCON * params.freq[0] * (UN + zj1) - 3.0 * e2;
let zxxm = XCON * params.freq[1] * (UN + zj2) - 3.0 * e2;
let zxx = (UN - dlt) * zxx0 + dlt * zxxm;
let combid = e2 / params.dlnfr[iji - 2] + (UN - dlt) * zxx;
let comaid = -e2 / params.dlnfr[iji - 2] + dlt * zxx;
(combid, comaid)
};
state.b[0][0] = combid;
state.b[0][1] = comaid;
// 注:需要 rad 数组来计算 vecl,这里简化
}
/// 上边界条件(ID = 1)。
fn brte_upper_boundary(
params: &mut BrteParams,
state: &mut BrteState,
id_idx: usize,
ij1: usize,
nhe: usize,
nre: usize,
npc: usize,
nmp: usize,
nse: usize,
gn: f64,
gp: f64,
isplin: i32,
) {
let ddp = (params.dm[1] - params.dm[0]) * HALF;
for ij in ij1..=params.nfreqe {
let ij_idx = ij - 1;
let ijt = params.ijfr[ij_idx];
let omeg0 = params.abso0[ij_idx] / params.dens[id_idx];
let omegp = params.absop[ij_idx] / params.dens[id_idx + 1];
let dzp = omeg0 + omegp;
let dtaup = dzp * ddp;
let alf1 = (params.fk0[ij_idx] * params.rad0[ij_idx]
- params.fkp[ij_idx] * params.radp[ij_idx])
/ dtaup;
let chiel0 = params.scat0[ij_idx];
let chielp = params.scatp[ij_idx];
let mut s0 = (params.emis0[ij_idx] + chiel0 * params.rad0[ij_idx]) / params.abso0[ij_idx];
let mut bs = HALF * dtaup;
let mut cs = 0.0;
let mut c2 = 0.0;
let mut gam2 = 0.0;
let mut sp = 0.0;
// Compton 项
if params.icompt > 0 {
let (_cma, _cmb, _cmc, _cme, cms, _cmd) =
compt0_brte(ijt, params.id, params.abso0[ij_idx], params.nfreq, params.kijt, params.elec, params.sige);
s0 += cms;
}
// Spline/Hermitian 方法
if isplin % 3 > 0 {
bs = dtaup * THIRD;
cs = HALF * bs;
sp = (params.emisp[ij_idx] + chielp * params.radp[ij_idx]) / params.absop[ij_idx];
c2 = cs / params.absop[ij_idx];
gam2 = cs * (params.radp[ij_idx] - sp);
}
// 辅助量
let alf2 = bs * (params.rad0[ij_idx] - s0);
let bet2 = alf2 + gam2;
let x1 = (alf1 - bet2) / dzp;
let b2 = (bs + params.q0[ijt - 1]) / params.abso0[ij_idx];
let mut b1 = x1 / params.dens[0] + params.uu0[ijt - 1] * s0 * params.dm[0] * HALF / params.dens[0];
let mut c1 = x1 / params.dens[1];
// 矩阵 B 元素
let rtn = omeg0 * params.wmm[0] * b1;
state.b[ij_idx][nhe - 1] = -gn * rtn;
b1 -= b2 * s0;
let rtnc = omegp * params.wmm[1] * c1;
state.c[ij_idx][nhe - 1] = -gn * rtnc;
c1 -= c2 * sp;
// 温度、电子密度、质量列
state.b[ij_idx][nre - 1] =
b1 * params.dabt0[ij_idx] + b2 * (params.demt0[ij_idx] + params.dst * params.rad0[ij_idx]);
state.c[ij_idx][nre - 1] =
c1 * params.dabtp[ij_idx] + c2 * (params.demtp[ij_idx] + params.dst * params.radp[ij_idx]);
let sigec_val = if ijt > 0 && ijt <= params.sigec.len() { params.sigec[ijt - 1] } else { 0.0 };
state.b[ij_idx][npc - 1] = b1 * params.dabn0[ij_idx]
+ b2 * (params.demn0[ij_idx] + (params.dsn + sigec_val) * params.rad0[ij_idx])
+ gn * rtn;
state.c[ij_idx][npc - 1] = c1 * params.dabnp[ij_idx]
+ c2 * (params.demnp[ij_idx] + (params.dsn + sigec_val) * params.radp[ij_idx])
+ gn * rtnc;
state.b[ij_idx][nmp - 1] = b1 * params.dabm0[ij_idx] + b2 * params.demm0[ij_idx] - gp * rtn;
state.c[ij_idx][nmp - 1] = c1 * params.dabmp[ij_idx] + c2 * params.demmp[ij_idx] - gp * rtnc;
// 能级列
for ii in 0..params.nlvexp {
if nse + ii < state.b[ij_idx].len() {
state.b[ij_idx][nse + ii] +=
b1 * params.drch0[ii][ij_idx] + b2 * params.dret0[ii][ij_idx];
}
if nse + ii < state.c[ij_idx].len() {
state.c[ij_idx][nse + ii] +=
c1 * params.drchp[ii][ij_idx] + c2 * params.dretp[ii][ij_idx];
}
}
// 对角元素
state.b[ij_idx][params.nfreqe - 1] = 0.0;
state.b[ij_idx][ij - 1] = -params.fk0[ij_idx] / dtaup
- params.fh[ijt - 1]
- bs * (UN - chiel0 / params.abso0[ij_idx])
+ params.q0[ijt - 1] * chiel0 / params.abso0[ij_idx];
state.c[ij_idx][params.nfreqe - 1] = 0.0;
state.c[ij_idx][ij - 1] = params.fkp[ij_idx] / dtaup - cs * (UN - chielp / params.absop[ij_idx]);
// RHS
state.vecl[ij_idx] = alf1 + bet2 + params.fh[ijt - 1] * params.rad0[ij_idx] - s0 * params.q0[ijt - 1];
if params.iwinbl < 0 {
state.vecl[ij_idx] -= params.hextrd[ijt - 1];
}
// Compton 附加项
if params.icompt > 4 {
brte_compton_terms(params, state, ij_idx, ijt, bs, nre, npc);
}
}
}
/// 内部深度点(1 < ID < ND)。
fn brte_internal(
params: &mut BrteParams,
state: &mut BrteState,
id: usize,
id_idx: usize,
ij1: usize,
nhe: usize,
nre: usize,
npc: usize,
nmp: usize,
nse: usize,
gn: f64,
gp: f64,
isplin: i32,
) {
let ddm = (params.dm[id_idx] - params.dm[id_idx - 1]) * HALF;
// 检查是否是下边界
if id == params.nd {
brte_lower_boundary(params, state, id, id_idx, ij1, nhe, nre, npc, nmp, nse, gn, gp, ddm);
return;
}
let ddp = (params.dm[id_idx + 1] - params.dm[id_idx]) * HALF;
for ij in ij1..=params.nfreqe {
let ij_idx = ij - 1;
let ijt = params.ijfr[ij_idx];
let omeg0 = params.abso0[ij_idx] / params.dens[id_idx];
let omegp = params.absop[ij_idx] / params.dens[id_idx + 1];
let omegm = params.absom[ij_idx] / params.dens[id_idx - 1];
let dzp = omeg0 + omegp;
let dzm = omeg0 + omegm;
let dtaup = dzp * ddp;
let dtaum = dzm * ddm;
let dtau0 = HALF * (dtaup + dtaum);
let frd = params.fk0[ij_idx] * params.rad0[ij_idx];
let alf1 = (frd - params.fkp[ij_idx] * params.radp[ij_idx]) / dtaup / dtau0;
let gam1 = (frd - params.fkm[ij_idx] * params.radm[ij_idx]) / dtaum / dtau0;
let bet1 = alf1 + gam1;
let x1 = HALF * bet1 / dtau0;
let mut a1 = (gam1 + x1 * dtaum) / dzm;
let mut c1 = (alf1 + x1 * dtaup) / dzp;
let mut b1 = (a1 + c1) / params.dens[id_idx];
a1 /= params.dens[id_idx - 1];
c1 /= params.dens[id_idx + 1];
let mut bs = UN;
let chielm = params.scatm[ij_idx];
let chiel0 = params.scat0[ij_idx];
let chielp = params.scatp[ij_idx];
let mut s0 = (params.emis0[ij_idx] + chiel0 * params.rad0[ij_idx]) / params.abso0[ij_idx];
let mut as_s = 0.0;
let mut cs_s = 0.0;
let mut a2 = 0.0;
let mut c2 = 0.0;
let mut bet2 = 0.0;
let mut sm = 0.0;
let mut sp = 0.0;
// Compton 项
if params.icompt > 0 {
let (_cma, _cmb, _cmc, _cme, cms, _cmd) =
compt0_brte(ijt, id, params.abso0[ij_idx], params.nfreq, params.kijt, params.elec, params.sige);
s0 += cms;
}
// Spline/Hermitian 方法
if isplin % 3 > 0 {
sm = (params.emism[ij_idx] + params.radm[ij_idx] * chielm) / params.absom[ij_idx];
sp = (params.emisp[ij_idx] + params.radp[ij_idx] * chielp) / params.absop[ij_idx];
if isplin == 1 {
// Spline collocation
as_s = dtaum / dtau0 * SIXTH;
cs_s = dtaup / dtau0 * SIXTH;
bs = 0.666666666666667;
let alf2 = as_s * (params.radm[ij_idx] - sm);
let gam2_s = cs_s * (params.radp[ij_idx] - sp);
bet2 = alf2 + gam2_s;
let x = HALF * bet2 / dtau0;
a2 = (gam2_s - x * dtaum) / dzm;
c2 = (alf2 - x * dtaup) / dzp;
} else {
// Hermitian
let as_h = dtaup * dtaup / dtaum / dtau0;
let cs_h = dtaum * dtaum / dtaup / dtau0;
let al3 = (params.radp[ij_idx] - sp - params.rad0[ij_idx] + s0) * SIXTH;
let ga3 = (params.radm[ij_idx] - sm - params.rad0[ij_idx] + s0) * SIXTH;
let av = al3 * cs_h;
let cv = ga3 * as_h;
as_s = (UN - HALF * as_h) * SIXTH;
cs_s = (UN - HALF * cs_h) * SIXTH;
bs = UN - as_s - cs_s;
let x = (av + cv) / dtau0 / 4.0;
a2 = (x * dtaum + HALF * cv - av) / dzm;
c2 = (x * dtaup + HALF * av - cv) / dzp;
bet2 = as_s * (params.radm[ij_idx] - sm) + cs_s * (params.radp[ij_idx] - sp);
}
b1 -= (a2 + c2) / params.dens[id_idx];
a1 -= a2 / params.dens[id_idx - 1];
c1 -= c2 / params.dens[id_idx + 1];
}
let a2_abs = as_s / params.absom[ij_idx];
let c2_abs = cs_s / params.absop[ij_idx];
let a3 = a2_abs * sm;
let c3 = c2_abs * sp;
let b2 = bs / params.abso0[ij_idx];
let b3 = b2 * s0;
// 矩阵元素
let rtna = omegm * params.wmm[id_idx - 1] * a1;
state.a[ij_idx][nhe - 1] = -gn * rtna;
let a1_adj = a1 - a3;
let rtn = omeg0 * params.wmm[id_idx] * b1;
state.b[ij_idx][nhe - 1] = -gn * rtn;
let b1_adj = b1 - b3;
let rtnc = omegp * params.wmm[id_idx + 1] * c1;
state.c[ij_idx][nhe - 1] = -gn * rtnc;
let c1_adj = c1 - c3;
let sigec_val = if ijt > 0 && ijt <= params.sigec.len() { params.sigec[ijt - 1] } else { 0.0 };
state.a[ij_idx][nre - 1] =
a1_adj * params.dabtm[ij_idx] + a2_abs * (params.demtm[ij_idx] + params.dst * params.radm[ij_idx]);
state.b[ij_idx][nre - 1] =
b1_adj * params.dabt0[ij_idx] + b2 * (params.demt0[ij_idx] + params.dst * params.rad0[ij_idx]);
state.c[ij_idx][nre - 1] =
c1_adj * params.dabtp[ij_idx] + c2_abs * (params.demtp[ij_idx] + params.dst * params.radp[ij_idx]);
state.a[ij_idx][npc - 1] = a1_adj * params.dabnm[ij_idx]
+ a2_abs * (params.demnm[ij_idx] + (params.dsn + sigec_val) * params.radm[ij_idx])
+ gn * rtna;
state.b[ij_idx][npc - 1] = b1_adj * params.dabn0[ij_idx]
+ b2 * (params.demn0[ij_idx] + (params.dsn + sigec_val) * params.rad0[ij_idx])
+ gn * rtn;
state.c[ij_idx][npc - 1] = c1_adj * params.dabnp[ij_idx]
+ c2_abs * (params.demnp[ij_idx] + (params.dsn + sigec_val) * params.radp[ij_idx])
+ gn * rtnc;
state.a[ij_idx][nmp - 1] = a1_adj * params.dabmm[ij_idx] + a2_abs * params.demmm[ij_idx] - gp * rtna;
state.b[ij_idx][nmp - 1] = b1_adj * params.dabm0[ij_idx] + b2 * params.demm0[ij_idx] - gp * rtn;
state.c[ij_idx][nmp - 1] = c1_adj * params.dabmp[ij_idx] + c2_abs * params.demmp[ij_idx] - gp * rtnc;
// 能级列
for ii in 0..params.nlvexp {
if nse + ii < state.a[ij_idx].len() {
state.a[ij_idx][nse + ii] += a1_adj * params.drchm[ii][ij_idx] + a2_abs * params.dretm[ii][ij_idx];
}
if nse + ii < state.b[ij_idx].len() {
state.b[ij_idx][nse + ii] += b1_adj * params.drch0[ii][ij_idx] + b2 * params.dret0[ii][ij_idx];
}
if nse + ii < state.c[ij_idx].len() {
state.c[ij_idx][nse + ii] += c1_adj * params.drchp[ii][ij_idx] + c2_abs * params.dretp[ii][ij_idx];
}
}
// 对角元素
state.a[ij_idx][params.nfreqe - 1] = 0.0;
state.a[ij_idx][ij - 1] = params.fkm[ij_idx] / dtaum / dtau0 - as_s * (UN - chielm / params.absom[ij_idx]);
state.b[ij_idx][params.nfreqe - 1] = 0.0;
state.b[ij_idx][ij - 1] = -params.fk0[ij_idx] / dtau0 * (UN / dtaup + UN / dtaum)
- bs * (UN - chiel0 / params.abso0[ij_idx]);
state.c[ij_idx][params.nfreqe - 1] = 0.0;
state.c[ij_idx][ij - 1] = params.fkp[ij_idx] / dtaup / dtau0 - cs_s * (UN - chielp / params.absop[ij_idx]);
// RHS
state.vecl[ij_idx] = bet1 + bet2 + bs * (params.rad0[ij_idx] - s0);
// Compton 附加项
if params.icompt > 4 {
brte_compton_terms(params, state, ij_idx, ijt, bs, nre, npc);
}
}
}
/// 下边界条件(ID = ND)。
fn brte_lower_boundary(
params: &mut BrteParams,
state: &mut BrteState,
id: usize,
id_idx: usize,
ij1: usize,
nhe: usize,
nre: usize,
npc: usize,
nmp: usize,
nse: usize,
gn: f64,
gp: f64,
ddm: f64,
) {
// 盘模型特殊处理
if params.idisk != 0 && params.ifz0 >= 0 {
brte_lower_disk(params, state, id, id_idx, ij1, nhe, nre, npc, nmp, nse, gn, gp, ddm);
return;
}
let t = if params.tempbd != 0.0 { params.tempbd } else { params.temp[id_idx] };
let tm = if params.tempbd != 0.0 { params.tempbd } else { params.temp[id_idx - 1] };
let hkt = params.hk / t;
let hktm = params.hk / tm;
for ij in ij1..=params.nfreqe {
let ij_idx = ij - 1;
let ijt = params.ijfr[ij_idx];
let chielm = params.scatm[ij_idx];
let chiel0 = params.scat0[ij_idx];
let omegm = params.absom[ij_idx] / params.dens[id_idx - 1];
let omeg0 = params.abso0[ij_idx] / params.dens[id_idx];
let dzm = omeg0 + omegm;
let dtaum = dzm * ddm;
let frd = params.fk0[ij_idx] * params.rad0[ij_idx] - params.fkm[ij_idx] * params.radm[ij_idx];
let mut gam1 = frd / dtaum;
let mut a1 = gam1 / dzm;
let mut as_s = 0.0;
let mut bs = 0.0;
let mut a2 = 0.0;
let mut b2 = 0.0;
let mut a3 = 0.0;
let mut b3 = 0.0;
let mut bet2 = 0.0;
let mut s0 = 0.0;
// 二阶边界条件
if params.ibc > 0 && params.ibc < 4 {
bs = dtaum * HALF;
s0 = (params.emis0[ij_idx] + chiel0 * params.rad0[ij_idx]) / params.abso0[ij_idx];
if params.icompt > 0 {
let (_cma, _cmb, _cmc, _cme, cms, _cmd) =
compt0_brte(ijt, id, params.abso0[ij_idx], params.nfreq, params.kijt, params.elec, params.sige);
s0 += cms;
}
let gam2 = bs * (params.rad0[ij_idx] - s0);
bet2 = gam2;
let x1 = bet2 / dzm;
a1 -= x1;
b2 = bs / params.abso0[ij_idx];
b3 = b2 * s0;
}
// Planck 函数
let fr = params.freq[ijt - 1];
let fr15 = fr * 1e-15;
let x = hkt * fr;
let ex = x.exp();
let xm = hktm * fr;
let exm = xm.exp();
let plan = params.bn * fr15 * fr15 * fr15 / (ex - UN) * params.rrdil;
// planm: 在 Fortran 中是循环内保持的变量,需要在此处初始化
let mut planm = params.bn * fr15 * fr15 * fr15 / (exm - UN) * params.rrdil;
if params.inre_idx == 0 || id >= params.ndre {
planm = params.bn * fr15 * fr15 * fr15 / (exm - UN) * params.rrdil;
let gam3 = (plan - planm) / dtaum * THIRD;
a1 -= gam3 / dzm;
gam1 -= gam3;
}
let c1 = a1;
let a1_dens = c1 / params.dens[id_idx - 1];
let b1 = c1 / params.dens[id_idx];
// 矩阵元素
let rtna = omegm * params.wmm[id_idx - 1] * a1_dens;
state.a[ij_idx][nhe - 1] = -gn * rtna;
let rtn = omeg0 * params.wmm[id_idx] * b1;
state.b[ij_idx][nhe - 1] = -gn * rtn;
let b1_adj = b1 - b3;
let sigec_val = if ijt > 0 && ijt <= params.sigec.len() { params.sigec[ijt - 1] } else { 0.0 };
let dplanm = planm * xm / tm / (UN - UN / exm);
state.a[ij_idx][nre - 1] = (a1_dens - a3) * params.dabtm[ij_idx]
+ a2 * (params.demtm[ij_idx] + params.dst * params.radm[ij_idx])
- dplanm / dtaum * THIRD;
let bb = HALF + THIRD / dtaum;
let dplan = plan * x / t / (UN - UN / ex);
state.b[ij_idx][nre - 1] = b1_adj * params.dabt0[ij_idx]
+ b2 * (params.demt0[ij_idx] + params.dst * params.rad0[ij_idx])
+ bb * dplan;
state.a[ij_idx][npc - 1] = (a1_dens - a3) * params.dabnm[ij_idx]
+ a2 * (params.demnm[ij_idx] + (params.dsn + sigec_val) * params.radm[ij_idx])
+ gn * rtna;
state.b[ij_idx][npc - 1] = b1_adj * params.dabn0[ij_idx]
+ b2 * (params.demn0[ij_idx] + (params.dsn + sigec_val) * params.rad0[ij_idx])
+ gn * rtn;
state.a[ij_idx][nmp - 1] = (a1_dens - a3) * params.dabmm[ij_idx] + a2 * params.demmm[ij_idx] - gp * rtna;
state.b[ij_idx][nmp - 1] = b1_adj * params.dabm0[ij_idx] + b2 * params.demm0[ij_idx] - gp * rtn;
// 能级列
for ii in 0..params.nlvexp {
if nse + ii < state.a[ij_idx].len() {
state.a[ij_idx][nse + ii] += (a1_dens - a3) * params.drchm[ii][ij_idx] + a2 * params.dretm[ii][ij_idx];
}
if nse + ii < state.b[ij_idx].len() {
state.b[ij_idx][nse + ii] += b1_adj * params.drch0[ii][ij_idx] + b2 * params.dret0[ii][ij_idx];
}
}
// 对角元素和 RHS
state.a[ij_idx][params.nfreqe - 1] = 0.0;
state.a[ij_idx][ij - 1] = params.fkm[ij_idx] / dtaum - as_s * (UN - chielm / params.absom[ij_idx]);
state.b[ij_idx][params.nfreqe - 1] = 0.0;
if params.ibc == 0 || params.ibc == 4 {
state.b[ij_idx][ij - 1] = -params.fk0[ij_idx] / dtaum
- bs * (UN - chiel0 / params.abso0[ij_idx])
- HALF;
state.vecl[ij_idx] = gam1 + bet2 - HALF * (plan - params.rad0[ij_idx]);
} else {
state.b[ij_idx][ij - 1] = -params.fk0[ij_idx] / dtaum
- bs * (UN - chiel0 / params.abso0[ij_idx])
- params.fhd[ijt - 1];
state.vecl[ij_idx] = gam1 + bet2 - HALF * plan + params.fhd[ijt - 1] * params.rad0[ij_idx];
}
// Compton 附加项
if params.icompt > 4 {
brte_compton_terms(params, state, ij_idx, ijt, bs, nre, npc);
}
}
}
/// 下边界条件(盘模型)。
fn brte_lower_disk(
params: &mut BrteParams,
state: &mut BrteState,
id: usize,
id_idx: usize,
ij1: usize,
nhe: usize,
nre: usize,
npc: usize,
nmp: usize,
nse: usize,
gn: f64,
gp: f64,
ddm: f64,
) {
for ij in ij1..=params.nfreqe {
let ij_idx = ij - 1;
let ijt = params.ijfr[ij_idx];
let chielm = params.scatm[ij_idx];
let chiel0 = params.scat0[ij_idx];
let omegm = params.absom[ij_idx] / params.dens[id_idx - 1];
let omeg0 = params.abso0[ij_idx] / params.dens[id_idx];
let dzm = omeg0 + omegm;
let dtaum = dzm * ddm;
let frd = params.fk0[ij_idx] * params.rad0[ij_idx] - params.fkm[ij_idx] * params.radm[ij_idx];
let gam1 = frd / dtaum;
let mut a1 = gam1 / dzm;
let bs = dtaum * HALF;
let s0 = (params.emis0[ij_idx] + chiel0 * params.rad0[ij_idx]) / params.abso0[ij_idx];
let gam2 = bs * (params.rad0[ij_idx] - s0);
let bet2 = gam2;
let x1 = bet2 / dzm;
a1 -= x1;
let b2 = bs / params.abso0[ij_idx];
let b3 = b2 * s0;
let c1 = a1;
let a1_dens = c1 / params.dens[id_idx - 1];
let b1 = c1 / params.dens[id_idx];
// 矩阵元素
let rtn = omegm * params.wmm[id_idx] * a1_dens;
state.a[ij_idx][nhe - 1] = -gn * rtn;
state.a[ij_idx][nmp - 1] = -gp * rtn;
let rtn_b = omeg0 * params.wmm[id_idx] * b1;
state.b[ij_idx][nhe - 1] = -gn * rtn_b;
state.b[ij_idx][nmp - 1] = -gp * rtn_b;
let sigec_val = if ijt > 0 && ijt <= params.sigec.len() { params.sigec[ijt - 1] } else { 0.0 };
state.a[ij_idx][nre - 1] = a1_dens * params.dabtm[ij_idx];
state.a[ij_idx][npc - 1] = a1_dens * params.dabnm[ij_idx]
+ gn * rtn;
state.b[ij_idx][nre - 1] = (b1 - b3) * params.dabt0[ij_idx] + b2 * (params.demt0[ij_idx] + params.dst * params.rad0[ij_idx]);
state.b[ij_idx][npc - 1] = (b1 - b3) * params.dabn0[ij_idx]
+ b2 * (params.demn0[ij_idx] + (params.dsn + sigec_val) * params.rad0[ij_idx])
+ gn * rtn_b;
state.a[ij_idx][nmp - 1] = a1_dens * params.dabmm[ij_idx] - gp * rtn;
state.b[ij_idx][nmp - 1] = (b1 - b3) * params.dabm0[ij_idx] + b2 * params.demm0[ij_idx] - gp * rtn_b;
// 对角元素
state.a[ij_idx][params.nfreqe - 1] = 0.0;
state.a[ij_idx][ij - 1] = params.fkm[ij_idx] / dtaum;
state.b[ij_idx][params.nfreqe - 1] = 0.0;
state.b[ij_idx][ij - 1] = -params.fk0[ij_idx] / dtaum - bs * (UN - chiel0 / params.abso0[ij_idx]);
state.vecl[ij_idx] = gam1 + bet2;
// Compton 附加项
if params.icompt > 4 {
brte_compton_terms(params, state, ij_idx, ijt, bs, nre, npc);
}
}
}
/// Compton 附加项。
fn brte_compton_terms(
params: &BrteParams,
state: &mut BrteState,
ij_idx: usize,
ijt: usize,
bs: f64,
nre: usize,
npc: usize,
) {
let (cma, cmb, cmc, _cme, cms, cmd) =
compt0_brte(ijt, params.id, params.abso0[ij_idx], params.nfreq, params.kijt, params.elec, params.sige);
state.b[ij_idx][ij_idx] += bs * (cmb + cms);
let iji = params.nfreq - params.kijt[ijt - 1] + 1;
if iji > 1 {
let ijm = params.ijex[params.ijorig[iji - 2]] as usize;
if ijm > 0 && ijm - 1 < state.b[ij_idx].len() {
state.b[ij_idx][ijm - 1] += bs * cma;
}
}
if iji < params.nfreq {
let ijp = params.ijex[params.ijorig[iji]] as usize;
if ijp > 0 && ijp - 1 < state.b[ij_idx].len() {
state.b[ij_idx][ijp - 1] += bs * cmc;
}
}
if params.inre_idx > 0 {
state.b[ij_idx][nre - 1] += cmd * bs;
}
if params.inpc_idx > 0 {
state.b[ij_idx][npc - 1] += cms * bs / params.elec[params.id - 1];
}
}
/// 低强度辐射场置零。
fn brte_zero_low_intensity(
params: &BrteParams,
state: &mut BrteState,
id_idx: usize,
ij1: usize,
) {
// 找到 nu*rad_nu 的峰值
let mut radsum = 0.0;
for ij in ij1..=params.nfreqe {
let ij_idx = ij - 1;
let val = params.freq[ij_idx] * params.radex[ij_idx][id_idx];
if val > radsum {
radsum = val;
}
}
// 如果远小于峰值,则置零
for ij in ij1..=params.nfreqe {
let ij_idx = ij - 1;
if params.freq[ij_idx] * params.radex[ij_idx][id_idx] < params.radzer * radsum {
for ii in 0..state.a[ij_idx].len() {
state.a[ij_idx][ii] = 0.0;
state.b[ij_idx][ii] = 0.0;
state.c[ij_idx][ii] = 0.0;
}
state.vecl[ij_idx] = 0.0;
state.b[ij_idx][ij_idx] = UN;
}
}
}
// ============================================================================
// 测试
// ============================================================================
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_brte_compile() {
assert!(true);
}
#[test]
fn test_brte_constants() {
assert!((XCON - 8.0935e-21).abs() < 1e-35);
assert!((YCON - 1.68638e-10).abs() < 1e-20);
assert!((SIXTH - 1.0 / 6.0).abs() < 1e-15);
assert!((THIRD - 1.0 / 3.0).abs() < 1e-15);
}
#[test]
fn test_compt0_brte_iji_1() {
let kijt = vec![100];
let elec = vec![1e12];
let result = compt0_brte(1, 1, 1e-8, 100, &kijt, &elec, 6.6516e-25);
assert!((result.0).abs() < 1e-15);
assert!((result.1).abs() < 1e-15);
assert!((result.2).abs() < 1e-15);
}
#[test]
fn test_brte_matrix_indices() {
let nfreqe = 10;
let inhe = 1;
let inre = 2;
let inpc = 3;
let inmp = 4;
let inse = 1;
assert_eq!(nfreqe + inhe, 11); // NHE
assert_eq!(nfreqe + inre, 12); // NRE
assert_eq!(nfreqe + inpc, 13); // NPC
assert_eq!(nfreqe + inmp, 14); // NMP
assert_eq!(nfreqe + inse - 1, 10); // NSE
}
}
+1005
View File
File diff suppressed because it is too large Load Diff
+632
View File
@@ -0,0 +1,632 @@
//! 氢碰撞速率计算。
//!
//! 重构自 TLUSTY `colh.f`
//!
//! 计算氢的碰撞电离和碰撞激发速率。
//! 标准表达式来自 Mihalas, Heasley, and Auer (1975)。
use super::butler::butler;
use super::ceh12::ceh12;
use super::cspec::cspec;
use super::irc::irc;
use crate::data::{COLH_CCOOL, COLH_CHOT};
use crate::state::constants::{EH, HK, TWO, UN};
// 物理常量
const CC0: f64 = 5.465e-11;
const CEX1: f64 = -30.20581;
const CEX2: f64 = 3.8608704;
const CEX3: f64 = 305.63574;
const ALF0: f64 = 1.8;
const ALF1: f64 = 0.4;
const BET0: f64 = 3.0;
const BET1: f64 = 1.2;
const O148: f64 = 0.148;
const CHMI: f64 = 5.59e-15;
// 指数积分系数
const EXPIA1: f64 = -0.57721566;
const EXPIA2: f64 = 0.99999193;
const EXPIA3: f64 = -0.24991055;
const EXPIA4: f64 = 0.05519968;
const EXPIA5: f64 = -0.00976004;
const EXPIA6: f64 = 0.00107857;
const EXPIB1: f64 = 0.2677734343;
const EXPIB2: f64 = 8.6347608925;
const EXPIB3: f64 = 18.059016973;
const EXPIB4: f64 = 8.5733287401;
const EXPIC1: f64 = 3.9584969228;
const EXPIC2: f64 = 21.0996530827;
const EXPIC3: f64 = 25.6329561486;
const EXPIC4: f64 = 9.5733223454;
/// COLH 输入参数
pub struct ColhParams {
/// 深度索引 (1-indexed)
pub id: usize,
/// 温度 (K)
pub t: f64,
/// 碰撞速率选择标志
pub icolhn: i32,
}
/// COLH 原子数据
pub struct ColhAtomicData<'a> {
/// 氢元素索引
pub ielh: usize,
/// H- 元素索引 (0 表示没有 H-)
pub ielhm: usize,
/// 氢原子索引
pub iath: usize,
/// 能级的第一个索引 (nelem)
pub nfirst: &'a [i32],
/// 能级的最后一个索引 (nelem)
pub nlast: &'a [i32],
/// 电离能级索引 (natom)
pub nka: &'a [i32],
/// 主量子数 (nlevel)
pub nquant: &'a [i32],
/// 上能级截止 (nelem)
pub icup: &'a [i32],
/// 跃迁索引数组 (nlevel × nlevel)
pub itra: &'a [i32],
/// 碰撞速率标志 (ntrans)
pub icol: &'a [i32],
/// 跃迁频率 (ntrans)
pub fr0: &'a [f64],
/// 振子强度 (ntrans)
pub osc0: &'a [f64],
/// 碰撞参数 (ntrans)
pub cpar: &'a [f64],
/// 电离能 (nlevel)
pub enion: &'a [f64],
/// 能级宽度选项 (nlevel)
pub ifwop: &'a [i32],
/// WNHINT 数组 (nlmx × nd)
pub wnhint: &'a [f64],
/// OSH 振子强度 (10 × 20)
pub osh: &'a [f64],
/// 最大谱线数
pub nlmx: usize,
}
/// COLH 输出
pub struct ColhOutput<'a> {
/// 碰撞速率数组 (ntrans)
pub col: &'a mut [f64],
}
// A 系数数组(用于高 n 碰撞电离)
static A: [[f64; 10]; 6] = [
[-86.7633398, 2632.8369, 7478.9556, -4202.8442, -47995.930,
-120942.89, -202300.81, -261373.03, -266337.91, -192293.20],
[100.919188, -2738.7485, -8495.4590, 1937.3763, 45825.371,
122209.39, 211928.67, 285044.75, 309455.47, 258802.22],
[-45.7813807, 1121.3976, 3794.6826, 340.35764, -16617.055,
-47390.313, -84973.688, -117833.95, -133243.61, -120363.95],
[10.1978559, -224.30670, -822.83636, -290.10489, 2905.7393,
8944.6025, 16556.992, 23544.543, 27419.742, 26002.143],
[-1.11223557, 21.923729, 86.619110, 48.840523, -246.99014,
-828.41028, -1581.2722, -2297.9321, -2738.1743, -2686.4087],
[0.0474198818, -0.83974838, -3.5534720, -2.6097214, 8.1972208,
30.267115, 59.521984, 88.178680, 107.05288, 107.73775],
];
/// 计算氢的碰撞速率。
///
/// # 参数
///
/// * `params` - 输入参数
/// * `atomic` - 原子数据
/// * `output` - 输出碰撞速率
pub fn colh(params: &ColhParams, atomic: &ColhAtomicData, output: &mut ColhOutput) {
let t = params.t;
let id = params.id;
let id_idx = id - 1;
let hkt = HK / t;
let ct = CC0 * t.sqrt();
let tk = hkt / EH;
let x = t.log10();
let x2 = x * x;
let x3 = x * x2;
let x4 = x2 * x2;
let x5 = x3 * x2;
let sqt = t.sqrt();
let xtt = [1.0, t, t * t, t * t * t];
let n0hn = atomic.nfirst[atomic.ielh] as usize - 1;
let n1h = atomic.nlast[atomic.ielh] as usize - 1;
let nkh = atomic.nka[atomic.iath] as usize - 1;
let n1q = if atomic.icup[atomic.ielh] > 0 {
atomic.icup[atomic.ielh] as usize
} else {
0
};
let n1hc = if atomic.ifwop[n1h] < 0 { n1h - 1 } else { n1h };
let nhl = if n1q > 0 {
n1q
} else {
n1h - n0hn + 1
};
for ii in n0hn..=n1h {
let i = ii - n0hn + 1;
let it = atomic.itra[ii * (n1h + 1) + nkh] as usize;
if it == 0 {
continue;
}
// 碰撞电离
if t > 1e6 {
// 高温使用 XSTAR 公式
let rno = 16.0;
let izc = 1;
let cs = irc(i as i32, t, izc, rno);
output.col[it - 1] = cs;
} else {
let ic = atomic.icol[it - 1];
let u0 = atomic.fr0[it - 1] * hkt;
if ic < 0 {
// 非标准公式
output.col[it - 1] = cspec(
ii as i32, nkh as i32, ic,
atomic.osc0[it - 1], atomic.cpar[it - 1], u0, t
);
} else if atomic.ifwop[ii] < 0 {
// 合并态电离
let ehk = EH / tk;
let n00q = atomic.nquant[n1h - 1] as usize + 1;
let mut sum1 = 0.0;
let mut sum2 = 0.0;
for img in n00q..=atomic.nlmx {
let xi = img as f64;
let xii = xi * xi;
sum1 += xii * xii * xi * atomic.wnhint[img * id + id_idx];
sum2 += xii * atomic.wnhint[img * id + id_idx] * (ehk / xii).exp();
}
output.col[it - 1] = ct * sum1 / sum2;
} else {
// 标准公式
let gam = if i > 10 {
(i * i * i) as f64
} else {
A[0][i - 1] + A[1][i - 1] * x + A[2][i - 1] * x2
+ A[3][i - 1] * x3 + A[4][i - 1] * x4 + A[5][i - 1] * x5
};
output.col[it - 1] = ct * (-u0).exp() * gam;
}
}
// 碰撞激发
let i1 = i + 1;
if i1 > nhl {
continue;
}
let xi = i as f64;
let vi = xi * xi;
let alf = ALF0 - ALF1 / vi;
let bet = BET0 - BET1 / xi;
let csca = 8.63e-6 / 2.0 / vi / sqt;
let mut csum = 0.0;
for j in i1..nhl {
let xj = j as f64;
let vj = xj * xj;
let jj = j + n0hn;
let (cs, is_lumped) = if jj > n1hc {
// 非显式能级,归入碰撞电离
let e = UN / vi - UN / vj;
let u0 = EH * e * tk;
let c1 = if j <= 20 {
atomic.osh[i * 20 + j]
} else {
atomic.osh[i * 20 + 19] * ((400.0 - vi) / 20.0 * xj / (vj - vi)).powi(3)
};
(compute_collision_rate(params.icolhn, i, j, t, u0, c1, csca, alf, bet, xi, xj, &xtt), true)
} else {
let ict = atomic.itra[ii * (n1h + 1) + jj] as usize;
if ict == 0 {
continue;
}
let ic = atomic.icol[ict - 1];
let u0 = atomic.fr0[ict - 1] * hkt;
let c1 = atomic.osc0[ict - 1];
let cs = if ic < 0 {
cspec(ii as i32, jj as i32, ic, c1, atomic.cpar[ict - 1], u0, t)
} else if ic == 0 {
compute_collision_rate(params.icolhn, i, j, t, u0, c1, csca, alf, bet, xi, xj, &xtt)
} else if ic == 1 {
ct * (-u0).exp() * (CEX1 + CEX2 * x + CEX3 / x / x)
} else {
ceh12(t)
};
(cs, false)
};
if is_lumped {
csum += cs;
} else {
let ict = atomic.itra[ii * (n1h + 1) + jj] as usize;
if ict > 0 {
output.col[ict - 1] = cs;
}
}
}
// 添加累积碰撞速率
let it_main = atomic.itra[ii * (n1h + 1) + nkh] as usize;
if it_main > 0 && n1q > 0 {
output.col[it_main - 1] += csum;
}
let ith = atomic.itra[ii * (n1h + 1) + n1h] as usize;
if atomic.ifwop[n1h] < 0 && ith > 0 {
output.col[ith - 1] = csum;
}
}
// H- 碰撞电离
if atomic.ielhm > 0 {
let it_hm = atomic.itra[atomic.nfirst[atomic.ielhm] as usize * (n1h + 1) + n0hn] as usize;
if it_hm > 0 {
let ic = atomic.icol[it_hm - 1];
if ic >= 0 {
output.col[it_hm - 1] = CHMI * t * t.sqrt();
} else {
let u0 = atomic.enion[atomic.nfirst[atomic.ielhm] as usize - 1] * tk;
output.col[it_hm - 1] = cspec(
atomic.nfirst[atomic.ielhm] as i32,
n0hn as i32,
ic,
atomic.osc0[it_hm - 1],
atomic.cpar[it_hm - 1],
u0,
t,
);
}
}
}
}
/// 计算碰撞激发速率
fn compute_collision_rate(
icolhn: i32,
i: usize,
j: usize,
t: f64,
u0: f64,
c1: f64,
csca: f64,
alf: f64,
bet: f64,
xi: f64,
xj: f64,
xtt: &[f64],
) -> f64 {
// 使用 Butler 新计算
if icolhn == 1 && j <= 7 {
let (cs, ierr) = butler(i as i32, j as i32, t, u0);
if ierr == 0 {
return cs;
}
}
// Giovanardi 公式
if icolhn == 2 && j <= 15 {
let cs = if t <= 60000.0 {
get_ccool(i, j, xtt)
} else {
get_chot(i, j, xtt)
};
return csca * cs * (-u0).exp();
}
// 标准公式 (Mihalas et al 1975)
let ct = CC0 * t.sqrt();
let e = u0 / (HK / t) / EH;
let cs = 4.0 * ct * c1 / (e * e);
let ex = (-u0).exp();
let e1 = if u0 <= UN {
-u0.ln() + EXPIA1 + u0 * (EXPIA2 + u0 * (EXPIA3 + u0 * (EXPIA4 + u0 * (EXPIA5 + u0 * EXPIA6))))
} else {
let u0_inv = 1.0 / u0;
(-u0).exp() * ((EXPIB1 + u0 * (EXPIB2 + u0 * (EXPIB3 + u0 * EXPIB4)))
/ (EXPIC1 + u0 * (EXPIC2 + u0 * (EXPIC3 + u0 * EXPIC4))))
* u0_inv
};
let mut e5 = e1;
for ix in 1..=4 {
e5 = (ex - u0 * e5) / ix as f64;
}
let mut cs = cs * u0 * (e1 + O148 * u0 * e5);
if j - i != 1 {
cs *= bet + TWO * (alf - bet) / (xj - xi);
}
cs
}
/// 获取 CCOOL 系数(从 data.rs 获取)
/// CCOOL(I, J, K) - I=1..4, J=1..14, K=1..15
/// Fortran 列优先: idx = (I-1) + 4*((J-1) + 14*(K-1))
fn get_ccool(i: usize, j: usize, xtt: &[f64]) -> f64 {
// Giovanardi 公式: CS = sum_{ica=1}^{4} CCOOL(ica, i, j) * XTT(ica)
let mut cs = 0.0;
for ica in 0..4 {
// Fortran: CCOOL(ica+1, i, j)
// 列优先索引: (ica) + 4*((i-1) + 14*(j-1))
let idx = ica + 4 * ((i - 1) + 14 * (j - 1));
if idx < COLH_CCOOL.len() {
cs += COLH_CCOOL[idx] * xtt[ica];
}
}
cs
}
/// 获取 CHOT 系数(从 data.rs 获取)
/// CHOT(I, J, K) - I=1..4, J=1..14, K=1..15
/// Fortran 列优先: idx = (I-1) + 4*((J-1) + 14*(K-1))
fn get_chot(i: usize, j: usize, xtt: &[f64]) -> f64 {
// Giovanardi 公式: CS = sum_{ica=1}^{4} CHOT(ica, i, j) * XTT(ica)
let mut cs = 0.0;
for ica in 0..4 {
// Fortran: CHOT(ica+1, i, j)
// 列优先索引: (ica) + 4*((i-1) + 14*(j-1))
let idx = ica + 4 * ((i - 1) + 14 * (j - 1));
if idx < COLH_CHOT.len() {
cs += COLH_CHOT[idx] * xtt[ica];
}
}
cs
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_colh_constants() {
assert!((CC0 - 5.465e-11).abs() < 1e-20);
assert!((ALF0 - 1.8).abs() < 1e-10);
assert!((BET0 - 3.0).abs() < 1e-10);
}
#[test]
fn test_get_ccool_returns_data() {
// 测试 get_ccool 函数是否正确返回数据
let xtt = [1.0, 10000.0, 1e8, 1e12]; // T=10000K 的幂次
// 测试几个有效索引
for i in 1..=14 {
for j in 1..=15 {
let cs = get_ccool(i, j, &xtt);
// 返回值应该是有限数
assert!(cs.is_finite(), "get_ccool({}, {}) returned non-finite: {}", i, j, cs);
}
}
}
#[test]
fn test_get_chot_returns_data() {
// 测试 get_chot 函数是否正确返回数据
let xtt = [1.0, 100000.0, 1e10, 1e15]; // T=100000K 的幂次
// 测试几个有效索引
for i in 1..=14 {
for j in 1..=15 {
let cs = get_chot(i, j, &xtt);
// 返回值应该是有限数
assert!(cs.is_finite(), "get_chot({}, {}) returned non-finite: {}", i, j, cs);
}
}
}
#[test]
fn test_get_ccool_temperature_scaling() {
// 验证 CCOOL 系数随温度的缩放
let xtt_low = [1.0, 5000.0, 2.5e7, 1.25e11];
let xtt_high = [1.0, 50000.0, 2.5e9, 1.25e14];
for i in 1..=5 {
for j in 1..=5 {
let cs_low = get_ccool(i, j, &xtt_low);
let cs_high = get_ccool(i, j, &xtt_high);
// 高温应该有更大的系数(大部分情况)
// 这里只验证都是有限数
assert!(cs_low.is_finite() && cs_high.is_finite());
}
}
}
#[test]
fn test_compute_collision_rate_standard() {
// 测试标准 Mihalas 公式
let t: f64 = 10000.0;
let i = 1;
let j = 2;
let u0: f64 = 0.5; // 激发能量 / kT
let c1: f64 = 0.4162; // Lyman-alpha 振子强度
let csca = 8.63e-6 / 2.0 / 1.0 / t.sqrt();
let alf = ALF0 - ALF1 / 1.0;
let bet = BET0 - BET1 / 1.0;
let xi = 1.0;
let xj = 2.0;
let xtt = [1.0, t, t*t, t*t*t];
// 使用 icolhn=0 强制使用标准公式
let cs = compute_collision_rate(0, i, j, t, u0, c1, csca, alf, bet, xi, xj, &xtt);
// 碰撞速率应该是有限的(可能很小但应该有效)
assert!(cs.is_finite(), "Collision rate should be finite");
}
#[test]
fn test_compute_collision_rate_butler() {
// 测试 Butler 公式 (icolhn=1, j<=7)
let t: f64 = 10000.0;
let i = 1;
let j = 2;
let u0: f64 = 0.5;
let c1: f64 = 0.4162;
let csca = 8.63e-6 / 2.0 / 1.0 / t.sqrt();
let alf = ALF0 - ALF1 / 1.0;
let bet = BET0 - BET1 / 1.0;
let xi = 1.0;
let xj = 2.0;
let xtt = [1.0, t, t*t, t*t*t];
let cs = compute_collision_rate(1, i, j, t, u0, c1, csca, alf, bet, xi, xj, &xtt);
assert!(cs > 0.0, "Butler collision rate should be positive");
assert!(cs.is_finite(), "Butler collision rate should be finite");
}
#[test]
fn test_compute_collision_rate_giovanardi_cool() {
// 测试 Giovanardi 冷公式 (icolhn=2, j<=15, T<=60000K)
let t: f64 = 10000.0;
let i = 1;
let j = 2;
let u0: f64 = 0.5;
let c1: f64 = 0.4162;
let csca = 8.63e-6 / 2.0 / 1.0 / t.sqrt();
let alf = ALF0 - ALF1 / 1.0;
let bet = BET0 - BET1 / 1.0;
let xi = 1.0;
let xj = 2.0;
let xtt = [1.0, t, t*t, t*t*t];
let cs = compute_collision_rate(2, i, j, t, u0, c1, csca, alf, bet, xi, xj, &xtt);
assert!(cs > 0.0, "Giovanardi cool collision rate should be positive");
assert!(cs.is_finite(), "Giovanardi cool collision rate should be finite");
}
#[test]
fn test_compute_collision_rate_giovanardi_hot() {
// 测试 Giovanardi 热公式 (icolhn=2, j<=15, T>60000K)
let t: f64 = 100000.0;
let i = 1;
let j = 2;
let u0: f64 = 0.1; // 更小的激发能量比
let c1: f64 = 0.4162;
let csca = 8.63e-6 / 2.0 / 1.0 / t.sqrt();
let alf = ALF0 - ALF1 / 1.0;
let bet = BET0 - BET1 / 1.0;
let xi = 1.0;
let xj = 2.0;
let xtt = [1.0, t, t*t, t*t*t];
let cs = compute_collision_rate(2, i, j, t, u0, c1, csca, alf, bet, xi, xj, &xtt);
assert!(cs > 0.0, "Giovanardi hot collision rate should be positive");
assert!(cs.is_finite(), "Giovanardi hot collision rate should be finite");
}
#[test]
fn test_collision_rate_temperature_dependence() {
// 验证碰撞速率随温度的变化
let i = 1;
let j = 2;
let u0_base: f64 = 10.0; // 跃迁能量
let c1: f64 = 0.4162;
let alf = ALF0 - ALF1 / 1.0;
let bet = BET0 - BET1 / 1.0;
let xi = 1.0;
let xj = 2.0;
let temperatures: [f64; 4] = [5000.0, 10000.0, 20000.0, 50000.0];
let mut rates = Vec::new();
for t in &temperatures {
let hkt = HK / t;
let u0 = u0_base * hkt / (HK / 10000.0) * 0.5; // 归一化激发能量
let csca = 8.63e-6 / 2.0 / 1.0 / t.sqrt();
let xtt = [1.0, *t, (*t)*(*t), (*t)*(*t)*(*t)];
let cs = compute_collision_rate(0, i, j, *t, u0, c1, csca, alf, bet, xi, xj, &xtt);
rates.push(cs);
}
// 碰撞速率应该是有限的
for rate in &rates {
assert!(rate.is_finite(), "Collision rate should be finite");
}
}
#[test]
fn test_collision_rate_level_dependence() {
// 验证碰撞速率随能级的变化
let t: f64 = 10000.0;
let c1: f64 = 0.1;
let xtt = [1.0, t, t*t, t*t*t];
let mut rates = Vec::new();
for j in 2..=10 {
let i = j - 1;
let xi = i as f64;
let xj = j as f64;
let u0 = 0.5; // 简化
let csca = 8.63e-6 / 2.0 / ((i*i) as f64) / t.sqrt();
let alf = ALF0 - ALF1 / ((i*i) as f64);
let bet = BET0 - BET1 / (i as f64);
let cs = compute_collision_rate(0, i, j, t, u0, c1, csca, alf, bet, xi, xj, &xtt);
rates.push(cs);
}
// 所有速率应该是正的有限数
for (j, rate) in rates.iter().enumerate() {
assert!(rate.is_finite() && rate > &0.0,
"Rate for transition {}->{} should be positive finite", j+1, j+2);
}
}
#[test]
fn test_exponential_integral_approximation() {
// 验证指数积分近似
// 小 u0 使用级数展开,大 u0 使用连分式
let small_u0 = 0.5_f64;
let large_u0 = 5.0_f64;
// 小 u0 的级数展开
let e1_small = -small_u0.ln() + EXPIA1
+ small_u0 * (EXPIA2 + small_u0 * (EXPIA3
+ small_u0 * (EXPIA4 + small_u0 * (EXPIA5 + small_u0 * EXPIA6))));
assert!(e1_small > 0.0, "E1 for small u0 should be positive");
// 大 u0 的连分式
let u0_inv = 1.0 / large_u0;
let e1_large = (-large_u0).exp()
* ((EXPIB1 + large_u0 * (EXPIB2 + large_u0 * (EXPIB3 + large_u0 * EXPIB4)))
/ (EXPIC1 + large_u0 * (EXPIC2 + large_u0 * (EXPIC3 + large_u0 * EXPIC4))))
* u0_inv;
assert!(e1_large > 0.0 && e1_large < 1.0, "E1 for large u0 should be between 0 and 1");
}
#[test]
fn test_hm_ionization_rate() {
// 验证 H- 碰撞电离速率公式
let t = 10000.0_f64;
let rate = CHMI * t * t.sqrt();
// H- 电离速率应该在合理范围内
assert!(rate > 0.0, "H- ionization rate should be positive");
assert!(rate.is_finite(), "H- ionization rate should be finite");
// CHMI = 5.59e-15, 所以 rate ≈ 5.59e-15 * 10000 * 100 = 5.59e-9
assert!(rate > 1e-12 && rate < 1e-6,
"H- ionization rate {:.2e} seems unreasonable", rate);
}
}
+156
View File
@@ -0,0 +1,156 @@
//! 原子和离子的内能和熵计算。
//!
//! 重构自 TLUSTY `entene.f`。
use crate::math::mpartf::mpartf;
const EV2ERG: f64 = 1.6018e-12;
const ENTCON: f64 = 103.973;
/// ENTENE 参数结构体
#[allow(non_snake_case)]
pub struct EnteneParams<'a> {
pub t: f64,
pub ah: f64,
pub anh: f64,
pub anpr: f64,
pub ane: f64,
pub rr: &'a [[f64; 2]], // 30 x 2 (元素 x 电离态)
pub enev: &'a [[f64; 2]], // 30 x 2 (元素 x 电离态)
pub amas: &'a [f64], // 30
pub natoms: usize,
pub bolk: f64,
pub un: f64, // 配分函数下限
}
/// ENTENE 输出结构体
#[allow(non_snake_case)]
pub struct EnteneOutput {
pub energ: f64,
pub entrop: f64,
}
/// 计算原子和离子的内能和熵。
#[allow(non_snake_case)]
pub fn entene(params: &EnteneParams) -> EnteneOutput {
let t = params.t;
let ah = params.ah;
let anh = params.anh;
let anpr = params.anpr;
let ane = params.ane;
let tk = params.bolk * t;
let tkk = tk * t;
let tkln15 = 1.5 * tk.ln();
let natoms = params.natoms.min(30);
let mut energ = 0.0;
let mut entrop = 0.0;
// 氢
let result = mpartf(1, 1, 0, t);
let mut u = result.u.max(2.0);
let mut dulog = result.dulog.max(0.0);
let alm = 1.5 * params.amas[0].ln();
energ = tk * dulog * anh;
entrop = (tkln15 - anh.ln() + u.ln() + alm + tkk * dulog + ENTCON) * anh;
energ = energ + params.enev[0][0] * EV2ERG * anpr;
entrop = entrop + (tkln15 - anpr.ln() + alm + ENTCON) * anpr;
// 其他物种
let xmax = 2.154e4 * (t / ane).sqrt().sqrt();
for i in 2..=natoms {
let i_idx = i - 1;
let mut chip = 0.0;
for j in 1..=2 {
let j_idx = j - 1;
if params.rr[i_idx][j_idx] > 1e-15 {
let mut aden = params.rr[i_idx][j_idx] * ah;
if aden < 1e-20 {
aden = 1e-20;
}
let result = mpartf(i, j, 0, t);
u = result.u.max(params.un);
dulog = result.dulog.max(0.0);
energ = energ + (chip * EV2ERG + tk * dulog) * aden;
entrop = entrop + (tkln15 - aden.ln() + u.ln()
+ 1.5 * params.amas[i_idx].ln()
+ tkk * dulog + ENTCON) * aden;
}
chip = chip + params.enev[i_idx][j_idx];
}
}
// 电子熵
let entel = tkln15 - ane.ln() - 11.2622 + ENTCON;
entrop = entrop + entel * ane;
entrop = entrop * params.bolk;
EnteneOutput { energ, entrop }
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
#[test]
fn test_entene_simple() {
let rr = [[0.0; 2]; 30];
let enev = [[0.0; 2]; 30];
let amas = [1.0; 30];
let params = EnteneParams {
t: 10000.0,
ah: 1e12,
anh: 1e11,
anpr: 1e11,
ane: 1e10,
rr: &rr,
enev: &enev,
amas: &amas,
natoms: 1,
bolk: 1.380649e-16,
un: 2.0,
};
let result = entene(&params);
// 验证结果为有限值
assert!(result.energ.is_finite());
assert!(result.entrop.is_finite());
}
#[test]
fn test_entene_hydrogen() {
let rr = [[0.0; 2]; 30];
let enev = [[0.0; 2]; 30];
let amas = [1.008; 30]; // H 原子质量
let params = EnteneParams {
t: 8000.0,
ah: 1e14,
anh: 1e13,
anpr: 1e13,
ane: 1e11,
rr: &rr,
enev: &enev,
amas: &amas,
natoms: 1,
bolk: 1.380649e-16,
un: 2.0,
};
let result = entene(&params);
// 结果应该是有限值
assert!(result.energ.is_finite());
assert!(result.entrop.is_finite());
}
}
+558
View File
@@ -0,0 +1,558 @@
//! 耦合系统求解器:流体静力学平衡 + 状态方程 + z-m 关系。
//!
//! 重构自 TLUSTY `hesol6.f`。
//!
//! 使用 Newton-Raphson 方法求解以下耦合系统:
//! 1. 流体静力学平衡方程
//! 2. Ptotal, Pgas 和 rho 的定义
//! 3. 状态方程
//! 4. z-m 关系
use crate::math::erfcx::erfcx;
use crate::math::matinv::matinv;
const UN: f64 = 1.0;
const HALF: f64 = 0.5;
const TWO: f64 = 2.0;
const ERROR: f64 = 1e-4;
const NITERH: usize = 15;
const MP: usize = 6;
const NP: usize = 6;
const IP: usize = 1; // PTOTAL index
const IG: usize = 2; // PGS index
const IR: usize = 3; // DENS index
const IN: usize = 4; // ANTT index
const IE: usize = 5; // ELEC index
const IZ: usize = 6; // ZD index
/// HESOL6 参数结构体
#[allow(non_snake_case)]
pub struct Hesol6Params<'a> {
// 标量
pub nd: usize,
pub iter: usize,
pub nfreqe: usize,
pub teff: f64,
pub sig4p: f64,
pub pck: f64,
pub qgrav: f64,
pub bolk: f64,
pub abrosd: &'a [f64],
pub fprd: &'a [f64],
pub grd: f64,
// 数组
pub temp: &'a mut [f64],
pub dm: &'a mut [f64],
pub zd: &'a mut [f64],
pub dens: &'a mut [f64],
pub elec: &'a mut [f64],
pub wmm: &'a [f64],
pub pradt: &'a [f64],
pub pgs: &'a mut [f64],
pub ptotal: &'a mut [f64],
// 辐射场相关 (可选)
pub ijfr: &'a [usize],
pub lskip: &'a [bool], // nd x nfreq
pub w: &'a [f64],
pub fh: &'a [f64],
pub radex: &'a [f64],
pub hextrd: &'a [f64],
pub absoe1: &'a [f64],
}
/// HESOL6 输出结构体
#[allow(non_snake_case)]
pub struct Hesol6Output {
pub iterh: usize,
pub chanm: f64,
pub err: [[f64; NP]; 100], // 简化:固定大小
pub chng: Vec<f64>,
}
/// HESOL6 辅助 COMMON 块
#[allow(non_snake_case)]
pub struct Hesol6Aux {
pub vsnd2: Vec<f64>,
pub hg1: f64,
pub hr1: f64,
pub rr1: f64,
}
/// 求解耦合系统。
#[allow(non_snake_case)]
pub fn hesol6(params: &mut Hesol6Params) -> Hesol6Output {
let nd = params.nd;
// 辐射压尺度高度
let mut hr1 = if params.iter == 0 {
params.sig4p * params.teff.powi(4) * params.pck * params.abrosd[0] / params.qgrav
} else {
let mut grd = params.grd;
if params.nfreqe > 0 {
for ij in 1..=params.nfreqe {
let ijt = params.ijfr[ij - 1] - 1;
let lskip_val = params.lskip[0 * params.nd + ijt]; // ID=1
if !lskip_val {
let fluxw = params.w[ijt]
* (params.fh[ijt] * params.radex[(ij - 1) * nd] - params.hextrd[ijt]);
grd += fluxw * params.absoe1[ij - 1];
}
}
}
params.pck / params.qgrav * (grd + params.fprd[0]) / params.dens[0]
};
// 初始化
let mut antt = vec![0.0; nd];
for id in 0..nd {
antt[id] = params.dens[id] / params.wmm[id] + params.elec[id];
params.pgs[id] = antt[id] * params.bolk * params.temp[id];
params.ptotal[id] = params.pradt[id] + params.pgs[id];
}
let mut p = params.ptotal.to_vec();
// 工作数组
let mut vec = [[0.0; 100]; NP]; // NP x ND
let mut vec1 = [[0.0; 100]; NP];
let mut vec2 = [[0.0; 100]; NP];
let mut vec3 = [[0.0; 100]; NP];
let mut anu = [[0.0; 100]; NP];
let mut d_arr = [[[0.0; 100]; NP]; NP]; // NP x NP x ND
let mut anerl = vec![0.0; nd];
let mut chng = vec![0.0; nd];
let mut err = [[0.0; NP]; 100];
// Newton-Raphson 迭代
let mut iterh = 0;
let mut lac2h = false;
let iach = 6;
let iacdh = 4;
let mut iach0 = iach - 3;
loop {
iterh += 1;
// ================== 前向消元 ==================
for id in 0..nd {
let id_1 = id + 1; // 1-indexed
antt[id] = params.dens[id] / params.wmm[id] + params.elec[id];
params.pgs[id] = antt[id] * params.bolk * params.temp[id];
params.ptotal[id] = params.pradt[id] + params.pgs[id];
p[id] = params.ptotal[id];
vec[IP - 1][id] = params.ptotal[id];
vec[IG - 1][id] = params.pgs[id];
vec[IR - 1][id] = params.dens[id];
vec[IN - 1][id] = antt[id];
vec[IE - 1][id] = params.elec[id];
vec[IZ - 1][id] = params.zd[id];
let mut a = [[0.0; NP]; NP];
let mut b = [[0.0; NP]; NP];
let mut c = [[0.0; NP]; NP];
let mut vl = [0.0; NP];
if id == 0 {
// 上边界条件
let hg1_val = (TWO * params.pgs[id] / params.dens[id] / params.qgrav).sqrt();
let mut x = (params.zd[id] - hr1) / hg1_val;
let f1 = if x < 3.0 {
if x < 0.0 { x = 0.0; }
8.86226925e-1 * (x * x).exp() * erfcx(x)
} else {
HALF * (UN - HALF / x / x) / x
};
let x1 = x * 1.01;
let f1d = if x1 < 3.0 {
8.86226925e-1 * (x1 * x1).exp() * erfcx(x1)
} else {
HALF * (UN - HALF / x1 / x1) / x1
};
let f1d_deriv = if x > 0.0 { (f1d - f1) * 100.0 / x } else { 0.0 };
b[IG - 1][IG - 1] = params.dens[id] * hg1_val * HALF / params.pgs[id] * (f1 - f1d_deriv * x);
b[IG - 1][IR - 1] = hg1_val * HALF * (f1 + f1d_deriv * x);
b[IG - 1][IZ - 1] = params.dens[id] * f1d_deriv;
vl[IG - 1] = params.dm[id] - params.dens[id] * hg1_val * f1;
b[IP - 1][IP - 1] = UN;
b[IP - 1][IG - 1] = -UN;
vl[IP - 1] = params.pradt[id] + params.pgs[id] - p[id];
} else {
// 内部点
let dmm = UN / (params.dm[id] - params.dm[id - 1]);
let qg = HALF * params.qgrav;
b[IP - 1][IP - 1] = dmm;
b[IP - 1][IZ - 1] = -qg;
a[IP - 1][IP - 1] = dmm;
a[IP - 1][IZ - 1] = qg;
vl[IP - 1] = qg * (params.zd[id] + params.zd[id - 1]) - (p[id] - p[id - 1]) * dmm;
b[IG - 1][IP - 1] = UN;
b[IG - 1][IG - 1] = -UN;
vl[IG - 1] = params.pradt[id] + params.pgs[id] - p[id];
}
// 密度方程
b[IR - 1][IR - 1] = UN;
b[IR - 1][IN - 1] = -params.wmm[id];
b[IR - 1][IE - 1] = params.wmm[id];
vl[IR - 1] = params.wmm[id] * (antt[id] - params.elec[id]) - params.dens[id];
// 状态方程
b[IN - 1][IG - 1] = UN;
b[IN - 1][IN - 1] = -params.bolk * params.temp[id];
vl[IN - 1] = antt[id] * params.bolk * params.temp[id] - params.pgs[id];
// 电子密度
anerl[id] = params.elec[id] / antt[id];
b[IE - 1][IE - 1] = UN;
b[IE - 1][IN - 1] = -anerl[id];
vl[IE - 1] = anerl[id] * antt[id] - params.elec[id];
// z-m 关系
if id < nd - 1 {
let dmp = (params.dm[id + 1] - params.dm[id]) * HALF;
b[IZ - 1][IR - 1] = dmp / params.dens[id].powi(2);
b[IZ - 1][IZ - 1] = UN;
c[IZ - 1][IR - 1] = -dmp / params.dens[id + 1].powi(2);
c[IZ - 1][IZ - 1] = UN;
vl[IZ - 1] = params.zd[id + 1] - params.zd[id]
+ dmp / params.dens[id] + dmp / params.dens[id + 1];
} else {
b[IZ - 1][IZ - 1] = UN;
vl[IZ - 1] = 0.0;
}
// 前向消元更新
if id > 0 {
b[IP - 1][IR - 1] -= a[IP - 1][IP - 1] * d_arr[IP - 1][IR - 1][id - 1]
+ a[IP - 1][IZ - 1] * d_arr[IZ - 1][IR - 1][id - 1];
b[IP - 1][IZ - 1] -= a[IP - 1][IP - 1] * d_arr[IP - 1][IZ - 1][id - 1]
+ a[IP - 1][IZ - 1] * d_arr[IZ - 1][IZ - 1][id - 1];
vl[IP - 1] += a[IP - 1][IP - 1] * anu[IP - 1][id - 1]
+ a[IP - 1][IZ - 1] * anu[IZ - 1][id - 1];
}
// 矩阵求逆
let mut b_vec: Vec<f64> = b.iter().flat_map(|row| row.iter().copied()).collect();
matinv(&mut b_vec, NP);
for i in 0..NP {
for j in 0..NP {
b[i][j] = b_vec[i * MP + j];
}
}
// 计算 ANU
for i in 0..NP {
let mut sum = 0.0;
for j in 0..NP {
sum += b[i][j] * vl[j];
}
anu[i][id] = sum;
}
// 存储 D 数组
if id < nd - 1 {
for i in 0..NP {
d_arr[i][IR - 1][id] = b[i][IZ - 1] * c[IZ - 1][IR - 1];
d_arr[i][IZ - 1][id] = b[i][IZ - 1] * c[IZ - 1][IZ - 1];
}
}
}
// ================== 回代 ==================
let mut chanm = 0.0;
for id in (0..nd).rev() {
chng[id] = 0.0;
let sol = if id == nd - 1 {
let mut s = [0.0; NP];
for i in 0..NP {
s[i] = anu[i][id];
}
s
} else {
for i in 0..NP {
anu[i][id] = anu[i][id] + d_arr[i][IR - 1][id] * anu[IR - 1][id + 1]
+ d_arr[i][IZ - 1][id] * anu[IZ - 1][id + 1];
}
let mut s = [0.0; NP];
for i in 0..NP {
s[i] = anu[i][id];
}
s
};
for i in 0..NP {
let chan = if vec[i][id] != 0.0 { sol[i] / vec[i][id] } else { 0.0 };
let chan_clamped = chan.max(-0.99).min(99.00);
if chan.abs() > chanm {
chanm = chan.abs();
}
if chan.abs() > chng[id] {
chng[id] = chan.abs();
}
vec[i][id] = vec[i][id] * (UN + chan_clamped);
}
params.ptotal[id] = vec[IP - 1][id];
p[id] = params.ptotal[id];
params.pgs[id] = vec[IG - 1][id];
params.dens[id] = vec[IR - 1][id];
antt[id] = vec[IN - 1][id];
params.elec[id] = vec[IE - 1][id];
params.zd[id] = vec[IZ - 1][id];
// 计算误差
if id == 0 {
let hg1_val = (TWO * params.pgs[id] / params.dens[id] / params.qgrav).sqrt();
let mut x = (params.zd[id] - hr1) / hg1_val;
let f1 = if x < 3.0 {
if x < 0.0 { x = 0.0; }
8.86226925e-1 * (x * x).exp() * erfcx(x)
} else {
HALF * (UN - HALF / x / x) / x
};
err[IP - 1][id] = (params.dens[id] * hg1_val * f1 - params.dm[id]) / params.dm[id];
} else if id < nd - 1 {
err[IP - 1][id] = (p[id + 1] - p[id]) / (params.dm[id + 1] - params.dm[id])
* TWO / params.qgrav / (params.zd[id + 1] + params.zd[id]) - UN;
} else {
err[IP - 1][id] = 0.0;
}
if p[id] != 0.0 {
err[IG - 1][id] = (params.pradt[id] + params.pgs[id] - p[id]) / p[id];
}
if params.dens[id] != 0.0 {
err[IR - 1][id] = (params.dens[id] - (antt[id] - params.elec[id]) * params.wmm[id])
/ params.dens[id];
}
if params.pgs[id] != 0.0 {
err[IN - 1][id] = (params.pgs[id] - antt[id] * params.bolk * params.temp[id])
/ params.pgs[id];
}
if params.elec[id] != 0.0 {
err[IE - 1][id] = (params.elec[id] - anerl[id] * antt[id]) / params.elec[id];
}
if id < nd - 1 {
err[IZ - 1][id] = (params.zd[id] - params.zd[id + 1]) * TWO
/ (params.dm[id + 1] - params.dm[id])
/ (UN / params.dens[id] + UN / params.dens[id + 1]) - UN;
} else {
err[IZ - 1][id] = params.zd[id];
}
}
// ================== 加速 ==================
if iterh >= iach0 && iterh < iach {
// 存储历史
let ipt = iterh % 3;
let ipt1 = (iach + 1) % 3;
let ipt2 = (iach + 2) % 3;
if iterh == iach0 {
for id in 0..nd {
for ix in 0..NP {
vec3[ix][id] = vec[ix][id];
}
}
} else if ipt == ipt1 {
for id in 0..nd {
for ix in 0..NP {
vec2[ix][id] = vec[ix][id];
}
}
} else if ipt == ipt2 {
for id in 0..nd {
for ix in 0..NP {
vec1[ix][id] = vec[ix][id];
}
}
}
} else if iterh >= iach && !lac2h {
// Ng 加速
let mut a1 = 0.0;
let mut b1 = 0.0;
let mut b2 = 0.0;
let mut c1 = 0.0;
let mut c2 = 0.0;
for ix in 0..NP {
for id in 0..nd {
let wt = if vec[ix][id] != 0.0 { 1.0 / vec[ix][id].abs() } else { 0.0 };
let d0 = vec[ix][id] - vec1[ix][id];
let d1 = d0 - vec1[ix][id] + vec2[ix][id];
let d2 = d0 - vec2[ix][id] + vec3[ix][id];
a1 += wt * d1 * d1;
b1 += wt * d1 * d2;
b2 += wt * d2 * d2;
c1 += wt * d0 * d1;
c2 += wt * d0 * d2;
}
}
let ab = b2 * a1 - b1 * b1;
if ab != 0.0 {
let aa = (b2 * c1 - b1 * c2) / ab;
let bb = (a1 * c2 - b1 * c1) / ab;
for id in 0..nd {
for ix in 0..NP {
vec[ix][id] = (UN - aa - bb) * vec[ix][id]
+ aa * vec1[ix][id]
+ bb * vec2[ix][id];
}
}
lac2h = true;
}
}
// 更新变量
for id in 0..nd {
params.ptotal[id] = vec[IP - 1][id];
p[id] = params.ptotal[id];
params.pgs[id] = vec[IG - 1][id];
params.dens[id] = vec[IR - 1][id];
antt[id] = vec[IN - 1][id];
params.elec[id] = vec[IE - 1][id];
params.zd[id] = vec[IZ - 1][id];
// 更新误差
if id == 0 {
let hg1_val = (TWO * params.pgs[id] / params.dens[id] / params.qgrav).sqrt();
let mut x = (params.zd[id] - hr1) / hg1_val;
let f1 = if x < 3.0 {
if x < 0.0 { x = 0.0; }
8.86226925e-1 * (x * x).exp() * erfcx(x)
} else {
HALF * (UN - HALF / x / x) / x
};
err[IP - 1][id] = (params.dens[id] * hg1_val * f1 - params.dm[id]) / params.dm[id];
} else if id < nd - 1 {
err[IP - 1][id] = (p[id] - p[id - 1]) / (params.dm[id] - params.dm[id - 1])
* TWO / params.qgrav / (params.zd[id] + params.zd[id - 1]) - UN;
}
if p[id] != 0.0 {
err[IG - 1][id] = (params.pradt[id] + params.pgs[id] - p[id]) / p[id];
}
if params.dens[id] != 0.0 {
err[IR - 1][id] = (params.dens[id] - (antt[id] - params.elec[id]) * params.wmm[id])
/ params.dens[id];
}
if params.pgs[id] != 0.0 {
err[IN - 1][id] = (params.pgs[id] - antt[id] * params.bolk * params.temp[id])
/ params.pgs[id];
}
if params.elec[id] != 0.0 {
err[IE - 1][id] = (params.elec[id] - anerl[id] * antt[id]) / params.elec[id];
}
if id > 0 {
err[IZ - 1][id] = (params.zd[id - 1] - params.zd[id]) * TWO
/ (params.dm[id] - params.dm[id - 1])
/ (UN / params.dens[id - 1] + UN / params.dens[id]) - UN;
} else {
err[IZ - 1][id] = 0.0;
}
}
// 收敛检查
if chanm <= ERROR || iterh >= NITERH {
break;
}
}
Hesol6Output {
iterh,
chanm: 0.0, // 最后的 chanm
err,
chng,
}
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
#[test]
fn test_hesol6_constants() {
assert_relative_eq!(UN, 1.0, epsilon = 1e-15);
assert_relative_eq!(HALF, 0.5, epsilon = 1e-15);
assert_relative_eq!(TWO, 2.0, epsilon = 1e-15);
assert_relative_eq!(ERROR, 1e-4, epsilon = 1e-10);
assert_eq!(NP, 6);
assert_eq!(MP, 6);
}
#[test]
fn test_hesol6_simple() {
// 简单测试:验证函数可以运行
let nd = 3;
let mut temp = vec![10000.0; nd];
let mut dm = vec![1e-4, 1e-3, 1e-2];
let mut zd = vec![0.0, 1e8, 2e8];
let mut dens = vec![1e-7; nd];
let mut elec = vec![1e-8; nd];
let mut pgs = vec![0.0; nd];
let mut ptotal = vec![0.0; nd];
let wmm = vec![1.0; nd];
let pradt = vec![0.0; nd];
let abrosd = vec![1.0; nd];
let fprd = vec![0.0; nd];
let mut params = Hesol6Params {
nd,
iter: 0,
nfreqe: 0,
teff: 10000.0,
sig4p: 5.670373e-5 * 4.0, // sigma * 4
pck: 2.99792458e10,
qgrav: 2.739e4,
bolk: 1.380649e-16,
abrosd: &abrosd,
fprd: &fprd,
grd: 0.0,
temp: &mut temp,
dm: &mut dm,
zd: &mut zd,
dens: &mut dens,
elec: &mut elec,
wmm: &wmm,
pradt: &pradt,
pgs: &mut pgs,
ptotal: &mut ptotal,
ijfr: &[],
lskip: &[],
w: &[],
fh: &[],
radex: &[],
hextrd: &[],
absoe1: &[],
};
let output = hesol6(&mut params);
// 验证迭代次数在合理范围内
assert!(output.iterh <= NITERH);
}
}
+338
View File
@@ -0,0 +1,338 @@
//! Rosseland 和 Planck 平均不透明度计算(使用 OPCTAB)。
//!
//! 重构自 TLUSTY `meanopt.f`
//!
//! 通过调用 OPCTAB 动态计算每个频率点的吸收和散射系数,
//! 然后计算 Rosseland 和 Planck 平均不透明度。
use super::opctab::{opctab, OpctabOutput};
use crate::state::constants::{HK, UN};
/// MEANOPT 输入参数
pub struct MeanoptParams<'a> {
/// 温度 (K)
pub t: f64,
/// 深度索引 (1-indexed)
pub id: usize,
/// 密度 (g cm^-3)
pub rho: f64,
/// 迭代次数
pub iter: i32,
/// Rayleigh 散射标志 (<0: 简单公式, >0: 调用 rayleigh, =0: 关闭)
pub ifrayl: i32,
/// 选项表标志 (<0: 添加电子散射)
pub ioptab: i32,
/// Rayleigh 参数 (当 ifrayl > 0 时需要)
pub rayleigh_params: Option<&'a super::rayleigh::RayleighParams<'a>>,
}
/// MEANOPT 模型状态
pub struct MeanoptModelState<'a> {
/// 频率数组 (nfreq)
pub freq: &'a [f64],
/// Planck 函数 (nfreq)
pub bnue: &'a [f64],
/// 权重 (nfreq)
pub w: &'a [f64],
}
/// MEANOPT 输出
pub struct MeanoptOutput {
/// Rosseland 平均不透明度 (per gram)
pub opros: f64,
/// Planck 平均不透明度 (per gram)
pub oppla: f64,
}
/// 计算 Rosseland 和 Planck 平均不透明度。
///
/// 通过调用 OPCTAB 为每个频率点动态计算吸收和散射系数。
///
/// # 参数
///
/// * `params` - 输入参数
/// * `model` - 模型状态(频率、Planck 函数、权重)
/// * `opctab_params_base` - OPCTAB 基础参数(不包括 fr, ij
/// * `opctab_table` - 不透明度表数据
/// * `opctab_model` - OPCTAB 模型状态
///
/// # 返回值
///
/// 包含 Rosseland 和 Planck 平均不透明度的结构体
pub fn meanopt(
params: &MeanoptParams,
model: &MeanoptModelState,
opctab_table: &super::opctab::OpctabTableData,
opctab_model: &mut super::opctab::OpctabModelState,
) -> MeanoptOutput {
let t = params.t;
let hkt = HK / t;
let mut abr = 0.0;
let mut sumdb = 0.0;
let mut abp = 0.0;
let mut sumb = 0.0;
let nfreq = model.freq.len();
for ij in 0..nfreq {
let fr = model.freq[ij];
let ex = (hkt * fr).exp();
let e1 = UN / (ex - UN);
let plan = model.bnue[ij] * e1 * model.w[ij];
let dplan = plan * hkt * fr * ex * e1;
// 调用 OPCTAB 获取吸收和散射系数
let opctab_params = super::opctab::OpctabParams {
fr,
ij: ij + 1, // 1-indexed
id: params.id,
t,
rho: params.rho,
igram: 1, // 返回每克
iter: params.iter,
ifrayl: params.ifrayl,
ioptab: params.ioptab,
rayleigh_params: params.rayleigh_params,
};
let OpctabOutput { ab, sct, .. } = opctab(&opctab_params, opctab_table, opctab_model);
abr += dplan / (ab + sct);
abp += plan * ab;
sumdb += dplan;
sumb += plan;
}
let opros = sumdb / abr;
let oppla = abp / sumb;
MeanoptOutput { opros, oppla }
}
#[cfg(test)]
mod tests {
use super::*;
use super::super::opctab::{OpctabTableData, OpctabModelState};
use approx::assert_relative_eq;
#[test]
fn test_meanopt_function_basic() {
let numtemp = 1;
let nd = 1;
let nfreq = 3;
let max_numrh = 1;
let tempvec = vec![9.2103];
let numrh = vec![1];
let rhomat = vec![-16.1181];
let absopac = vec![0.0, 0.0, 0.0];
let raysc = vec![0.0];
let freq = vec![1e14, 3e14, 1e15];
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: 0.0,
};
let t = 10000.0;
let rho = 1e-7;
// 简化的 Planck 函数
let bnue: Vec<f64> = freq.iter().map(|f| 1.4743e-2 * f.powi(3)).collect();
let w: Vec<f64> = vec![0.3, 0.4, 0.3]; // 权重
let params = MeanoptParams {
t,
id: 1,
rho,
iter: 1,
ifrayl: 0,
ioptab: 0,
rayleigh_params: None,
};
let model = MeanoptModelState {
freq: &freq,
bnue: &bnue,
w: &w,
};
let elec = vec![1e-10];
let dens = vec![rho];
let mut opctab_model = OpctabModelState {
elec: &elec,
dens: &dens,
raysct: None,
eospar: None,
};
// 调用函数
let result = meanopt(&params, &model, &table, &mut opctab_model);
// 验证输出
assert!(result.opros > 0.0, "Rosseland mean should be positive");
assert!(result.oppla > 0.0, "Planck mean should be positive");
assert!(result.opros.is_finite(), "Rosseland mean should be finite");
assert!(result.oppla.is_finite(), "Planck mean should be finite");
}
#[test]
fn test_meanopt_function_uniform_opacity() {
let numtemp = 1;
let nd = 1;
let nfreq = 3;
let max_numrh = 1;
let tempvec = vec![9.2103];
let numrh = vec![1];
let rhomat = vec![-16.1181];
let absopac = vec![0.0, 0.0, 0.0];
let raysc = vec![0.0];
let freq = vec![1e14, 3e14, 1e15];
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: 0.0,
};
let t = 10000.0;
let rho = 1e-7;
// 使用相同的权重
let bnue: Vec<f64> = vec![1.0, 1.0, 1.0];
let w: Vec<f64> = vec![1.0/3.0, 1.0/3.0, 1.0/3.0];
let params = MeanoptParams {
t,
id: 1,
rho,
iter: 1,
ifrayl: 0,
ioptab: 0,
rayleigh_params: None,
};
let model = MeanoptModelState {
freq: &freq,
bnue: &bnue,
w: &w,
};
let elec = vec![1e-10];
let dens = vec![rho];
let mut opctab_model = OpctabModelState {
elec: &elec,
dens: &dens,
raysct: None,
eospar: None,
};
let result = meanopt(&params, &model, &table, &mut opctab_model);
// 验证结果合理
assert!(result.opros > 0.0 && result.opros < 10.0,
"Rosseland mean should be reasonable: got {}", result.opros);
assert!(result.oppla > 0.0 && result.oppla < 10.0,
"Planck mean should be reasonable: got {}", result.oppla);
}
#[test]
fn test_meanopt_function_temperature_dependence() {
let numtemp = 1;
let nd = 1;
let nfreq = 3;
let max_numrh = 1;
let tempvec = vec![9.2103];
let numrh = vec![1];
let rhomat = vec![-16.1181];
let absopac = vec![0.0, 0.0, 0.0];
let raysc = vec![0.0];
let freq = vec![1e14, 3e14, 1e15];
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: 0.0,
};
let bnue: Vec<f64> = freq.iter().map(|f| 1.4743e-2 * f.powi(3)).collect();
let w: Vec<f64> = vec![0.3, 0.4, 0.3];
let rho = 1e-7;
let elec = vec![1e-10];
let dens = vec![rho];
// 测试不同温度
let temperatures = [5000.0, 10000.0, 20000.0];
let mut results = Vec::new();
for t in &temperatures {
let params = MeanoptParams {
t: *t,
id: 1,
rho,
iter: 1,
ifrayl: 0,
ioptab: 0,
rayleigh_params: None,
};
let model = MeanoptModelState {
freq: &freq,
bnue: &bnue,
w: &w,
};
let mut opctab_model = OpctabModelState {
elec: &elec,
dens: &dens,
raysct: None,
eospar: None,
};
results.push(meanopt(&params, &model, &table, &mut opctab_model));
}
// 所有结果应该是有效的
for (i, r) in results.iter().enumerate() {
assert!(r.opros > 0.0 && r.opros.is_finite(),
"Rosseland mean at T={} should be valid", temperatures[i]);
assert!(r.oppla > 0.0 && r.oppla.is_finite(),
"Planck mean at T={} should be valid", temperatures[i]);
}
}
}
+47
View File
@@ -8,13 +8,19 @@ mod allardt;
mod angset;
mod betah;
mod bkhsgo;
mod bre;
mod brez;
mod brte;
mod brtez;
mod bhe;
mod bpopf;
mod bpope;
mod butler;
mod carbon;
mod ceh12;
mod cion;
mod ckoest;
mod colh;
mod collhe;
mod compt0;
mod comset;
@@ -30,6 +36,7 @@ mod dwnfr;
mod dwnfr0;
mod dwnfr1;
mod emat;
mod entene;
mod erfcx;
mod expo;
mod expint;
@@ -45,6 +52,7 @@ mod gntk;
mod gridp;
mod grcor;
mod hephot;
mod hesol6;
mod hidalg;
mod indexx;
mod inicom;
@@ -63,13 +71,20 @@ mod linspl;
mod locate;
mod matinv;
mod meanop;
mod meanopt;
mod minv3;
mod mpartf;
mod odfhst;
mod odfhyd;
mod odfmer;
mod odffr;
mod opadd0;
mod opact1;
mod opactd;
mod opctab;
mod pfcno;
mod pffe;
mod prd;
mod prdini;
mod profil;
mod quartc;
@@ -87,6 +102,7 @@ mod rayleigh;
mod rybmat;
mod rayset;
mod reiman;
mod rteang;
mod rte_sc;
mod rtefe2;
mod rtedf1;
@@ -101,6 +117,7 @@ mod sbfoh;
mod setdrt;
mod sghe12;
mod sgmer;
mod sigmar;
mod sffhmi;
mod tabint;
mod taufr1;
@@ -110,6 +127,7 @@ mod stark0;
mod starka;
mod szirc;
mod tiopf;
mod tlocal;
mod tdpini;
mod traini;
mod tridag;
@@ -135,13 +153,22 @@ pub use allardt::{allardt, AllardData};
pub use angset::angset;
pub use betah::betah;
pub use bkhsgo::bkhsgo;
pub use bre::{bre, BreParams, BreState};
pub use brez::{brez, BrezParams, BrezState};
pub use brte::{brte, BrteParams, BrteState};
pub use brtez::{brtez, BrtezParams, BrtezState};
pub use bhe::{bhe, bhed, bhez, BheParams, BheState, MatKey};
pub use bpopf::{bpopf, BpopfParams};
pub use bpope::{
bpope, BpopeAtomicData, BpopeConfig, BpopeFreqData, BpopeMatrixData, BpopeModelState,
BpopeOutput, BpopeParams,
};
pub use butler::butler;
pub use carbon::carbon;
pub use ceh12::ceh12;
pub use cion::cion;
pub use ckoest::ckoest;
pub use colh::{colh, ColhAtomicData, ColhOutput, ColhParams};
pub use collhe::collhe;
pub use comset::{comset, ComsetParams, ComsetResult};
pub use cross::{cross, crossd};
@@ -156,6 +183,7 @@ pub use dwnfr::dwnfr;
pub use dwnfr0::dwnfr0;
pub use dwnfr1::dwnfr1;
pub use emat::emat;
pub use entene::{entene, EnteneOutput, EnteneParams};
pub use erfcx::{erfcin, erfcx};
pub use expo::expo;
pub use expint::{eint, expinx};
@@ -171,6 +199,7 @@ pub use gntk::gntk;
pub use gridp::gridp;
pub use grcor::grcor;
pub use hephot::hephot;
pub use hesol6::{hesol6, Hesol6Aux, Hesol6Output, Hesol6Params};
pub use hidalg::hidalg;
pub use indexx::indexx;
pub use inicom::inicom;
@@ -189,13 +218,26 @@ pub use linspl::{linspl, LinsplParams};
pub use locate::locate;
pub use matinv::matinv;
pub use meanop::meanop;
pub use meanopt::{meanopt, MeanoptModelState, MeanoptOutput, MeanoptParams};
pub use minv3::minv3;
pub use mpartf::{mpartf, MpartfResult};
pub use opadd0::{opadd0, Opadd0Params, Opadd0FreqData, Opadd0OutputState};
pub use opact1::{
opact1, Opact1ModelState, Opact1OutputState, Opact1Params,
};
pub use opactd::{
opactd, OpactdExpData, OpactdModelState, OpactdOutputState, OpactdParams,
};
pub use opctab::{opctab, OpctabParams, OpctabTableData, OpctabModelState, OpctabOutput};
pub use odfhst::odfhst;
pub use odfhyd::{
odfhyd, OdfhydAtomicData, OdfhydConfig, OdfhydModelState, OdfhydOdfData, OdfhydParams,
};
pub use odfmer::{odfmer, OdfmerAtomicData, OdfmerModelState, OdfmerParams};
pub use odffr::{odffr, OdffrParams, OdffrAtomicData, OdffrModelData, OdffrOutputState};
pub use pfcno::pfcno;
pub use pffe::pffe;
pub use prd::prd;
pub use prdini::prdini;
pub use profil::{profil, ProfilParams};
pub use pfni::pfni;
@@ -215,6 +257,7 @@ pub use rayleigh::{
};
pub use rayset::rayset;
pub use reiman::reiman;
pub use rteang::{rteang, RteangOutput, RteangParams};
pub use rte_sc::rte_sc;
pub use rtefe2::rtefe2;
pub use rtedf1::{rtedf1, Rtedf1AliState, Rtedf1ModelState, Rtedf1Params};
@@ -230,6 +273,7 @@ pub use sbfoh::sbfoh;
pub use setdrt::setdrt;
pub use sghe12::sghe12;
pub use sgmer::{sgmer0, sgmer1, sgmerd};
pub use sigmar::sigmar;
pub use sffhmi::sffhmi;
pub use sffhmi_add::sffhmi_add;
pub use spsigk::spsigk;
@@ -239,6 +283,9 @@ pub use stark0::stark0;
pub use starka::starka;
pub use szirc::szirc;
pub use tiopf::tiopf;
pub use tlocal::{
tlocal, TlocalConfig, TlocalFactrs, TlocalFlxaux, TlocalModelState, TlocalParams,
};
pub use tdpini::tdpini;
pub use traini::traini;
pub use tridag::tridag;
+191
View File
@@ -0,0 +1,191 @@
//! 配分函数计算器(Irwin 多项式数据)。
//!
//! 重构自 TLUSTY `mpartf.f`。
//!
//! 使用 Irwin (1981, ApJ Suppl. 45, 621) 的多项式数据计算配分函数。
//!
//! 公式: ln u(T) = Σ a(i) * (ln T)^(i-1), i=1..6
/// 配分函数结果
#[derive(Debug, Clone, Copy)]
pub struct MpartfResult {
/// 配分函数 u(线性标度)
pub u: f64,
/// d ln(u) / d ln(T)
pub dulog: f64,
}
/// 原子配分函数数据(简化版,仅包含常用元素)
///
/// 数据格式: a(1)..a(6) 对于 ion=1,2,3
/// 来自 Irwin 1981
static ATOMIC_DATA: &[(usize, [[f64; 6]; 3])] = &[
// (jatom, [a1..a6 for ion=1, ion=2, ion=3])
// H (1)
(1, [
[-0.440174, 0.704848, -0.132879, 0.0, 0.0, 0.0], // H I
[-1.07213, 0.0, 0.0, 0.0, 0.0, 0.0], // H II
[0.0; 6], // 无 H III
]),
// He (2)
(2, [
[-0.406163, 0.0, 0.0, 0.0, 0.0, 0.0], // He I
[-1.07213, 0.0, 0.0, 0.0, 0.0, 0.0], // He II
[0.0; 6], // 无 He III (实际存在但简化)
]),
// C (6)
(6, [
[1.52679, 0.561367, -0.0393169, 0.0, 0.0, 0.0], // C I
[0.719738, 0.269122, -0.0277346, 0.0, 0.0, 0.0], // C II
[0.0278285, 0.341071, -0.0448487, 0.0, 0.0, 0.0], // C III
]),
// N (7)
(7, [
[0.889762, 0.366812, -0.0223886, 0.0, 0.0, 0.0], // N I
[0.552811, 0.186518, -0.0119729, 0.0, 0.0, 0.0], // N II
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], // N III (简化)
]),
// O (8)
(8, [
[1.17507, 0.242649, -0.0109539, 0.0, 0.0, 0.0], // O I
[0.675772, 0.123212, -0.00544146, 0.0, 0.0, 0.0], // O II
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], // O III (简化)
]),
// Si (14)
(14, [
[1.39128, 0.628384, -0.0572266, 0.0, 0.0, 0.0], // Si I
[0.674652, 0.262046, -0.0224006, 0.0, 0.0, 0.0], // Si II
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], // Si III (简化)
]),
// Fe (26)
(26, [
[1.73902, 1.44676, -0.219931, 0.0, 0.0, 0.0], // Fe I
[0.845661, 0.640025, -0.0900389, 0.0, 0.0, 0.0], // Fe II
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], // Fe III (简化)
]),
];
/// 统计权重数组(用于 T > 16000K 时的外推)
/// IGLE 数组来自 Fortran
static IGLE: &[usize] = &[2, 1, 2, 1, 6, 9, 4, 9, 6, 1, 2, 1, 6, 9, 4, 9, 6, 1, 10, 21, 28, 25, 6, 25, 28, 21, 10, 21];
/// 计算配分函数。
///
/// # 参数
/// * `jatom` - 元素周期表中的原子序数(1-92)
/// * `ion` - 电离态(1=中性,2=一次电离,3=二次电离)
/// * `indmol` - 分子物种索引(Tsuji 索引,0 表示原子)
/// * `t` - 温度(K
///
/// # 返回值
/// (u, dulog) - 配分函数及其对温度的对数导数
#[allow(non_snake_case)]
pub fn mpartf(jatom: usize, ion: usize, indmol: usize, t: f64) -> MpartfResult {
let mut u = 1.0;
let mut dulog = 0.0;
// 温度范围检查
if t < 1000.0 {
panic!("mpartf: temperature {} < 1000 K", t);
}
if t > 16000.0 {
// 高温外推:使用统计权重
if indmol == 0 && jatom > 0 && jatom <= 28 && ion > 0 {
let igle_idx = jatom - ion + 1;
if igle_idx > 0 && igle_idx <= IGLE.len() {
u = IGLE[igle_idx - 1] as f64;
}
}
return MpartfResult { u, dulog };
}
let tl = t.ln();
// 原子物种
if jatom > 0 && ion > 0 && ion <= 3 && indmol == 0 {
// 查找原子数据
if let Some((_, data)) = ATOMIC_DATA.iter().find(|(z, _)| *z == jatom) {
let a = &data[ion - 1];
// 检查是否有效数据
if a[0] != 0.0 || a[1] != 0.0 {
let ulog = a[0] + tl * (a[1] + tl * (a[2] + tl * (a[3] + tl * (a[4] + tl * a[5]))));
// B III 特殊处理
let ulog = if jatom == 5 && ion == 3 { 1.0 } else { ulog };
u = ulog.exp();
dulog = a[1] + tl * (a[2] * 2.0 + tl * (a[3] * 3.0 + tl * (a[4] * 4.0 + tl * a[5] * 5.0)));
}
} else {
// 未找到数据,使用简单估计
u = 1.0;
dulog = 0.0;
}
}
// 分子物种(简化:不实现分子数据)
if indmol > 0 {
// 分子配分函数需要更复杂的数据表
// 这里返回简单估计
u = 1.0;
dulog = 0.0;
}
MpartfResult { u, dulog }
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
#[test]
fn test_mpartf_hydrogen() {
// H I at T = 10000 K
let result = mpartf(1, 1, 0, 10000.0);
assert!(result.u > 0.0);
assert!(result.dulog.is_finite());
}
#[test]
fn test_mpartf_helium() {
// He I at T = 10000 K
let result = mpartf(2, 1, 0, 10000.0);
assert!(result.u > 0.0);
}
#[test]
fn test_mpartf_high_temp() {
// 高温外推
let result = mpartf(1, 1, 0, 20000.0);
// H I: IGLE[1-1+1-1] = IGLE[0] = 2
assert_relative_eq!(result.u, 2.0, epsilon = 1e-10);
}
#[test]
#[should_panic(expected = "temperature")]
fn test_mpartf_low_temp() {
// 低温应该 panic
let _ = mpartf(1, 1, 0, 500.0);
}
#[test]
fn test_mpartf_carbon() {
// C I at T = 8000 K
let result = mpartf(6, 1, 0, 8000.0);
assert!(result.u > 1.0); // C I 配分函数应该 > 1
}
#[test]
fn test_mpartf_iron() {
// Fe I at T = 6000 K
// 注意:简化数据可能不包含完整 Fe 数据
let result = mpartf(26, 1, 0, 6000.0);
// 结果应该是正数
assert!(result.u > 0.0);
assert!(result.dulog.is_finite());
}
}
+298
View File
@@ -0,0 +1,298 @@
//! 氢线系列的 ODFOpacity Distribution Function)计算。
//!
//! 重构自 TLUSTY `odfhyd.f`
//!
//! 计算氢线系列的不透明度分布函数。
use super::divstr::divstr;
use super::indexx::indexx;
use super::odfhst::odfhst;
use crate::state::constants::{HALF, TWO, UN};
use crate::state::model::StrAux;
// 物理常量
/// 玻尔兹曼常数 / 氢质量
const CDOP: f64 = 2.0 * 1.38054e-16 / 1.6726e-24;
/// 光速 (Angstrom/s)
const CA: f64 = 2.997925e18;
/// 转换因子
const CCM: f64 = CA / 1.0e8;
/// 氢的 Rydberg 频率
const FRH: f64 = 3.28805e15;
/// Rydberg 波长 (Angstrom)
const RYDEL: f64 = 911.764;
/// 2/3
const TTW: f64 = 2.0 / 3.0;
/// Stark 加宽常量
const C00: f64 = 1.25e-9;
/// 仪器加宽常量
const CID: f64 = 0.02654;
/// ODFHYD 输入参数
pub struct OdfhydParams {
/// 深度索引 (1-indexed)
pub id: usize,
/// 跃迁索引 (1-indexed)
pub itr: usize,
}
/// ODFHYD 配置参数
pub struct OdfhydConfig {
/// ODF 采样标志 (0: 标准 ODF, >0: 采样 ODF)
pub ispodf: i32,
/// 最大谱线数
pub nlmx: usize,
/// 光速 (Angstrom/s)
pub cas: f64,
/// 湍流速度 (cm/s)
pub vtb: f64,
}
/// ODFHYD 原子数据
pub struct OdfhydAtomicData<'a> {
/// 下能级索引 (ntrans)
pub ilow: &'a [i32],
/// 上能级索引 (ntrans)
pub iup: &'a [i32],
/// 主量子数 (nlevel)
pub nquant: &'a [i32],
/// ODF 的下量子数 (nlevel)
pub nqlodf: &'a [i32],
/// 跃迁起始频率索引 (ntrans)
pub ifr0: &'a [i32],
/// 跃迁结束频率索引 (ntrans)
pub ifr1: &'a [i32],
/// XKIJ 系数 (ntrans × nlmx)
pub xkij: &'a [f64],
/// FIJ 系数 (ntrans × nlmx)
pub fij: &'a [f64],
/// ODF 索引 (ntrans)
pub jndodf: &'a [i32],
}
/// ODFHYD 模型状态
pub struct OdfhydModelState<'a> {
/// 温度 (nd)
pub temp: &'a [f64],
/// 电子密度 (nd)
pub elec: &'a [f64],
/// WNHINT 数组 (nlmx × nd)
pub wnhint: &'a [f64],
/// Stark 展宽参数
pub straux: StrAux,
}
/// ODFHYD ODF 数据
pub struct OdfhydOdfData<'a> {
/// ODF 频率数 (ntrans)
pub nfrodf: &'a [i32],
/// ODF 频率 (mfro × ntrans)
pub fros: &'a [f64],
/// ODF 频率宽度 (mfro × ntrans)
pub wnus: &'a [f64],
/// 谱线轮廓 (nd × nfreq 或其他)
pub prflin: &'a mut [f64],
/// 频率数组 (nfreq)
pub freq: &'a [f64],
/// KFR0 索引 (ntrans)
pub kfr0: &'a [i32],
/// INDEXP 标志 (ntrans)
pub indexp: &'a [i32],
}
/// XI2 函数:计算 1/n^2
fn xi2(n: i32) -> f64 {
1.0 / (n as f64 * n as f64)
}
/// 计算氢线系列的 ODF。
///
/// # 参数
///
/// * `params` - 输入参数(id, itr
/// * `config` - 配置参数
/// * `atomic` - 原子数据
/// * `model` - 模型状态
/// * `odf_data` - ODF 数据
///
/// # 返回值
///
/// 更新后的 PRFLIN 数组
pub fn odfhyd(
params: &OdfhydParams,
config: &OdfhydConfig,
atomic: &OdfhydAtomicData,
model: &mut OdfhydModelState,
odf_data: &mut OdfhydOdfData,
) {
let id = params.id;
let id_idx = id - 1;
let itr = params.itr;
let itr_idx = itr - 1;
let jo = atomic.jndodf[itr_idx] as usize - 1;
let ispodf = config.ispodf;
// 确定频率点数和初始化数组
let nf = if ispodf == 0 {
odf_data.nfrodf[jo] as usize
} else {
(atomic.ifr1[itr_idx] - atomic.ifr0[itr_idx] + 1) as usize
};
let mut sig = vec![0.0; nf];
let mut sgt = vec![0.0; nf];
let mut odf = vec![0.0; nf];
let mut iodf = vec![0; nf];
let mut ynus = vec![0.0; nf];
let mut alam = vec![0.0; nf];
// 初始化频率和波长
if ispodf == 0 {
for ij in 0..nf {
iodf[ij] = 0;
sig[ij] = 0.0;
odf[ij] = 0.0;
ynus[ij] = odf_data.fros[ij + jo * 1000]; // 假设 MFRO = 1000
alam[ij] = config.cas / ynus[ij];
}
} else {
let ifr0 = atomic.ifr0[itr_idx] as usize - 1;
for ij in 0..nf {
sig[ij] = 0.0;
ynus[ij] = odf_data.freq[ifr0 + ij];
alam[ij] = config.cas / ynus[ij];
}
}
// 获取能级索引
let ii = atomic.ilow[itr_idx] as usize - 1;
let jj = atomic.iup[itr_idx] as usize - 1;
// 计算电子密度因子
let anes = (TTW * model.elec[id_idx].ln()).exp();
let f00 = C00 * anes;
// 计算跃迁频率
let nquant_ii = atomic.nquant[ii];
let nquant_jj = atomic.nquant[jj];
let fra = FRH * (xi2(nquant_ii) - xi2(nquant_jj));
// 计算 Doppler 宽度
let dopo = fra / CCM * (CDOP * model.temp[id_idx] + config.vtb * config.vtb).sqrt();
// 遍历谱线系列
let nqlodf_ii = atomic.nqlodf[ii] as usize;
for j in nqlodf_ii..=config.nlmx {
let wl = RYDEL / (xi2(nquant_ii) - xi2(j as i32));
let fxk = f00 * atomic.xkij[jo * config.nlmx + j];
let dbeta = wl * wl / CA / fxk;
let betad = dbeta * dopo;
let fid = CID * atomic.fij[jo * config.nlmx + j] * dbeta;
// 调用 DIVSTR
let (adh, divh) = divstr(betad, 1);
// 获取 Stark 宽度
let wprob = model.wnhint[j * id + id_idx];
// 更新 straux 中的参数
model.straux.betad = betad;
model.straux.adh = adh;
model.straux.divh = divh;
// 调用 ODFHST
odfhst(nf, fxk, fid, wprob, wl, &alam, &model.straux, &mut sgt);
// 累加截面
for ij in 0..nf {
sig[ij] += sgt[ij];
}
}
// 后处理(仅对标准 ODF
if ispodf == 0 {
// 排序
iodf = indexx(&sig);
// 重新排列 ODF
for ij in 0..nf {
odf[ij] = sig[iodf[ij]];
}
// 计算频率网格
let i0 = atomic.ifr0[itr_idx] as usize;
let i1 = atomic.ifr1[itr_idx] as usize;
if odf_data.indexp[itr_idx].abs() == 2 {
ynus[0] = odf_data.freq[i0 - 1];
}
let mut iw1 = iodf[0];
for ij in 1..nf {
let iw2 = iodf[ij];
if ij > 1 && ij < nf - 1 {
ynus[ij] = ynus[ij - 1]
- HALF * (odf_data.wnus[iw1 + jo * 1000] + odf_data.wnus[iw2 + jo * 1000]);
} else if ij == 1 {
ynus[ij] = ynus[ij - 1]
- HALF * (TWO * odf_data.wnus[iw1 + jo * 1000] + odf_data.wnus[iw2 + jo * 1000]);
} else if ij == nf - 1 {
ynus[ij] = ynus[ij - 1]
- HALF * (odf_data.wnus[iw1 + jo * 1000] + TWO * odf_data.wnus[iw2 + jo * 1000]);
}
iw1 = iw2;
}
// 插值到频率网格
odf_data.prflin[id_idx * 100000 + i1 - 1] = 1e-35;
for ij0 in i0..i1 {
let mut ji = 1;
for ij in 1..nf {
ji = ij;
if ynus[ij] <= odf_data.freq[ij0 - 1] {
break;
}
}
let prfln = if ji > 0 && ji < nf {
odf[ji - 1]
+ (odf[ji] - odf[ji - 1]) * (odf_data.freq[ij0 - 1] - ynus[ji - 1])
/ (ynus[ji] - ynus[ji - 1])
} else {
0.0
};
odf_data.prflin[id_idx * 100000 + ij0 - 1] = prfln;
}
} else {
// 采样 ODF 情况
let kfr0 = odf_data.kfr0[itr_idx] as usize;
for ij in 0..nf {
odf_data.prflin[id_idx * 100000 + kfr0 + ij - 1] = sig[ij];
}
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_xi2() {
assert!((xi2(1) - 1.0).abs() < 1e-15);
assert!((xi2(2) - 0.25).abs() < 1e-15);
assert!((xi2(3) - 1.0 / 9.0).abs() < 1e-15);
}
#[test]
fn test_odfhyd_constants() {
// CDOP = 2 * BOLK / HMASS = 2 * 1.38054e-16 / 1.6726e-24 ≈ 1.65e8
// 注意:Fortran 中的 BOLK 可能是 1.38054e-16HMASS 是 1.6726e-24
// CDOP ≈ 1.65e8 cm/s/K^0.5
assert!(CDOP > 1e7 && CDOP < 1e9);
assert!((CA - 2.997925e18).abs() < 1e13);
assert!((FRH - 3.28805e15).abs() < 1e10);
}
}
+219
View File
@@ -0,0 +1,219 @@
//! 合并态超线的不透明度分布函数。
//!
//! 重构自 TLUSTY `odfmer.f`
//!
//! 计算合并态超线的 ODF(仅当温度变化足够大时更新)。
use super::odfhyd::{odfhyd, OdfhydAtomicData, OdfhydConfig, OdfhydModelState, OdfhydOdfData, OdfhydParams};
/// 温度变化阈值
const CHTL: f64 = 1e-3;
/// ODFMER 输入参数
pub struct OdfmerParams {
/// 初始化标志 (1: 初始化)
pub init: i32,
}
/// ODFMER 原子数据
pub struct OdfmerAtomicData<'a> {
/// 是否是谱线跃迁 (ntrans)
pub line: &'a [bool],
/// INDEXP 标志 (ntrans)
pub indexp: &'a [i32],
}
/// ODFMER 模型状态
pub struct OdfmerModelState<'a> {
/// 温度变化 (nd)
pub chant: &'a [f64],
}
/// 计算合并态超线的 ODF。
///
/// 仅当温度变化超过阈值 CHTL 时才重新计算。
///
/// # 参数
///
/// * `params` - 输入参数
/// * `atomic` - 原子数据
/// * `model` - 模型状态
/// * `odfhyd_config` - ODFHYD 配置
/// * `odfhyd_atomic` - ODFHYD 原子数据
/// * `odfhyd_model` - ODFHYD 模型状态
/// * `odfhyd_odf` - ODFHYD ODF 数据
pub fn odfmer(
params: &OdfmerParams,
atomic: &OdfmerAtomicData,
model: &OdfmerModelState,
odfhyd_config: &OdfhydConfig,
odfhyd_atomic: &OdfhydAtomicData,
odfhyd_model: &mut OdfhydModelState,
odfhyd_odf: &mut OdfhydOdfData,
) {
let ntrans = atomic.line.len();
let nd = model.chant.len();
for itr in 0..ntrans {
// 只处理谱线跃迁且 INDEXP == 2
if !atomic.line[itr] || atomic.indexp[itr].abs() != 2 {
continue;
}
for id in 0..nd {
// 仅在初始化或温度变化超过阈值时更新
if params.init == 1 || model.chant[id].abs() >= CHTL {
let odfhyd_params = OdfhydParams {
id: id + 1, // 1-indexed
itr: itr + 1, // 1-indexed
};
odfhyd(&odfhyd_params, odfhyd_config, odfhyd_atomic, odfhyd_model, odfhyd_odf);
}
}
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_odfmer_constant() {
assert!((CHTL - 1e-3).abs() < 1e-15);
}
#[test]
fn test_odfmer_temperature_threshold() {
// 验证温度变化阈值逻辑
let chtl = CHTL;
// 当温度变化小于阈值时,不应更新 ODF
let chant_small: f64 = 0.5e-3; // < CHTL
assert!(chant_small.abs() < chtl,
"Small temperature change should not trigger ODF update");
// 当温度变化大于等于阈值时,应更新 ODF
let chant_large: f64 = 1.5e-3; // > CHTL
assert!(chant_large.abs() >= chtl,
"Large temperature change should trigger ODF update");
}
#[test]
fn test_odfmer_initialization_flag() {
// 验证初始化标志逻辑
let init: i32 = 1;
// 初始化时应始终更新
assert!(init == 1, "init=1 should trigger update regardless of temperature change");
let init_normal: i32 = 0;
if init_normal != 1 {
// 正常迭代时,根据温度变化决定
let chant: f64 = 2e-3;
assert!(chant.abs() >= CHTL,
"Normal iteration: update only if temperature change exceeds threshold");
}
}
#[test]
fn test_odfmer_indexp_filter() {
// 验证 INDEXP 筛选逻辑
// 只处理 INDEXP == 2 或 INDEXP == -2 的跃迁
let indexp_values: [i32; 6] = [-2, -1, 0, 1, 2, 3];
let should_process: Vec<bool> = indexp_values.iter()
.map(|&idx| idx.abs() == 2)
.collect();
assert_eq!(should_process, vec![true, false, false, false, true, false],
"Only INDEXP = ±2 should be processed for merged states");
}
#[test]
fn test_odfmer_line_filter() {
// 验证谱线跃迁筛选
let line = true;
let indexp: i32 = 2;
let should_process = line && indexp.abs() == 2;
assert!(should_process, "Should process line transitions with INDEXP = ±2");
// 连续谱跃迁不应处理
let line_cont = false;
let indexp_neg: i32 = -2;
let should_skip = line_cont && indexp_neg.abs() == 2;
assert!(!should_skip || !line_cont, "Should skip continuum transitions");
}
#[test]
fn test_odfmer_depth_loop() {
// 验证深度点循环逻辑
let nd = 10;
let chant: Vec<f64> = vec![0.0; nd]; // 所有深度温度变化为零
let init: i32 = 0;
// 当温度变化为零且非初始化时,不应调用 ODFHYD
for id in 0..nd {
let should_update = init == 1 || chant[id].abs() >= CHTL;
assert!(!should_update,
"No update when temperature change is zero and not initializing");
}
// 部分深度有温度变化
let mut chant_partial = vec![0.0_f64; nd];
chant_partial[5] = 2e-3; // 第 6 个深度有变化
let update_count = chant_partial.iter()
.filter(|&&c| init == 1 || c.abs() >= CHTL)
.count();
assert_eq!(update_count, 1,
"Only one depth should trigger update");
}
#[test]
fn test_odfmer_transition_loop() {
// 验证跃迁循环逻辑
let ntrans = 5;
let line = vec![true, false, true, true, false];
let indexp: Vec<i32> = vec![2, 2, 1, -2, 2];
let processed: Vec<bool> = (0..ntrans)
.map(|itr| line[itr] && indexp[itr].abs() == 2)
.collect();
// 跃迁 0, 3 应该被处理
assert_eq!(processed, vec![true, false, false, true, false],
"Only transitions with line=true and INDEXP=±2 should be processed");
}
#[test]
fn test_odfmer_combined_logic() {
// 综合测试:验证跃迁和深度的组合筛选
let ntrans = 4;
let nd = 3;
let line = vec![true, true, false, true];
let indexp: Vec<i32> = vec![2, 1, 2, -2];
let chant: Vec<f64> = vec![0.0, 0.002, 0.001]; // 深度 1, 2 有变化
let init: i32 = 0;
// 计算应该被处理的 (跃迁, 深度) 对
let mut expected_updates = 0;
for itr in 0..ntrans {
if line[itr] && indexp[itr].abs() == 2 {
for id in 0..nd {
if init == 1 || chant[id].abs() >= CHTL {
expected_updates += 1;
}
}
}
}
// 跃迁 0, 3 应该被处理(满足 line=true 且 indexp=±2
// 深度 1, 2 应该被处理(温度变化 >= CHTL)
// 所以预期更新次数 = 2 跃迁 × 2 深度 = 4
assert_eq!(expected_updates, 4,
"Expected 4 ODF updates for the given configuration");
}
}
+405
View File
@@ -0,0 +1,405 @@
//! 所有深度点的吸收、发射和散射系数计算。
//!
//! 重构自 TLUSTY `opact1.f`
//!
//! 对于给定频率点,计算所有深度点的吸收、发射和散射系数。
use super::opctab::{opctab, OpctabParams, OpctabTableData, OpctabModelState, OpctabOutput};
use crate::state::constants::{HK, UN};
/// OPACT1 输入参数
pub struct Opact1Params<'a> {
/// 频率索引 (1-indexed)
pub ij: usize,
/// 迭代次数
pub iter: i32,
/// Rayleigh 散射标志 (<0: 简单公式, >0: 调用 rayleigh, =0: 关闭)
pub ifrayl: i32,
/// 选项表标志 (<0: 添加电子散射, >0: 使用选项表)
pub ioptab: i32,
/// Rayleigh 参数 (当 ifrayl > 0 时需要)
pub rayleigh_params: Option<&'a super::rayleigh::RayleighParams<'a>>,
}
/// OPACT1 模型状态
pub struct Opact1ModelState<'a> {
/// 温度 (nd)
pub temp: &'a [f64],
/// 密度 (nd)
pub dens: &'a [f64],
/// 频率数组 (nfreq)
pub freq: &'a [f64],
/// Planck 函数 (nfreq)
pub bnue: &'a [f64],
/// HKT1 数组 (nd) - HK/T
pub hkt1: &'a mut [f64],
/// XKF 数组 (nd)
pub xkf: &'a mut [f64],
/// XKF1 数组 (nd)
pub xkf1: &'a mut [f64],
/// XKFB 数组 (nd)
pub xkfb: &'a mut [f64],
}
/// OPACT1 输出状态
pub struct Opact1OutputState<'a> {
/// 吸收系数 (nd)
pub abso1: &'a mut [f64],
/// 发射系数 (nd)
pub emis1: &'a mut [f64],
/// 散射系数 (nd)
pub scat1: &'a mut [f64],
/// 累积吸收系数 (nd) - 用于选项表
pub absot: &'a mut [f64],
}
/// 计算所有深度点的吸收、发射和散射系数。
///
/// # 参数
///
/// * `params` - 输入参数
/// * `model` - 模型状态
/// * `output` - 输出状态
/// * `table` - 不透明度表数据
/// * `opctab_model` - OPCTAB 模型状态
pub fn opact1(
params: &Opact1Params,
model: &mut Opact1ModelState,
output: &mut Opact1OutputState,
table: &OpctabTableData,
opctab_model: &mut OpctabModelState,
) {
let ij = params.ij;
let ij_idx = ij - 1;
let fr = model.freq[ij_idx];
let nd = model.temp.len();
for id in 0..nd {
let t = model.temp[id];
let rho = model.dens[id];
// 更新 HKT1, XKF, XKF1, XKFB
model.hkt1[id] = HK / t;
model.xkf[id] = (-model.hkt1[id] * fr).exp();
model.xkf1[id] = UN - model.xkf[id];
model.xkfb[id] = model.xkf[id] * model.bnue[ij_idx];
let plan = model.xkfb[id] / model.xkf1[id];
// 调用 OPCTAB
let opctab_params = OpctabParams {
fr,
ij,
id: id + 1, // 1-indexed
t,
rho,
igram: 0, // 返回每体积
iter: params.iter,
ifrayl: params.ifrayl,
ioptab: params.ioptab,
rayleigh_params: params.rayleigh_params,
};
let OpctabOutput { ab, sc: _, sct } = opctab(&opctab_params, table, opctab_model);
if params.ioptab < 0 {
output.abso1[id] = ab + sct;
output.scat1[id] = sct;
output.absot[id] += output.abso1[id] / model.dens[id];
} else if params.ioptab > 0 {
output.abso1[id] += ab;
}
output.emis1[id] += ab * plan;
}
}
#[cfg(test)]
mod tests {
use super::*;
use super::super::opctab::{OpctabTableData, OpctabModelState};
use approx::assert_relative_eq;
#[test]
fn test_opact1_function_basic() {
// numtemp == nd 时使用直接访问路径
let numtemp = 2;
let nd = 2;
let nfreq = 2;
let max_numrh = 1;
let tempvec = vec![9.2103, 9.3927];
let numrh = vec![1, 1];
let rhomat = vec![-16.1181, -15.9];
// absopac: (numtemp × max_numrh × nfreq) = 2 × 1 × 2 = 4 个元素
let absopac = vec![-5.0, -4.5, -5.2, -4.7];
let raysc = vec![1e-20, 1.2e-20];
let freq = vec![1e15, 2e15];
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 temp: Vec<f64> = vec![10000.0, 12000.0];
let dens: Vec<f64> = vec![1e-7, 1.5e-7];
let bnue: Vec<f64> = vec![1e10, 5e10];
let mut hkt1 = vec![0.0; nd];
let mut xkf = vec![0.0; nd];
let mut xkf1 = vec![0.0; nd];
let mut xkfb = vec![0.0; nd];
let mut abso1 = vec![0.0; nd];
let mut emis1 = vec![0.0; nd];
let mut scat1 = vec![0.0; nd];
let mut absot = vec![0.0; nd];
let elec = vec![1e-10, 1.2e-10];
let params = Opact1Params {
ij: 1,
iter: 1,
ifrayl: -1,
ioptab: -1,
rayleigh_params: None,
};
let mut model = Opact1ModelState {
temp: &temp,
dens: &dens,
freq: &freq,
bnue: &bnue,
hkt1: &mut hkt1,
xkf: &mut xkf,
xkf1: &mut xkf1,
xkfb: &mut xkfb,
};
let mut output = Opact1OutputState {
abso1: &mut abso1,
emis1: &mut emis1,
scat1: &mut scat1,
absot: &mut absot,
};
let mut opctab_model = OpctabModelState {
elec: &elec,
dens: &dens,
raysct: None,
eospar: None,
};
// 调用函数
opact1(&params, &mut model, &mut output, &table, &mut opctab_model);
// 验证输出
for id in 0..nd {
// hkt1 应该被正确计算
assert_relative_eq!(model.hkt1[id], HK / temp[id], epsilon = 1e-15);
// xkf = exp(-hkt1 * freq)
let expected_xkf = (-model.hkt1[id] * freq[0]).exp();
assert_relative_eq!(model.xkf[id], expected_xkf, epsilon = 1e-15);
// 吸收系数应该为正
assert!(output.abso1[id] > 0.0, "Absorption coefficient should be positive");
assert!(output.emis1[id] >= 0.0, "Emission coefficient should be non-negative");
}
}
#[test]
fn test_opact1_function_planck_relation() {
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![-5.0];
let raysc = vec![1e-20];
let freq = vec![1e15];
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 temp: Vec<f64> = vec![10000.0];
let dens: Vec<f64> = vec![1e-7];
let bnue: Vec<f64> = vec![1e10];
let mut hkt1 = vec![0.0; nd];
let mut xkf = vec![0.0; nd];
let mut xkf1 = vec![0.0; nd];
let mut xkfb = vec![0.0; nd];
let mut abso1 = vec![0.0; nd];
let mut emis1 = vec![0.0; nd];
let mut scat1 = vec![0.0; nd];
let mut absot = vec![0.0; nd];
let elec = vec![1e-10];
let params = Opact1Params {
ij: 1,
iter: 1,
ifrayl: -1,
ioptab: -1,
rayleigh_params: None,
};
let mut model = Opact1ModelState {
temp: &temp,
dens: &dens,
freq: &freq,
bnue: &bnue,
hkt1: &mut hkt1,
xkf: &mut xkf,
xkf1: &mut xkf1,
xkfb: &mut xkfb,
};
let mut output = Opact1OutputState {
abso1: &mut abso1,
emis1: &mut emis1,
scat1: &mut scat1,
absot: &mut absot,
};
let mut opctab_model = OpctabModelState {
elec: &elec,
dens: &dens,
raysct: None,
eospar: None,
};
opact1(&params, &mut model, &mut output, &table, &mut opctab_model);
// 验证 Planck 源函数关系: emis1 = ab * plan = ab * xkfb / xkf1
let ab = (-5.0f64).exp() * dens[0];
let plan = model.xkfb[0] / model.xkf1[0];
let expected_emis = ab * plan;
assert_relative_eq!(output.emis1[0], expected_emis, epsilon = 1e-10);
}
#[test]
fn test_opact1_function_multiple_depths() {
// numtemp == nd 时使用直接访问路径
let numtemp = 3;
let nd = 3;
let nfreq = 1;
let max_numrh = 1;
let tempvec = vec![8.987, 9.2103, 9.615];
let numrh = vec![1, 1, 1];
let rhomat = vec![-18.42, -16.1181, -13.8];
// absopac: (numtemp × max_numrh × nfreq) = 3 × 1 × 1 = 3 个元素
let absopac = vec![-5.0, -4.8, -4.5];
let raysc = vec![1e-20, 1.1e-20, 1.2e-20];
let freq = vec![1e15];
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 temp: Vec<f64> = vec![8000.0, 10000.0, 15000.0];
let dens: Vec<f64> = vec![1e-8, 1e-7, 1e-6];
let bnue: Vec<f64> = vec![1e10];
let mut hkt1 = vec![0.0; nd];
let mut xkf = vec![0.0; nd];
let mut xkf1 = vec![0.0; nd];
let mut xkfb = vec![0.0; nd];
let mut abso1 = vec![0.0; nd];
let mut emis1 = vec![0.0; nd];
let mut scat1 = vec![0.0; nd];
let mut absot = vec![0.0; nd];
let elec = vec![1e-12, 1e-10, 1e-8];
let params = Opact1Params {
ij: 1,
iter: 1,
ifrayl: -1,
ioptab: -1,
rayleigh_params: None,
};
let mut model = Opact1ModelState {
temp: &temp,
dens: &dens,
freq: &freq,
bnue: &bnue,
hkt1: &mut hkt1,
xkf: &mut xkf,
xkf1: &mut xkf1,
xkfb: &mut xkfb,
};
let mut output = Opact1OutputState {
abso1: &mut abso1,
emis1: &mut emis1,
scat1: &mut scat1,
absot: &mut absot,
};
let mut opctab_model = OpctabModelState {
elec: &elec,
dens: &dens,
raysct: None,
eospar: None,
};
opact1(&params, &mut model, &mut output, &table, &mut opctab_model);
// 验证所有深度点都被处理
for id in 0..nd {
assert!(output.abso1[id] > 0.0, "Depth {} should have positive absorption", id);
assert!(output.emis1[id] >= 0.0, "Depth {} should have non-negative emission", id);
assert!(output.absot[id] > 0.0, "Depth {} should have positive cumulative absorption", id);
}
// 验证温度越高,hkt1 越小
assert!(model.hkt1[0] > model.hkt1[1], "hkt1 should decrease with temperature");
assert!(model.hkt1[1] > model.hkt1[2], "hkt1 should decrease with temperature");
}
}
+672
View File
@@ -0,0 +1,672 @@
//! 吸收、发射和散射系数及其导数计算。
//!
//! 重构自 TLUSTY `opactd.f`
//!
//! 与 OPACT1 类似,但额外计算温度和密度导数。
use super::opctab::{opctab, OpctabParams, OpctabTableData, OpctabModelState, OpctabOutput};
use crate::state::constants::{HK, UN};
/// 微分步长
const DELT: f64 = 1e-3;
const DELR: f64 = 1e-3;
/// OPACTD 输入参数
pub struct OpactdParams<'a> {
/// 频率索引 (1-indexed)
pub ij: usize,
/// Rayleigh/B 矩阵标志 (0: 普通模式, >0: 计算密度导数)
pub ifryb: i32,
/// 能量方程标志 (<=0: 密度不是状态参数)
pub inhe: i32,
/// 迭代次数
pub iter: i32,
/// Rayleigh 散射标志 (<0: 简单公式, >0: 调用 rayleigh, =0: 关闭)
pub ifrayl: i32,
/// 选项表标志
pub ioptab: i32,
/// Rayleigh 参数
pub rayleigh_params: Option<&'a super::rayleigh::RayleighParams<'a>>,
}
/// OPACTD 模型状态
pub struct OpactdModelState<'a> {
/// 温度 (nd)
pub temp: &'a [f64],
/// 密度 (nd)
pub dens: &'a [f64],
/// 频率数组 (nfreq)
pub freq: &'a [f64],
/// Planck 函数 (nfreq)
pub bnue: &'a [f64],
/// HKT1 数组 (nd)
pub hkt1: &'a [f64],
/// 密度对温度的导数 (nd)
pub drhodt: &'a [f64],
/// XKF 数组 (nd)
pub xkf: &'a mut [f64],
/// XKF1 数组 (nd)
pub xkf1: &'a mut [f64],
/// XKFB 数组 (nd)
pub xkfb: &'a mut [f64],
}
/// OPACTD 输出状态
pub struct OpactdOutputState<'a> {
/// 吸收系数 (nd)
pub abso1: &'a mut [f64],
/// 发射系数 (nd)
pub emis1: &'a mut [f64],
/// 散射系数 (nd)
pub scat1: &'a mut [f64],
/// 累积吸收系数 (nd)
pub absot: &'a mut [f64],
/// 吸收系数温度导数 (nd)
pub dabt1: &'a mut [f64],
/// 发射系数温度导数 (nd)
pub demt1: &'a mut [f64],
/// 散射系数温度导数 (nd)
pub dsct1: &'a mut [f64],
/// 吸收系数密度导数 (nd)
pub dabn1: &'a mut [f64],
/// 发射系数密度导数 (nd)
pub demn1: &'a mut [f64],
/// 散射系数密度导数 (nd)
pub dscn1: &'a mut [f64],
}
/// OPACTD 显式频率数据
pub struct OpactdExpData<'a> {
/// 显式频率索引 (nfreq)
pub ijex: &'a [i32],
/// 显式吸收系数 (nfreqe × nd)
pub absoex: &'a mut [f64],
/// 显式散射系数 (nfreqe × nd)
pub scatex: &'a mut [f64],
/// 显式发射系数 (nfreqe × nd)
pub emisex: &'a mut [f64],
/// 显式吸收温度导数 (nfreqe × nd)
pub dabtex: &'a mut [f64],
/// 显式发射温度导数 (nfreqe × nd)
pub demtex: &'a mut [f64],
/// 显式吸收密度导数 (nfreqe × nd)
pub dabnex: &'a mut [f64],
/// 显式发射密度导数 (nfreqe × nd)
pub demnex: &'a mut [f64],
/// 显式频率点数
pub nfreqe: usize,
/// 深度点数
pub nd: usize,
}
/// 计算吸收、发射和散射系数及其导数。
///
/// # 参数
///
/// * `params` - 输入参数
/// * `model` - 模型状态
/// * `output` - 输出状态
/// * `exp_data` - 显式频率数据(可选)
/// * `table` - 不透明度表数据
/// * `opctab_model` - OPCTAB 模型状态
pub fn opactd(
params: &OpactdParams,
model: &mut OpactdModelState,
output: &mut OpactdOutputState,
exp_data: Option<&mut OpactdExpData>,
table: &OpctabTableData,
opctab_model: &mut OpctabModelState,
) {
let ij = params.ij;
let ij_idx = ij - 1;
let fr = model.freq[ij_idx];
let nd = model.temp.len();
// 确定模式
let imodf = if params.ifryb > 0 { 1 } else { 0 };
for id in 0..nd {
let t = model.temp[id];
let t1 = t * (UN + DELT);
let rho = model.dens[id];
let rho1 = rho * (UN + DELR);
// 更新 XKF, XKF1, XKFB
model.xkf[id] = (-model.hkt1[id] * fr).exp();
model.xkf1[id] = UN - model.xkf[id];
model.xkfb[id] = model.xkf[id] * model.bnue[ij_idx];
let plan = model.xkfb[id] / model.xkf1[id];
let dplan = plan / model.xkf1[id] * model.hkt1[id] * fr / t;
// 调用 OPCTAB 三次:原始、温度扰动、密度扰动
let opctab_params_base = OpctabParams {
fr,
ij,
id: id + 1,
t,
rho,
igram: imodf,
iter: params.iter,
ifrayl: params.ifrayl,
ioptab: params.ioptab,
rayleigh_params: params.rayleigh_params,
};
// 原始值
let OpctabOutput { ab, sc: _, sct } = opctab(&opctab_params_base, table, opctab_model);
// 温度扰动
let opctab_params_t = OpctabParams {
t: t1,
..opctab_params_base
};
let OpctabOutput { ab: ab1, sc: _, sct: sct1 } =
opctab(&opctab_params_t, table, opctab_model);
// 密度扰动
let opctab_params_rho = OpctabParams {
rho: rho1,
..opctab_params_base
};
let OpctabOutput { ab: ab2, sc: _, sct: sct2 } =
opctab(&opctab_params_rho, table, opctab_model);
// 存储系数
output.abso1[id] = ab + sct;
output.scat1[id] = sct;
output.emis1[id] = ab * plan;
output.absot[id] = if imodf == 0 {
output.abso1[id] / model.dens[id]
} else {
output.abso1[id]
};
// 温度导数
output.dabt1[id] = (ab1 - ab) / t / DELT;
output.demt1[id] = ab * dplan + output.dabt1[id] * plan;
output.dsct1[id] = (sct1 - sct) / t / DELT;
output.dabt1[id] += output.dsct1[id];
if params.ifryb > 0 {
// 密度导数
output.dabn1[id] = (ab2 - ab) / rho / DELR;
output.demn1[id] = output.dabn1[id] * plan;
output.dscn1[id] = (sct2 - sct) / rho / DELR;
output.dabn1[id] += output.dscn1[id];
// 如果密度不是状态参数,修改温度导数
if params.inhe <= 0 {
output.dabt1[id] += output.dabn1[id] * model.drhodt[id];
output.demt1[id] += output.demn1[id] * model.drhodt[id];
output.dsct1[id] += output.dscn1[id] * model.drhodt[id];
output.dabn1[id] = 0.0;
output.demn1[id] = 0.0;
output.dscn1[id] = 0.0;
}
}
}
// 存储显式频率数据
if let Some(exp) = exp_data {
let ijex_val = if ij_idx < exp.ijex.len() { exp.ijex[ij_idx] } else { 0 };
if ijex_val > 0 && params.ifryb <= 0 {
let ije = (ijex_val - 1) as usize;
for id in 0..nd {
exp.absoex[ije * exp.nd + id] = output.abso1[id];
exp.scatex[ije * exp.nd + id] = output.scat1[id];
exp.emisex[ije * exp.nd + id] = output.emis1[id];
exp.dabtex[ije * exp.nd + id] = output.dabt1[id];
exp.demtex[ije * exp.nd + id] = output.demt1[id];
exp.dabnex[ije * exp.nd + id] = output.dabn1[id];
exp.demnex[ije * exp.nd + id] = output.demn1[id];
}
}
}
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
#[test]
fn test_opactd_function_basic() {
// 创建测试数据 - numtemp == nd 使用直接访问路径
let numtemp = 2;
let nd = 2;
let nfreq = 2;
let max_numrh = 1;
let tempvec = vec![9.2103, 9.3927];
let numrh = vec![1, 1];
let rhomat = vec![-16.1181, -15.9];
let absopac = vec![-5.0, -4.5, -5.2, -4.7];
let raysc = vec![1e-20, 1.2e-20];
let freq = vec![1e15, 2e15];
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 temp: Vec<f64> = vec![10000.0, 12000.0];
let dens: Vec<f64> = vec![1e-7, 1.5e-7];
let bnue: Vec<f64> = vec![1e10, 5e10];
let hkt1: Vec<f64> = temp.iter().map(|t| HK / t).collect();
let drhodt: Vec<f64> = vec![0.0; nd];
let mut xkf = vec![0.0; nd];
let mut xkf1 = vec![0.0; nd];
let mut xkfb = vec![0.0; nd];
let mut abso1 = vec![0.0; nd];
let mut emis1 = vec![0.0; nd];
let mut scat1 = vec![0.0; nd];
let mut absot = vec![0.0; nd];
let mut dabt1 = vec![0.0; nd];
let mut demt1 = vec![0.0; nd];
let mut dsct1 = vec![0.0; nd];
let mut dabn1 = vec![0.0; nd];
let mut demn1 = vec![0.0; nd];
let mut dscn1 = vec![0.0; nd];
let elec = vec![1e-10, 1.2e-10];
let params = OpactdParams {
ij: 1,
ifryb: 0, // 不计算密度导数
inhe: 0,
iter: 1,
ifrayl: -1,
ioptab: -1,
rayleigh_params: None,
};
let mut model = OpactdModelState {
temp: &temp,
dens: &dens,
freq: &freq,
bnue: &bnue,
hkt1: &hkt1,
drhodt: &drhodt,
xkf: &mut xkf,
xkf1: &mut xkf1,
xkfb: &mut xkfb,
};
let mut output = OpactdOutputState {
abso1: &mut abso1,
emis1: &mut emis1,
scat1: &mut scat1,
absot: &mut absot,
dabt1: &mut dabt1,
demt1: &mut demt1,
dsct1: &mut dsct1,
dabn1: &mut dabn1,
demn1: &mut demn1,
dscn1: &mut dscn1,
};
let mut opctab_model = OpctabModelState {
elec: &elec,
dens: &dens,
raysct: None,
eospar: None,
};
// 调用函数
opactd(&params, &mut model, &mut output, None, &table, &mut opctab_model);
// 验证输出
for id in 0..nd {
// 吸收系数应该为正
assert!(output.abso1[id] > 0.0, "Absorption coefficient should be positive");
assert!(output.emis1[id] >= 0.0, "Emission coefficient should be non-negative");
// 温度导数应该是有限的
assert!(output.dabt1[id].is_finite(), "Temperature derivative should be finite");
assert!(output.demt1[id].is_finite(), "Emission temperature derivative should be finite");
// xkf 应该被正确计算
let expected_xkf = (-hkt1[id] * freq[0]).exp();
assert_relative_eq!(model.xkf[id], expected_xkf, epsilon = 1e-15);
}
}
#[test]
fn test_opactd_function_with_density_derivatives() {
// 测试密度导数计算 (ifryb > 0)
let numtemp = 2;
let nd = 2;
let nfreq = 2;
let max_numrh = 1;
let tempvec = vec![9.2103, 9.3927];
let numrh = vec![1, 1];
let rhomat = vec![-16.1181, -15.9];
let absopac = vec![-5.0, -4.5, -5.2, -4.7];
let raysc = vec![1e-20, 1.2e-20];
let freq = vec![1e15, 2e15];
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 temp: Vec<f64> = vec![10000.0, 12000.0];
let dens: Vec<f64> = vec![1e-7, 1.5e-7];
let bnue: Vec<f64> = vec![1e10, 5e10];
let hkt1: Vec<f64> = temp.iter().map(|t| HK / t).collect();
let drhodt: Vec<f64> = vec![1e-12, 1e-12]; // 密度对温度的导数
let mut xkf = vec![0.0; nd];
let mut xkf1 = vec![0.0; nd];
let mut xkfb = vec![0.0; nd];
let mut abso1 = vec![0.0; nd];
let mut emis1 = vec![0.0; nd];
let mut scat1 = vec![0.0; nd];
let mut absot = vec![0.0; nd];
let mut dabt1 = vec![0.0; nd];
let mut demt1 = vec![0.0; nd];
let mut dsct1 = vec![0.0; nd];
let mut dabn1 = vec![0.0; nd];
let mut demn1 = vec![0.0; nd];
let mut dscn1 = vec![0.0; nd];
let elec = vec![1e-10, 1.2e-10];
let params = OpactdParams {
ij: 1,
ifryb: 1, // 计算密度导数
inhe: 1, // 密度是状态参数
iter: 1,
ifrayl: -1,
ioptab: -1,
rayleigh_params: None,
};
let mut model = OpactdModelState {
temp: &temp,
dens: &dens,
freq: &freq,
bnue: &bnue,
hkt1: &hkt1,
drhodt: &drhodt,
xkf: &mut xkf,
xkf1: &mut xkf1,
xkfb: &mut xkfb,
};
let mut output = OpactdOutputState {
abso1: &mut abso1,
emis1: &mut emis1,
scat1: &mut scat1,
absot: &mut absot,
dabt1: &mut dabt1,
demt1: &mut demt1,
dsct1: &mut dsct1,
dabn1: &mut dabn1,
demn1: &mut demn1,
dscn1: &mut dscn1,
};
let mut opctab_model = OpctabModelState {
elec: &elec,
dens: &dens,
raysct: None,
eospar: None,
};
opactd(&params, &mut model, &mut output, None, &table, &mut opctab_model);
// 验证密度导数被计算
for id in 0..nd {
assert!(output.dabn1[id].is_finite(), "Density derivative should be finite");
assert!(output.demn1[id].is_finite(), "Emission density derivative should be finite");
// 当 ifryb > 0 且 inhe > 0 时,密度导数不应该被清零
}
}
#[test]
fn test_opactd_function_density_temperature_coupling() {
// 测试密度-温度耦合 (inhe <= 0)
let numtemp = 2;
let nd = 2;
let nfreq = 2;
let max_numrh = 1;
let tempvec = vec![9.2103, 9.3927];
let numrh = vec![1, 1];
let rhomat = vec![-16.1181, -15.9];
let absopac = vec![-5.0, -4.5, -5.2, -4.7];
let raysc = vec![1e-20, 1.2e-20];
let freq = vec![1e15, 2e15];
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 temp: Vec<f64> = vec![10000.0, 12000.0];
let dens: Vec<f64> = vec![1e-7, 1.5e-7];
let bnue: Vec<f64> = vec![1e10, 5e10];
let hkt1: Vec<f64> = temp.iter().map(|t| HK / t).collect();
let drhodt: Vec<f64> = vec![-1e-12, -1e-12]; // 负的 d(rho)/dT
let mut xkf = vec![0.0; nd];
let mut xkf1 = vec![0.0; nd];
let mut xkfb = vec![0.0; nd];
let mut abso1 = vec![0.0; nd];
let mut emis1 = vec![0.0; nd];
let mut scat1 = vec![0.0; nd];
let mut absot = vec![0.0; nd];
let mut dabt1 = vec![0.0; nd];
let mut demt1 = vec![0.0; nd];
let mut dsct1 = vec![0.0; nd];
let mut dabn1 = vec![0.0; nd];
let mut demn1 = vec![0.0; nd];
let mut dscn1 = vec![0.0; nd];
let elec = vec![1e-10, 1.2e-10];
let params = OpactdParams {
ij: 1,
ifryb: 1, // 计算密度导数
inhe: 0, // 密度不是状态参数 - 应触发耦合
iter: 1,
ifrayl: -1,
ioptab: -1,
rayleigh_params: None,
};
let mut model = OpactdModelState {
temp: &temp,
dens: &dens,
freq: &freq,
bnue: &bnue,
hkt1: &hkt1,
drhodt: &drhodt,
xkf: &mut xkf,
xkf1: &mut xkf1,
xkfb: &mut xkfb,
};
let mut output = OpactdOutputState {
abso1: &mut abso1,
emis1: &mut emis1,
scat1: &mut scat1,
absot: &mut absot,
dabt1: &mut dabt1,
demt1: &mut demt1,
dsct1: &mut dsct1,
dabn1: &mut dabn1,
demn1: &mut demn1,
dscn1: &mut dscn1,
};
let mut opctab_model = OpctabModelState {
elec: &elec,
dens: &dens,
raysct: None,
eospar: None,
};
opactd(&params, &mut model, &mut output, None, &table, &mut opctab_model);
// 验证密度导数被清零(因为 inhe <= 0
for id in 0..nd {
assert_relative_eq!(output.dabn1[id], 0.0, epsilon = 1e-20);
assert_relative_eq!(output.demn1[id], 0.0, epsilon = 1e-20);
assert_relative_eq!(output.dscn1[id], 0.0, epsilon = 1e-20);
}
}
#[test]
fn test_opactd_function_planck_relation() {
// 验证 Planck 源函数关系
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![-5.0];
let raysc = vec![1e-20];
let freq = vec![1e15];
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 temp: Vec<f64> = vec![10000.0];
let dens: Vec<f64> = vec![1e-7];
let bnue: Vec<f64> = vec![1e10];
let hkt1: Vec<f64> = temp.iter().map(|t| HK / t).collect();
let drhodt: Vec<f64> = vec![0.0];
let mut xkf = vec![0.0; nd];
let mut xkf1 = vec![0.0; nd];
let mut xkfb = vec![0.0; nd];
let mut abso1 = vec![0.0; nd];
let mut emis1 = vec![0.0; nd];
let mut scat1 = vec![0.0; nd];
let mut absot = vec![0.0; nd];
let mut dabt1 = vec![0.0; nd];
let mut demt1 = vec![0.0; nd];
let mut dsct1 = vec![0.0; nd];
let mut dabn1 = vec![0.0; nd];
let mut demn1 = vec![0.0; nd];
let mut dscn1 = vec![0.0; nd];
let elec = vec![1e-10];
let params = OpactdParams {
ij: 1,
ifryb: 0,
inhe: 0,
iter: 1,
ifrayl: -1,
ioptab: -1,
rayleigh_params: None,
};
let mut model = OpactdModelState {
temp: &temp,
dens: &dens,
freq: &freq,
bnue: &bnue,
hkt1: &hkt1,
drhodt: &drhodt,
xkf: &mut xkf,
xkf1: &mut xkf1,
xkfb: &mut xkfb,
};
let mut output = OpactdOutputState {
abso1: &mut abso1,
emis1: &mut emis1,
scat1: &mut scat1,
absot: &mut absot,
dabt1: &mut dabt1,
demt1: &mut demt1,
dsct1: &mut dsct1,
dabn1: &mut dabn1,
demn1: &mut demn1,
dscn1: &mut dscn1,
};
let mut opctab_model = OpctabModelState {
elec: &elec,
dens: &dens,
raysct: None,
eospar: None,
};
opactd(&params, &mut model, &mut output, None, &table, &mut opctab_model);
// 验证发射系数 = ab * plan
// ab 从 opctab 返回 = exp(absopac) * rho (当 igram=0)
let ab_per_vol = (-5.0f64).exp() * dens[0];
let plan = model.xkfb[0] / model.xkf1[0];
let expected_emis = ab_per_vol * plan;
assert_relative_eq!(output.emis1[0], expected_emis, epsilon = 1e-10);
}
#[test]
fn test_opactd_constants() {
assert!((DELT - 1e-3).abs() < 1e-15);
assert!((DELR - 1e-3).abs() < 1e-15);
}
}
+374
View File
@@ -0,0 +1,374 @@
//! 部分重分布(PRD)线发射和散射系数修正。
//!
//! 重构自 TLUSTY `prd.f`
//!
//! 在 PRD 情况下修改线发射系数和散射系数。
use super::gami::gami;
use crate::state::constants::{TWO, UN};
// 物理常量
/// Einstein A21 系数
const A21: f64 = 4.699e8;
/// 2π
const PI2: f64 = 6.28318531;
/// 辐射阻尼常量
const GR: f64 = 2.0 * 4.8e-8;
/// PRD 输入参数
pub struct PrdParams {
/// 频率索引 (1-indexed)
pub ij: usize,
}
/// PRD 配置参数
pub struct PrdConfig {
/// ODF 采样标志
pub ispodf: i32,
/// 深度点数
pub nd: usize,
/// PRD 分离阈值
pub xpdiv: f64,
}
/// PRD 原子数据
pub struct PrdAtomicData<'a> {
/// 跃迁的谱线索引 (nfreq)
pub ijlin: &'a [i32],
/// 跃迁的 PRD 索引 (ntrans)
pub iprd: &'a [i32],
/// 跃迁频率 (ntrans)
pub fr0: &'a [f64],
/// 下能级索引 (ntrans)
pub ilow: &'a [i32],
/// 跃迁起始频率索引 (ntrans)
pub ifr0: &'a [i32],
/// 跃迁结束频率索引 (ntrans)
pub ifr1: &'a [i32],
/// KFR0 索引 (ntrans)
pub kfr0: &'a [i32],
/// INDEXP 标志 (ntrans)
pub indexp: &'a [i32],
/// H- 的第一个能级索引
pub nfirst_elh: usize,
}
/// PRD 模型状态
pub struct PrdModelState<'a> {
/// 温度 (nd)
pub temp: &'a [f64],
/// 电子密度 (nd)
pub elec: &'a [f64],
/// 吸收系数 (ntrans × nd)
pub abtra: &'a [f64],
/// 发射系数 (ntrans × nd)
pub emtra: &'a [f64],
/// 谱线轮廓 (nd × nfreq)
pub prflin: &'a [f64],
/// Doppler 宽度 (ntrans_prd × nd)
pub doptr: &'a [f64],
/// 相干性因子 (ntrans_prd × nd)
pub coher: &'a mut [f64],
/// 占据数 (nlevel × nd)
pub popul: &'a [f64],
/// XKFB 数组 (nd)
pub xkfb: &'a [f64],
/// 散射系数 (nd)
pub scat1: &'a mut [f64],
/// 发射系数 (nd)
pub emis1: &'a mut [f64],
}
/// PRD 频率数据
pub struct PrdFreqData<'a> {
/// 频率数组 (nfreq)
pub freq: &'a [f64],
/// 主谱线索引 (nfreq)
pub ijlin: &'a [i32],
/// 重叠谱线数 (nfreq)
pub nlines: &'a [i32],
/// 重叠谱线索引 (nliness × nfreq)
pub itrlin: &'a [i32],
}
/// 在 PRD 情况下修改线发射系数和散射系数。
///
/// # 参数
///
/// * `params` - 输入参数
/// * `config` - 配置参数
/// * `atomic` - 原子数据
/// * `model` - 模型状态
/// * `freq_data` - 频率数据
pub fn prd(
params: &PrdParams,
config: &PrdConfig,
atomic: &PrdAtomicData,
model: &mut PrdModelState,
freq_data: &PrdFreqData,
) {
let ij = params.ij;
if ij == 0 {
return;
}
let ij_idx = ij - 1;
let fr = freq_data.freq[ij_idx];
let nd = config.nd;
if config.ispodf == 0 {
// 标准 ODF 情况
// 处理主谱线
if freq_data.ijlin[ij_idx] > 0 {
let itr = (freq_data.ijlin[ij_idx] - 1) as usize;
let itrprd = atomic.iprd[itr];
if itrprd > 0 {
let dfr = (fr - atomic.fr0[itr]).abs();
// 检查是否是 H- 的跃迁
if atomic.ilow[itr] as usize == atomic.nfirst_elh {
let omeg = dfr * PI2;
let gra = A21 + GR * model.popul[atomic.nfirst_elh * nd];
for id in 0..nd {
model.coher[(itrprd as usize - 1) * nd + id] =
A21 / (gra + gami(2, "elec", omeg, model.temp[id], model.elec[id]));
}
}
for id in 0..nd {
let sg = model.prflin[id * 100000 + ij_idx];
let sg_final =
if dfr / model.doptr[(itrprd as usize - 1) * nd + id] <= config.xpdiv {
0.0
} else {
sg
};
let scalin = sg_final
* model.abtra[itr * nd + id]
* model.coher[(itrprd as usize - 1) * nd + id];
model.scat1[id] += scalin;
let scem = sg_final
* model.emtra[itr * nd + id]
* model.coher[(itrprd as usize - 1) * nd + id]
* model.xkfb[id];
model.emis1[id] -= scem;
}
}
}
// 处理重叠谱线
if freq_data.nlines[ij_idx] > 0 {
let nlines = freq_data.nlines[ij_idx] as usize;
for ilint in 0..nlines {
let itr = (freq_data.itrlin[ilint * 100000 + ij_idx] - 1) as usize;
let itrprd = atomic.iprd[itr];
if itrprd == 0 {
continue;
}
// 找到频率范围
let mut ij0 = atomic.ifr0[itr] as usize;
let ifr1 = atomic.ifr1[itr] as usize;
for ijt in ij0..=ifr1 {
ij0 = ijt;
if freq_data.freq[ijt - 1] <= fr {
break;
}
}
let ij1 = ij0 - 1;
let a1 = (fr - freq_data.freq[ij0 - 1])
/ (freq_data.freq[ij1] - freq_data.freq[ij0 - 1]);
let a2 = UN - a1;
let dfr = (fr - atomic.fr0[itr]).abs();
// 检查是否是 H- 的跃迁
if atomic.ilow[itr] as usize == atomic.nfirst_elh {
let omeg = dfr * PI2;
let gra = A21 + GR * model.popul[atomic.nfirst_elh * nd];
for id in 0..nd {
model.coher[(itrprd as usize - 1) * nd + id] =
A21 / (gra + gami(2, "elec", omeg, model.temp[id], model.elec[id]));
}
}
for id in 0..nd {
let sg = a1 * model.prflin[id * 100000 + ij1]
+ a2 * model.prflin[id * 100000 + ij0 - 1];
let sg_final =
if dfr / model.doptr[(itrprd as usize - 1) * nd + id] <= config.xpdiv {
0.0
} else {
sg
};
let scalin = sg_final
* model.abtra[itr * nd + id]
* model.coher[(itrprd as usize - 1) * nd + id];
let scem = sg_final
* model.emtra[itr * nd + id]
* model.coher[(itrprd as usize - 1) * nd + id]
* model.xkfb[id];
model.scat1[id] += scalin;
model.emis1[id] -= scem;
}
}
}
} else {
// ODF 采样选项
if freq_data.nlines[ij_idx] > 0 {
let nlines = freq_data.nlines[ij_idx] as usize;
for ilint in 0..nlines {
let itr = (freq_data.itrlin[ilint * 100000 + ij_idx] - 1) as usize;
let itrprd = atomic.iprd[itr];
if itrprd == 0 {
continue;
}
let kj = ij - atomic.ifr0[itr] as usize + atomic.kfr0[itr] as usize;
let indxpa = atomic.indexp[itr].abs();
if indxpa != 3 && indxpa != 4 {
let dfr = (fr - atomic.fr0[itr]).abs();
// 检查是否是 H- 的跃迁
if atomic.ilow[itr] as usize == atomic.nfirst_elh {
let omeg = dfr * PI2;
let gra = A21 + GR * model.popul[atomic.nfirst_elh * nd];
for id in 0..nd {
model.coher[(itrprd as usize - 1) * nd + id] =
A21 / (gra + gami(2, "elec", omeg, model.temp[id], model.elec[id]));
}
}
for id in 0..nd {
let sg = model.prflin[id * 100000 + kj - 1];
let sg_final =
if dfr / model.doptr[(itrprd as usize - 1) * nd + id] <= config.xpdiv
{
0.0
} else {
sg
};
let scalin = sg_final
* model.abtra[itr * nd + id]
* model.coher[(itrprd as usize - 1) * nd + id];
model.scat1[id] += scalin;
model.emis1[id] -= 0.0; // SCEM 在这个分支中未定义
}
}
}
}
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_prd_constants() {
assert!((A21 - 4.699e8).abs() < 1e3);
assert!((PI2 - 6.28318531).abs() < 1e-8);
assert!((GR - 9.6e-8).abs() < 1e-10);
}
#[test]
fn test_prd_zero_ij() {
// 当 ij = 0 时,应该直接返回
let params = PrdParams { ij: 0 };
let config = PrdConfig {
ispodf: 0,
nd: 10,
xpdiv: 10.0,
};
// 创建空数据
let ijlin = vec![0; 100];
let iprd = vec![0; 100];
let fr0 = vec![0.0; 100];
let ilow = vec![0; 100];
let ifr0 = vec![0; 100];
let ifr1 = vec![0; 100];
let kfr0 = vec![0; 100];
let indexp = vec![0; 100];
let atomic = PrdAtomicData {
ijlin: &ijlin,
iprd: &iprd,
fr0: &fr0,
ilow: &ilow,
ifr0: &ifr0,
ifr1: &ifr1,
kfr0: &kfr0,
indexp: &indexp,
nfirst_elh: 0,
};
let temp = vec![10000.0; 10];
let elec = vec![1e12; 10];
let abtra = vec![1e-10; 1000];
let emtra = vec![1e-10; 1000];
let prflin = vec![1.0; 1000000];
let doptr = vec![1.0; 100];
let mut coher = vec![1.0; 100];
let popul = vec![1e10; 1000];
let xkfb = vec![1.0; 10];
let mut scat1 = vec![0.0; 10];
let mut emis1 = vec![0.0; 10];
let mut model = PrdModelState {
temp: &temp,
elec: &elec,
abtra: &abtra,
emtra: &emtra,
prflin: &prflin,
doptr: &doptr,
coher: &mut coher,
popul: &popul,
xkfb: &xkfb,
scat1: &mut scat1,
emis1: &mut emis1,
};
let freq = vec![1e15; 100];
let ijlin = vec![0; 100];
let nlines = vec![0; 100];
let itrlin = vec![0; 10000];
let freq_data = PrdFreqData {
freq: &freq,
ijlin: &ijlin,
nlines: &nlines,
itrlin: &itrlin,
};
prd(&params, &config, &atomic, &mut model, &freq_data);
// ij = 0 时,scat1 和 emis1 应该保持不变
for id in 0..10 {
assert!((model.scat1[id] - 0.0).abs() < 1e-15);
assert!((model.emis1[id] - 0.0).abs() < 1e-15);
}
}
}
+15 -14
View File
@@ -68,7 +68,7 @@ pub fn ratmat(
let ntrans = config.basnum.ntrans as usize;
let natom = config.basnum.natom as usize;
let lte = config.inppar.lte;
let ipslte = config.inppar.ipslte;
let ipslte = config.basnum.ipslte;
// 如果使用 Slater 迭代,清零辐射速率
if ipslte != 0 {
@@ -105,8 +105,8 @@ pub fn ratmat(
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;
llte[i] = llt || lte || atomic.ionpar.ilte[i] >= 1 || ilt >= 2;
llte[i] = llte[i] || id >= config.basnum.idlte as usize;
for j in 0..nlevel {
a[j][i] = 0.0;
@@ -141,20 +141,20 @@ pub fn ratmat(
// 计算跃迁速率
for itr in 0..ntrans {
let i = atomic.trapar.ilow[itr] as usize;
if atomic.atopar.iifix[atomic.trapar.iatm[i] as usize] == 1 {
if atomic.atopar.iifix[atomic.levpar.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])
aij[itr] = (model.crates.colrat[itr][id_idx] + model.rrrates.rru[itr][id_idx])
* model.wmcomp.wop[j][id_idx];
// 向下总速率
if atomic.trapar.line[itr] {
if atomic.trapar.line[itr] != 0 {
// 束缚-束缚跃迁
aji[itr] = (model.rrrates.coltar[itr][id_idx]
aji[itr] = (model.crates.coltar[itr][id_idx]
+ model.rrrates.rrd[itr][id_idx]
* atomic.levpar.g[i] / atomic.levpar.g[j]
* (hkt * atomic.trapar.fr0[itr]).exp())
@@ -167,7 +167,7 @@ pub fn ratmat(
} else {
UN
};
aji[itr] = model.rrrates.coltar[itr][id_idx] * model.wmcomp.wop[i][id_idx]
aji[itr] = model.crates.coltar[itr][id_idx] * model.wmcomp.wop[i][id_idx]
+ model.rrrates.rrd[itr][id_idx] * sbw[i] * corr;
}
@@ -181,10 +181,10 @@ pub fn ratmat(
// 填充速率矩阵
for itr in 0..ntrans {
let i = atomic.trapar.ilow[itr] as usize;
if atomic.atopar.iifix[atomic.trapar.iatm[i] as usize] == 1 {
if atomic.atopar.iifix[atomic.levpar.iatm[i] as usize] == 1 {
continue;
}
let nrefi = atomic.atopar.nrefs[atomic.trapar.iatm[i] as usize][id_idx] as usize;
let nrefi = atomic.atopar.nrefs[atomic.levpar.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;
@@ -273,8 +273,8 @@ pub fn ratmat(
}
b[nrefii] += model.modpar.dens[id_idx]
/ model.modpar.wmm[id_idx]
/ model.modpar.ytot[id_idx]
/ config.inppar.wmm[id_idx]
/ config.inppar.ytot[id_idx]
* atomic.atopar.abund[iat][id_idx];
}
@@ -284,6 +284,7 @@ pub fn ratmat(
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
#[test]
fn test_ratmat_ioptab_negative() {
@@ -349,8 +350,8 @@ mod tests {
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;
config.inppar.wmm[0] = 1.0;
config.inppar.ytot[0] = 1.0;
atomic.atopar.abund[0][0] = 1.0;
let mut iical = vec![1, 2, 3, 4, 5];
+200
View File
@@ -0,0 +1,200 @@
//! 辐射转移方程的角度积分点初始化。
//!
//! 重构自 TLUSTY `RTEANG.f`。
use super::gauleg;
/// RTEANG 的输入参数。
pub struct RteangParams {
/// 外部辐照的角度半宽(以弧度为单位)
/// 如果 WANGLE <= 0,则使用标准高斯积分
pub wangle: f64,
/// 外部辐照强度数组(长度 = nfreq)
pub extin: Vec<f64>,
}
/// RTEANG 的输出状态。
pub struct RteangOutput {
/// 角度点(cos(theta)
pub amu: Vec<f64>,
/// 角度权重
pub wtmu: Vec<f64>,
/// 辐照角度权重
pub fmu: Vec<f64>,
/// 角度点数
pub nmu: usize,
/// 表面 J 积分的贡献
pub extj: Vec<f64>,
/// 表面 H 积分的贡献
pub exth: Vec<f64>,
}
/// 初始化辐射转移方程的角度积分点。
///
/// # 参数
/// - `params`: 输入参数(wangle, extin
/// - `nmu_standard`: 标准角度点数(用于无辐照情况)
///
/// # 返回
/// - `RteangOutput`: 包含角度点和权重
///
/// # 算法说明
/// - 如果 `wangle <= 0`:使用标准 NMU 点高斯积分
/// - 如果 `wangle > 0`:使用 5 点积分方案,考虑外部辐照
pub fn rteang(params: &RteangParams, nmu_standard: usize) -> RteangOutput {
let nfreq = params.extin.len();
let half = 0.5_f64;
let one = 1.0_f64;
let pi = std::f64::consts::PI;
// 5 点积分的常量
const NMU5: usize = 5;
const NMU3: usize = 3;
// 1/sqrt(3)
const INV_SQRT3: f64 = 0.577350269189626;
let x = params.wangle * half;
// 初始化输出数组(最大可能的尺寸)
let max_nmu = NMU5.max(nmu_standard);
let mut amu = vec![0.0; max_nmu];
let mut wtmu = vec![0.0; max_nmu];
let mut fmu = vec![0.0; max_nmu];
let mut nmu: usize;
let mut xj = 0.0;
let mut xh = 0.0;
if x <= 0.0 {
// 无外部辐照:使用标准高斯积分
let (amu0, wtmu0) = gauleg(0.0, one, nmu_standard);
nmu = nmu_standard;
for i in 0..nmu {
amu[i] = amu0[i];
wtmu[i] = wtmu0[i];
fmu[i] = 0.0;
}
} else {
// 有外部辐照:使用特殊的 5 点积分方案
let x0 = half - x;
let x1 = half + x;
// 3 点高斯积分在 [-1, 1] 上
let (amu0, wtmu0) = gauleg(-one, one, NMU3);
for i in 0..NMU3 {
amu[i] = x0 * amu0[i] + x1;
wtmu[i] = x0 * wtmu0[i];
fmu[i] = 0.0;
}
nmu = NMU5;
// 第 4 和第 5 个角度点
let i4 = NMU3; // 0-indexed
let i5 = NMU3 + 1;
amu[i4] = x * (one + INV_SQRT3);
amu[i5] = x * (one - INV_SQRT3);
for i in NMU3..NMU5 {
wtmu[i] = x;
// 计算辐照角度权重
let amu_sq = amu[i] * amu[i];
let arg = (params.wangle * params.wangle - amu_sq) / (one - amu_sq);
if arg > 0.0 {
fmu[i] = (arg.sqrt().asin()) / pi;
} else {
fmu[i] = 0.0;
}
xj += wtmu[i] * fmu[i];
xh += wtmu[i] * amu[i] * fmu[i];
}
}
// 计算表面积分贡献
let mut extj = vec![0.0; nfreq];
let mut exth = vec![0.0; nfreq];
for ij in 0..nfreq {
extj[ij] = xj * params.extin[ij] * half;
exth[ij] = xh * params.extin[ij] * half;
}
// 截断数组到实际大小
amu.truncate(nmu);
wtmu.truncate(nmu);
fmu.truncate(nmu);
RteangOutput {
amu,
wtmu,
fmu,
nmu,
extj,
exth,
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_rteang_no_irradiation() {
// 无外部辐照的情况
let params = RteangParams {
wangle: 0.0,
extin: vec![1.0; 100],
};
let result = rteang(&params, 3);
// 应该使用 3 点高斯积分
assert_eq!(result.nmu, 3);
// 检查权重和(应该接近 1
let wsum: f64 = result.wtmu.iter().sum();
assert!((wsum - 1.0).abs() < 1e-10);
// 检查 fmu 全为零
for &f in &result.fmu {
assert_eq!(f, 0.0);
}
// 检查 extj 和 exth 为零(因为 xj = xh = 0
for &e in &result.extj {
assert_eq!(e, 0.0);
}
for &e in &result.exth {
assert_eq!(e, 0.0);
}
}
#[test]
fn test_rteang_with_irradiation() {
// 有外部辐照的情况
let params = RteangParams {
wangle: 0.1, // 小角度
extin: vec![1.0; 100],
};
let result = rteang(&params, 3);
// 应该使用 5 点积分
assert_eq!(result.nmu, 5);
// 检查 extj 和 exth 非零
assert!(result.extj[0] > 0.0);
assert!(result.exth[0] > 0.0);
}
#[test]
fn test_gauss_legendre_weights() {
// 验证高斯积分权重的归一性
let (x, w) = gauleg(0.0, 1.0, 5);
let sum: f64 = w.iter().sum();
assert!((sum - 1.0).abs() < 1e-14);
// 验证积分 x^2 在 [0,1] 上等于 1/3
let integral: f64 = x.iter().zip(w.iter()).map(|(xi, wi)| xi * xi * wi).sum();
assert!((integral - 1.0 / 3.0).abs() < 1e-14);
}
}
+130
View File
@@ -0,0 +1,130 @@
//! 表面质量密度计算(假设电子散射主导的不透明度)。
//!
//! 重构自 TLUSTY `sigmar.f`
//!
//! 模型假设:1-zone 模型(rho, T_g, mu 随高度常数),
//! t_r,phi = -alpha P,耗散单位光学深度常数。
//!
//! 参考: Krolik, Chapter 7
use super::laguer::laguer;
use num_complex::Complex64;
/// 计算表面质量密度。
///
/// 假设不透明度是电子散射主导的(或 kappa 与密度无关)。
///
/// # 参数
/// * `alpha` - 薄片参数 P (有效光学厚度)
/// * `xmdt` - M_dot * R (质量 × 旋转频率)
/// * `tef` - 有效温度 T_eff (K)
/// * `omega` - 角速度 Ω (rad/s)
/// * `relr` - 径向相对论因子
/// * `relt` - 横向相对论因子
/// * `relz` - 垂直相对论因子
///
/// # 返回值
/// 表面质量密度 Σ (g/cm²)
///
/// # 算法说明
/// 求解 10 阶方程找 x^4 项的根,使用 Laguerre 方法。
pub fn sigmar(alpha: f64, xmdt: f64, tef: f64, omega: f64, relr: f64, relt: f64, relz: f64) -> f64 {
const ZERO: f64 = 0.0;
const ONE: f64 = 1.0;
const TRES: f64 = 3.0;
const FOUR: f64 = 4.0;
const HALF: f64 = 0.5;
const FOURTH: f64 = 0.25;
const EPS: f64 = 1e-5;
// 物理常数
const C: f64 = 2.9979e10; // 光速 (cm/s)
const SIGMAB: f64 = 5.6703e-5; // 汤姆逊散射截面 (cm²)
const BK: f64 = 1.3807e-16; // 玻尔兹曼常数 (erg/K)
const PI: f64 = std::f64::consts::PI;
// 假设完全电离的纯氢
let kappa: f64 = 0.39;
let mu: f64 = 0.5 * 1.6726e-24; // 平均分子质量 (g)
let fac1 = relz * (HALF * TRES * C * omega / alpha / kappa / SIGMAB / tef.powi(4)).powi(2);
let fac2 = (HALF * kappa).powf(FOURTH) * BK * tef / mu;
let fac3 = xmdt * omega * relt / PI;
// 构造 11 次方程的系数 (索引 0-10)
let mut coeff: [Complex64; 11] = [Complex64::new(0.0, 0.0); 11];
coeff[0] = Complex64::new(fac1 * (HALF * fac3).powi(2), ZERO);
coeff[1] = Complex64::new(ZERO, ZERO);
coeff[2] = Complex64::new(ZERO, ZERO);
coeff[3] = Complex64::new(ZERO, ZERO);
coeff[4] = Complex64::new(-(TRES * fac3) / (8.0 * alpha), ZERO);
coeff[5] = Complex64::new(-fac1 * fac3 * alpha * fac2, ZERO);
coeff[6] = Complex64::new(ZERO, ZERO);
coeff[7] = Complex64::new(ZERO, ZERO);
coeff[8] = Complex64::new(ZERO, ZERO);
coeff[9] = Complex64::new(FOURTH * fac2, ZERO);
coeff[10] = Complex64::new(fac1 * (alpha * fac2).powi(2), ZERO);
// 计算辐射压主导时的 sigma
let sigrad = FOUR * omega * C * C * relt * relz / alpha / kappa.powi(2) / SIGMAB / tef.powi(4) / relr;
// 计算气压主导时的 sigma
let siggas = ((mu * xmdt * omega * relt / PI / alpha / BK / tef).powi(4) / 8.0 / kappa).powf(0.2);
// 初始猜测
let mut x_guess = Complex64::new(ONE / (ONE / sigrad + ONE / siggas).powf(FOURTH), ZERO);
// 使用 Laguerre 方法找根
laguer(&coeff, &mut x_guess);
// 检查结果有效性
if x_guess.im.abs() < EPS && x_guess.re > ZERO {
x_guess.re.powi(4)
} else {
ONE / (ONE / sigrad + ONE / siggas)
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_sigmar_basic() {
// 基本测试:确保函数能运行
let result = sigmar(
1.0, // alpha
1e-10, // xmdt
10000.0, // tef
1e-3, // omega
1.0, // relr
1.0, // relt
1.0, // relz
);
assert!(result > 0.0, "sigmar should return positive value, got {}", result);
}
#[test]
fn test_sigmar_typical_disk() {
// 典型吸积盘参数测试
let alpha = 0.1;
let xmdt = 1e-8;
let tef = 5000.0;
let omega = 1e-4;
let relr = 1.0;
let relt = 1.0;
let relz = 1.0;
let result = sigmar(alpha, xmdt, tef, omega, relr, relt, relz);
assert!(result > 0.0, "sigmar should return positive value, got {}", result);
assert!(result.is_finite());
}
#[test]
fn test_sigmar_high_temperature() {
// 高温测试
let result = sigmar(1.0, 1e-10, 50000.0, 1e-3, 1.0, 1.0, 1.0);
assert!(result > 0.0, "sigmar should return positive value, got {}", result);
}
}
+221
View File
@@ -0,0 +1,221 @@
//! 灰模型的局部温度计算。
//!
//! 重构自 TLUSTY `tlocal.f`
//!
//! 根据光学深度计算灰模型的局部温度。
use super::quartc::quartc;
use crate::state::constants::{MDEPTH, UN};
/// TLOCAL 输入参数
pub struct TlocalParams {
/// 深度索引 (1-indexed)
pub id: usize,
/// 通量平均不透明度的当前估计
pub tauf: f64,
}
/// TLOCAL 模型状态(来自 MODELQ.FOR
pub struct TlocalModelState<'a> {
/// 温度数组 (nd)
pub temp: &'a [f64],
/// 粘性系数数组 (nd)
pub viscd: &'a [f64],
/// 热深度数组 (nd)
pub tauthe: &'a [f64],
/// 普朗克平均不透明度数组 (nd)
pub abplad: &'a [f64],
/// Rosseland 平均不透明度数组 (nd)
pub abrosd: &'a [f64],
}
/// TLOCAL 配置参数(来自 BASICS.FOR
pub struct TlocalConfig {
/// 深度点数
pub nd: usize,
/// 有效温度
pub teff: f64,
/// 恒星温度
pub tstar: f64,
/// 稀释因子
pub wdil: f64,
/// 分数
pub fractv: f64,
/// Compton 标志
pub icompt: i32,
/// Compton 灰色标志
pub icomgr: i32,
/// 盘温度(如果 > 0 则使用恒定温度)
pub tdisk: f64,
/// 质量列向量 (nd)
pub dm: Vec<f64>,
}
/// TLOCAL 辅助通量变量(来自 COMMON/FLXAUX/
pub struct TlocalFlxaux {
/// T^4 因子
pub t4: f64,
}
/// TLOCAL 因子变量(来自 COMMON/FACTRS/
pub struct TlocalFactrs<'a> {
/// GAMJ 数组
pub gamj: &'a [f64],
/// GAMH
pub gamh: f64,
/// FAK0
pub fak0: f64,
}
/// 计算灰模型的局部温度。
///
/// # 参数
///
/// * `params` - 输入参数(id, tauf
/// * `config` - 配置参数
/// * `model` - 模型状态
/// * `flxaux` - 辅助通量变量
/// * `factrs` - 因子变量
///
/// # 返回值
///
/// 局部温度 T (K)
pub fn tlocal(
params: &TlocalParams,
config: &TlocalConfig,
model: &TlocalModelState,
flxaux: &TlocalFlxaux,
factrs: &TlocalFactrs,
) -> f64 {
let id = params.id;
let id_idx = id - 1;
let nd = config.nd;
let nd_idx = nd - 1;
// 如果是盘模型且 tdisk > 0,使用恒定温度
if config.tdisk > 0.0 {
return config.tdisk;
}
// 计算粘性项
let vis = model.viscd[id_idx] / (3.0 * config.dm[nd_idx]);
// 计算额外辐照项
let extra = 4.0 * factrs.fak0 * config.wdil * (config.tstar / config.teff).powi(4);
// 获取 GAMJ 和 GAMH
let gj = factrs.gamj[id_idx];
let gh = factrs.gamh * 0.57735; // 5.7735e-1
// 计算 gg
let gg = (params.tauf - model.tauthe[id_idx]) * config.fractv + gh + extra;
// 简化情况:icompt == 0 或 icomgr == 0
if config.icompt == 0 || config.icomgr == 0 {
let t = (0.75 * flxaux.t4 * (gj * gg + vis / model.abplad[id_idx])).powf(0.25);
return t;
}
// 完整计算
// 常量
const C1: f64 = 0.8112;
const C3: f64 = 6.745e-10;
const C4: f64 = 0.96;
const C34: f64 = C3 * C4;
let epsbar = model.abplad[id_idx] / model.abrosd[id_idx];
let mut tfor = C1 * config.teff * epsbar.powf(-0.125);
let tf0 = tfor;
let b: f64;
if (params.tauf > UN && tfor < model.temp[id_idx]) || params.tauf >= 100.0 {
tfor = 0.0;
// 注意:Fortran 中先计算 b = gg*(c3-c34),然后 b = 0.
// 最终 b = 0
b = 0.0;
} else {
b = gg * C3;
}
let a = epsbar / (0.75 * flxaux.t4);
let c = gg * (epsbar * gj + C34 * tfor) + vis / model.abrosd[id_idx];
// 解四次方程
let t1 = quartc(a, b, c);
t1
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_tlocal_tdisk_positive() {
// 当 tdisk > 0 时,应返回 tdisk
let params = TlocalParams { id: 1, tauf: 1.0 };
let config = TlocalConfig {
nd: 10,
teff: 10000.0,
tstar: 10000.0,
wdil: 0.5,
fractv: 1.0,
icompt: 0,
icomgr: 0,
tdisk: 5000.0,
dm: vec![1.0; 10],
};
let temp = vec![10000.0; 10];
let model = TlocalModelState {
temp: &temp,
viscd: &[0.0; 10],
tauthe: &[0.0; 10],
abplad: &[1.0; 10],
abrosd: &[1.0; 10],
};
let flxaux = TlocalFlxaux { t4: 1.0 };
let gamj = vec![1.0; MDEPTH];
let factrs = TlocalFactrs {
gamj: &gamj,
gamh: 1.0,
fak0: 1.0,
};
let result = tlocal(&params, &config, &model, &flxaux, &factrs);
assert!((result - 5000.0).abs() < 1e-10);
}
#[test]
fn test_tlocal_simple_case() {
// icompt = 0 的简化情况
let params = TlocalParams { id: 1, tauf: 1.0 };
let config = TlocalConfig {
nd: 10,
teff: 10000.0,
tstar: 10000.0,
wdil: 0.5,
fractv: 1.0,
icompt: 0,
icomgr: 0,
tdisk: 0.0,
dm: vec![1.0; 10],
};
let temp = vec![10000.0; 10];
let model = TlocalModelState {
temp: &temp,
viscd: &[0.0; 10],
tauthe: &[0.0; 10],
abplad: &[1.0; 10],
abrosd: &[1.0; 10],
};
let flxaux = TlocalFlxaux { t4: 1.0 };
let gamj = vec![1.0; MDEPTH];
let factrs = TlocalFactrs {
gamj: &gamj,
gamh: 1.0,
fak0: 1.0,
};
let result = tlocal(&params, &config, &model, &flxaux, &factrs);
assert!(result > 0.0);
}
}