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

feat: 添加 9 个 SYNSPEC 数学模块 (第8批)

Browse Source 
新增模块:
- 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 <noreply@anthropic.com>
This commit is contained in:
Asfmq 2026-03-25 06:52:44 +08:00
parent c0e70ef895
commit 0674b4f174
11 changed files with 1796 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
.claude/skills/fortran-analyzer/scripts/analyze_fortran.py
View File

@ -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:
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 ""

93
src/math/count_words.rs Normal file
View File

@ -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<char> = 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);
}
}

132
src/math/divhe2.rs Normal file
View File

@ -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);
}
}

151
src/math/extprf.rs Normal file
View File

@ -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);
}
}

236
src/math/feautr.rs Normal file
View File

@ -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, &params);
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, &params);
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, &params);
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, &params_f20);
let r40 = feautr(freq, &params_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, &params);
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, &params);
// 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);
}
}

263
src/math/gamhe.rs Normal file
View File

@ -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, &params, &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, &params, &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, &params, &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, &params, &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, &params, &data);
let gam2 = gamhe(1, 10000.0, 2e13, 2e12, &params, &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, &params, &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, &params, &data);
let gam2 = gamhe(2, 10000.0, 1e13, 0.0, &params, &data);
// GAM ∝ 1/波长^2，较长波长应该有较小的 GAM
// (假设其他参数相似)
assert!(gam1.is_finite() && gam2.is_finite());
}
}

177
src/math/griem.rs Normal file
View File

@ -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, &params);
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, &params);
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, &params);
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, &params);
let gam_ion = griem(10000.0, 1e13, 2, 1.0, &wgr, &params);
// 离子的 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, &params);
let gam2 = griem(10000.0, 2e13, 1, 1.0, &wgr, &params);
// 密度翻倍，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, &params);
let gam2 = griem(10000.0, 1e13, 1, 2.0, &wgr, &params);
// 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, &params);
// 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, &params);
// 负值应该被钳制为 0
assert_relative_eq!(gam, 0.0);
}
}

262
src/math/intrp.rs Normal file
View File

@ -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<f64> {
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<f64> = vec![];
let absop: Vec<f64> = 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<f64> = (0..n).map(|i| i as f64).collect();
let absop: Vec<f64> = (0..n).map(|i| (i as f64).sin()).collect();
let m = 1000;
let wlgrid: Vec<f64> = (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);
}
}
}

18
src/math/mod.rs
View File

@ -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};

290
src/math/partdv.rs Normal file
View File

@ -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(&params);
/// 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(&params);
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);
}
}

136
src/math/sffhmi_old.rs Normal file
View File

@ -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);
}
}