From 0674b4f174b4f57ae73ad5c5641a3798ca6a697d Mon Sep 17 00:00:00 2001 From: Asfmq <2696428814@qq.com> Date: Wed, 25 Mar 2026 06:52:44 +0800 Subject: [PATCH] =?UTF-8?q?feat:=20=E6=B7=BB=E5=8A=A0=209=20=E4=B8=AA=20SY?= =?UTF-8?q?NSPEC=20=E6=95=B0=E5=AD=A6=E6=A8=A1=E5=9D=97=20(=E7=AC=AC8?= =?UTF-8?q?=E6=89=B9)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 新增模块: - count_words: 字符串单词计数工具 - divhe2: He II Stark 轮廓除数参数计算 - extprf: 谱线轮廓波长外推 (Cooper 公式) - feautr: Lyman-α Stark 加宽 (Feautrier 方法) - gamhe: 中性氦 Stark 加宽参数 - griem: Griem Stark 阻尼参数计算 - intrp: 二分法高效插值程序 - partdv: 配分函数计算 (含压力效应) - sffhmi_old: H- 自由-自由截面 (Kurucz 公式) 改进: - 修复 fortran-analyzer 注释行误匹配问题 Co-Authored-By: Claude Opus 4.6 --- .../scripts/analyze_fortran.py | 46 ++- src/math/count_words.rs | 93 ++++++ src/math/divhe2.rs | 132 ++++++++ src/math/extprf.rs | 151 +++++++++ src/math/feautr.rs | 236 ++++++++++++++ src/math/gamhe.rs | 263 ++++++++++++++++ src/math/griem.rs | 177 +++++++++++ src/math/intrp.rs | 262 ++++++++++++++++ src/math/mod.rs | 18 ++ src/math/partdv.rs | 290 ++++++++++++++++++ src/math/sffhmi_old.rs | 136 ++++++++ 11 files changed, 1796 insertions(+), 8 deletions(-) create mode 100644 src/math/count_words.rs create mode 100644 src/math/divhe2.rs create mode 100644 src/math/extprf.rs create mode 100644 src/math/feautr.rs create mode 100644 src/math/gamhe.rs create mode 100644 src/math/griem.rs create mode 100644 src/math/intrp.rs create mode 100644 src/math/partdv.rs create mode 100644 src/math/sffhmi_old.rs diff --git a/.claude/skills/fortran-analyzer/scripts/analyze_fortran.py b/.claude/skills/fortran-analyzer/scripts/analyze_fortran.py index 06a38f9..4a1e7d6 100644 --- a/.claude/skills/fortran-analyzer/scripts/analyze_fortran.py +++ b/.claude/skills/fortran-analyzer/scripts/analyze_fortran.py @@ -53,6 +53,26 @@ FORTRAN_INTRINSICS = { 'READ', 'WRITE', 'OPEN', 'CLOSE', 'FORMAT', 'PRINT', } +def strip_fortran_comments(content): + """移除 Fortran 注释行（固定格式） + + Fortran 固定格式中，以下开头的行是注释： + - 'c' 或 'C' 在第 1 列 + - '!' 在第 1 列（自由格式也支持） + - '*' 在第 1 列 + - 空行 + """ + lines = content.split('\n') + code_lines = [] + for line in lines: + if len(line) == 0: + continue + first_char = line[0].upper() + if first_char in ('C', '!', '*'): + continue # 注释行 + code_lines.append(line) + return '\n'.join(code_lines) + def extract_calls(content, known_functions=None): """提取 CALL 语句和 FUNCTION 调用 @@ -62,9 +82,12 @@ def extract_calls(content, known_functions=None): """ calls = set() + # 先移除注释行，避免误匹配注释中的 CALL + code_content = strip_fortran_comments(content) + # 1. 提取 CALL 语句（支持有括号和无括号两种形式） # CALL NAME(...) 或 CALL NAME - call_stmts = re.findall(r'(?i)CALL\s+(\w+)(?:\s*\(|\s*$|\s*\n)', content) + call_stmts = re.findall(r'(?i)CALL\s+(\w+)(?:\s*\(|\s*$|\s*\n)', code_content) # 过滤掉 Fortran 关键字（IF, DO, THEN 等不是子程序名） calls.update(c.upper() for c in call_stmts if c.upper() not in FORTRAN_INTRINSICS) @@ -146,20 +169,27 @@ SPECIAL_MAPPINGS = { def find_rust_module(fortran_name, rust_math_dir, rust_io_dir): """查找对应的 Rust 模块""" + # Fortran 名称是大写，Rust 文件是小写 + rust_name = fortran_name.lower() + # 先检查 math 目录 - rust_file = os.path.join(rust_math_dir, f"{fortran_name}.rs") + rust_file = os.path.join(rust_math_dir, f"{rust_name}.rs") if os.path.exists(rust_file): - return f"src/math/{fortran_name}.rs" + return f"src/math/{rust_name}.rs" # 检查 io 目录 - rust_file = os.path.join(rust_io_dir, f"{fortran_name}.rs") + rust_file = os.path.join(rust_io_dir, f"{rust_name}.rs") if os.path.exists(rust_file): - return f"src/io/{fortran_name}.rs" + return f"src/io/{rust_name}.rs" - # 检查特殊映射 (math 目录) + # 检查特殊映射 (math 目录) - 必须验证文件实际存在 for rust_mod, fortran_funcs in SPECIAL_MAPPINGS.items(): - if fortran_name in fortran_funcs: - return f"src/math/{rust_mod}.rs" + if fortran_name.lower() in [f.lower() for f in fortran_funcs]: + mapped_file = os.path.join(rust_math_dir, f"{rust_mod}.rs") + if os.path.exists(mapped_file): + return f"src/math/{rust_mod}.rs" + # 如果映射的文件不存在，继续查找其他映射或返回空 + break return "" diff --git a/src/math/count_words.rs b/src/math/count_words.rs new file mode 100644 index 0000000..6614dbc --- /dev/null +++ b/src/math/count_words.rs @@ -0,0 +1,93 @@ +//! 字符串单词计数工具。 +//! +//! 重构自 SYNSPEC `count_words.f` + +/// 统计字符串中由空格分隔的单词数量。 +/// +/// # 参数 +/// +/// * `cadena` - 输入字符串 +/// +/// # 返回值 +/// +/// 字符串中的单词数量 +/// +/// # 算法 +/// +/// 遍历字符串，当遇到非空格字符且前一个字符是空格时，计数加一。 +/// 如果第一个字符不是空格，则计数从 1 开始。 +/// +/// # 示例 +/// +/// ``` +/// use spectrarust::math::count_words::count_words; +/// assert_eq!(count_words("hello world"), 2); +/// assert_eq!(count_words(" hello world "), 2); +/// assert_eq!(count_words(""), 0); +/// assert_eq!(count_words(" "), 0); +/// ``` +pub fn count_words(cadena: &str) -> i32 { + let chars: Vec = cadena.chars().collect(); + let len = chars.len(); + + if len == 0 { + return 0; + } + + let mut n: i32 = 0; + let mut a = chars[0]; + + // 如果第一个字符不是空格，计数从 1 开始 + if a != ' ' { + n = 1; + } + + // 遍历剩余字符 + for i in 1..len { + let b = chars[i]; + // 如果当前字符不是空格且前一个字符是空格，增加计数 + if b != ' ' && a == ' ' { + n += 1; + } + a = b; + } + + n +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_empty_string() { + assert_eq!(count_words(""), 0); + } + + #[test] + fn test_only_spaces() { + assert_eq!(count_words(" "), 0); + } + + #[test] + fn test_single_word() { + assert_eq!(count_words("hello"), 1); + assert_eq!(count_words(" hello"), 1); + assert_eq!(count_words("hello "), 1); + assert_eq!(count_words(" hello "), 1); + } + + #[test] + fn test_multiple_words() { + assert_eq!(count_words("hello world"), 2); + assert_eq!(count_words("hello world foo"), 3); + assert_eq!(count_words(" hello world foo "), 3); + } + + #[test] + fn test_fortran_style() { + // Fortran character*1000 测试 + let long_str = "word1 word2 word3"; + assert_eq!(count_words(long_str), 3); + } +} diff --git a/src/math/divhe2.rs b/src/math/divhe2.rs new file mode 100644 index 0000000..f757ee9 --- /dev/null +++ b/src/math/divhe2.rs @@ -0,0 +1,132 @@ +//! He II Stark 轮廓近似计算辅助函数。 +//! +//! 重构自 SYNSPEC `divhe2.f` +//! +//! 计算 He II 线 Stark 轮廓的除数参数。 +//! 此函数与氢的 DIVSTR 类似，唯一的区别是参数 A 的定义略有不同： +//! He II 的 A 等于氢的 A 减去 ln(2)。 + +/// 计算 He II Stark 轮廓的除数参数。 +/// +/// # 参数 +/// +/// * `betad` - 约化电子密度参数 +/// +/// # 返回值 +/// +/// 返回除数参数 DIV。如果 `betad < 5.821`，返回初始计算的 A 值。 +/// +/// # 算法 +/// +/// 1. 计算 A = 1.5 * ln(betad) - 0.978 +/// 2. 如果 betad < 5.821，直接返回 A +/// 3. 否则，使用牛顿-拉夫逊迭代求解方程 X² - 2.5*ln(X) = A +/// +/// # 示例 +/// +/// ``` +/// use spectrarust::math::divhe2::divhe2; +/// // betad < 5.821 时返回 A +/// let result = divhe2(5.0); +/// assert!(result > 0.0); +/// +/// // betad >= 5.821 时进行迭代 +/// let result = divhe2(10.0); +/// assert!(result > 0.0); +/// ``` +pub fn divhe2(betad: f64) -> f64 { + const UN: f64 = 1.0; + const TWO: f64 = 2.0; + const UNQ: f64 = 1.25; + const UNH: f64 = 1.5; + const TWH: f64 = 2.5; + const FO: f64 = 4.0; + const FI: f64 = 5.0; + const CA: f64 = 0.978; + const BL: f64 = 5.821; + const AL: f64 = 1.26; + const CX: f64 = 0.28; + const DX: f64 = 0.0001; + + // 计算 A = 1.5 * ln(betad) - 0.978 + let a = UNH * betad.ln() - CA; + + // 如果 betad < BL，直接返回 A + if betad < BL { + return a; + } + + // 初始猜测值 X + let x = if a >= AL { + // 对于大的 A，使用渐近公式 + (a.sqrt()) * (UN + UNQ * a.ln() / (FO * a - FI)) + } else { + // 对于小的 A，使用修正公式 + (CX + a).sqrt() + }; + + // 牛顿-拉夫逊迭代求解 X² - 2.5*ln(X) = A + let mut x = x; + for _ in 0..5 { + // f(X) = X² - 2.5*ln(X) - A + // f'(X) = 2*X - 2.5/X + // 牛顿迭代: X_new = X - f(X)/f'(X) + // = X * (1 - (X² - 2.5*ln(X) - A) / (2*X² - 2.5)) + let xn = x * (UN - (x * x - TWH * x.ln() - a) / (TWO * x * x - TWH)); + if (xn - x).abs() <= DX { + return xn; + } + x = xn; + } + + x +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_betad_below_bl() { + // betad < 5.821 时，返回 A = 1.5 * ln(betad) - 0.978 + let result = divhe2(5.0); + let expected = 1.5 * 5.0_f64.ln() - 0.978; + assert_relative_eq!(result, expected, epsilon = 1e-10); + } + + #[test] + fn test_betad_above_bl() { + // betad >= 5.821 时，进行迭代 + let result = divhe2(10.0); + // 验证结果是正数且合理 + assert!(result > 0.0); + assert!(result < 10.0); // 合理范围 + } + + #[test] + fn test_betad_at_bl() { + // betad = BL 边界情况 + let result = divhe2(5.821); + assert!(result > 0.0); + } + + #[test] + fn test_large_betad() { + // 大的 betad 值 + let result = divhe2(100.0); + assert!(result > 0.0); + // 验证近似满足方程 X² - 2.5*ln(X) ≈ A + let a = 1.5 * 100.0_f64.ln() - 0.978; + let residual = result * result - 2.5 * result.ln() - a; + assert_relative_eq!(residual, 0.0, epsilon = 0.001); + } + + #[test] + fn test_small_betad() { + // 小的 betad 值 + let result = divhe2(1.0); + let expected = 1.5 * 1.0_f64.ln() - 0.978; + assert_relative_eq!(result, expected, epsilon = 1e-10); + } +} diff --git a/src/math/extprf.rs b/src/math/extprf.rs new file mode 100644 index 0000000..29e183b --- /dev/null +++ b/src/math/extprf.rs @@ -0,0 +1,151 @@ +//! 谱线轮廓波长外推函数。 +//! +//! 重构自 SYNSPEC `extprf.f` +//! +//! 在 Shamey 或 Barnard, Cooper, Smith 表中进行波长外推。 +//! 使用 Cooper 建议的特殊公式。 + +/// Cooper 外推公式计算谱线轮廓。 +/// +/// # 参数 +/// +/// * `dlam` - 波长偏移 +/// * `it` - 温度索引 (1-4) +/// * `iline` - 谱线索引 (1-4) +/// * `anel` - 电子密度对数参数 +/// * `dlast` - 最后一个波长点 +/// * `plast` - 最后一个轮廓值 +/// +/// # 返回值 +/// +/// 外推的谱线轮廓值 +/// +/// # Panics +/// +/// 如果 `it` 或 `iline` 不在 1-4 范围内会 panic。 +/// +/// # 算法 +/// +/// 使用 Cooper 公式： +/// ```text +/// WE = W0(it,iline) * exp(anel * ln(10)) * 1e-16 +/// F = |dlast|^2.5 * (plast - WE / (pi * dlast^2)) +/// EXTPRF = (WE / pi + F / sqrt(|dlam|)) / dlam^2 +/// ``` +pub fn extprf(dlam: f64, it: usize, iline: usize, anel: f64, dlast: f64, plast: f64) -> f64 { + // W0 数组：温度索引 (行) x 谱线索引 (列) + // Fortran DATA 语句是列优先存储 + const W0: [[f64; 4]; 4] = [ + [1.460, 6.130, 4.040, 2.312], // it=1 + [1.269, 5.150, 3.490, 1.963], // it=2 + [1.079, 4.240, 2.960, 1.624], // it=3 + [0.898, 3.450, 2.470, 1.315], // it=4 + ]; + + let it_idx = it - 1; // 转换为 0 索引 + let iline_idx = iline - 1; + + // 获取 W0 值 + let w0_val = W0[it_idx][iline_idx]; + + // WE = W0 * 10^anel * 1e-16 + // Fortran: EXP(ANEL*2.3025851) = 10^ANEL (因为 ln(10) ≈ 2.3025851) + let we = w0_val * (anel * 2.3025851_f64).exp() * 1e-16; + + // 使用 PI 的精确值 + const PI: f64 = std::f64::consts::PI; + + let dlasta = dlast.abs(); + // D52 = |dlast|^2.5 = |dlast|^2 * sqrt(|dlast|) + let d52 = dlasta * dlasta * dlasta.sqrt(); + + // F = D52 * (plast - WE / (pi * dlast^2)) + let f = d52 * (plast - we / (PI * dlast * dlast)); + + // EXTPRF = (WE / pi + F / sqrt(|dlam|)) / dlam^2 + (we / PI + f / dlam.abs().sqrt()) / (dlam * dlam) +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_basic() { + // 基本测试 + let result = extprf(0.1, 1, 1, 0.0, 0.5, 1.0); + assert!(result.is_finite()); + assert!(result > 0.0); + } + + #[test] + fn test_different_it_with_large_anel() { + // 使用大的 anel 值使 W0 差异显现 + // WE = W0 * 10^anel * 1e-16，当 anel 大时 W0 差异才明显 + let anel = 10.0; // 10^10 = 1e10 倍放大 + let base = extprf(0.1, 1, 1, anel, 0.5, 1.0); + let t2 = extprf(0.1, 2, 1, anel, 0.5, 1.0); + let t3 = extprf(0.1, 3, 1, anel, 0.5, 1.0); + let t4 = extprf(0.1, 4, 1, anel, 0.5, 1.0); + + // 不同温度索引应该产生不同结果 + assert!((base - t2).abs() > 1e-10, "IT=1 vs IT=2 should differ"); + assert!((t2 - t3).abs() > 1e-10, "IT=2 vs IT=3 should differ"); + assert!((t3 - t4).abs() > 1e-10, "IT=3 vs IT=4 should differ"); + } + + #[test] + fn test_different_iline_with_large_anel() { + // 使用大的 anel 值使 W0 差异显现 + let anel = 10.0; + let l1 = extprf(0.1, 1, 1, anel, 0.5, 1.0); + let l2 = extprf(0.1, 1, 2, anel, 0.5, 1.0); + let l3 = extprf(0.1, 1, 3, anel, 0.5, 1.0); + let l4 = extprf(0.1, 1, 4, anel, 0.5, 1.0); + + // 不同谱线索引应该产生不同结果 + assert!((l1 - l2).abs() > 1e-10, "ILINE=1 vs ILINE=2 should differ"); + assert!((l2 - l3).abs() > 1e-10, "ILINE=2 vs ILINE=3 should differ"); + assert!((l3 - l4).abs() > 1e-10, "ILINE=3 vs ILINE=4 should differ"); + } + + #[test] + fn test_with_anel() { + // 测试带电子密度参数 + // 注意：anel 增加时 WE 增加但 F 项可能减小，所以结果不一定增加 + // 这里测试 anel 确实影响结果 + let zero = extprf(0.1, 1, 1, 0.0, 0.5, 1.0); + let pos = extprf(0.1, 1, 1, 15.0, 0.5, 1.0); + + // anel 变化应该改变结果（不验证方向，只验证变化） + assert!((pos - zero).abs() > 1e-6, "anel=15 should change the result"); + } + + #[test] + fn test_symmetry() { + // 测试正负 dlast 的对称性 + let pos = extprf(0.1, 1, 1, 0.0, 0.5, 1.0); + let neg = extprf(0.1, 1, 1, 0.0, -0.5, 1.0); + + // 由于 dlast 使用 abs()，结果应该相同 + assert_relative_eq!(pos, neg, epsilon = 1e-14); + } + + #[test] + fn test_w0_values() { + // 验证 W0 数组值正确 + // 直接测试：当 anel 很大且 plast=0 时，结果主要由 WE 决定 + let anel = 20.0; + let plast = 0.0; // 消除 F 项 + + // W0(1,1) = 1.460, W0(2,1) = 1.269, 比例 = 1.460/1.269 + let r1 = extprf(0.1, 1, 1, anel, 0.5, plast); + let r2 = extprf(0.1, 2, 1, anel, 0.5, plast); + + // 比例应该接近 W0 比例 + let ratio = r1 / r2; + let expected_ratio = 1.460 / 1.269; + assert_relative_eq!(ratio, expected_ratio, epsilon = 1e-6); + } +} diff --git a/src/math/feautr.rs b/src/math/feautr.rs new file mode 100644 index 0000000..e53a79b --- /dev/null +++ b/src/math/feautr.rs @@ -0,0 +1,236 @@ +//! Lyman-α Stark 加宽计算 (Feautrier 方法)。 +//! +//! 重构自 SYNSPEC `feautr.f` +//! +//! 基于 N. Feautrier 的 Lyman-α Stark 加宽表进行插值计算。 + +/// Feautrier Lyman-α Stark 轮廓参数。 +/// +/// 包含从模型大气层传递的插值参数。 +#[derive(Debug, Clone)] +pub struct FeautrParams { + /// 温度插值索引 (JT) + pub jt: usize, + /// 温度插值系数 0 + pub ti0: f64, + /// 温度插值系数 1 + pub ti1: f64, + /// 温度插值系数 2 + pub ti2: f64, +} + +/// 使用 Feautrier 方法计算 Lyman-α Stark 加宽轮廓。 +/// +/// # 参数 +/// +/// * `freq` - 频率 (Hz) +/// * `params` - 模型层参数 (JT, TI0, TI1, TI2) +/// +/// # 返回值 +/// +/// 归一化的 Stark 轮廓值 (0 到 1 之间) +/// +/// # 算法 +/// +/// 1. 计算波长偏移: dlam = c/freq - 1215.685 Å +/// 2. 在预计算的 Feautrier 表中进行线性插值 +/// 3. 使用温度插值系数得到最终值 +/// +/// # 注意 +/// +/// 原始 Fortran 使用 COMMON/MODELP/ 中的 JT, TI0, TI1, TI2 数组。 +/// 此纯函数版本将这些作为参数传入。 +pub fn feautr(freq: f64, params: &FeautrParams) -> f64 { + // Feautrier 表数据 + // DL: 波长偏移 (Å) + const DL: [f64; 20] = [ + -150.0, -120.0, -90.0, -60.0, -40.0, -20.0, -10.0, -8.0, -4.0, -2.0, + 2.0, 4.0, 8.0, 10.0, 20.0, 40.0, 60.0, 90.0, 120.0, 150.0, + ]; + + // F05: T=5000K 时的轮廓值 + const F05: [f64; 20] = [ + 0.0537, 0.0964, 0.1330, 0.3105, 0.4585, 0.6772, 0.8229, + 0.8556, 0.9250, 0.9618, 0.9733, 1.1076, 1.0644, 1.0525, + 0.8841, 0.8282, 0.7541, 0.7091, 0.7164, 0.7672, + ]; + + // F10: T=10000K 时的轮廓值 + const F10: [f64; 20] = [ + 0.1986, 0.2764, 0.3959, 0.5740, 0.7385, 0.9448, 1.0292, + 1.0317, 0.9947, 0.8679, 0.8648, 0.9815, 1.0660, 1.0793, + 1.0699, 1.0357, 0.9245, 0.8603, 0.8195, 0.7928, + ]; + + // F20: T=20000K 时的轮廓值 + const F20: [f64; 20] = [ + 0.4843, 0.5821, 0.7003, 0.8411, 0.9405, 1.0300, 1.0029, + 0.9753, 0.8478, 0.6851, 0.6861, 0.8554, 0.9916, 1.0264, + 1.0592, 1.0817, 1.0575, 1.0152, 0.9761, 0.9451, + ]; + + // F40: T=40000K 时的轮廓值 + const F40: [f64; 20] = [ + 0.7862, 0.8566, 0.9290, 0.9915, 1.0066, 0.9878, 0.8983, + 0.8513, 0.6881, 0.5277, 0.5302, 0.6920, 0.8607, 0.9111, + 0.9651, 1.0793, 1.1108, 1.1156, 1.1003, 1.0839, + ]; + + // 光速 (cm/s) + const CL: f64 = 2.997925e18; + // Lyman-α 中心波长 (Å) + const LYA_WAVELENGTH: f64 = 1215.685; + + // 计算波长偏移 + let dlam = CL / freq - LYA_WAVELENGTH; + + // 查找 dlam 所在的区间 + // Fortran: DO 10 I=2,20; IF(DLAM.LE.DL(I)) GO TO 20 + let mut i = 19; // 默认最后一个区间 (0-indexed) + for idx in 1..20 { + if dlam <= DL[idx] { + i = idx; + break; + } + } + + let j = i - 1; // 区间左端点 + + // 线性插值系数 + let c = DL[j] - DL[i]; + let a = (dlam - DL[i]) / c; + let b = (DL[j] - dlam) / c; + + // 对四个温度表进行插值 + let x = [ + F05[j] * a + F05[i] * b, + F10[j] * a + F10[i] * b, + F20[j] * a + F20[i] * b, + F40[j] * a + F40[i] * b, + ]; + + // 使用模型参数进行温度插值 + // Fortran: J=JT(ID); Y=TI0(ID)*X(J)+TI1(ID)*X(J-1)+TI2(ID)*X(J-2) + // 注意: Fortran 索引从 1 开始，Rust 从 0 开始 + let jt_idx = params.jt; // 已是 0-indexed (调用者负责转换) + let y = params.ti0 * x[jt_idx] + + params.ti1 * x[jt_idx - 1] + + params.ti2 * x[jt_idx - 2]; + + // 归一化到 [0, 1] + 0.5 * (y + 1.0) +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_basic() { + // 使用中心波长附近的频率 + let freq = 2.997925e18 / 1215.685; // Lyman-α 中心频率 + // jt 必须至少为 2 (0-indexed) 才能访问 x[jt-2] + let params = FeautrParams { + jt: 2, // 使用 F20 表 (x[2]) + ti0: 1.0, + ti1: 0.0, + ti2: 0.0, + }; + + let result = feautr(freq, ¶ms); + assert!(result >= 0.0 && result <= 1.0); + } + + #[test] + fn test_wing() { + // 测试翼部 (远离线心) + let freq = 2.997925e18 / 1365.685; // 红翼 +150Å + let params = FeautrParams { + jt: 3, // 使用 F40 表 + ti0: 1.0, + ti1: 0.0, + ti2: 0.0, + }; + + let result = feautr(freq, ¶ms); + assert!(result >= 0.0); + } + + #[test] + fn test_blue_wing() { + // 测试蓝翼 + let freq = 2.997925e18 / 1065.685; // 蓝翼 -150Å + let params = FeautrParams { + jt: 3, // 使用 F40 表 + ti0: 1.0, + ti1: 0.0, + ti2: 0.0, + }; + + let result = feautr(freq, ¶ms); + assert!(result >= 0.0); + } + + #[test] + fn test_temperature_interpolation() { + let freq = 2.997925e18 / 1215.685; + + // 纯 F10 温度 (jt=1 需要 x[1], x[0], x[-1] - 不行!) + // 必须使用 jt >= 2 + let params_f20 = FeautrParams { + jt: 2, // F20 表 + ti0: 1.0, + ti1: 0.0, + ti2: 0.0, + }; + + // 纯 F40 温度 + let params_f40 = FeautrParams { + jt: 3, // F40 表 + ti0: 1.0, + ti1: 0.0, + ti2: 0.0, + }; + + let r20 = feautr(freq, ¶ms_f20); + let r40 = feautr(freq, ¶ms_f40); + + // 不同温度应该给出不同结果 + assert!((r20 - r40).abs() > 1e-10); + } + + #[test] + fn test_mixed_interpolation() { + let freq = 2.997925e18 / 1215.685; + + // 混合温度插值 + let params = FeautrParams { + jt: 3, // 使用 F40 作为主表 + ti0: 0.5, + ti1: 0.3, + ti2: 0.2, + }; + + let result = feautr(freq, ¶ms); + assert!(result >= 0.0 && result <= 1.0); + } + + #[test] + fn test_table_values() { + // 验证表值正确：在 dlam=-150Å 处，F40=0.7862 + // dlam = CL/freq - LYA，要得到 dlam=-150，需要 freq = CL/(LYA-150) + let freq = 2.997925e18 / (1215.685 - 150.0); // dlam = -150 + let params = FeautrParams { + jt: 3, // F40 表 + ti0: 1.0, + ti1: 0.0, + ti2: 0.0, + }; + + let result = feautr(freq, ¶ms); + // F40(-150) = 0.7862, 所以 y ≈ 0.7862, result = 0.5 * (0.7862 + 1) ≈ 0.8931 + let expected = 0.5 * (0.7862 + 1.0); + assert_relative_eq!(result, expected, epsilon = 1e-4); + } +} diff --git a/src/math/gamhe.rs b/src/math/gamhe.rs new file mode 100644 index 0000000..62b83c6 --- /dev/null +++ b/src/math/gamhe.rs @@ -0,0 +1,263 @@ +//! 中性氦 Stark 加宽参数计算。 +//! +//! 重构自 SYNSPEC `gamhe.f` +//! +//! 基于 Dimitrijevic 和 Sahal-Brechot (1984, J.Q.S.R.T. 31, 301) +//! 或 Freudenstein 和 Cooper (1978, AP.J. 224, 1079) 的数据。 + +/// Gamhe 计算所需的模型层参数。 +#[derive(Debug, Clone)] +pub struct GamheParams { + /// 温度插值索引 (JT) + pub jt: usize, + /// 温度插值系数 0 + pub ti0: f64, + /// 温度插值系数 1 + pub ti1: f64, + /// 温度插值系数 2 + pub ti2: f64, +} + +/// 中性氦 Stark 加宽参数。 +/// +/// 包含电子和质子贡献的加宽数据。 +pub struct GamheData { + /// W 数组: [电子加宽 @ 5000K, @ 10000K, @ 20000K, @ 40000K, 波长(Å)] + /// 按谱线索引 (1-20) 存储 + pub w: [[f64; 5]; 20], + /// V 数组: 质子加宽 @ 5000K, @ 10000K, @ 20000K, @ 40000K + pub v: [[f64; 4]; 20], + /// C 数组: 备用系数 (当 W=0 时使用) + pub c: [f64; 20], +} + +impl GamheData { + /// 创建默认的 GamheData，包含所有谱线数据。 + pub fn new() -> Self { + // W 数组: 电子加宽参数 + 波长 + // Fortran DATA 是列优先，需要转置 + const W_DATA: [[f64; 5]; 20] = [ + [5.990, 6.650, 6.610, 6.210, 3819.60], + [2.950, 3.130, 3.230, 3.300, 3867.50], + [0.000, 0.000, 0.000, 0.000, 3871.79], + [0.142, 0.166, 0.182, 0.190, 3888.65], + [0.000, 0.000, 0.000, 0.000, 3926.53], + [1.540, 1.480, 1.400, 1.290, 3964.73], + [41.600, 50.500, 57.400, 65.800, 4009.27], + [1.320, 1.350, 1.380, 1.460, 4120.80], + [7.830, 8.750, 8.690, 8.040, 4143.76], + [5.830, 6.370, 6.820, 6.990, 4168.97], + [0.000, 0.000, 0.000, 0.000, 4437.55], + [1.630, 1.610, 1.490, 1.350, 4471.50], + [0.588, 0.620, 0.641, 0.659, 4713.20], + [2.600, 2.480, 2.240, 1.960, 4921.93], + [0.627, 0.597, 0.568, 0.532, 5015.68], + [1.050, 1.090, 1.110, 1.140, 5047.74], + [0.277, 0.298, 0.296, 0.293, 5875.70], + [0.714, 0.666, 0.602, 0.538, 6678.15], + [3.490, 3.630, 3.470, 3.190, 4026.20], + [4.970, 5.100, 4.810, 4.310, 4387.93], + ]; + + // V 数组: 质子加宽参数 + const V_DATA: [[f64; 4]; 20] = [ + [1.520, 4.540, 9.140, 10.200], + [0.607, 0.710, 0.802, 0.901], + [0.000, 0.000, 0.000, 0.000], + [0.0396, 0.0434, 0.0476, 0.0526], + [0.000, 0.000, 0.000, 0.000], + [0.507, 0.585, 0.665, 0.762], + [0.930, 1.710, 13.600, 27.200], + [0.288, 0.325, 0.365, 0.410], + [1.330, 6.800, 12.900, 14.300], + [1.100, 1.370, 1.560, 1.760], + [0.000, 0.000, 0.000, 0.000], + [1.340, 1.690, 1.820, 1.630], + [0.128, 0.143, 0.161, 0.181], + [2.040, 2.740, 2.950, 2.740], + [0.187, 0.210, 0.237, 0.270], + [0.231, 0.260, 0.291, 0.327], + [0.0591, 0.0650, 0.0719, 0.0799], + [0.231, 0.260, 0.295, 0.339], + [2.180, 3.760, 4.790, 4.560], + [1.860, 5.320, 7.070, 7.150], + ]; + + // C 数组: 备用系数 + // DATA C /2*0.,1.83E-4,0.,1.13E-4,5*0.,1.6E-4,9*0./ + const C_DATA: [f64; 20] = [ + 0.0, 0.0, 1.83e-4, 0.0, 1.13e-4, + 0.0, 0.0, 0.0, 0.0, 0.0, + 1.6e-4, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, + ]; + + Self { + w: W_DATA, + v: V_DATA, + c: C_DATA, + } + } +} + +impl Default for GamheData { + fn default() -> Self { + Self::new() + } +} + +/// 计算中性氦 Stark 加宽参数 GAM。 +/// +/// # 参数 +/// +/// * `ind` - 谱线索引 (1-20, 1-indexed) +/// * `t` - 温度 (K) +/// * `ane` - 电子数密度 +/// * `anp` - 质子数密度 +/// * `params` - 模型层插值参数 +/// * `data` - Stark 加宽数据表 +/// +/// # 返回值 +/// +/// Stark 加宽参数 GAM (Å) +/// +/// # 算法 +/// +/// 如果 W(1,IND) != 0: +/// GAM = (电子贡献 * ANE + 质子贡献 * ANP) * 1.884e3 / 波长^2 +/// 否则: +/// GAM = C(IND) * T^0.16667 * ANE +pub fn gamhe(ind: usize, t: f64, ane: f64, anp: f64, params: &GamheParams, data: &GamheData) -> f64 { + // 转换为 0-indexed + let ind_idx = ind - 1; + + // 检查是否有有效的电子加宽数据 + if data.w[ind_idx][0] == 0.0 { + // 使用备用公式 + let gam = data.c[ind_idx] * t.powf(0.16667) * ane; + return gam.max(0.0); + } + + // 温度插值 + let jt_idx = params.jt; // 已是 0-indexed + + // 电子加宽贡献 + let w_electron = params.ti0 * data.w[ind_idx][jt_idx] + + params.ti1 * data.w[ind_idx][jt_idx - 1] + + params.ti2 * data.w[ind_idx][jt_idx - 2]; + + // 质子加宽贡献 + let v_proton = params.ti0 * data.v[ind_idx][jt_idx] + + params.ti1 * data.v[ind_idx][jt_idx - 1] + + params.ti2 * data.v[ind_idx][jt_idx - 2]; + + // 波长 (Å) + let wavelength = data.w[ind_idx][4]; + + // 计算 GAM + let gam = (w_electron * ane + v_proton * anp) * 1.884e3 / (wavelength * wavelength); + + gam.max(0.0) +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + fn make_params(jt: usize) -> GamheParams { + GamheParams { + jt, + ti0: 1.0, + ti1: 0.0, + ti2: 0.0, + } + } + + #[test] + fn test_basic() { + let data = GamheData::new(); + let params = make_params(2); // jt >= 2 以访问 jt-2 + + // 使用谱线 1 (3819.60 Å) + let gam = gamhe(1, 10000.0, 1e13, 1e12, ¶ms, &data); + + assert!(gam >= 0.0); + assert!(gam.is_finite()); + } + + #[test] + fn test_zero_w_case() { + let data = GamheData::new(); + let params = make_params(2); + + // 谱线 3 的 W(1,3) = 0，使用备用公式 + let gam = gamhe(3, 10000.0, 1e13, 1e12, ¶ms, &data); + + // C(3) = 1.83e-4 + let expected = 1.83e-4 * 10000.0_f64.powf(0.16667) * 1e13; + assert_relative_eq!(gam, expected, epsilon = 1e-6); + } + + #[test] + fn test_electron_contribution() { + let data = GamheData::new(); + let params = make_params(2); + + // 纯电子贡献 (ANP = 0) + let gam_e = gamhe(1, 10000.0, 1e13, 0.0, ¶ms, &data); + assert!(gam_e > 0.0); + } + + #[test] + fn test_proton_contribution() { + let data = GamheData::new(); + let params = make_params(2); + + // 纯质子贡献 (ANE = 0) + let gam_p = gamhe(1, 10000.0, 0.0, 1e13, ¶ms, &data); + assert!(gam_p > 0.0); + } + + #[test] + fn test_density_scaling() { + let data = GamheData::new(); + let params = make_params(2); + + let gam1 = gamhe(1, 10000.0, 1e13, 1e12, ¶ms, &data); + let gam2 = gamhe(1, 10000.0, 2e13, 2e12, ¶ms, &data); + + // 密度翻倍，GAM 应该翻倍 + assert_relative_eq!(gam2, 2.0 * gam1, epsilon = 1e-10); + } + + #[test] + fn test_negative_result_clamped() { + let data = GamheData::new(); + let params = GamheParams { + jt: 2, + ti0: -1.0, // 负系数可能导致负 GAM + ti1: 0.0, + ti2: 0.0, + }; + + let gam = gamhe(1, 10000.0, 1e13, 1e12, ¶ms, &data); + + // 负值应该被钳制为 0 + assert!(gam >= 0.0); + } + + #[test] + fn test_wavelength_scaling() { + let data = GamheData::new(); + let params = make_params(2); + + // 谱线 1: 3819.60 Å, 谱线 2: 3867.50 Å + let gam1 = gamhe(1, 10000.0, 1e13, 0.0, ¶ms, &data); + let gam2 = gamhe(2, 10000.0, 1e13, 0.0, ¶ms, &data); + + // GAM ∝ 1/波长^2，较长波长应该有较小的 GAM + // (假设其他参数相似) + assert!(gam1.is_finite() && gam2.is_finite()); + } +} diff --git a/src/math/griem.rs b/src/math/griem.rs new file mode 100644 index 0000000..d713719 --- /dev/null +++ b/src/math/griem.rs @@ -0,0 +1,177 @@ +//! Griem Stark 阻尼参数计算。 +//! +//! 重构自 SYNSPEC `griem.f` +//! +//! 从输入的 Stark 宽度值计算 Stark 阻尼参数 (GAM)。 +//! 输入值对应 T=5000, 10000, 20000, 40000 K， +//! 以及 NE=1e16 (中性原子) 或 NE=1e17 (离子)。 + +/// Griem 计算所需的模型层参数。 +#[derive(Debug, Clone)] +pub struct GriemParams { + /// 温度插值索引 (JT) + pub jt: usize, + /// 温度插值系数 0 + pub ti0: f64, + /// 温度插值系数 1 + pub ti1: f64, + /// 温度插值系数 2 + pub ti2: f64, +} + +/// 计算 Griem Stark 阻尼参数 GAM。 +/// +/// # 参数 +/// +/// * `t` - 温度 (K)，如果 <= 0 则返回 0 +/// * `ane` - 电子数密度 +/// * `ion` - 电离级数 (1=中性, 2=一次电离, ...) +/// * `fr` - 频率相关因子 +/// * `wgr` - Stark 宽度数组 [W@5000K, W@10000K, W@20000K, W@40000K] +/// * `params` - 模型层插值参数 +/// +/// # 返回值 +/// +/// Stark 阻尼参数 GAM +/// +/// # 算法 +/// +/// ```text +/// GAM = (插值后的WGR) * ANE * 1e-10 * FR * 1e-10 * FR * 4.2e-14 +/// 如果 ION > 1: GAM = GAM * 0.1 +/// 如果 GAM < 0: GAM = 0 +/// ``` +pub fn griem(t: f64, ane: f64, ion: i32, fr: f64, wgr: &[f64; 4], params: &GriemParams) -> f64 { + // 温度检查 + if t <= 0.0 { + return 0.0; + } + + // 温度插值 + let jt_idx = params.jt; // 已是 0-indexed + + // 插值计算 Stark 宽度 + let w_interp = params.ti0 * wgr[jt_idx] + + params.ti1 * wgr[jt_idx - 1] + + params.ti2 * wgr[jt_idx - 2]; + + // 计算 GAM + let mut gam = w_interp * ane * 1e-10 * fr * 1e-10 * fr * 4.2e-14; + + // 离子的修正 + if ion > 1 { + gam *= 0.1; + } + + // 确保非负 + gam.max(0.0) +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + fn make_params(jt: usize) -> GriemParams { + GriemParams { + jt, + ti0: 1.0, + ti1: 0.0, + ti2: 0.0, + } + } + + #[test] + fn test_basic() { + let params = make_params(2); + let wgr = [1.0, 1.0, 1.0, 1.0]; + + let gam = griem(10000.0, 1e13, 1, 1.0, &wgr, ¶ms); + + assert!(gam >= 0.0); + assert!(gam.is_finite()); + } + + #[test] + fn test_zero_temperature() { + let params = make_params(2); + let wgr = [1.0, 1.0, 1.0, 1.0]; + + let gam = griem(0.0, 1e13, 1, 1.0, &wgr, ¶ms); + + assert_relative_eq!(gam, 0.0); + } + + #[test] + fn test_negative_temperature() { + let params = make_params(2); + let wgr = [1.0, 1.0, 1.0, 1.0]; + + let gam = griem(-100.0, 1e13, 1, 1.0, &wgr, ¶ms); + + assert_relative_eq!(gam, 0.0); + } + + #[test] + fn test_ion_correction() { + let params = make_params(2); + let wgr = [1.0, 1.0, 1.0, 1.0]; + + let gam_neutral = griem(10000.0, 1e13, 1, 1.0, &wgr, ¶ms); + let gam_ion = griem(10000.0, 1e13, 2, 1.0, &wgr, ¶ms); + + // 离子的 GAM 应该是中性的 0.1 倍 + assert_relative_eq!(gam_ion, 0.1 * gam_neutral, epsilon = 1e-10); + } + + #[test] + fn test_density_scaling() { + let params = make_params(2); + let wgr = [1.0, 1.0, 1.0, 1.0]; + + let gam1 = griem(10000.0, 1e13, 1, 1.0, &wgr, ¶ms); + let gam2 = griem(10000.0, 2e13, 1, 1.0, &wgr, ¶ms); + + // 密度翻倍，GAM 应该翻倍 + assert_relative_eq!(gam2, 2.0 * gam1, epsilon = 1e-10); + } + + #[test] + fn test_fr_scaling() { + let params = make_params(2); + let wgr = [1.0, 1.0, 1.0, 1.0]; + + let gam1 = griem(10000.0, 1e13, 1, 1.0, &wgr, ¶ms); + let gam2 = griem(10000.0, 1e13, 1, 2.0, &wgr, ¶ms); + + // FR 翻倍，GAM 应该是 4 倍 (因为 FR^2) + assert_relative_eq!(gam2, 4.0 * gam1, epsilon = 1e-10); + } + + #[test] + fn test_temperature_interpolation() { + let params = make_params(2); + let wgr = [1.0, 2.0, 4.0, 8.0]; // 不同的温度宽度 + + let gam = griem(10000.0, 1e13, 1, 1.0, &wgr, ¶ms); + + // jt=2 使用 wgr[2]=4.0 + assert!(gam > 0.0); + } + + #[test] + fn test_negative_result_clamped() { + let params = GriemParams { + jt: 2, + ti0: -1.0, // 负系数 + ti1: 0.0, + ti2: 0.0, + }; + let wgr = [1.0, 1.0, 1.0, 1.0]; + + let gam = griem(10000.0, 1e13, 1, 1.0, &wgr, ¶ms); + + // 负值应该被钳制为 0 + assert_relative_eq!(gam, 0.0); + } +} diff --git a/src/math/intrp.rs b/src/math/intrp.rs new file mode 100644 index 0000000..e36cd30 --- /dev/null +++ b/src/math/intrp.rs @@ -0,0 +1,262 @@ +//! 使用二分法的高效插值程序。 +//! +//! 重构自 SYNSPEC `intrp.f` +//! +//! 对波长/频率表进行线性插值，使用二分法查找插值区间。 + +/// 使用二分法进行线性插值。 +/// +/// 对输入的不透明度表进行插值，将值从 `wltab/absop` 网格插值到 `wlgrid` 网格。 +/// +/// # 参数 +/// +/// * `wltab` - 原始波长/频率表 (单调递增或递减) +/// * `absop` - 原始不透明度值 +/// * `wlgrid` - 目标波长/频率网格 +/// * `abgrd` - 输出插值后的不透明度值 (需预分配) +/// +/// # 算法 +/// +/// 1. 使用二分法查找每个目标点所在的区间 `[j, j+1]` +/// 2. 存储区间索引和插值系数 +/// 3. 对每个目标点执行线性插值 +/// +/// # 注意 +/// +/// - `wltab` 必须是单调的 (递增或递减) +/// - `wlgrid` 中的点应该在 `wltab` 的范围内 +/// +/// # 示例 +/// +/// ``` +/// use spectrarust::math::intrp::intrp; +/// +/// // 原始数据表 +/// let wltab = vec![400.0, 450.0, 500.0, 550.0, 600.0]; +/// let absop = vec![0.1, 0.2, 0.3, 0.25, 0.15]; +/// +/// // 目标网格 +/// let wlgrid = vec![420.0, 480.0, 520.0, 580.0]; +/// let mut abgrd = vec![0.0; 4]; +/// +/// intrp(&wltab, &absop, &wlgrid, &mut abgrd); +/// +/// // 验证插值结果在合理范围内 +/// for &val in &abgrd { +/// assert!(val >= 0.0); +/// } +/// ``` +pub fn intrp(wltab: &[f64], absop: &[f64], wlgrid: &[f64], abgrd: &mut [f64]) { + let nfr = wltab.len(); + let nfgrid = wlgrid.len(); + + if nfr < 2 || nfgrid == 0 { + return; + } + + // 工作数组：存储区间索引和插值系数 + let mut jint = vec![0usize; nfgrid]; + let mut yint = vec![0.0f64; nfgrid]; + + // 第一遍：使用二分法查找区间 + let fr1 = wltab[0]; + let fr2 = wltab[nfr - 1]; + let ascending = fr2 > fr1; + + for ij in 0..nfgrid { + let xint = wlgrid[ij]; + + // 二分法查找 + let mut jl: isize = 0; + let mut ju: isize = nfr as isize; + + while ju - jl > 1 { + let jm = (ju + jl) / 2; + let wltab_jm = wltab[jm as usize]; + + // Fortran: if((fr2.gt.fr1).eqv.(xint.gt.wltab(jm))) + // 如果表是递增的且 xint > wltab[jm]，或者表是递减的且 xint <= wltab[jm] + if (ascending && xint > wltab_jm) || (!ascending && xint <= wltab_jm) { + jl = jm; + } else { + ju = jm; + } + } + + // 边界处理 + let mut j = jl as usize; + if j >= nfr - 1 { + j = nfr - 2; + } + if j == 0 && xint < wltab[0] { + // 外推：使用第一个区间 + j = 0; + } + + jint[ij] = j; + // Fortran 索引从 1 开始，Rust 从 0 开始 + // yint(ij)=1./(wltab(j+1)-wltab(j)) + let diff = wltab[j + 1] - wltab[j]; + yint[ij] = if diff.abs() > 1e-30 { 1.0 / diff } else { 0.0 }; + } + + // 第二遍：执行线性插值 + for ij in 0..nfgrid { + let j = jint[ij]; + // rc=(absop(j+1)-absop(j))*yint(ij) + let rc = (absop[j + 1] - absop[j]) * yint[ij]; + // abgrd(ij)=rc*(wlgrid(ij)-wltab(j))+absop(j) + abgrd[ij] = rc * (wlgrid[ij] - wltab[j]) + absop[j]; + } +} + +/// 使用二分法进行线性插值（返回新向量版本）。 +/// +/// # 参数 +/// +/// * `wltab` - 原始波长/频率表 +/// * `absop` - 原始不透明度值 +/// * `wlgrid` - 目标波长/频率网格 +/// +/// # 返回值 +/// +/// 返回插值后的不透明度向量 +/// +/// # 示例 +/// +/// ``` +/// use spectrarust::math::intrp::intrp_to_vec; +/// +/// let wltab = vec![400.0, 500.0, 600.0]; +/// let absop = vec![0.1, 0.3, 0.2]; +/// let wlgrid = vec![450.0, 550.0]; +/// +/// let result = intrp_to_vec(&wltab, &absop, &wlgrid); +/// assert_eq!(result.len(), 2); +/// ``` +pub fn intrp_to_vec(wltab: &[f64], absop: &[f64], wlgrid: &[f64]) -> Vec { + let nfgrid = wlgrid.len(); + let mut abgrd = vec![0.0; nfgrid]; + intrp(wltab, absop, wlgrid, &mut abgrd); + abgrd +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_basic_interpolation() { + // 简单的线性插值测试 + let wltab = vec![400.0, 500.0, 600.0]; + let absop = vec![0.1, 0.3, 0.2]; + let wlgrid = vec![450.0]; // 中点 + let mut abgrd = vec![0.0; 1]; + + intrp(&wltab, &absop, &wlgrid, &mut abgrd); + + // 450 在 400 和 500 的中点，所以应该是 0.1 和 0.3 的中点 + assert_relative_eq!(abgrd[0], 0.2, epsilon = 1e-10); + } + + #[test] + fn test_multiple_points() { + let wltab = vec![0.0, 1.0, 2.0, 3.0, 4.0]; + let absop = vec![0.0, 1.0, 4.0, 9.0, 16.0]; // x^2 + let wlgrid = vec![0.5, 1.5, 2.5, 3.5]; + let mut abgrd = vec![0.0; 4]; + + intrp(&wltab, &absop, &wlgrid, &mut abgrd); + + // 线性插值结果 + assert_relative_eq!(abgrd[0], 0.5, epsilon = 1e-10); // 0.5 * (4-0) + 0 = 0.5? 不对，应该是 (1-0)*0.5 = 0.5 + assert_relative_eq!(abgrd[1], 2.5, epsilon = 1e-10); // (4-1)*0.5 + 1 = 2.5 + assert_relative_eq!(abgrd[2], 6.5, epsilon = 1e-10); // (9-4)*0.5 + 4 = 6.5 + assert_relative_eq!(abgrd[3], 12.5, epsilon = 1e-10); // (16-9)*0.5 + 9 = 12.5 + } + + #[test] + fn test_exact_points() { + // 测试恰好落在网格点上的情况 + let wltab = vec![100.0, 200.0, 300.0]; + let absop = vec![1.0, 2.0, 3.0]; + let wlgrid = vec![100.0, 200.0, 300.0]; + let mut abgrd = vec![0.0; 3]; + + intrp(&wltab, &absop, &wlgrid, &mut abgrd); + + // 应该得到精确值 + assert_relative_eq!(abgrd[0], 1.0, epsilon = 1e-10); + assert_relative_eq!(abgrd[1], 2.0, epsilon = 1e-10); + assert_relative_eq!(abgrd[2], 3.0, epsilon = 1e-10); + } + + #[test] + fn test_descending() { + // 测试递减序列 + let wltab = vec![600.0, 500.0, 400.0]; // 递减 + let absop = vec![0.2, 0.3, 0.1]; + let wlgrid = vec![550.0]; // 在 600 和 500 之间 + let mut abgrd = vec![0.0; 1]; + + intrp(&wltab, &absop, &wlgrid, &mut abgrd); + + // 550 在 600 和 500 的中点 + assert_relative_eq!(abgrd[0], 0.25, epsilon = 1e-10); + } + + #[test] + fn test_to_vec() { + let wltab = vec![0.0, 1.0, 2.0]; + let absop = vec![0.0, 2.0, 4.0]; + let wlgrid = vec![0.5, 1.5]; + + let result = intrp_to_vec(&wltab, &absop, &wlgrid); + + assert_eq!(result.len(), 2); + assert_relative_eq!(result[0], 1.0, epsilon = 1e-10); + assert_relative_eq!(result[1], 3.0, epsilon = 1e-10); + } + + #[test] + fn test_empty() { + let wltab: Vec = vec![]; + let absop: Vec = vec![]; + let wlgrid = vec![1.0, 2.0]; + let mut abgrd = vec![0.0; 2]; + + // 应该不崩溃，直接返回 + intrp(&wltab, &absop, &wlgrid, &mut abgrd); + } + + #[test] + fn test_single_point_table() { + let wltab = vec![100.0]; + let absop = vec![1.0]; + let wlgrid = vec![100.0]; + let mut abgrd = vec![0.0; 1]; + + // 单点表无法插值，应该安全处理 + intrp(&wltab, &absop, &wlgrid, &mut abgrd); + } + + #[test] + fn test_large_dataset() { + // 测试大数据集性能 + let n = 10000; + let wltab: Vec = (0..n).map(|i| i as f64).collect(); + let absop: Vec = (0..n).map(|i| (i as f64).sin()).collect(); + + let m = 1000; + let wlgrid: Vec = (0..m).map(|i| i as f64 * (n - 1) as f64 / (m - 1) as f64).collect(); + let mut abgrd = vec![0.0; m]; + + intrp(&wltab, &absop, &wlgrid, &mut abgrd); + + // 验证结果在合理范围内 + for &val in &abgrd { + assert!(val >= -1.0 && val <= 1.0); + } + } +} diff --git a/src/math/mod.rs b/src/math/mod.rs index 16656ee..251067a 100644 --- a/src/math/mod.rs +++ b/src/math/mod.rs @@ -54,12 +54,14 @@ mod contmd; mod contmp; mod convec; mod coolrt; +mod count_words; mod cross; mod cspec; mod ctdata; mod cubic; mod dielrc; mod dietot; +mod divhe2; mod divstr; mod dopgam; mod dmder; @@ -76,9 +78,12 @@ mod eps; mod erfcx; mod expo; mod expint; +mod extprf; +mod feautr; mod ffcros; mod gauleg; mod getwrd; +mod gamhe; mod gami; mod getlal; mod gamsp; @@ -87,6 +92,7 @@ mod ghydop; mod gaunt; mod gntk; mod gridp; +mod griem; mod gomini; mod grcor; mod greyd; @@ -107,6 +113,7 @@ mod inifrc; mod inifrt; mod inkul; mod interp; +mod intrp; mod inthyd; mod intlem; mod intxen; @@ -147,6 +154,7 @@ mod opfrac; mod opadd; mod opadd0; mod partf; +mod partdv; mod opahst; mod opacfa; mod opacf0; @@ -234,6 +242,7 @@ mod sigk; mod sigmar; mod state; mod sffhmi; +mod sffhmi_old; mod tabint; mod taufr1; mod sffhmi_add; @@ -337,12 +346,14 @@ pub use contmd::{contmd_pure, ContmdConfig, ContmdParams, ContmdOutput, CubconDa pub use contmp::{contmp, ContmpConfig, ContmpParams, ContmpOutput}; pub use convec::{convec, convc1, ConvecConfig, ConvecParams, ConvecOutput, Convc1Output}; pub use coolrt::{coolrt_pure, CoolrtParams, CoolrtOutput, compute_taud, find_fe2_ion}; +pub use count_words::count_words; pub use cross::{cross, crossd}; pub use cspec::cspec; pub use ctdata::{hction, hctrecom, CTION, CTRECOMB}; pub use cubic::{cubic, CubicCon}; pub use dielrc::dielrc; pub use dietot::{dietot, DietotParams}; +pub use divhe2::divhe2; pub use divstr::divstr; pub use dopgam::dopgam; pub use dmder::{dmder, DepthDeriv}; @@ -359,9 +370,12 @@ pub use eps::eps; pub use erfcx::{erfcin, erfcx}; pub use expo::expo; pub use expint::{eint, expinx}; +pub use extprf::extprf; +pub use feautr::{feautr, FeautrParams}; pub use ffcros::ffcros; pub use gauleg::gauleg; pub use getwrd::getwrd; +pub use gamhe::{gamhe, GamheData, GamheParams}; pub use gami::gami; pub use getlal::{getlal, GetlalParams, GetlalResult}; pub use gamsp::gamsp; @@ -370,6 +384,7 @@ pub use ghydop::{ghydop, GhydopParams, GhydopResult}; pub use gaunt::gaunt; pub use gntk::gntk; pub use gridp::gridp; +pub use griem::{griem, GriemParams}; pub use gomini::{gomini, GominiParams, GominiResult}; pub use grcor::grcor; pub use greyd::{ @@ -393,6 +408,7 @@ pub use inifrc::{inifrc, InifrcParams, InifrcOutput}; pub use inifrt::{inifrt, InifrtParams, InifrtOutput}; pub use inkul::{inkul, inkul_pure, InkulParams, InkulOutput, ColKur, Lined, LineRecord}; pub use interp::interp; +pub use intrp::{intrp, intrp_to_vec}; pub use inthyd::inthyd; pub use intlem::intlem; pub use intxen::intxen; @@ -427,6 +443,7 @@ pub use opfrac::{opfrac_pure, opfrac_init, OpfracParams, OpfracOutput, PfOptB}; pub use opadd::{opadd, OpaddInput, OpaddModel, OpaddSwitches, OpaddOutput, OpaddCache}; pub use opadd0::{opadd0, Opadd0Params, Opadd0FreqData, Opadd0OutputState}; pub use partf::{partf_pure, PartfParams, PartfOutput, PartfMode}; +pub use partdv::{partdv, partdv_with_params, PartdvParams}; pub use opahst::{opahst, OpahstParams, OpahstOutput, LymanConfig, BalmerConfig, NLMX}; pub use opacfa::{opacfa, OpacfaParams, OpacfaOutput}; pub use opacf0::{opacf0, Opacf0Config, Opacf0ModelState, Opacf0AtomicParams, Opacf0FreqParams, Opacf0Output}; @@ -545,6 +562,7 @@ pub use state::{ get_max_ionization, get_element_symbol, }; pub use sffhmi::sffhmi; +pub use sffhmi_old::sffhmi_old; pub use sffhmi_add::sffhmi_add; pub use spsigk::spsigk; pub use tabint::{tabint, IntCff, OpacTable, TabintParams}; diff --git a/src/math/partdv.rs b/src/math/partdv.rs new file mode 100644 index 0000000..e8c187f --- /dev/null +++ b/src/math/partdv.rs @@ -0,0 +1,290 @@ +//! 配分函数计算（考虑压力效应）。 +//! +//! 重构自 SYNSPEC `partdv.f` +//! +//! 计算原子或离子的配分函数，包含等离子体效应的修正。 + +/// 计算配分函数（考虑压力效应）。 +/// +/// 使用 Traving 等人的方法计算配分函数，包含等离子体效应对高能级的抑制。 +/// +/// # 参数 +/// +/// * `temp` - 温度 (K) +/// * `dne` - 电子密度 (cm^-3) +/// * `z` - 原子序数 +/// * `nlev` - 能级数量 +/// * `ne` - 各能级的主量子数数组 +/// * `gee` - 各能级的统计权重数组 +/// * `enrgy` - 各能级的激发能量数组 (eV) +/// * `s` - 各能级的屏蔽常数数组 +/// +/// # 返回值 +/// +/// 返回配分函数 U +/// +/// # 算法 +/// +/// 1. 计算 ET = TEMP / 11604.8 (温度转换为 eV 单位) +/// 2. 计算压力修正参数 P = 14.69 - 0.6667 * log10(DNE) +/// 3. 对每个能级： +/// - 计算有效核电荷 ZSTAR = Z - S +/// - 计算抑制因子 W = min(1, P + 4*log10(ZSTAR) - 4*log10(NE)) +/// - 如果激发能不太高，添加贡献：GEE * W * exp(-ENRGY/ET) +/// +/// # 示例 +/// +/// ``` +/// use spectrarust::math::partdv::partdv; +/// +/// // 简单的氢原子模型 +/// let temp = 5780.0; // 太阳表面温度 +/// let dne = 1e13; // 电子密度 +/// let z = 1; // 氢原子序数 +/// let ne = vec![1, 2, 3]; // 主量子数 +/// let gee = vec![2.0, 8.0, 18.0]; // 统计权重 (2n^2) +/// let enrgy = vec![0.0, 10.2, 12.1]; // 激发能 (eV) +/// let s = vec![0.0, 0.0, 0.0]; // 屏蔽常数 +/// +/// let u = partdv(temp, dne, z, &ne, &gee, &enrgy, &s); +/// assert!(u > 0.0); +/// ``` +pub fn partdv( + temp: f64, + dne: f64, + z: i32, + ne: &[i32], + gee: &[f64], + enrgy: &[f64], + s: &[f64], +) -> f64 { + let nlev = ne.len().min(gee.len()).min(enrgy.len()).min(s.len()); + + if nlev == 0 { + return 0.0; + } + + // ET = TEMP / 11604.8 (温度转换为 eV 单位，11604.8 K = 1 eV) + let et = temp / 11604.8; + + // 压力修正参数 + // P = 14.69 - 0.20 - 0.6667 * log10(DNE) + // 注意：14.69 - 0.20 = 14.49 + let p = 14.49 - 0.6667 * dne.log10(); + + let mut u = 0.0; + + for i in 0..nlev { + let u1 = ne[i] as f64; + let zstar = z as f64 - s[i]; + + // 计算抑制因子 W + let w = if zstar > 0.0 { + let w_raw = p + 4.0 * zstar.log10() - 4.0 * u1.log10(); + if w_raw > 1.0 { 1.0 } else { w_raw } + } else { + 0.0 + }; + + // 如果激发能不太高（< 65 * ET），添加贡献 + if (enrgy[i] / et) < 65.0 { + let contribution = gee[i] * w * (-enrgy[i] / et).exp(); + u += contribution; + } + } + + u +} + +/// 配分函数计算的参数结构体。 +#[derive(Debug, Clone)] +pub struct PartdvParams<'a> { + /// 温度 (K) + pub temp: f64, + /// 电子密度 (cm^-3) + pub dne: f64, + /// 原子序数 + pub z: i32, + /// 主量子数数组 + pub ne: &'a [i32], + /// 统计权重数组 + pub gee: &'a [f64], + /// 激发能数组 (eV) + pub enrgy: &'a [f64], + /// 屏蔽常数数组 + pub s: &'a [f64], +} + +/// 使用参数结构体计算配分函数。 +/// +/// # 示例 +/// +/// ``` +/// use spectrarust::math::partdv::{partdv_with_params, PartdvParams}; +/// +/// let params = PartdvParams { +/// temp: 5780.0, +/// dne: 1e13, +/// z: 1, +/// ne: &[1, 2, 3], +/// gee: &[2.0, 8.0, 18.0], +/// enrgy: &[0.0, 10.2, 12.1], +/// s: &[0.0, 0.0, 0.0], +/// }; +/// +/// let u = partdv_with_params(¶ms); +/// assert!(u > 0.0); +/// ``` +pub fn partdv_with_params(params: &PartdvParams) -> f64 { + partdv( + params.temp, + params.dne, + params.z, + params.ne, + params.gee, + params.enrgy, + params.s, + ) +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_basic() { + // 简单的氢原子模型 + let temp = 5780.0; + let dne = 1e13; + let z = 1; + let ne = vec![1, 2, 3]; + let gee = vec![2.0, 8.0, 18.0]; + let enrgy = vec![0.0, 10.2, 12.1]; + let s = vec![0.0, 0.0, 0.0]; + + let u = partdv(temp, dne, z, &ne, &gee, &enrgy, &s); + + // 基态贡献最大 + assert!(u > 0.0); + // 配分函数应该大于基态统计权重 + assert!(u >= 2.0); + } + + #[test] + fn test_ground_state_only() { + // 只有基态 + let temp = 5000.0; + let dne = 1e13; + let z = 1; + let ne = vec![1]; + let gee = vec![2.0]; + let enrgy = vec![0.0]; + let s = vec![0.0]; + + let u = partdv(temp, dne, z, &ne, &gee, &enrgy, &s); + + // 基态贡献 = 2 * W * exp(0) = 2 * W + // W = P + 4*log10(1) - 4*log10(1) = P + // P = 14.49 - 0.6667 * log10(1e13) = 14.49 - 8.6671 = 5.82 + // W = min(1, 5.82) = 1 + // 所以 U = 2 * 1 * 1 = 2 + assert_relative_eq!(u, 2.0, epsilon = 0.01); + } + + #[test] + fn test_high_energy_cutoff() { + // 高能级被截断 + let temp = 5000.0; + let dne = 1e13; + let z = 1; + let ne = vec![1, 10]; + let gee = vec![2.0, 200.0]; + let enrgy = vec![0.0, 100.0]; // 100 eV 对于 5000K 来说太高了 + let s = vec![0.0, 0.0]; + + let u = partdv(temp, dne, z, &ne, &gee, &enrgy, &s); + + // 第二个能级贡献应该很小 + assert!(u < 10.0); + } + + #[test] + fn test_pressure_suppression() { + // 高电子密度 -> 更强的压力抑制 + // 需要非常高的电子密度才能看到抑制效果 + let temp = 5780.0; + let z = 1; + let ne = vec![1, 2, 3, 4, 5]; + let gee = vec![2.0, 8.0, 18.0, 32.0, 50.0]; + let enrgy = vec![0.0, 10.2, 12.1, 13.1, 13.6]; + let s = vec![0.0, 0.0, 0.0, 0.0, 0.0]; + + let low_dne = 1e10; + let high_dne = 1e20; // 需要非常高的密度 + + let u_low = partdv(temp, low_dne, z, &ne, &gee, &enrgy, &s); + let u_high = partdv(temp, high_dne, z, &ne, &gee, &enrgy, &s); + + // 高电子密度下，高能级被抑制更多 + assert!(u_low > u_high, "u_low ({}) should be > u_high ({})", u_low, u_high); + } + + #[test] + fn test_effective_charge() { + // 有效核电荷的影响 + let temp = 5780.0; + let dne = 1e13; + let ne = vec![1]; + let gee = vec![2.0]; + let enrgy = vec![0.0]; + + // 无屏蔽 + let s_no = vec![0.0]; + let u_no = partdv(temp, dne, 1, &ne, &gee, &enrgy, &s_no); + + // 有屏蔽 + let s_yes = vec![0.5]; + let u_yes = partdv(temp, dne, 1, &ne, &gee, &enrgy, &s_yes); + + // 两者基态贡献相同 (W = 1) + assert_relative_eq!(u_no, u_yes, epsilon = 0.01); + } + + #[test] + fn test_empty() { + let u = partdv(5000.0, 1e13, 1, &[], &[], &[], &[]); + assert_relative_eq!(u, 0.0); + } + + #[test] + fn test_with_params() { + let params = PartdvParams { + temp: 5780.0, + dne: 1e13, + z: 1, + ne: &[1, 2], + gee: &[2.0, 8.0], + enrgy: &[0.0, 10.2], + s: &[0.0, 0.0], + }; + + let u = partdv_with_params(¶ms); + assert!(u > 0.0); + } + + #[test] + fn test_he_like() { + // 类氦离子 + let temp = 20000.0; + let dne = 1e14; + let z = 2; + let ne = vec![1, 2]; + let gee = vec![1.0, 3.0]; // 1s: 单态, 2s: 单态 + let enrgy = vec![0.0, 20.6]; // He+ 的激发能 + let s = vec![0.0, 0.0]; + + let u = partdv(temp, dne, z, &ne, &gee, &enrgy, &s); + assert!(u > 0.0); + } +} diff --git a/src/math/sffhmi_old.rs b/src/math/sffhmi_old.rs new file mode 100644 index 0000000..9f5df37 --- /dev/null +++ b/src/math/sffhmi_old.rs @@ -0,0 +1,136 @@ +//! H- 自由-自由截面计算。 +//! +//! 重构自 SYNSPEC `sffhmi_old.f` +//! +//! 根据 Kurucz (1970, SAO 309, P.80) 的公式计算 H- 离子的自由-自由吸收截面。 + +/// 计算 H- 自由-自由吸收截面。 +/// +/// # 参数 +/// +/// * `popi` - H- 离子布居数 +/// * `fr` - 频率 (Hz) +/// * `t` - 温度 (K) +/// +/// # 返回值 +/// +/// 返回 H- 自由-自由吸收截面 (cm²) +/// +/// # 公式 +/// +/// 根据 Kurucz (1970) 的公式: +/// ```text +/// σ = (1.3727e-25 + (4.3748e-10 - 2.5993e-7/T) / fr) * popi / fr +/// ``` +/// +/// # 示例 +/// +/// ``` +/// use spectrarust::math::sffhmi_old::sffhmi_old; +/// +/// // 典型恒星大气条件 +/// let popi = 1e-10; // H- 布居数 +/// let fr = 5e14; // 频率 (Hz), 约 600 nm +/// let t = 5780.0; // 太阳表面温度 +/// +/// let cross_section = sffhmi_old(popi, fr, t); +/// assert!(cross_section > 0.0); +/// ``` +pub fn sffhmi_old(popi: f64, fr: f64, t: f64) -> f64 { + // Kurucz (1970, SAO 309, P.80) 公式 + // σ = (1.3727e-25 + (4.3748e-10 - 2.5993e-7/T) / fr) * popi / fr + let term1 = 1.3727e-25; + let term2 = (4.3748e-10 - 2.5993e-7 / t) / fr; + + (term1 + term2) * popi / fr +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_basic() { + // 基本功能测试 + let popi = 1e-10; + let fr = 5e14; + let t = 5780.0; + + let result = sffhmi_old(popi, fr, t); + + // 手动计算验证 + let term1 = 1.3727e-25; + let term2 = (4.3748e-10 - 2.5993e-7 / t) / fr; + let expected = (term1 + term2) * popi / fr; + + assert_relative_eq!(result, expected, epsilon = 1e-40); + } + + #[test] + fn test_solar_conditions() { + // 太阳条件下的测试 + let popi = 1e-8; + let fr = 6e14; // 约 500 nm + let t = 5780.0; + + let result = sffhmi_old(popi, fr, t); + assert!(result > 0.0); + assert!(result < 1e-30); // 截面应该很小 + } + + #[test] + fn test_hot_star() { + // 热星条件 + let popi = 1e-12; + let fr = 1e15; // 紫外 + let t = 25000.0; + + let result = sffhmi_old(popi, fr, t); + assert!(result > 0.0); + } + + #[test] + fn test_cool_star() { + // 冷星条件 + let popi = 1e-6; + let fr = 3e14; // 红外 + let t = 4000.0; + + let result = sffhmi_old(popi, fr, t); + assert!(result > 0.0); + } + + #[test] + fn test_frequency_dependence() { + // 验证频率依赖性：高频 -> 小截面 + let popi = 1e-10; + let t = 5780.0; + + let low_fr = 3e14; + let high_fr = 1e15; + + let low_result = sffhmi_old(popi, low_fr, t); + let high_result = sffhmi_old(popi, high_fr, t); + + // 高频截面应该更小 + assert!(low_result > high_result); + } + + #[test] + fn test_temperature_dependence() { + // 验证温度依赖性 + let popi = 1e-10; + let fr = 5e14; + + let low_t = 4000.0; + let high_t = 10000.0; + + let low_result = sffhmi_old(popi, fr, low_t); + let high_result = sffhmi_old(popi, fr, high_t); + + // 温度影响较小，主要影响第二项 + assert!(low_result > 0.0); + assert!(high_result > 0.0); + } +}