diff --git a/src/math/gvdw.rs b/src/math/gvdw.rs new file mode 100644 index 0000000..e514c06 --- /dev/null +++ b/src/math/gvdw.rs @@ -0,0 +1,206 @@ +//! Van der Waals broadening parameter evaluation. +//! +//! Refactored from SYNSPEC `gvdw.f` +//! +//! Evaluates the Van der Waals broadening parameter using either: +//! - Standard expression (classical) +//! - EXOMOL data (broadening by H2 and He) + +/// Parameters needed for GVDW calculation. +/// +/// These are typically extracted from COMMON blocks in the original Fortran. +pub struct GvdwParams { + /// Mode selector for the line list + /// - 0: standard expression + /// - >0: EXOMOL data (H2 and He broadening) + pub ivdwli: i32, + + /// Depth index (1-based in Fortran, 0-based in Rust) + pub id: usize, + + /// Line index (1-based in Fortran, 0-based in Rust) + pub il: usize, + + /// Standard Van der Waals width parameter for the line + pub gwm: f64, + + /// Van der Waals coefficient at depth id + pub vdwc: f64, + + /// Temperature at depth id (K) + pub temp: f64, + + /// He number density ratio at depth id (rrr(id, 1, 2)) + /// Index: (depth, ion=1, atom=2 for He) + pub anhe: f64, + + /// H2 number density at depth id + pub anh2: f64, + + // EXOMOL parameters + /// H2 broadening parameter + pub gvdwh2: f64, + + /// H2 temperature exponent + pub gexph2: f64, + + /// He broadening parameter + pub gvdwhe: f64, + + /// He temperature exponent + pub gexphe: f64, +} + +/// Calculates the Van der Waals broadening parameter. +/// +/// # Arguments +/// * `params` - Parameters struct containing all required data +/// +/// # Returns +/// The Van der Waals broadening parameter (width) +/// +/// # Algorithm +/// Two modes based on `ivdwli`: +/// - `ivdwli == 0`: Standard expression: `gwm * vdwc` +/// - `ivdwli > 0`: EXOMOL formula with H2 and He broadening: +/// ``` +/// con * T * ((296/T)^gexph2 * gvdwh2 * anh2 + (296/T)^gexphe * gvdwhe * anhe) +/// ``` +/// where `con = 4.1388e-12` (derived from `1e-6 * c * k`) +pub fn gvdw(params: &GvdwParams) -> f64 { + // Classical, original expression + if params.ivdwli == 0 { + return params.gwm * params.vdwc; + } + + // EXOMOL form - broadening by H2 and He + // con = 1.e-6 * c * k (where c is speed of light, k is Boltzmann constant) + const CON: f64 = 4.1388e-12; + const REF_TEMP: f64 = 296.0; + + let t = params.temp; + let anhe = params.anhe; + + let temp_ratio = REF_TEMP / t; + let h2_term = temp_ratio.powf(params.gexph2) * params.gvdwh2 * params.anh2; + let he_term = temp_ratio.powf(params.gexphe) * params.gvdwhe * anhe; + + CON * t * (h2_term + he_term) +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_gvdw_standard_mode() { + // Standard mode (ivdwli = 0) + let params = GvdwParams { + ivdwli: 0, + id: 0, + il: 0, + gwm: 2.5, + vdwc: 1.0e-7, + temp: 10000.0, + anhe: 0.1, + anh2: 0.0, + gvdwh2: 0.0, + gexph2: 0.0, + gvdwhe: 0.0, + gexphe: 0.0, + }; + + let result = gvdw(¶ms); + let expected = 2.5 * 1.0e-7; + assert_relative_eq!(result, expected, epsilon = 1e-15); + } + + #[test] + fn test_gvdw_exomol_mode() { + // EXOMOL mode (ivdwli > 0) + let params = GvdwParams { + ivdwli: 1, + id: 0, + il: 0, + gwm: 0.0, // not used in EXOMOL mode + vdwc: 0.0, // not used in EXOMOL mode + temp: 5800.0, + anhe: 0.085, // He abundance + anh2: 0.001, // H2 abundance + gvdwh2: 1.5e-7, + gexph2: 0.7, + gvdwhe: 0.5e-7, + gexphe: 0.5, + }; + + let result = gvdw(¶ms); + + // Manual calculation + let con: f64 = 4.1388e-12; + let t: f64 = 5800.0; + let temp_ratio: f64 = 296.0 / t; + let h2_term: f64 = temp_ratio.powf(0.7) * 1.5e-7 * 0.001; + let he_term: f64 = temp_ratio.powf(0.5) * 0.5e-7 * 0.085; + let expected: f64 = con * t * (h2_term + he_term); + + assert_relative_eq!(result, expected, epsilon = 1e-20); + } + + #[test] + fn test_gvdw_exomol_zero_abundance() { + // EXOMOL mode with zero H2 abundance (should only use He term) + let params = GvdwParams { + ivdwli: 1, + id: 0, + il: 0, + gwm: 0.0, + vdwc: 0.0, + temp: 10000.0, + anhe: 0.1, + anh2: 0.0, + gvdwh2: 1.0e-7, + gexph2: 0.5, + gvdwhe: 2.0e-7, + gexphe: 0.3, + }; + + let result = gvdw(¶ms); + + // Only He term contributes + let con: f64 = 4.1388e-12; + let t: f64 = 10000.0; + let temp_ratio: f64 = 296.0 / t; + let he_term: f64 = temp_ratio.powf(0.3) * 2.0e-7 * 0.1; + let expected: f64 = con * t * he_term; + + assert_relative_eq!(result, expected, epsilon = 1e-20); + } + + #[test] + fn test_gvdw_reference_temperature() { + // At reference temperature (296 K), temperature factor = 1 + let params = GvdwParams { + ivdwli: 1, + id: 0, + il: 0, + gwm: 0.0, + vdwc: 0.0, + temp: 296.0, + anhe: 0.1, + anh2: 0.1, + gvdwh2: 1.0e-7, + gexph2: 0.5, + gvdwhe: 1.0e-7, + gexphe: 0.5, + }; + + let result = gvdw(¶ms); + + // At T=296K, (296/T)^exp = 1 for any exponent + let con: f64 = 4.1388e-12; + let expected: f64 = con * 296.0 * (1.0e-7 * 0.1 + 1.0e-7 * 0.1); + + assert_relative_eq!(result, expected, epsilon = 1e-20); + } +} diff --git a/src/math/he2sew.rs b/src/math/he2sew.rs new file mode 100644 index 0000000..3a1f9a1 --- /dev/null +++ b/src/math/he2sew.rs @@ -0,0 +1,306 @@ +//! He II 线不透明度窗口初始化。 +//! +//! 重构自 SYNSPEC `he2sew.f` +//! +//! 设置 He II 线在频率窗口中的处理参数。 + +/// He II 电离阈值频率 (Hz)。 +/// +/// 对应 Fortran DATA FRHE 数组，是 He II Lyman 系列各线的阈值频率。 +/// FRHE(n) = R_inf * c / n², 其中 n = 1..12 +const FRHE: [f64; 12] = [ + 1.3158153e16, 3.2895381e15, 1.4624854e15, + 8.2261878e14, 5.2647201e14, 3.6560459e14, + 2.6860713e14, 2.0565220e14, 1.6249055e14, + 1.3161730e14, 1.0877460e14, 9.1400851e13, +]; + +/// 光速 (Å/s)，用于波长转换。 +const CLIGHT_A: f64 = 2.997925e17; + +/// He II 线窗口参数。 +/// +/// 包含 He II 线在给定频率点的处理参数。 +#[derive(Debug, Clone, Default)] +pub struct He2WindowParams { + /// He II 线处理标志 (-1: 排除, 1: 包含) + pub ihe2lw: i32, + /// He II 线系列索引 (1-12) + pub ilwhew: i32, + /// 主量子数上限 1 + pub mhe10w: i32, + /// 主量子数上限 2 + pub mhe20w: i32, +} + +/// 初始化 He II 线窗口参数。 +/// +/// 根据频率和表面重力确定 He II 线是否被包含在计算中， +/// 并设置相应的处理参数。 +/// +/// # 参数 +/// * `freq` - 频率 (Hz) +/// * `grav` - 表面重力 log g (cgs) +/// * `ifhe2` - He II 线处理标志 (<=0: 不处理) +/// +/// # 返回值 +/// 返回 He II 窗口参数结构体 +/// +/// # Fortran 源码 +/// ```fortran +/// SUBROUTINE HE2SEW(IJ) +/// ``` +pub fn he2sew(freq: f64, grav: f64, ifhe2: i32) -> He2WindowParams { + // 默认值：He II 线被排除 + let mut result = He2WindowParams { + ihe2lw: -1, + ilwhew: 0, + mhe10w: 60, + mhe20w: 60, + }; + + // 如果 He II 线处理标志 <= 0，直接返回 + if ifhe2 <= 0 { + return result; + } + + // 计算波长 (Å) + let al0 = CLIGHT_A / freq; + let al1 = al0; // 在原代码中 al0 和 al1 相同 + + // 根据表面重力确定哪些波长范围被排除 + // grav < 6: 低重力 (巨星) + // grav < 7: 中等重力 + // grav >= 7: 高重力 (白矮星) + if grav < 6.0 { + // 低重力恒星的排除范围 + if is_excluded_low_gravity(al0, al1) { + return result; + } + } else if grav < 7.0 { + // 中等重力恒星的排除范围 + if is_excluded_mid_gravity(al0, al1) { + return result; + } + } else { + // 高重力恒星的排除范围 + if is_excluded_high_gravity(al0, al1) { + return result; + } + } + + // He II 线被包含 + result.ihe2lw = 1; + result.mhe10w = 60; + result.mhe20w = 60; + + // 确定线系列索引 + result.ilwhew = determine_line_index(al0, al1); + + // 计算主量子数上限 + let idx = (result.ilwhew - 1) as usize; + let frion = FRHE[idx]; + let fr1 = frion * (result.ilwhew as f64) * (result.ilwhew as f64); + + if frion > freq { + let ratio = fr1 / (frion - freq); + result.mhe10w = (ratio.sqrt() as i32); + } + + result +} + +/// 检查低重力恒星 (grav < 6) 的排除范围。 +fn is_excluded_low_gravity(al0: f64, al1: f64) -> bool { + // 这些波长范围对应特定的 He II 线系 + // 条件格式: AL0 > lower .AND. AL1 < upper + const RANGES: [(f64, f64); 14] = [ + (31.0, 91.1), + (26.1, 29.8), + (24.8, 25.1), + (122.1, 162.9), + (165.1, 204.9), + (109.0, 120.9), + (103.0, 107.9), + (99.7, 102.0), + (320.8, 364.4), + (273.8, 319.8), + (251.6, 272.8), + (239.0, 250.6), + (231.1, 238.0), + (225.8, 230.1), + ]; + + for &(lower, upper) in &RANGES { + if al0 > lower && al1 < upper { + return true; + } + } + false +} + +/// 检查中等重力恒星 (6 <= grav < 7) 的排除范围。 +fn is_excluded_mid_gravity(al0: f64, al1: f64) -> bool { + const RANGES: [(f64, f64); 9] = [ + (33.0, 91.1), + (124.1, 160.9), + (167.1, 202.9), + (111.0, 118.9), + (322.8, 364.4), + (275.8, 317.8), + (253.6, 270.8), + (241.0, 248.6), + (233.1, 236.0), + ]; + + for &(lower, upper) in &RANGES { + if al0 > lower && al1 < upper { + return true; + } + } + false +} + +/// 检查高重力恒星 (grav >= 7) 的排除范围。 +fn is_excluded_high_gravity(al0: f64, al1: f64) -> bool { + const RANGES: [(f64, f64); 3] = [ + (39.0, 91.1), + (134.1, 150.9), + (177.1, 202.9), + ]; + + for &(lower, upper) in &RANGES { + if al0 > lower && al1 < upper { + return true; + } + } + false +} + +/// 根据波长确定线系列索引。 +/// +/// 索引 1-12 对应 He II 的不同线系。 +fn determine_line_index(al0: f64, _al1: f64) -> i32 { + // 使用 al0 进行判断 (原代码中 AL1.LT.91. 使用 al1，其他使用 al0) + // 但由于 al0 == al1，这里统一使用 al0 + if _al1 < 91.0 { + 1 + } else if al0 < 204.0 { + 2 + } else if al0 < 364.0 { + 3 + } else if al0 < 569.0 { + 4 + } else if al0 < 819.0 { + 5 + } else if al0 < 1116.0 { + 6 + } else if al0 < 1457.0 { + 7 + } else if al0 < 1844.0 { + 8 + } else if al0 < 2277.0 { + 9 + } else if al0 < 2756.0 { + 10 + } else if al0 < 3279.0 { + 11 + } else { + 12 + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_he2sew_disabled() { + // ifhe2 <= 0 时应返回默认值（排除） + let result = he2sew(1e15, 4.0, 0); + assert_eq!(result.ihe2lw, -1); + assert_eq!(result.ilwhew, 0); + } + + #[test] + fn test_he2sew_low_gravity_excluded() { + // 低重力恒星，波长在排除范围内 (31-91.1 Å) + // 频率 = c / lambda = 2.997925e17 / 50 = 5.99585e15 Hz + let freq = CLIGHT_A / 50.0; + let result = he2sew(freq, 4.0, 1); + assert_eq!(result.ihe2lw, -1); // 被排除 + } + + #[test] + fn test_he2sew_low_gravity_included() { + // 低重力恒星，波长不在排除范围内 + // 使用波长 95 Å，检查所有排除范围后确认不在其中 + let freq = CLIGHT_A / 95.0; + let result = he2sew(freq, 4.0, 1); + assert_eq!(result.ihe2lw, 1); // 被包含 + assert_eq!(result.ilwhew, 2); // al0=95, al1=95, 95 > 91 所以 ilwhew=2 + } + + #[test] + fn test_he2sew_mid_gravity_excluded() { + // 中等重力恒星 (6 <= grav < 7) + // 波长 50 Å 在排除范围 33-91.1 Å 内 + let freq = CLIGHT_A / 50.0; + let result = he2sew(freq, 6.5, 1); + assert_eq!(result.ihe2lw, -1); + } + + #[test] + fn test_he2sew_high_gravity_excluded() { + // 高重力恒星 (grav >= 7) + // 波长 50 Å 在排除范围 39-91.1 Å 内 + let freq = CLIGHT_A / 50.0; + let result = he2sew(freq, 7.5, 1); + assert_eq!(result.ihe2lw, -1); + } + + #[test] + fn test_he2sew_high_gravity_included() { + // 高重力恒星，波长 100 Å 不在排除范围内 + let freq = CLIGHT_A / 100.0; + let result = he2sew(freq, 7.5, 1); + assert_eq!(result.ihe2lw, 1); + // al0=100, al1=100: al1 > 91, al0 < 204, 所以 ilwhew=2 + assert_eq!(result.ilwhew, 2); + } + + #[test] + fn test_determine_line_index() { + // 测试线系列索引确定 + assert_eq!(determine_line_index(100.0, 90.0), 1); // al1 < 91 + assert_eq!(determine_line_index(150.0, 100.0), 2); // al0 < 204 + assert_eq!(determine_line_index(300.0, 100.0), 3); // al0 < 364 + assert_eq!(determine_line_index(500.0, 100.0), 4); // al0 < 569 + assert_eq!(determine_line_index(700.0, 100.0), 5); // al0 < 819 + assert_eq!(determine_line_index(1000.0, 100.0), 6); // al0 < 1116 + assert_eq!(determine_line_index(1300.0, 100.0), 7); // al0 < 1457 + assert_eq!(determine_line_index(1700.0, 100.0), 8); // al0 < 1844 + assert_eq!(determine_line_index(2100.0, 100.0), 9); // al0 < 2277 + assert_eq!(determine_line_index(2600.0, 100.0), 10); // al0 < 2756 + assert_eq!(determine_line_index(3100.0, 100.0), 11); // al0 < 3279 + assert_eq!(determine_line_index(4000.0, 100.0), 12); // >= 3279 + } + + #[test] + fn test_frhe_values() { + // 验证 FRHE 数组值 + assert_relative_eq!(FRHE[0], 1.3158153e16, epsilon = 1e10); + assert_relative_eq!(FRHE[11], 9.1400851e13, epsilon = 1e7); + } + + #[test] + fn test_mhe10w_calculation() { + // 当 frion > freq 时，mhe10w 应该被重新计算 + // 使用 FRHE[0] = 1.3158153e16，频率略低于此值 + let freq = 1.3e16; + let result = he2sew(freq, 7.5, 1); + // 验证 mhe10w 被正确设置 + assert!(result.mhe10w <= 60); + } +} diff --git a/src/math/heset.rs b/src/math/heset.rs new file mode 100644 index 0000000..e71fab7 --- /dev/null +++ b/src/math/heset.rs @@ -0,0 +1,476 @@ +//! 设置氦线参数。 +//! +//! 重构自 SYNSPEC `heset.f`。 +//! +//! 为 INISET 辅助程序，设置以下量： +//! - `iprf0`: He I 线标准吸收轮廓系数评估过程的索引（见 GAMHE） +//! - `ilwn`, `iupn`: 仅当 NLTE 选项开启时使用；与给定谱线相关的下能级和上能级索引 + +// ============================================================================ +// 物理常数 +// ============================================================================ + +/// 普朗克常数 (erg·s) +const H: f64 = 6.6256e-27; + +/// 光速 (cm/s) +const CL: f64 = 2.997925e10; + +// ============================================================================ +// He I 线数据表 +// ============================================================================ + +/// He I 线的 IT 数组 (0 或 1) +/// 用于确定能级选择逻辑 +const IT: [i32; 24] = [1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0]; + +/// He I 线的 NU 数组 (主量子数) +const NU: [i32; 24] = [6, 6, 9, 3, 8, 4, 7, 5, 6, 6, 5, 4, 4, 4, 3, 4, 3, 3, 5, 5, 7, 8, 10, 2]; + +/// He I 线的 JU 数组 (统计权重相关) +const JU: [i32; 24] = [ + 15, 3, 5, 9, 5, 3, 5, 3, 5, 1, 1, 15, 3, 5, 3, 1, 15, 5, 15, 5, 1, 1, 1, 9, +]; + +/// He I 线波长表 (Å) +/// 用于匹配输入波长以确定 IPRF 索引 +const HE1_WAVELENGTHS: [f64; 24] = [ + 3819.60, 3867.50, 3871.79, 3888.65, 3926.53, 3964.73, 4009.27, 4120.80, 4143.76, 4168.97, + 4437.55, 4471.50, 4713.20, 4921.93, 5015.68, 5047.74, 5875.70, 6678.15, 4026.20, 4387.93, + 4023.97, 3935.91, 3833.55, 10830.0, +]; + +/// He II 电离能 (cm⁻¹) +const HE2_IONIZATION_ENERGY: f64 = 438916.146; + +// ============================================================================ +// 参数结构体 +// ============================================================================ + +/// HESET 输入参数。 +#[derive(Debug, Clone)] +pub struct HesetInput<'a> { + /// 谱线索引 + pub il: i32, + /// 谱线波长 (nm) + pub alm: f64, + /// 下能级激发势 (cm⁻¹) + pub excl: f64, + /// 上能级激发势 (cm⁻¹) + pub excu: f64, + /// 电离度 (1=中性, 2=一次电离, ...) + pub ion: i32, + /// NLTE 模式标志 (>5 表示 LTE) + pub inlte: i32, + /// He I 元素索引 (0 表示未使用) + pub ielhe1: i32, + /// He II 元素索引 (0 表示未使用) + pub ielhe2: i32, + /// 离子起始能级索引数组 + pub nfirst: &'a [i32], + /// 离子终止能级索引数组 + pub nlast: &'a [i32], + /// 能级主量子数数组 + pub nquant: &'a [i32], + /// 能级电离能数组 (cm⁻¹) + pub enion: &'a [f64], + /// 能级统计权重数组 + pub g: &'a [f64], +} + +/// HESET 输出结果。 +#[derive(Debug, Clone, Default)] +pub struct HesetOutput { + /// He I 线标准吸收轮廓系数评估过程的索引 + pub iprf0: i32, + /// 下能级索引 (NLTE) + pub ilwn: i32, + /// 上能级索引 (NLTE) + pub iupn: i32, +} + +// ============================================================================ +// HESET 函数 +// ============================================================================ + +/// 设置氦线参数。 +/// +/// 根据 He I 或 He II 谱线的波长和激发势，确定： +/// - 标准吸收轮廓索引 (IPRF0) +/// - NLTE 能级索引 (ILWN, IUPN) +/// +/// # 参数 +/// +/// * `input` - 输入参数结构体 +/// * `iprf0_in` - 初始 IPRF0 值（仅当找到匹配时更新） +/// +/// # 返回 +/// +/// 包含 `iprf0`, `ilwn`, `iupn` 的输出结构体 +/// +/// # Fortran 原始代码 +/// +/// ```fortran +/// SUBROUTINE HESET(IL,ALM,EXCL,EXCU,ION,IPRF0,ILWN,IUPN) +/// ``` +pub fn heset(input: &HesetInput, iprf0_in: i32) -> HesetOutput { + let mut result = HesetOutput { + iprf0: iprf0_in, + ilwn: 0, + iupn: 0, + }; + + // ========== He I 处理 ========== + if input.ion != 1 { + // 跳转到 He II 处理 + return handle_he2(input, result); + } + + // 将波长从 nm 转换为 Å + let alam = input.alm * 10.0; + + // 查找匹配的 He I 线 + let mut iprf: i32 = 0; + for (i, &wavelength) in HE1_WAVELENGTHS.iter().enumerate() { + if (alam - wavelength).abs() < 1.0 { + iprf = (i + 1) as i32; // Fortran 1-indexed + break; + } + } + + // 更新 IPRF0（仅当 IPRF 在 1-20 范围内） + if iprf > 0 && iprf <= 20 { + result.iprf0 = iprf; + } + + // NLTE 能级索引处理 + if input.inlte > 5 || input.ielhe1 == 0 { + return result; + } + + let n0i = input.nfirst[input.ielhe1 as usize - 1]; // Fortran 1-indexed -> Rust 0-indexed + let n1i = input.nlast[input.ielhe1 as usize - 1]; + let hc = CL * H; + let eion = input.enion[n0i as usize - 1] / hc; + + let mut ilw: i32 = 0; + let mut iun: i32 = 0; + let mut nql: i32 = 0; + + if iprf > 0 { + nql = NU[(iprf - 1) as usize]; + } + + let mut igl: i32 = 0; + + // 遍历所有 He I 能级 + for i in n0i..=n1i { + let idx = (i - 1) as usize; // 转换为 Rust 0-indexed + let nq = input.nquant[idx]; + let ex = eion - input.enion[idx] / hc; + + // 匹配下能级 + if (input.excl - ex).abs() < 100.0 { + ilw = i; + igl = (input.g[idx] + 0.001) as i32; + } + + // 匹配上能级（根据量子数和统计权重） + if nq == nql && iprf > 0 { + let ig = (input.g[idx] + 0.001) as i32; + let iprf_idx = (iprf - 1) as usize; + + if IT[iprf_idx] == 0 { + // IT = 0 分支 + if nq == 2 && ig == JU[iprf_idx] { + iun = i; + } + if nq == 3 { + if ig == JU[iprf_idx] { + if ig == 1 || ig == 5 { + iun = i; + } + if ig == 3 && igl == 1 { + iun = i; + } + } else if ig == 9 { + iun = i; + } + } + if nq == 4 { + if ig == JU[iprf_idx] { + if ig == 1 || ig == 5 || ig == 7 { + iun = i; + } + if ig == 3 && igl == 1 { + iun = i; + } + } else if ig == 16 { + iun = i; + } + } + if ig == 25 || ig == 36 { + iun = i; + } + if ig == 49 || ig == 64 || ig == 81 { + iun = i; + } + if ig == 100 || ig == 121 || ig == 144 { + iun = i; + } + } else { + // IT = 1 分支 + if nq == 3 { + if ig == JU[iprf_idx] { + if ig == 9 || ig == 15 { + iun = i; + } + if ig == 3 && igl == 9 { + iun = i; + } + } else if ig == 27 { + iun = i; + } + } + if nq == 4 { + if ig == JU[iprf_idx] { + if ig == 9 || ig == 15 || ig == 21 { + iun = i; + } + if ig == 3 && igl == 9 { + iun = i; + } + } else if ig == 48 { + iun = i; + } + } + if ig == 75 { + iun = i; + } + if ig == 108 || ig == 147 || ig == 192 { + iun = i; + } + if ig == 243 || ig == 300 || ig == 363 { + iun = i; + } + } + + // 高能级处理 (n² = g) + if nq == 2 && ig == 16 { + iun = i; + } + if nq == 3 && ig == 36 { + iun = i; + } + if nq == 4 && ig == 64 { + iun = i; + } + if nq == 5 && ig == 100 { + iun = i; + } + if nq == 6 && ig == 144 { + iun = i; + } + if nq == 7 && ig == 196 { + iun = i; + } + if nq == 8 && ig == 256 { + iun = i; + } + if nq == 9 && ig == 324 { + iun = i; + } + if nq == 10 && ig == 400 { + iun = i; + } + } + } + + result.ilwn = ilw; + result.iupn = iun; + result +} + +/// 处理 He II 线。 +fn handle_he2(input: &HesetInput, mut result: HesetOutput) -> HesetOutput { + if input.ion != 2 || input.ielhe2 <= 0 { + return result; + } + + let n0i = input.nfirst[input.ielhe2 as usize - 1]; + let nlhe2 = input.nlast[input.ielhe2 as usize - 1] - n0i + 1; + + // 计算下能级主量子数 + let xl = (1.0 / (1.0 - input.excl / HE2_IONIZATION_ENERGY)).sqrt(); + let mut ilw = xl as i32; + if (ilw as f64) < xl { + ilw += 1; + } + + // 计算上能级主量子数 + let xu = (1.0 / (1.0 - input.excu / HE2_IONIZATION_ENERGY)).sqrt(); + let mut iun = xu as i32; + if (iun as f64) < xu { + iun += 1; + } + + // 设置输出（如果不超过能级数） + if ilw <= nlhe2 { + result.ilwn = ilw + n0i - 1; + } + if iun <= nlhe2 { + result.iupn = iun + n0i - 1; + } + + result +} + +// ============================================================================ +// 测试 +// ============================================================================ + +#[cfg(test)] +mod tests { + use super::*; + + /// 创建测试用的输入参数 + fn create_test_input() -> HesetInput<'static> { + // 模拟 He I 和 He II 的能级数据 + static NFIRST: [i32; 10] = [1, 5, 10, 0, 0, 0, 0, 0, 0, 0]; + static NLAST: [i32; 10] = [4, 9, 20, 0, 0, 0, 0, 0, 0, 0]; + static NQUANT: [i32; 25] = [ + 1, 2, 2, 2, 1, 2, 2, 3, 3, 3, 1, 2, 2, 3, 3, 3, 4, 4, 4, 4, 0, 0, 0, 0, 0, + ]; + static ENION: [f64; 25] = [ + 198310.0, 166240.0, 166240.0, 166240.0, 194470.0, 164070.0, 164070.0, 162970.0, + 162970.0, 162970.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + ]; + static G: [f64; 25] = [ + 1.0, 3.0, 1.0, 9.0, 2.0, 6.0, 2.0, 12.0, 4.0, 36.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + ]; + + HesetInput { + il: 1, + alm: 587.57, // He I 5876 Å 线 + excl: 159850.0, + excu: 211200.0, + ion: 1, + inlte: 0, // NLTE 模式 + ielhe1: 2, // He I 元素索引 + ielhe2: 3, // He II 元素索引 + nfirst: &NFIRST, + nlast: &NLAST, + nquant: &NQUANT, + enion: &ENION, + g: &G, + } + } + + #[test] + fn test_he1_wavelength_matching() { + // 测试 He I 5876 Å 线 (IPRF=17) + let input = HesetInput { + alm: 587.57, // nm -> 5875.7 Å + ion: 1, + ..create_test_input() + }; + + let result = heset(&input, 0); + assert_eq!(result.iprf0, 17, "应该匹配 He I 5876 Å 线 (IPRF=17)"); + } + + #[test] + fn test_he1_4471_line() { + // 测试 He I 4471 Å 线 (IPRF=12) + let input = HesetInput { + alm: 447.15, // nm -> 4471.5 Å + ion: 1, + ..create_test_input() + }; + + let result = heset(&input, 0); + assert_eq!(result.iprf0, 12, "应该匹配 He I 4471 Å 线 (IPRF=12)"); + } + + #[test] + fn test_he1_10830_line() { + // 测试 He I 10830 Å 线 (IPRF=24, 但 IPRF0 只在 1-20 范围内更新) + let input = HesetInput { + alm: 1083.0, // nm -> 10830 Å + ion: 1, + ..create_test_input() + }; + + let result = heset(&input, 0); + // IPRF=24 > 20, 所以 IPRF0 不更新 + assert_eq!(result.iprf0, 0, "IPRF=24 超出 1-20 范围，IPRF0 不应更新"); + } + + #[test] + fn test_he2_level_calculation() { + // 测试 He II 能级计算 + let input = HesetInput { + alm: 1640.0, // He II 1640 Å 线 + excl: 0.0, // n=1 的激发势为 0 + excu: 328150.0, // n=2 的激发势 (约为 3/4 * 438916) + ion: 2, + ielhe2: 3, + ..create_test_input() + }; + + let result = heset(&input, 0); + // He II 的能级索引计算 + assert!(result.ilwn > 0 || result.iupn > 0, "He II 应该计算出能级索引"); + } + + #[test] + fn test_non_helium_ion() { + // 测试非氦离子（不匹配任何 He I/He II 线） + let input = HesetInput { + alm: 500.0, + ion: 3, // 不是 He I 或 He II + ..create_test_input() + }; + + let result = heset(&input, 0); + assert_eq!(result.iprf0, 0); + assert_eq!(result.ilwn, 0); + assert_eq!(result.iupn, 0); + } + + #[test] + fn test_lte_mode() { + // 测试 LTE 模式（inlte > 5） + let input = HesetInput { + alm: 587.57, + ion: 1, + inlte: 6, // LTE 模式 + ..create_test_input() + }; + + let result = heset(&input, 0); + assert_eq!(result.iprf0, 17, "LTE 模式下 IPRF0 仍应正确设置"); + // LTE 模式下不计算 NLTE 能级索引 + } + + #[test] + fn test_constants() { + // 验证常数与 Fortran 一致 + use approx::assert_relative_eq; + + assert_relative_eq!(H, 6.6256e-27, epsilon = 1e-32); + assert_relative_eq!(CL, 2.997925e10, epsilon = 1e-5); + assert_relative_eq!(HE2_IONIZATION_ENERGY, 438916.146, epsilon = 1e-3); + } + + #[test] + fn test_data_arrays() { + // 验证数据数组长度 + assert_eq!(IT.len(), 24); + assert_eq!(NU.len(), 24); + assert_eq!(JU.len(), 24); + assert_eq!(HE1_WAVELENGTHS.len(), 24); + } +} diff --git a/src/math/hylset.rs b/src/math/hylset.rs new file mode 100644 index 0000000..fef861e --- /dev/null +++ b/src/math/hylset.rs @@ -0,0 +1,423 @@ +//! 氢线不透明度初始化过程。 +//! +//! 重构自 SYNSPEC `hylset.f`。 +//! +//! 设置氢线在频率窗口中的处理参数。 + +// ============================================================================ +// 物理常数 +// ============================================================================ + +/// 光速 (Å/s)，用于波长转换 +const CLIGHT_A: f64 = 2.997925e17; + +/// H I 电离频率 (Hz) - Lyman 系限 +const FR_H_LYMAN: f64 = 3.28805e15; + +/// Balmer 系限频率 (不同版本) +const FR_BALMER_VAC: f64 = 8.2225e14; // 真空波长 +const FR_BALMER_AIR: f64 = 8.22013e14; // 空气波长 + +/// 参考频率用于量子数计算 +const FR_REF_VAC: f64 = 3.289017e15; +const FR_REF_STD: f64 = 3.28805e15; + +// ============================================================================ +// Balmer 线波长表 (Å) +// ============================================================================ + +/// Balmer 系列线波长表 (Hα 到 Hο) +/// Hα = 656.28 nm, Hβ = 486.13 nm, ... +const BALMER_WAVELENGTHS: [f64; 15] = [ + 656.28, 486.13, 434.05, 410.17, 397.01, + 388.91, 383.54, 379.79, 377.06, 375.02, + 373.44, 372.19, 371.20, 370.39, 369.72, +]; + +// ============================================================================ +// 参数结构体 +// ============================================================================ + +/// HYLSET 输入参数。 +#[derive(Debug, Clone)] +pub struct HylsetParams { + /// 氢原子处理标志 (<= 0 表示不处理氢线) + pub iath: i32, + + /// 频率范围下限 (Hz) - FREQ(1) + pub freq1: f64, + + /// 频率范围上限 (Hz) - FREQ(2) + pub freq2: f64, + + /// 计算模式 (1: LTE, 2: NLTE 等) + pub imode: i32, + + /// 氢线轮廓处理标志 + pub ihydpr: i32, + + /// 表面重力 log g (cgs) + pub grav: f64, + + /// 真空波长限制 (Å) + pub vaclim: f64, +} + +/// HYLSET 输出结果。 +#[derive(Debug, Clone, Default)] +pub struct HylsetOutput { + /// 氢线处理标志 + /// - -1: 氢线被排除 + /// - 0: 待定 + /// - 1: 氢线被包含 + pub ihyl: i32, + + /// 氢线系列索引 + /// - 1: Lyman 系列 (波长 < 364.6 nm) + /// - 2: Balmer 系列 (364.6 <= 波长 < 820 nm) + /// - 3: Paschen 及更高系列 (波长 >= 820 nm) + pub ilowh: i32, + + /// 主量子数上限 1 (用于线强度计算) + pub m10: i32, + + /// 主量子数上限 2 (用于线强度计算) + pub m20: i32, +} + +// ============================================================================ +// HYLSET 函数 +// ============================================================================ + +/// 初始化氢线不透明度参数。 +/// +/// 根据频率范围和其他参数确定氢线是否被包含在计算中， +/// 并设置相应的处理参数。 +/// +/// # 参数 +/// +/// * `params` - 输入参数结构体 +/// +/// # 返回 +/// +/// 包含 `ihyl`, `ilowh`, `m10`, `m20` 的输出结构体 +/// +/// # 算法 +/// +/// 1. 默认 IHYL = -1 (氢线被排除) +/// 2. 如果 IATH <= 0，直接返回 +/// 3. 如果频率上限 >= H I 电离频率，直接返回 +/// 4. 计算波长范围并检查排除区域 +/// 5. 根据波长范围设置 ILOWH 和 M10/M20 +/// 6. 对于 Balmer 系列，检查是否有 Balmer 线在范围内 +/// +/// # Fortran 源码 +/// +/// ```fortran +/// SUBROUTINE HYLSET +/// ``` +pub fn hylset(params: &HylsetParams) -> HylsetOutput { + // 默认值：氢线被排除 + let mut result = HylsetOutput { + ihyl: -1, + ilowh: 0, + m10: 0, + m20: 40, + }; + + // 如果氢原子处理标志 <= 0，直接返回 + if params.iath <= 0 { + return result; + } + + // 如果频率上限 >= H I 电离频率，直接返回 + if params.freq2 >= FR_H_LYMAN { + return result; + } + + // 计算波长范围 (Å) + let al0 = CLIGHT_A / params.freq1; + let al1 = CLIGHT_A / params.freq2; + + // 检查是否在排除区域 + // 区域 1: 200 < AL0 且 AL1 < 364.6 nm (Lyman 系限附近) + if al0 > 200.0 && al1 < 364.6 { + return result; + } + // 区域 2: 560 < AL0 且 AL1 < 580 nm + if al0 > 560.0 && al1 < 580.0 { + return result; + } + // 区域 3: 720 < AL0 且 AL1 < 820.3 nm + if al0 > 720.0 && al1 < 820.3 { + return result; + } + + // 氢线状态：待定 + result.ihyl = 0; + result.m20 = 40; + + // 根据波长范围确定系列 + if al1 < 364.6 { + // Lyman 系列 + result.ilowh = 1; + let frion = FR_H_LYMAN; + result.m10 = ((FR_H_LYMAN / (frion - params.freq2).abs()).sqrt()) as i32; + if frion > params.freq1 { + result.m20 = ((FR_H_LYMAN / (frion - params.freq1)).sqrt()) as i32; + } + result.ihyl = 1; + + // 额外排除检查 + if al0 > 123.0 { + result.ihyl = 0; + } + if al0 > 104.0 && al1 < 120.0 { + result.ihyl = 0; + } + if al0 > 98.5 && al1 < 102.0 { + result.ihyl = 0; + } + + // NLTE 或高重力情况下强制包含 + if params.imode == 2 || params.ihydpr != 0 || params.grav >= 6.0 { + result.ihyl = 1; + } + } else if al1 < 820.0 { + // Balmer 系列 + result.ilowh = 2; + + let frion = if params.vaclim < 3600.0 { + FR_BALMER_VAC + } else { + FR_BALMER_AIR + }; + + result.m10 = if params.vaclim < 3600.0 { + ((FR_REF_VAC / (frion - params.freq2).abs()).sqrt()) as i32 + } else { + ((FR_REF_STD / (frion - params.freq2).abs()).sqrt()) as i32 + }; + + if frion > params.freq1 { + result.m20 = ((FR_REF_VAC / (frion - params.freq1)).sqrt()) as i32; + } + + // 检查是否有 Balmer 线在范围内 + for &al in &BALMER_WAVELENGTHS { + if al < al0 - 1.0 || al > al1 + 1.0 { + continue; + } + result.ihyl = 1; + break; + } + + // NLTE 或高重力情况下强制包含 + if params.imode == 2 || params.ihydpr != 0 || params.grav >= 6.0 { + result.ihyl = 1; + } + } else { + // Paschen 及更高系列 + result.ilowh = 3; + result.ihyl = 1; + } + + // 最终设置：总是包含氢线 + result.ihyl = 1; + + result +} + +// ============================================================================ +// 测试 +// ============================================================================ + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + /// 创建默认测试参数 + fn create_test_params() -> HylsetParams { + HylsetParams { + iath: 1, + freq1: 4.0e14, // 750 nm + freq2: 8.0e14, // 375 nm + imode: 1, + ihydpr: 0, + grav: 4.0, + vaclim: 2000.0, + } + } + + #[test] + fn test_hylset_disabled() { + // IATH <= 0 时应返回排除状态 + let params = HylsetParams { + iath: 0, + ..create_test_params() + }; + let result = hylset(¶ms); + assert_eq!(result.ihyl, -1); + } + + #[test] + fn test_hylset_freq_too_high() { + // 频率上限 >= H I 电离频率时应返回排除状态 + let params = HylsetParams { + freq2: 4.0e15, // > 3.28805e15 + ..create_test_params() + }; + let result = hylset(¶ms); + assert_eq!(result.ihyl, -1); + } + + #[test] + fn test_hylset_exclusion_zone_lyman() { + // 200 < AL0 且 AL1 < 364.6 nm (Lyman 系限附近) + let params = HylsetParams { + freq1: CLIGHT_A / 300.0, // AL0 = 300 nm + freq2: CLIGHT_A / 350.0, // AL1 = 350 nm + ..create_test_params() + }; + let result = hylset(¶ms); + assert_eq!(result.ihyl, -1); + } + + #[test] + fn test_hylset_exclusion_zone_560_580() { + // 560 < AL0 且 AL1 < 580 nm + let params = HylsetParams { + freq1: CLIGHT_A / 570.0, // AL0 = 570 nm + freq2: CLIGHT_A / 575.0, // AL1 = 575 nm + ..create_test_params() + }; + let result = hylset(¶ms); + assert_eq!(result.ihyl, -1); + } + + #[test] + fn test_hylset_exclusion_zone_720_820() { + // 720 < AL0 且 AL1 < 820.3 nm + let params = HylsetParams { + freq1: CLIGHT_A / 750.0, // AL0 = 750 nm + freq2: CLIGHT_A / 800.0, // AL1 = 800 nm + ..create_test_params() + }; + let result = hylset(¶ms); + assert_eq!(result.ihyl, -1); + } + + #[test] + fn test_hylset_lyman_series() { + // Lyman 系列 (AL1 < 364.6 nm) + let params = HylsetParams { + freq1: CLIGHT_A / 120.0, // AL0 = 120 nm + freq2: CLIGHT_A / 200.0, // AL1 = 200 nm + imode: 2, // NLTE 模式 + ..create_test_params() + }; + let result = hylset(¶ms); + assert_eq!(result.ilowh, 1); + assert_eq!(result.ihyl, 1); + } + + #[test] + fn test_hylset_balmer_series() { + // Balmer 系列 (364.6 <= AL1 < 820 nm) + // 包含 Hβ (486.13 nm) + // 注意：避免所有排除区域 + let params = HylsetParams { + freq1: CLIGHT_A / 550.0, // AL0 = 550 nm (不在 560-580 排除区域) + freq2: CLIGHT_A / 450.0, // AL1 = 450 nm + ..create_test_params() + }; + let result = hylset(¶ms); + assert_eq!(result.ilowh, 2); + assert_eq!(result.ihyl, 1); // 包含 Hβ + } + + #[test] + fn test_hylset_balmer_series_forced() { + // Balmer 系列，通过 NLTE 或高重力强制包含 + let params = HylsetParams { + freq1: CLIGHT_A / 500.0, // AL0 = 500 nm + freq2: CLIGHT_A / 400.0, // AL1 = 400 nm + grav: 7.0, // 高重力 + ..create_test_params() + }; + let result = hylset(¶ms); + assert_eq!(result.ilowh, 2); + assert_eq!(result.ihyl, 1); + } + + #[test] + fn test_hylset_paschen_series() { + // Paschen 及更高系列 (AL1 >= 820 nm) + let params = HylsetParams { + freq1: CLIGHT_A / 1000.0, // AL0 = 1000 nm + freq2: CLIGHT_A / 900.0, // AL1 = 900 nm + ..create_test_params() + }; + let result = hylset(¶ms); + assert_eq!(result.ilowh, 3); + assert_eq!(result.ihyl, 1); + } + + #[test] + fn test_hylset_m10_m20_calculation() { + // 测试 M10 和 M20 的计算 + let params = HylsetParams { + freq1: CLIGHT_A / 150.0, // AL0 = 150 nm + freq2: CLIGHT_A / 300.0, // AL1 = 300 nm + imode: 2, + ..create_test_params() + }; + let result = hylset(¶ms); + assert_eq!(result.ilowh, 1); + assert!(result.m10 > 0); + assert!(result.m20 > 0); + } + + #[test] + fn test_balmer_wavelengths() { + // 验证 Balmer 波长表 + assert_relative_eq!(BALMER_WAVELENGTHS[0], 656.28, epsilon = 0.01); // Hα + assert_relative_eq!(BALMER_WAVELENGTHS[1], 486.13, epsilon = 0.01); // Hβ + assert_relative_eq!(BALMER_WAVELENGTHS[2], 434.05, epsilon = 0.01); // Hγ + } + + #[test] + fn test_constants() { + // 验证常数与 Fortran 一致 + use approx::assert_relative_eq; + + assert_relative_eq!(CLIGHT_A, 2.997925e17, epsilon = 1e12); + assert_relative_eq!(FR_H_LYMAN, 3.28805e15, epsilon = 1e9); + } + + #[test] + fn test_vaclim_effect() { + // 测试 vaclim 对 Balmer 系列计算的影响 + let params_low_vaclim = HylsetParams { + freq1: CLIGHT_A / 500.0, + freq2: CLIGHT_A / 400.0, + vaclim: 3000.0, // < 3600 + ..create_test_params() + }; + + let params_high_vaclim = HylsetParams { + freq1: CLIGHT_A / 500.0, + freq2: CLIGHT_A / 400.0, + vaclim: 4000.0, // >= 3600 + ..create_test_params() + }; + + let result_low = hylset(¶ms_low_vaclim); + let result_high = hylset(¶ms_high_vaclim); + + // 两种情况下 M10 可能有细微差异 + assert_eq!(result_low.ilowh, 2); + assert_eq!(result_high.ilowh, 2); + } +} diff --git a/src/math/inibla.rs b/src/math/inibla.rs new file mode 100644 index 0000000..bd0a4d7 --- /dev/null +++ b/src/math/inibla.rs @@ -0,0 +1,422 @@ +//! 线列表初始化驱动程序。 +//! +//! 重构自 SYNSPEC `inibla.f` +//! +//! 处理当前波长范围的部分线列表，设置 Doppler 宽度和 Van der Waals 宽度等参数。 + +// ============================================================================ +// 物理常数 +// ============================================================================ + +/// 光速 (cm/s) +pub const CL: f64 = 2.997925e10; + +/// 普朗克常数 (erg·s) +pub const H: f64 = 6.6256e-27; + +/// 玻尔兹曼常数 (erg/K) +pub const BOLK: f64 = 1.38054e-16; + +/// Planck 函数常数 BN = 2*h*c²/c³ = 2*h/c² +pub const BN: f64 = 1.4743e-2; + +/// h/k (s·K) +pub const HK: f64 = 4.79928144e-11; + +/// Doppler 宽度常数 DP0 = 1/c +pub const DP0: f64 = 3.33564e-11; + +/// Doppler 宽度常数 DP1 = 2*k/m_H * 1e8 +pub const DP1: f64 = 1.651e8; + +/// Van der Waals 宽度常数 VW1 +pub const VW1: f64 = 0.42; + +/// Van der Waals 宽度常数 VW2 +pub const VW2: f64 = 0.45; + +/// 温度归一化因子 +pub const TENM4: f64 = 1e-4; + +// ============================================================================ +// 参数结构体 +// ============================================================================ + +/// INIBLA 输入参数。 +#[derive(Debug, Clone)] +pub struct IniblaParams<'a> { + /// 线数量 (如果为 0，直接返回) + pub nlin: i32, + + /// 频率数组 (Hz) + pub freq: &'a [f64], + + /// 频率点数 + pub nfreq: i32, + + /// 频率窗口标志 (>0 表示使用窗口) + pub ifwin: i32, + + /// 中心频率数组 (Hz) + pub freqc: &'a [f64], + + /// 中心频率点数 + pub nfreqc: i32, + + /// 深度点数 + pub nd: usize, + + /// 温度数组 (K) + pub temp: &'a [f64], + + /// 电子密度数组 + pub elec: &'a [f64], + + /// 氢密度数组 (中性氢) + pub ah: &'a [f64], + + /// 氦密度数组 + pub ahe: &'a [f64], + + /// H2 分子密度数组 + pub anh2: &'a [f64], + + /// 原子质量数组 + pub amas: &'a [f64], + + /// 湍流速度数组 + pub vturb: &'a [f64], + + /// 氢原子处理标志 (>0 表示处理氢) + pub iath: i32, +} + +/// INIBLA 输出结果。 +#[derive(Debug, Clone)] +pub struct IniblaOutput { + /// h*k*频率 / 温度 (用于激发) - EXHK + pub exhk: Vec, + + /// Planck 函数值 - PLAN + pub plan: Vec, + + /// 受激因子 (1 - e^(-hν/kT)) - STIM + pub stim: Vec, + + /// Van der Waals 宽度 - VDWC + pub vdwc: Vec, + + /// Doppler 宽度数组 (原子数 × 深度点数) - DOPA1 + pub dopa1: Vec, +} + +// ============================================================================ +// INIBLA 函数 +// ============================================================================ + +/// 初始化线列表参数。 +/// +/// 设置 Doppler 宽度、Van der Waals 宽度、Planck 函数等参数， +/// 用于后续的线不透明度计算。 +/// +/// # 参数 +/// +/// * `params` - 输入参数结构体 +/// +/// # 返回 +/// +/// 包含各种预计算值的输出结构体 +/// +/// # Fortran 原始代码 +/// +/// ```fortran +/// SUBROUTINE INIBLA +/// ``` +pub fn inibla(params: &IniblaParams) -> IniblaOutput { + let nd = params.nd; + let n_atoms = params.amas.len(); + + // 初始化输出数组 + let mut output = IniblaOutput { + exhk: vec![0.0; nd], + plan: vec![0.0; nd], + stim: vec![0.0; nd], + vdwc: vec![0.0; nd], + dopa1: vec![0.0; nd * n_atoms], + }; + + // 如果没有线，直接返回 + if params.nlin == 0 { + return output; + } + + // 计算参考频率 XX + let xx = if params.nfreq >= 2 { + 0.5 * (params.freq[0] + params.freq[1]) + } else { + params.freq[0] + }; + + // 如果使用频率窗口，使用窗口中心频率 + let mut xx = if params.ifwin > 0 && params.nfreqc > 0 { + 0.5 * (params.freqc[0] + params.freqc[params.nfreqc as usize - 1]) + } else { + xx + }; + + // 计算 Planck 函数常数 + let bnu = BN * (xx * 1e-15).powi(3); + let hkf = HK * xx; + + // 如果使用频率窗口，归一化 xx + if params.ifwin > 0 { + xx = 1.0; + } + + // 遍历所有深度点 + for id in 0..nd { + let t = params.temp[id]; + let _ane = params.elec[id]; + + // 计算激发因子 + let exh = (hkf / t).exp(); + output.exhk[id] = 1.0 / exh; + + // 计算 Planck 函数 + output.plan[id] = bnu / (exh - 1.0); + + // 计算受激因子 + output.stim[id] = 1.0 - output.exhk[id]; + + // 计算 Van der Waals 宽度 + let ah = params.ah[id]; + let ahe = params.ahe[id]; + let anh2 = params.anh2[id]; + output.vdwc[id] = (ah + VW1 * ahe + 0.85 * anh2) * (t * TENM4).powf(VW2); + + // 计算每个原子的 Doppler 宽度 + for (iat, &amas) in params.amas.iter().enumerate() { + if amas > 0.0 { + let vturb_sq = params.vturb[id] * params.vturb[id]; + let dopa1_val = 1.0 / (xx * DP0 * (DP1 * t / amas + vturb_sq).sqrt()); + output.dopa1[iat * nd + id] = dopa1_val; + } + } + } + + output +} + +// ============================================================================ +// 辅助函数 +// ============================================================================ + +/// 计算单个深度点的 Doppler 宽度。 +/// +/// # 参数 +/// +/// * `xx` - 参考频率 (Hz) +/// * `temp` - 温度 (K) +/// * `amas` - 原子质量 (原子质量单位) +/// * `vturb` - 湍流速度 (cm/s) +/// +/// # 返回 +/// +/// Doppler 宽度参数 +pub fn compute_doppler_width(xx: f64, temp: f64, amas: f64, vturb: f64) -> f64 { + 1.0 / (xx * DP0 * (DP1 * temp / amas + vturb * vturb).sqrt()) +} + +/// 计算单个深度点的 Van der Waals 宽度。 +/// +/// # 参数 +/// +/// * `ah` - 中性氢密度 +/// * `ahe` - 氦密度 +/// * `anh2` - H2 分子密度 +/// * `temp` - 温度 (K) +/// +/// # 返回 +/// +/// Van der Waals 宽度参数 +pub fn compute_vdw_width(ah: f64, ahe: f64, anh2: f64, temp: f64) -> f64 { + (ah + VW1 * ahe + 0.85 * anh2) * (temp * TENM4).powf(VW2) +} + +/// 计算 Planck 函数值。 +/// +/// # 参数 +/// +/// * `freq` - 频率 (Hz) +/// * `temp` - 温度 (K) +/// +/// # 返回 +/// +/// Planck 函数值 (erg/s/cm²/Hz/sr) +pub fn compute_planck(freq: f64, temp: f64) -> f64 { + let hkf = HK * freq; + let exh = (hkf / temp).exp(); + let bnu = BN * (freq * 1e-15).powi(3); + bnu / (exh - 1.0) +} + +// ============================================================================ +// 测试 +// ============================================================================ + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_inibla_no_lines() { + // 当 nlin = 0 时应返回零数组 + let freq = [1e15]; + let temp = [10000.0]; + let elec = [1e12]; + let ah = [1e10]; + let ahe = [1e9]; + let anh2 = [1e8]; + let amas = [1.0, 4.0]; + let vturb = [1e5]; + + let params = IniblaParams { + nlin: 0, + freq: &freq, + nfreq: 1, + ifwin: 0, + freqc: &[], + nfreqc: 0, + nd: 1, + temp: &temp, + elec: &elec, + ah: &ah, + ahe: &ahe, + anh2: &anh2, + amas: &amas, + vturb: &vturb, + iath: 1, + }; + + let result = inibla(¶ms); + assert_eq!(result.exhk.len(), 1); + assert_eq!(result.plan.len(), 1); + assert_eq!(result.stim.len(), 1); + assert_eq!(result.vdwc.len(), 1); + } + + #[test] + fn test_inibla_with_lines() { + // 当 nlin > 0 时应计算正确的值 + let freq = [1e15, 1.1e15]; + let temp = [10000.0, 9000.0]; + let elec = [1e12, 9e11]; + let ah = [1e10, 9e9]; + let ahe = [1e9, 9e8]; + let anh2 = [1e8, 9e7]; + let amas = [1.0, 4.0]; + let vturb = [1e5, 9e4]; + + let params = IniblaParams { + nlin: 100, + freq: &freq, + nfreq: 2, + ifwin: 0, + freqc: &[], + nfreqc: 0, + nd: 2, + temp: &temp, + elec: &elec, + ah: &ah, + ahe: &ahe, + anh2: &anh2, + amas: &amas, + vturb: &vturb, + iath: 1, + }; + + let result = inibla(¶ms); + + // 验证输出数组长度 + assert_eq!(result.exhk.len(), 2); + assert_eq!(result.plan.len(), 2); + assert_eq!(result.stim.len(), 2); + assert_eq!(result.vdwc.len(), 2); + assert_eq!(result.dopa1.len(), 4); // 2 atoms * 2 depth points + + // 验证所有值都是正数 + for &v in &result.exhk { + assert!(v > 0.0 && v < 1.0); + } + for &v in &result.plan { + assert!(v > 0.0); + } + for &v in &result.stim { + assert!(v > 0.0 && v < 1.0); + } + for &v in &result.vdwc { + assert!(v > 0.0); + } + } + + #[test] + fn test_compute_doppler_width() { + let xx = 1e15; // Hz + let temp = 10000.0; // K + let amas = 1.0; // 氢原子质量 + let vturb = 1e5; // cm/s + + let result = compute_doppler_width(xx, temp, amas, vturb); + + // 验证结果为正数 + assert!(result > 0.0); + + // 手动计算 + let expected = 1.0 / (xx * DP0 * (DP1 * temp / amas + vturb * vturb).sqrt()); + assert_relative_eq!(result, expected, epsilon = 1e-20); + } + + #[test] + fn test_compute_vdw_width() { + let ah = 1e10; + let ahe = 1e9; + let anh2 = 1e8; + let temp = 10000.0; + + let result = compute_vdw_width(ah, ahe, anh2, temp); + + // 验证结果为正数 + assert!(result > 0.0); + + // 手动计算 + let expected = (ah + VW1 * ahe + 0.85 * anh2) * (temp * TENM4).powf(VW2); + assert_relative_eq!(result, expected, epsilon = 1e-20); + } + + #[test] + fn test_compute_planck() { + let freq = 1e15; // Hz + let temp = 10000.0; // K + + let result = compute_planck(freq, temp); + + // 验证结果为正数 + assert!(result > 0.0); + + // 验证 Planck 函数值在合理范围内 + // 对于 T=10000K, ν=1e15 Hz，Planck 函数约为 1e-5 erg/s/cm²/Hz/sr + assert!(result > 1e-10 && result < 1e10); + } + + #[test] + fn test_constants() { + // 验证常数与 Fortran 一致 + assert_relative_eq!(DP0, 3.33564e-11, epsilon = 1e-16); + assert_relative_eq!(DP1, 1.651e8, epsilon = 1e3); + assert_relative_eq!(VW1, 0.42, epsilon = 1e-10); + assert_relative_eq!(VW2, 0.45, epsilon = 1e-10); + assert_relative_eq!(BN, 1.4743e-2, epsilon = 1e-6); + assert_relative_eq!(HK, 4.79928144e-11, epsilon = 1e-18); + } +} diff --git a/src/math/iniblm.rs b/src/math/iniblm.rs new file mode 100644 index 0000000..285f2e7 --- /dev/null +++ b/src/math/iniblm.rs @@ -0,0 +1,292 @@ +//! 分子线列表初始化驱动程序。 +//! +//! 重构自 SYNSPEC `iniblm.f` +//! +//! 处理当前波长范围的部分分子线列表，设置 Doppler 宽度等参数。 + +// ============================================================================ +// 物理常数 +// ============================================================================ + +/// 光速 (cm/s) +pub const CL: f64 = 2.997925e10; + +/// 普朗克常数 (erg·s) +pub const H: f64 = 6.6256e-27; + +/// 玻尔兹曼常数 (erg/K) +pub const BOLK: f64 = 1.38054e-16; + +/// Planck 函数常数 BN = 2*h*c²/c³ = 2*h/c² +pub const BN_MOL: f64 = 1.4743e-2; + +/// h/k (s·K) +pub const HK_MOL: f64 = 4.79928144e-11; + +/// Doppler 宽度常数 DP0 = 1/c +pub const DP0_MOL: f64 = 3.33564e-11; + +/// Doppler 宽度常数 DP1 = 2*k/m_H * 1e8 +pub const DP1_MOL: f64 = 1.651e8; + +// ============================================================================ +// 参数结构体 +// ============================================================================ + +/// INIBLM 输入参数。 +#[derive(Debug, Clone)] +pub struct IniblmParams<'a> { + /// 分子线数量 (如果为 0，直接返回) + pub nmol: i32, + + /// 频率数组 (Hz) + pub freq: &'a [f64], + + /// 频率点数 + pub nfreq: i32, + + /// 深度点数 + pub nd: usize, + + /// 温度数组 (K) + pub temp: &'a [f64], + + /// 分子质量数组 + pub ammol: &'a [f64], + + /// 湍流速度数组 + pub vturb: &'a [f64], +} + +/// INIBLM 输出结果。 +#[derive(Debug, Clone)] +pub struct IniblmOutput { + /// h*k*频率 / 温度 (用于激发) - EXHK + pub exhk: Vec, + + /// Planck 函数值 - PLAN + pub plan: Vec, + + /// 受激因子 (1 - e^(-hν/kT)) - STIM + pub stim: Vec, + + /// Doppler 宽度数组 (分子数 × 深度点数) - DOPMOL + pub dopmol: Vec, +} + +// ============================================================================ +// INIBLM 函数 +// ============================================================================ + +/// 初始化分子线列表参数。 +/// +/// 设置 Doppler 宽度、Planck 函数等参数， +/// 用于后续的分子线不透明度计算。 +/// +/// # 参数 +/// +/// * `params` - 输入参数结构体 +/// +/// # 返回 +/// +/// 包含各种预计算值的输出结构体 +/// +/// # Fortran 原始代码 +/// +/// ```fortran +/// SUBROUTINE INIBLM +/// ``` +pub fn iniblm(params: &IniblmParams) -> IniblmOutput { + let nd = params.nd; + let nmol = params.nmol as usize; + + // 初始化输出数组 + let mut output = IniblmOutput { + exhk: vec![0.0; nd], + plan: vec![0.0; nd], + stim: vec![0.0; nd], + dopmol: vec![0.0; nd * nmol], + }; + + // 如果没有分子线，直接返回 + if params.nmol == 0 { + return output; + } + + // 计算参考频率 + let xx = if params.nfreq >= 2 { + 0.5 * (params.freq[0] + params.freq[1]) + } else { + params.freq[0] + }; + + // 计算 Planck 函数常数 + let bnu = BN_MOL * (xx * 1e-15).powi(3); + let hkf = HK_MOL * xx; + + // 遍历所有深度点 + for id in 0..nd { + let t = params.temp[id]; + + // 计算激发因子 + let exh = (hkf / t).exp(); + output.exhk[id] = 1.0 / exh; + + // 计算 Planck 函数 + output.plan[id] = bnu / (exh - 1.0); + + // 计算受激因子 + output.stim[id] = 1.0 - output.exhk[id]; + + // 计算每个分子的 Doppler 宽度 + for (imol, &ammol) in params.ammol.iter().enumerate() { + if ammol > 0.0 { + let vturb_sq = params.vturb[id] * params.vturb[id]; + let dopmol_val = 1.0 / (xx * DP0_MOL * (DP1_MOL * t / ammol + vturb_sq).sqrt()); + output.dopmol[imol * nd + id] = dopmol_val; + } + } + } + + output +} + +// ============================================================================ +// 辅助函数 +// ============================================================================ + +/// 计算单个深度点的分子 Doppler 宽度。 +/// +/// # 参数 +/// +/// * `xx` - 参考频率 (Hz) +/// * `temp` - 温度 (K) +/// * `ammol` - 分子质量 (原子质量单位) +/// * `vturb` - 湍流速度 (cm/s) +/// +/// # 返回 +/// +/// Doppler 宽度参数 +pub fn compute_molecular_doppler_width(xx: f64, temp: f64, ammol: f64, vturb: f64) -> f64 { + 1.0 / (xx * DP0_MOL * (DP1_MOL * temp / ammol + vturb * vturb).sqrt()) +} + +/// 计算分子 Planck 函数值。 +/// +/// # 参数 +/// +/// * `freq` - 频率 (Hz) +/// * `temp` - 温度 (K) +/// +/// # 返回 +/// +/// Planck 函数值 (erg/s/cm²/Hz/sr) +pub fn compute_molecular_planck(freq: f64, temp: f64) -> f64 { + let hkf = HK_MOL * freq; + let exh = (hkf / temp).exp(); + let bnu = BN_MOL * (freq * 1e-15).powi(3); + bnu / (exh - 1.0) +} + +// ============================================================================ +// 测试 +// ============================================================================ + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_iniblm_no_molecules() { + // 当 nmol = 0 时应返回零数组 + let freq = [1e15]; + let temp = [10000.0]; + let ammol = [16.0, 28.0]; // H2, CO + let vturb = [1e5]; + + let params = IniblmParams { + nmol: 0, + freq: &freq, + nfreq: 1, + nd: 1, + temp: &temp, + ammol: &ammol, + vturb: &vturb, + }; + + let result = iniblm(¶ms); + assert_eq!(result.exhk.len(), 1); + assert_eq!(result.plan.len(), 1); + assert_eq!(result.stim.len(), 1); + } + #[test] + fn test_iniblm_with_molecules() { + // 当 nmol > 0 时应计算正确的值 + let freq = [1e15, 1.1e15]; + let temp = [10000.0, 9000.0]; + let ammol = [16.0, 28.0]; // H2, CO + let vturb = [1e5, 9e4]; + + let params = IniblmParams { + nmol: 2, + freq: &freq, + nfreq: 2, + nd: 2, + temp: &temp, + ammol: &ammol, + vturb: &vturb, + }; + + let result = iniblm(¶ms); + + // 验证输出数组长度 + assert_eq!(result.exhk.len(), 2); + assert_eq!(result.plan.len(), 2); + assert_eq!(result.stim.len(), 2); + assert_eq!(result.dopmol.len(), 4); // 2 molecules * 2 depth points + + // 验证所有值都是正数 + for &v in &result.exhk { + assert!(v > 0.0 && v < 1.0); + } + for &v in &result.plan { + assert!(v > 0.0); + } + for &v in &result.stim { + assert!(v > 0.0 && v < 1.0); + } + } + #[test] + fn test_compute_molecular_doppler_width() { + let xx = 1e15; // Hz + let temp = 10000.0; // K + let ammol = 16.0; // H2 分子质量 + let vturb = 1e5; // cm/s + + let result = compute_molecular_doppler_width(xx, temp, ammol, vturb); + // 验证结果为正数 + assert!(result > 0.0); + // 手动计算 + let expected = 1.0 / (xx * DP0_MOL * (DP1_MOL * temp / ammol + vturb * vturb).sqrt()); + assert_relative_eq!(result, expected, epsilon = 1e-20); + } + #[test] + fn test_compute_molecular_planck() { + let freq = 1e15; // Hz + let temp = 10000.0; // K + let result = compute_molecular_planck(freq, temp); + // 验证结果为正数 + assert!(result > 0.0); + // 验证 Planck 函数值在合理范围内 + assert!(result > 1e-10 && result < 1e10); + } + #[test] + fn test_constants() { + // 验证常数与 Fortran 一致 + assert_relative_eq!(DP0_MOL, 3.33564e-11, epsilon = 1e-16); + assert_relative_eq!(DP1_MOL, 1.651e8, epsilon = 1e3); + assert_relative_eq!(BN_MOL, 1.4743e-2, epsilon = 1e-6); + assert_relative_eq!(HK_MOL, 4.79928144e-11, epsilon = 1e-18); + } +} diff --git a/src/math/mod.rs b/src/math/mod.rs index 09011d4..985f520 100644 --- a/src/math/mod.rs +++ b/src/math/mod.rs @@ -96,16 +96,22 @@ mod gridp; mod griem; mod gomini; mod grcor; +mod gvdw; +mod he2sew; mod greyd; mod h2minus; mod hephot; mod hedif; +mod heset; mod hesol6; mod hesolv; mod hidalg; pub mod indexx; mod ijali2; mod inifrs; +mod hylset; +mod inibla; +mod iniblm; mod inilam; mod inpdis; mod ijalis; @@ -142,6 +148,7 @@ mod meanopt; mod minv3; mod mpartf; mod moleq; +mod molop; mod newdm; mod newdmt; mod newpop; @@ -187,6 +194,9 @@ mod pfni; mod pzert; mod pzevld; mod pfspec; +mod phe2; +mod phtion; +mod phtx; mod pgset; mod psolve; pub mod quit; @@ -239,6 +249,7 @@ mod setdrt; mod sghe12; mod sgmer; mod sgmer1; +mod sgmerg; mod sigave; mod sigk; mod sigmar; @@ -393,19 +404,25 @@ pub use gntk::gntk; pub use gridp::gridp; pub use griem::{griem, GriemParams}; pub use gomini::{gomini, GominiParams, GominiResult}; +pub use gvdw::{gvdw, GvdwParams}; pub use grcor::grcor; pub use greyd::{ greyd_pure, format_greyd_iter, GreydConfig, GreydState, GreydOutput, GreydIterOutput, }; pub use h2minus::h2minus; +pub use he2sew::{he2sew, He2WindowParams}; pub use hephot::hephot; pub use hedif::{hedif, hedif_io, HedifParams, HedifResult}; +pub use heset::{heset, HesetInput, HesetOutput}; +pub use hylset::{hylset, HylsetParams, HylsetOutput}; pub use hesol6::{hesol6, Hesol6Aux, Hesol6Output, Hesol6Params}; pub use hesolv::{hesolv_pure, HesolvAux, HesolvConfig, HesolvModelState, HesolvAtomicParams, HesolvParams, HesolvOutput}; pub use hidalg::hidalg; pub use indexx::indexx; pub use inifrs::{inifrs, InifrsConfig, InifrsFreqControl, InifrsOutput}; +pub use inibla::{inibla, IniblaParams, IniblaOutput, compute_doppler_width, compute_vdw_width, compute_planck}; +pub use iniblm::{iniblm, IniblmParams, IniblmOutput, compute_molecular_doppler_width, compute_molecular_planck}; pub use inilam::{inilam_pure, InilamConfig, InilamModelState, InilamAtomicParams, InilamFreqParams, InilamOutput}; pub use inpdis::{inpdis_pure, inpdis_io, InpDisParams, InpDisResult}; pub use ijalis::{ijalis, ijalis_io, IjalisParams, IjalisOutput}; @@ -444,6 +461,7 @@ pub use meanopt::{meanopt, MeanoptModelState, MeanoptOutput, MeanoptParams}; pub use minv3::minv3; pub use mpartf::{mpartf, MpartfResult}; pub use moleq::{moleq_pure, MoleqParams, MoleqOutput, MoleculeEqData, parse_molecule_data}; +pub use molop::{molop_pure, MolopConfig, MolopModelState, MolopFreqParams, MolLineData, MolModelState, MolopOutput}; pub use newdmt::{newdmt_pure, NewdmtConfig, NewdmtModelState, NewdmtPrsAux, NewdmtFactrs}; pub use newpop::{newpop, NewpopParams, NewpopResult}; pub use osccor::{osccor, OsccorParams, OsccorOutput, format_oscillation_message}; @@ -495,6 +513,9 @@ pub use pfni::pfni; pub use pzert::pzert; pub use pzevld::pzevld; pub use pfspec::pfspec; +pub use phe2::{phe2, Phe2Params, Phe2Output}; +pub use phtion::{phtion, phtion_pure, PhtionParams, PhtionOutput, PhotcsData, MFRQ}; +pub use phtx::{phtx, PhtxParams, PhtxOutput, PhtxState, LevelPhotoData, HhePhotoData}; pub use pgset::{pgset, PgsetParams, PgsetOutput, MDEPTH}; pub use psolve::psolve; pub use quartc::quartc; @@ -561,6 +582,7 @@ pub use sbfoh::sbfoh; pub use setdrt::setdrt; pub use sghe12::sghe12; pub use sgmer::{sgmer0, sgmer1, sgmerd}; +pub use sgmerg::{sgmerg, sgmerg_pure, SgmergParams}; pub use sigmar::sigmar; pub use sigave::{sigave_from_data, sigave_pure, SigaveParams, SigaveOutput}; pub use sigk::{sigk, SigkParams}; diff --git a/src/math/molop.rs b/src/math/molop.rs new file mode 100644 index 0000000..4c84221 --- /dev/null +++ b/src/math/molop.rs @@ -0,0 +1,468 @@ +//! 分子谱线不透明度和发射率计算。 +//! +//! 重构自 SYNSPEC `molop.f` +//! +//! 计算给定深度点的总分子谱线不透明度 (ABLIN) 和发射率 (EMLIN)。 + +use super::voigtk::{voigtk, MVOI}; + +/// 常量 +const UN: f64 = 1.0; +const EXT0: f64 = 3.17; +const TEN: f64 = 10.0; + +/// 分子谱线数据参数 +/// +/// 对应 LINDAT.FOR 中的 MOLTOT COMMON 块 +pub struct MolLineData<'a> { + /// 分子谱线频率 [Hz] + pub freqm: &'a [f64], + /// 下能级激发能 [K] + pub exclm: &'a [f64], + /// gf 值 (log10) + pub gfm: &'a [f64], + /// 辐射加宽参数 + pub grm: &'a [f64], + /// Stark 加宽参数 + pub gsm: &'a [f64], + /// van der Waals 加宽参数 (按深度) + pub gvdw: &'a [f64], + /// 分子数据索引 + pub indatm: &'a [i32], + /// 分子线索引 + pub inmlin: &'a [i32], + /// 中心频率索引 + pub ijcmtr: &'a [i32], + /// 分子谱线数 + pub nlinml: i32, + /// 分子列表激活标志 + pub inactm: i32, +} + +/// 分子状态参数 +/// +/// 对应 MODELP.FOR 中的 MOLPAR COMMON 块 +pub struct MolModelState<'a> { + /// 分子数密度 [cm^-3] + pub rrmol: &'a [f64], + /// 分子多普勒宽度 + pub dopmol: &'a [f64], +} + +/// 深度点状态参数 +/// +/// 对应 MODELP.FOR 和 LINDAT.FOR 中的相关变量 +pub struct MolopModelState<'a> { + /// 温度 [K] + pub temp: f64, + /// 电子密度 [cm^-3] + pub elec: f64, + /// 分子状态 + pub mol: MolModelState<'a>, + /// Planck 函数 + pub plan: f64, + /// 受激因子 + pub stim: f64, +} + +/// 频率网格参数 +/// +/// 对应 SYNTHP.FOR 中的 FREQSY COMMON 块 +pub struct MolopFreqParams<'a> { + /// 频率数组 [Hz] + pub freq: &'a [f64], + /// 频率点数 + pub nfreq: usize, + /// 参考频率点数（用于输出范围） + pub nfreqs: usize, +} + +/// 配置参数 +/// +/// 对应 PARAMS.FOR 和 LINDAT.FOR 中的相关变量 +pub struct MolopConfig { + /// 分子温度限制 [K] + pub tmolim: f64, + /// 频率转换因子 + pub dfrcon: f64, +} + +/// 输出结果 +pub struct MolopOutput { + /// 不透明度数组 [cm^-1] + pub ablin: Vec, + /// 发射率数组 + pub emlin: Vec, +} + +/// 分子谱线不透明度和发射率计算（纯函数版本）。 +/// +/// # 参数 +/// +/// * `config` - 配置参数 +/// * `model` - 模型状态 +/// * `line_data` - 分子谱线数据 +/// * `freq` - 频率网格参数 +/// * `avab` - 平均不透明度（用于确定线宽范围） +/// * `voigt_tables` - Voigt 函数预计算表格 (h0, h1, h2) +/// +/// # 返回值 +/// +/// 返回不透明度和发射率数组 +pub fn molop_pure( + config: &MolopConfig, + model: &MolopModelState, + line_data: &MolLineData, + freq: &MolopFreqParams, + avab: f64, + voigt_tables: (&[f64; MVOI], &[f64; MVOI], &[f64; MVOI]), +) -> MolopOutput { + let nfreq = freq.nfreq; + let mut ablin = vec![0.0; nfreq]; + let mut emlin = vec![0.0; nfreq]; + + // 检查温度是否超过分子温度限制 + if model.temp > config.tmolim { + return MolopOutput { ablin, emlin }; + } + + // 检查是否有分子谱线 + if line_data.nlinml == 0 { + return MolopOutput { ablin, emlin }; + } + + // 检查分子列表是否激活 + if line_data.inactm != 0 { + return MolopOutput { ablin, emlin }; + } + + let (h0tab, h1tab, h2tab) = voigt_tables; + let tem1 = UN / model.temp; + let ane = model.elec; + + // 遍历所有贡献的分子谱线 + for i in 0..line_data.nlinml as usize { + let il = line_data.inmlin[i] as usize; + let imol = line_data.indatm[il] as usize; + + let dop1 = model.mol.dopmol[imol]; + let agam = (line_data.grm[il] + line_data.gsm[il] * ane + line_data.gvdw[il]) * dop1; + + let fr0 = line_data.freqm[il]; + let ab0 = (line_data.gfm[il] - line_data.exclm[il] * tem1).exp() + * model.mol.rrmol[imol] + * dop1 + * model.stim; + + // 设置谱线贡献的频率范围限制 + let ex0 = ab0 / avab * agam; + let mut ext = EXT0; + if ex0 > TEN { + ext = ex0.sqrt(); + } + ext = ext / dop1; + let xijext = config.dfrcon * ext + 1.5; + + let ij1 = ((line_data.ijcmtr[i] as f64 - xijext).max(3.0)) as usize; + let ij2 = ((line_data.ijcmtr[i] as f64 + xijext).min(freq.nfreqs as f64)) as usize; + + // 只有当范围有效时才计算 + if ij1 < nfreq && ij2 > 2 { + for ij in ij1..=ij2.min(nfreq - 1) { + let xf = (freq.freq[ij] - fr0).abs() * dop1; + ablin[ij] = ablin[ij] + ab0 * voigtk(agam, xf, h0tab, h1tab, h2tab); + } + } + } + + // 计算发射率（从第3个频率点开始） + for ij in 3..nfreq { + emlin[ij] = emlin[ij] + ablin[ij] * model.plan; + } + + MolopOutput { ablin, emlin } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + fn create_test_voigt_tables() -> ([f64; MVOI], [f64; MVOI], [f64; MVOI]) { + let mut h0 = [0.0; MVOI]; + let mut h1 = [0.0; MVOI]; + let mut h2 = [0.0; MVOI]; + + // H(0,v) = exp(-v^2) / sqrt(pi) + let inv_sqrt_pi = 0.5641895835477563; + for i in 0..MVOI { + let v = (i as f64 - 1.5) / 200.0; + h0[i] = (-v * v).exp() * inv_sqrt_pi; + h1[i] = -2.0 * v * h0[i]; + h2[i] = (4.0 * v * v - 2.0) * h0[i]; + } + + (h0, h1, h2) + } + + #[test] + fn test_molop_high_temperature() { + let (h0, h1, h2) = create_test_voigt_tables(); + + let config = MolopConfig { + tmolim: 5000.0, + dfrcon: 1.0, + }; + + let mol_data = vec![1e10, 1e10, 1e10]; + let dop_data = vec![1e-10, 1e-10, 1e-10]; + + let model = MolopModelState { + temp: 6000.0, // 超过 tmolim + elec: 1e12, + mol: MolModelState { + rrmol: &mol_data, + dopmol: &dop_data, + }, + plan: 1.0, + stim: 1.0, + }; + + let freq_data = vec![0.0; 100]; + let freq = MolopFreqParams { + freq: &freq_data, + nfreq: 100, + nfreqs: 100, + }; + + let gfm = vec![]; + let exclm = vec![]; + let freqm = vec![]; + let grm = vec![]; + let gsm = vec![]; + let gvdw = vec![]; + let indatm = vec![]; + let inmlin = vec![]; + let ijcmtr = vec![]; + + let line_data = MolLineData { + freqm: &freqm, + exclm: &exclm, + gfm: &gfm, + grm: &grm, + gsm: &gsm, + gvdw: &gvdw, + indatm: &indatm, + inmlin: &inmlin, + ijcmtr: &ijcmtr, + nlinml: 10, + inactm: 0, + }; + + let result = molop_pure(&config, &model, &line_data, &freq, 1e-10, (&h0, &h1, &h2)); + + // 高温时应该返回零 + for i in 0..result.ablin.len() { + assert_relative_eq!(result.ablin[i], 0.0); + assert_relative_eq!(result.emlin[i], 0.0); + } + } + + #[test] + fn test_molop_no_lines() { + let (h0, h1, h2) = create_test_voigt_tables(); + + let config = MolopConfig { + tmolim: 10000.0, + dfrcon: 1.0, + }; + + let mol_data = vec![1e10, 1e10, 1e10]; + let dop_data = vec![1e-10, 1e-10, 1e-10]; + + let model = MolopModelState { + temp: 5000.0, + elec: 1e12, + mol: MolModelState { + rrmol: &mol_data, + dopmol: &dop_data, + }, + plan: 1.0, + stim: 1.0, + }; + + let freq_data = vec![0.0; 100]; + let freq = MolopFreqParams { + freq: &freq_data, + nfreq: 100, + nfreqs: 100, + }; + + let gfm = vec![]; + let exclm = vec![]; + let freqm = vec![]; + let grm = vec![]; + let gsm = vec![]; + let gvdw = vec![]; + let indatm = vec![]; + let inmlin = vec![]; + let ijcmtr = vec![]; + + let line_data = MolLineData { + freqm: &freqm, + exclm: &exclm, + gfm: &gfm, + grm: &grm, + gsm: &gsm, + gvdw: &gvdw, + indatm: &indatm, + inmlin: &inmlin, + ijcmtr: &ijcmtr, + nlinml: 0, // 没有谱线 + inactm: 0, + }; + + let result = molop_pure(&config, &model, &line_data, &freq, 1e-10, (&h0, &h1, &h2)); + + // 没有谱线时应该返回零 + for i in 0..result.ablin.len() { + assert_relative_eq!(result.ablin[i], 0.0); + assert_relative_eq!(result.emlin[i], 0.0); + } + } + + #[test] + fn test_molop_inactive_list() { + let (h0, h1, h2) = create_test_voigt_tables(); + + let config = MolopConfig { + tmolim: 10000.0, + dfrcon: 1.0, + }; + + let mol_data = vec![1e10, 1e10, 1e10]; + let dop_data = vec![1e-10, 1e-10, 1e-10]; + + let model = MolopModelState { + temp: 5000.0, + elec: 1e12, + mol: MolModelState { + rrmol: &mol_data, + dopmol: &dop_data, + }, + plan: 1.0, + stim: 1.0, + }; + + let freq_data = vec![0.0; 100]; + let freq = MolopFreqParams { + freq: &freq_data, + nfreq: 100, + nfreqs: 100, + }; + + let gfm = vec![]; + let exclm = vec![]; + let freqm = vec![]; + let grm = vec![]; + let gsm = vec![]; + let gvdw = vec![]; + let indatm = vec![]; + let inmlin = vec![]; + let ijcmtr = vec![]; + + let line_data = MolLineData { + freqm: &freqm, + exclm: &exclm, + gfm: &gfm, + grm: &grm, + gsm: &gsm, + gvdw: &gvdw, + indatm: &indatm, + inmlin: &inmlin, + ijcmtr: &ijcmtr, + nlinml: 10, + inactm: 1, // 不激活 + }; + + let result = molop_pure(&config, &model, &line_data, &freq, 1e-10, (&h0, &h1, &h2)); + + // 不激活时应该返回零 + for i in 0..result.ablin.len() { + assert_relative_eq!(result.ablin[i], 0.0); + assert_relative_eq!(result.emlin[i], 0.0); + } + } + + #[test] + fn test_molop_basic_calculation() { + let (h0, h1, h2) = create_test_voigt_tables(); + + let config = MolopConfig { + tmolim: 10000.0, + dfrcon: 1.0, + }; + + // 设置一条简单的分子谱线 + let mol_data = vec![1e10]; // 分子数密度 + let dop_data = vec![0.1]; // 多普勒宽度 + + let model = MolopModelState { + temp: 5000.0, + elec: 1e12, + mol: MolModelState { + rrmol: &mol_data, + dopmol: &dop_data, + }, + plan: 2.0, + stim: 1.0, + }; + + // 频率网格：中心在 1e15 Hz + let nfreq = 100; + let freq_data: Vec = (0..nfreq).map(|i| 1e15 + (i as f64 - 50.0) * 1e11).collect(); + let freq = MolopFreqParams { + freq: &freq_data, + nfreq, + nfreqs: nfreq, + }; + + // 一条分子谱线 + let gfm = vec![0.0]; // log gf = 0 + let exclm = vec![0.0]; // 激发能 = 0 + let freqm = vec![1e15]; // 中心频率 + let grm = vec![0.1]; // 辐射加宽 + let gsm = vec![0.01]; // Stark 加宽 + let gvdw = vec![0.05]; // van der Waals 加宽 + let indatm = vec![0]; // 分子索引 + let inmlin = vec![0]; // 谱线索引 + let ijcmtr = vec![50]; // 中心频率索引 + + let line_data = MolLineData { + freqm: &freqm, + exclm: &exclm, + gfm: &gfm, + grm: &grm, + gsm: &gsm, + gvdw: &gvdw, + indatm: &indatm, + inmlin: &inmlin, + ijcmtr: &ijcmtr, + nlinml: 1, + inactm: 0, + }; + + let result = molop_pure(&config, &model, &line_data, &freq, 1e-10, (&h0, &h1, &h2)); + + // 检查中心附近有不透明度 + let center_idx = 50; + assert!(result.ablin[center_idx] > 0.0, "中心应该有不透明度"); + + // 检查发射率与不透明度的关系 + for i in 3..nfreq { + if result.ablin[i] > 0.0 { + assert_relative_eq!(result.emlin[i], result.ablin[i] * model.plan, epsilon = 1e-10); + } + } + } +} diff --git a/src/math/phe2.rs b/src/math/phe2.rs new file mode 100644 index 0000000..c3518b9 --- /dev/null +++ b/src/math/phe2.rs @@ -0,0 +1,457 @@ +//! He II 线不透明度和发射率计算。 +//! +//! 重构自 SYNSPEC `phe2.f`。 +//! +//! 使用 Schoening 和 Butler 计算的轮廓系数评估给定 He II 线的不透明度和发射率。 + +// ============================================================================ +// 物理常数 +// ============================================================================ + +/// 普朗克常数 × 光速 (erg·cm) +const HC: f64 = 1.986477e-16; + +/// He II 电离能 (cm⁻¹) +const HE2_IONIZATION_ENERGY: f64 = 438916.146; + +/// He II 线波长因子 (Å·cm⁻¹) +/// WLIN = 227.838 / (1/II - 1/JJ) for n <= 2 +/// WLIN = 227.7776 / (1/II - 1/JJ) for n > 2 +const WLIN_FACTOR_1: f64 = 227.838; +const WLIN_FACTOR_2: f64 = 227.7776; + +/// 玻尔兹曼常数相关的激发能因子 +/// 631479 = HC * HE2_IONIZATION_ENERGY / k (简化形式) +const EXCITATION_FACTOR: f64 = 631479.0; + +/// 萨哈因子常数 +/// 4.1412E-16 / (k * sqrt(T)) 的前置因子 +const SAHA_PREFACTOR: f64 = 4.1412e-16; + +/// 转换因子 (2.3025851 = ln(10)) +const LN10: f64 = 2.3025851; + +/// 发射率转换因子 (1.4747E-2) +const EMIS_FACTOR: f64 = 1.4747e-2; + +/// 频率转换因子 (1E-15) +const FREQ_SCALE: f64 = 1e-15; + +/// 频率到波长转换因子 (Å·Hz) +const FREQ_TO_WAVE: f64 = 2.997925e17; + +// ============================================================================ +// He II 振子强度数据 +// ============================================================================ + +/// He II 线的振子强度 (OSCHE2) +/// +/// 对应 Fortran DATA 语句中的 19 个值。 +/// 这些是 Schoening 和 Butler 计算的振子强度。 +const OSCHE2: [f64; 19] = [ + 6.407e-1, 1.506e-1, 5.584e-2, 2.768e-2, + 1.604e-2, 1.023e-2, 6.980e-3, + 8.421e-1, 3.230e-2, 1.870e-2, 1.196e-2, 8.187e-3, + 5.886e-3, 4.393e-3, 3.375e-3, 2.656e-3, + 1.038, 1.793e-1, 6.549e-2, +]; + +// ============================================================================ +// 参数结构体 +// ============================================================================ + +/// PHE2 输入参数。 +/// +/// 包含计算 He II 线不透明度所需的所有数据。 +#[derive(Debug, Clone)] +pub struct Phe2Params<'a> { + /// 谱线索引 (ISPEC) + /// 定义在 HE2INI 中，He II 线从 ISPEC=6 开始 + pub ispec: i32, + + /// 深度索引 (ID) + pub id: i32, + + /// He II 元素索引 (IELHE2) + /// 0 表示不处理 He II NLTE + pub ielhe2: i32, + + /// NLTE 模式标志 (INLTE) + /// > 0 表示使用 NLTE 布居数 + pub inlte: i32, + + /// 频率点数 (NFREQ) + pub nfreq: i32, + + /// 频率数组 (FREQ, Hz) + pub freq: &'a [f64], + + /// 波长数组 (WLAM, Å) + pub wlam: &'a [f64], + + /// 温度 (TEMP, K) + pub temp: f64, + + /// 电子密度 (ELEC) + pub elec: f64, + + /// He III 布居数 (来自 POPUL 或 RRR) + pub he3_pop: f64, + + /// He II NLTE 能级数 (NLHE2 = NLAST - NFIRST + 1) + pub nlhe2: i32, + + /// He II 能级起始索引 (NFIRST) + pub nfirst_he2: i32, + + /// 各能级布居数 (POPUL) + pub popul: &'a [f64], + + /// He II 线轮廓数据 (PRFHE2) + /// PRFHE2(ILINE, ID, IWL) - 线索引 × 深度 × 波长点 + pub prfhe2: &'a [f64], + + /// He II 线波长表 (WLHE2) + /// WLHE2(ILINE, IWL) - log10(|波长 - 线心波长|) + pub wlhe2: &'a [f64], + + /// He II 线波长点数 (NWLHE2) + pub nwlhe2: i32, + + /// He II 线下能级主量子数 (ILHE2) + pub ilhe2: i32, + + /// He II 线上能级主量子数 (IUHE2) + pub iuhe2: i32, + + /// 激光 delta 标志 (lasdel) + /// 当 lasdel=true 且 popi <= popj 时跳过该频率点 + pub lasdel: bool, +} + +/// PHE2 输出结果。 +#[derive(Debug, Clone)] +pub struct Phe2Output { + /// 吸收系数数组 (ABLIN) + pub ablin: Vec, + + /// 发射系数数组 (EMLIN) + pub emlin: Vec, +} + +// ============================================================================ +// PHE2 函数 +// ============================================================================ + +/// 计算 He II 线不透明度和发射率。 +/// +/// 使用 Schoening 和 Butler 预计算的轮廓系数表， +/// 通过插值获得给定频率点的不透明度。 +/// +/// # 参数 +/// +/// * `params` - 输入参数结构体 +/// +/// # 返回 +/// +/// 返回包含吸收系数和发射系数数组的输出结构体 +/// +/// # Fortran 原始代码 +/// +/// ```fortran +/// SUBROUTINE PHE2(ISPEC,ID,ABLIN,EMLIN) +/// ``` +pub fn phe2(params: &Phe2Params) -> Phe2Output { + let nfreq = params.nfreq as usize; + + // 初始化输出数组 + let mut ablin = vec![0.0; nfreq]; + let mut emlin = vec![0.0; nfreq]; + + // 计算线索引 (He II 线从 ISPEC=6 开始) + let iline = (params.ispec - 5) as usize; + + // 检查线索引有效性 + if iline == 0 || iline > 19 { + return Phe2Output { ablin, emlin }; + } + + // 提取轮廓数据到局部数组 + let nwl = params.nwlhe2 as usize; + let mut prf0 = vec![0.0; nwl]; + let mut wll = vec![0.0; nwl]; + + for iwl in 0..nwl { + // PRFHE2(ILINE, ID, IWL) 存储为 [iline][id][iwl] 或线性化 + // 假设存储顺序为 [iline * nd * nwl + id * nwl + iwl] + let id_idx = (params.id - 1) as usize; // Fortran 1-indexed + let prf_idx = (iline - 1) * 100 * 36 + id_idx * 36 + iwl; // 假设 MDEPTH=100, 36 wavelength points max + let wl_idx = (iline - 1) * 36 + iwl; + + if prf_idx < params.prfhe2.len() { + prf0[iwl] = params.prfhe2[prf_idx]; + } + if wl_idx < params.wlhe2.len() { + wll[iwl] = params.wlhe2[wl_idx]; + } + } + + // 计算线心波长 + let i = params.ilhe2; + let j = params.iuhe2; + let ii = (i * i) as f64; + let jj = (j * j) as f64; + + let wlin = if i <= 2 { + WLIN_FACTOR_1 / (1.0 / ii - 1.0 / jj) + } else { + WLIN_FACTOR_2 / (1.0 / ii - 1.0 / jj) + }; + + let t = params.temp; + + // He III 布居数 (LTE 或 NLTE) + let (pp, nlhe2) = if params.ielhe2 > 0 && params.inlte > 0 { + (params.he3_pop, params.nlhe2) + } else { + (params.he3_pop, 0) + }; + + // 计算下能级布居数 + let pp_scaled = pp * params.elec * SAHA_PREFACTOR / t / t.sqrt() * ii; + + let popi = if i <= nlhe2 && params.inlte > 0 { + // NLTE 布居数 + let level_idx = (params.nfirst_he2 + i - 2) as usize; // 转换为 0-indexed + if level_idx < params.popul.len() { + params.popul[level_idx] + } else { + 0.0 + } + } else { + // LTE 布居数 + pp_scaled * (EXCITATION_FACTOR / t / ii).exp() + }; + + // 计算上能级布居数 + let popj = if j <= nlhe2 { + // NLTE 布居数 (归一化到下能级) + let level_idx = (params.nfirst_he2 + j - 2) as usize; + if level_idx < params.popul.len() { + params.popul[level_idx] * ii / jj + } else { + 0.0 + } + } else { + // LTE 布居数 + pp_scaled * (EXCITATION_FACTOR / t / jj).exp() + }; + + // 振子强度因子 + let fid = 0.02654 * OSCHE2[iline - 1]; + + // 遍历频率点 + for ij in 2..nfreq { + // 计算与线心的波长差 + let al = (params.wlam[ij] - wlin).abs(); + let al = if al < 1e-4 { 1e-4 } else { al }; + let al = al.log10(); + + // 在轮廓表中插值 + let mut iw0 = 0; + for iwl in 0..(nwl - 1) { + iw0 = iwl; + if al <= wll[iwl + 1] { + break; + } + } + + let iw1 = iw0 + 1; + if iw1 >= nwl { + continue; + } + + // 线性插值 + let denom = wll[iw1] - wll[iw0]; + let prh = if denom.abs() > 1e-10 { + (prf0[iw0] * (wll[iw1] - al) + prf0[iw1] * (al - wll[iw0])) / denom + } else { + prf0[iw0] + }; + + // 计算线强度 + let sg = (prh * LN10).exp() * fid; + + // 检查布居数反转 + if (popi - popj) <= 0.0 && params.lasdel { + continue; + } + + // 累加不透明度和发射率 + ablin[ij] += sg * (popi - popj); + emlin[ij] += sg * popj * EMIS_FACTOR * (params.freq[ij] * FREQ_SCALE).powi(3); + } + + Phe2Output { ablin, emlin } +} + +// ============================================================================ +// 测试 +// ============================================================================ + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + /// 创建测试用的参数 + fn create_test_params() -> Phe2Params<'static> { + // 创建简单的测试数据 + static FREQ: [f64; 10] = [ + 1e15, 1.1e15, 1.2e15, 1.3e15, 1.4e15, + 1.5e15, 1.6e15, 1.7e15, 1.8e15, 1.9e15, + ]; + static WLAM: [f64; 10] = [ + 2997.9, 2725.4, 2498.3, 2306.1, 2141.4, + 1998.6, 1873.7, 1763.5, 1665.5, 1577.9, + ]; + static POPUL: [f64; 20] = [ + 1e10, 5e9, 2e9, 1e9, 5e8, + 2e8, 1e8, 5e7, 2e7, 1e7, + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, + ]; + static PRFHE2: [f64; 100 * 36 * 19] = [0.0; 100 * 36 * 19]; + static WLHE2: [f64; 36 * 19] = [0.0; 36 * 19]; + + Phe2Params { + ispec: 6, // 第一条 He II 线 + id: 1, + ielhe2: 3, + inlte: 0, + nfreq: 10, + freq: &FREQ, + wlam: &WLAM, + temp: 20000.0, + elec: 1e13, + he3_pop: 1e10, + nlhe2: 0, + nfirst_he2: 11, + popul: &POPUL, + prfhe2: &PRFHE2, + wlhe2: &WLHE2, + nwlhe2: 10, + ilhe2: 2, + iuhe2: 3, + lasdel: false, + } + } + + #[test] + fn test_phe2_basic() { + let params = create_test_params(); + let result = phe2(¶ms); + + // 验证输出数组长度 + assert_eq!(result.ablin.len(), 10); + assert_eq!(result.emlin.len(), 10); + } + + #[test] + fn test_phe2_invalid_line_index() { + let mut params = create_test_params(); + params.ispec = 1; // 不是 He II 线 (ISPEC < 6) + + let result = phe2(¶ms); + + // 应该返回全零数组 + assert!(result.ablin.iter().all(|&x| x == 0.0)); + assert!(result.emlin.iter().all(|&x| x == 0.0)); + } + + #[test] + fn test_phe2_line_index_calculation() { + // 验证线索引计算 + // ISPEC=6 -> iline=1 + // ISPEC=7 -> iline=2 + assert_eq!(6 - 5, 1); + assert_eq!(7 - 5, 2); + } + + #[test] + fn test_wlin_calculation() { + // 测试线心波长计算 + // He II Lyman-alpha: n=1 -> n=2 (lower -> upper) + // I = lower level, J = upper level + let i = 1; // lower level + let j = 2; // upper level + let ii = (i * i) as f64; + let jj = (j * j) as f64; + let wlin = WLIN_FACTOR_1 / (1.0 / ii - 1.0 / jj); + + // He II Lyman-alpha 约为 303.8 Å + assert_relative_eq!(wlin, 303.784, epsilon = 0.1); + } + + #[test] + fn test_constants() { + // 验证常数 + assert_relative_eq!(LN10, 2.3025851, epsilon = 1e-7); + assert_relative_eq!(EMIS_FACTOR, 1.4747e-2, epsilon = 1e-6); + assert_relative_eq!(EXCITATION_FACTOR, 631479.0, epsilon = 1.0); + } + + #[test] + fn test_osche2_data() { + // 验证振子强度数组 + assert_eq!(OSCHE2.len(), 19); + assert_relative_eq!(OSCHE2[0], 6.407e-1, epsilon = 1e-4); + assert_relative_eq!(OSCHE2[18], 6.549e-2, epsilon = 1e-5); + } + + #[test] + fn test_population_inversion_skip() { + // 当 lasdel=true 且 popi <= popj 时应该跳过 + let mut params = create_test_params(); + params.lasdel = true; + + // 设置使得 popi <= popj 的情况 + // 这需要特定的温度和布居数设置 + let result = phe2(¶ms); + + // 由于轮廓数据为零，结果应该很小或为零 + // 这是一个边界情况测试 + assert!(result.ablin.iter().all(|&x| x >= 0.0)); + } + + #[test] + fn test_lte_population() { + // 测试 LTE 布居数计算 + let params = create_test_params(); + + // 使用 LTE 模式 (inlte = 0) + let t = params.temp; + let pp = params.he3_pop; + let elec = params.elec; + let ii = 4.0; // n=2 + + let pp_scaled = pp * elec * SAHA_PREFACTOR / t / t.sqrt() * ii; + let popi_lte = pp_scaled * (EXCITATION_FACTOR / t / ii).exp(); + + // 验证 LTE 布居数为正 + assert!(popi_lte > 0.0); + } + + #[test] + fn test_nlte_population() { + // 测试 NLTE 布居数模式 + let mut params = create_test_params(); + params.inlte = 1; // NLTE 模式 + params.nlhe2 = 5; // 有 5 个 NLTE 能级 + + let result = phe2(¶ms); + + // 验证输出 + assert_eq!(result.ablin.len(), 10); + } +} diff --git a/src/math/phtion.rs b/src/math/phtion.rs new file mode 100644 index 0000000..cb634ef --- /dev/null +++ b/src/math/phtion.rs @@ -0,0 +1,416 @@ +//! 光致电离不透明度计算。 +//! +//! 重构自 SYNSPEC `phtion.f` +//! +//! 计算详细光致电离的吸收和发射系数（从 READPH 读取的表格数据）。 + +use crate::state::constants::{HK, BN}; + +/// 光致电离截面数据（对应 COMMON/PHOTCS/） +#[derive(Debug, Clone)] +pub struct PhotcsData { + /// 光致电离截面表 PHOT(MFRQ, MPHOT) + /// [频率索引][光致电离索引] + pub phot: Vec>, + /// 波长下限 WPHT0 + pub wpht0: f64, + /// 波长上限 WPHT1 + pub wpht1: f64, + /// 原子序数/离子信息 APHT(MPHOT) + pub apht: Vec, + /// 激发能 EPHT(MPHOT) [cm^-1] + pub epht: Vec, + /// 统计权重 GPHT(MPHOT) + pub gpht: Vec, + /// 能级索引 JPHT(MPHOT)（负值表示使用 LTE 数） + pub jpht: Vec, + /// 光致电离数 NPHT + pub npht: usize, +} + +impl Default for PhotcsData { + fn default() -> Self { + Self { + phot: Vec::new(), + wpht0: 0.0, + wpht1: 0.0, + apht: Vec::new(), + epht: Vec::new(), + gpht: Vec::new(), + jpht: Vec::new(), + npht: 0, + } + } +} + +/// PHTION 输入参数 +pub struct PhtionParams<'a> { + /// 深度点温度 + pub temp: f64, + /// 频率数组 FRE(MFRQ) + pub fre: &'a [f64], + /// 频率数 + pub nfre: usize, + /// 原子数密度 RRR(ND, MION, MATOM) + pub rrr: &'a [f64], + /// 能级布居数 POPUL(MLEVEL, ND) + pub popul: &'a [f64], + /// 光致电离截面数据 + pub photcs: &'a PhotcsData, +} + +/// PHTION 输入/输出参数（吸收和发射系数） +pub struct PhtionOutput<'a> { + /// 吸收系数 ABSO(MFRQ) + pub abso: &'a mut [f64], + /// 发射系数 EMIS(MFRQ) + pub emis: &'a mut [f64], +} + +/// 光致电离截面数据频率点数 +pub const MFRQ: usize = 135000; + +/// c3 常数：hc/k [cm K] +const C3: f64 = 1.4387886; + +/// 计算光致电离不透明度。 +/// +/// # 参数 +/// * `params` - 输入参数 +/// * `output` - 输出数组（吸收和发射系数会被累加） +/// +/// # 说明 +/// 如果 NPHT <= 0，函数直接返回不做任何计算。 +/// +/// 对于每个光致电离跃迁： +/// - 如果 JPHT(I) <= 0：使用 LTE 布居数（从原子数密度计算） +/// - 如果 JPHT(I) > 0：使用 NLTE 布居数 +pub fn phtion(params: &PhtionParams, output: &mut PhtionOutput) { + let photcs = params.photcs; + + // 如果没有光致电离数据，直接返回 + if photcs.npht == 0 { + return; + } + + let t = params.temp; + let nfre = params.nfre; + + // 预计算 Planck 函数和受激因子 + let mut planf = vec![0.0; nfre]; + let mut stimu = vec![0.0; nfre]; + + for ij in 0..nfre { + let xx = params.fre[ij]; + let x15 = xx * 1.0e-15; + let bnu = BN * x15 * x15 * x15; + let hkf = HK * xx; + let exh = (hkf / t).exp(); + planf[ij] = bnu / (exh - 1.0); + stimu[ij] = 1.0 - 1.0 / exh; + } + + // 遍历所有光致电离跃迁 + for i in 0..photcs.npht { + let pop: f64; + + if photcs.jpht[i] <= 0 { + // 使用 LTE 布居数 + // APHT(I) 编码了原子和离子信息 + // IAT = int(APHT(I)) 是原子序数 + // X = (APHT(I) - IAT + 1e-4) * 1e2 + // ION = int(X) + 1 是离子索引 + let iat = photcs.apht[i] as i32; + let x = (photcs.apht[i] - iat as f64 + 1.0e-4) * 1.0e2; + let ion = x as i32 as usize + 1; + + // 计算布居数 + // 需要索引 rrr 数组：rrr(id, ion, iat) + // 假设 rrr 是按 [深度][离子][原子] 排列的 3D 数组 + // 这里简化为使用 params.rrr 的线性索引 + let rrr_idx = ion; // 简化索引 + if rrr_idx < params.rrr.len() { + pop = params.rrr[rrr_idx] + * photcs.gpht[i] + * (-photcs.epht[i] * C3 / t).exp(); + } else { + continue; + } + } else { + // 使用 NLTE 布居数 + let jj = (photcs.jpht[i] - 1) as usize; // Fortran 1-indexed -> Rust 0-indexed + if jj < params.popul.len() { + pop = params.popul[jj]; + } else { + continue; + } + } + + // 累加吸收和发射系数 + for ij in 0..nfre { + if i < photcs.phot.len() && ij < photcs.phot[i].len() { + let ab = photcs.phot[i][ij] * pop * stimu[ij]; + output.abso[ij] += ab; + output.emis[ij] += ab * planf[ij]; + } + } + } +} + +/// 纯计算版本的 PHTION（返回新数组） +/// +/// # 参数 +/// * `temp` - 温度 +/// * `fre` - 频率数组 +/// * `rrr` - 原子数密度 +/// * `popul` - 能级布居数 +/// * `photcs` - 光致电离截面数据 +/// * `abso_init` - 初始吸收系数 +/// * `emis_init` - 初始发射系数 +/// +/// # 返回 +/// (吸收系数, 发射系数) +pub fn phtion_pure( + temp: f64, + fre: &[f64], + rrr: &[f64], + popul: &[f64], + photcs: &PhotcsData, + abso_init: &[f64], + emis_init: &[f64], +) -> (Vec, Vec) { + let nfre = fre.len(); + let mut abso = abso_init.to_vec(); + let mut emis = emis_init.to_vec(); + + let params = PhtionParams { + temp, + fre, + nfre, + rrr, + popul, + photcs, + }; + + let mut output = PhtionOutput { + abso: &mut abso, + emis: &mut emis, + }; + + phtion(¶ms, &mut output); + + (abso, emis) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn create_test_photcs() -> PhotcsData { + // 创建简单的测试数据：1 个光致电离跃迁，4 个频率点 + let phot = vec![vec![1.0e-18, 2.0e-18, 3.0e-18, 4.0e-18]]; + PhotcsData { + phot, + wpht0: 1000.0, + wpht1: 2000.0, + apht: vec![1.01], // H I + epht: vec![0.0], // 基态 + gpht: vec![2.0], // 统计权重 + jpht: vec![0], // 使用 LTE + npht: 1, + } + } + + #[test] + fn test_phtion_no_data() { + // 测试无数据时直接返回 + let photcs = PhotcsData::default(); + let fre = vec![1.0e15, 2.0e15, 3.0e15]; + let rrr = vec![1.0e10]; + let popul = vec![]; + + let params = PhtionParams { + temp: 10000.0, + fre: &fre, + nfre: 3, + rrr: &rrr, + popul: &popul, + photcs: &photcs, + }; + + let mut abso = vec![0.0; 3]; + let mut emis = vec![0.0; 3]; + let mut output = PhtionOutput { + abso: &mut abso, + emis: &mut emis, + }; + + phtion(¶ms, &mut output); + + // 应该没有变化 + assert_eq!(abso, vec![0.0, 0.0, 0.0]); + assert_eq!(emis, vec![0.0, 0.0, 0.0]); + } + + #[test] + fn test_phtion_lte_population() { + // 测试 LTE 布居数计算 + let photcs = create_test_photcs(); + let fre = vec![1.0e15, 2.0e15, 3.0e15, 4.0e15]; + // APHT = 1.01 -> IAT=1, X=1.01, ION=2 -> rrr[2] + let rrr = vec![0.0, 0.0, 1.0e10]; // ion index 2 -> rrr[2] + let popul = vec![]; + + let params = PhtionParams { + temp: 10000.0, + fre: &fre, + nfre: 4, + rrr: &rrr, + popul: &popul, + photcs: &photcs, + }; + + let mut abso = vec![0.0; 4]; + let mut emis = vec![0.0; 4]; + let mut output = PhtionOutput { + abso: &mut abso, + emis: &mut emis, + }; + + phtion(¶ms, &mut output); + + // 吸收系数应该为正 + for &a in &abso { + assert!(a >= 0.0); + } + // 发射系数应该为正 + for &e in &emis { + assert!(e >= 0.0); + } + } + + #[test] + fn test_phtion_nlte_population() { + // 测试 NLTE 布居数 + let mut photcs = create_test_photcs(); + photcs.jpht = vec![1]; // 使用能级 1 的布居数 + + let fre = vec![1.0e15, 2.0e15, 3.0e15, 4.0e15]; + let rrr = vec![]; + let popul = vec![1.0e8]; // 能级 1 的布居数 + + let params = PhtionParams { + temp: 10000.0, + fre: &fre, + nfre: 4, + rrr: &rrr, + popul: &popul, + photcs: &photcs, + }; + + let mut abso = vec![0.0; 4]; + let mut emis = vec![0.0; 4]; + let mut output = PhtionOutput { + abso: &mut abso, + emis: &mut emis, + }; + + phtion(¶ms, &mut output); + + // 吸收系数应该为正 + for &a in &abso { + assert!(a > 0.0); + } + } + + #[test] + fn test_phtion_accumulates() { + // 测试累加效果 + let photcs = create_test_photcs(); + let fre = vec![1.0e15, 2.0e15, 3.0e15, 4.0e15]; + // APHT = 1.01 -> IAT=1, X=1.01, ION=2 -> rrr[2] + let rrr = vec![0.0, 0.0, 1.0e10]; + let popul = vec![]; + + let params = PhtionParams { + temp: 10000.0, + fre: &fre, + nfre: 4, + rrr: &rrr, + popul: &popul, + photcs: &photcs, + }; + + // 初始非零值 + let mut abso = vec![1.0e-20, 2.0e-20, 3.0e-20, 4.0e-20]; + let mut emis = vec![1.0e-30, 2.0e-30, 3.0e-30, 4.0e-30]; + let mut output = PhtionOutput { + abso: &mut abso, + emis: &mut emis, + }; + + phtion(¶ms, &mut output); + + // 值应该增加（累加） + assert!(abso[0] > 1.0e-20); + assert!(abso[1] > 2.0e-20); + } + + #[test] + fn test_phtion_pure() { + let photcs = create_test_photcs(); + let fre = vec![1.0e15, 2.0e15, 3.0e15, 4.0e15]; + // APHT = 1.01 -> IAT=1, X=1.01, ION=2 -> rrr[2] + let rrr = vec![0.0, 0.0, 1.0e10]; + let popul = vec![]; + + let abso_init = vec![0.0; 4]; + let emis_init = vec![0.0; 4]; + + let (abso, emis) = phtion_pure(10000.0, &fre, &rrr, &popul, &photcs, &abso_init, &emis_init); + + // 吸收和发射系数应该为正 + for a in &abso { + assert!(*a >= 0.0); + } + for e in &emis { + assert!(*e >= 0.0); + } + } + + #[test] + fn test_planck_function_values() { + // 测试 Planck 函数计算 + // 在高频/低温时，受激因子接近 1 + let photcs = create_test_photcs(); + let fre = vec![1.0e15]; // 低频 + // APHT = 1.01 -> IAT=1, X=1.01, ION=2 -> rrr[2] + let rrr = vec![0.0, 0.0, 1.0e10]; + let popul = vec![]; + + let params = PhtionParams { + temp: 10000.0, + fre: &fre, + nfre: 1, + rrr: &rrr, + popul: &popul, + photcs: &photcs, + }; + + let mut abso = vec![0.0]; + let mut emis = vec![0.0]; + let mut output = PhtionOutput { + abso: &mut abso, + emis: &mut emis, + }; + + phtion(¶ms, &mut output); + + // 发射/吸收比应该接近 Planck 函数 + if abso[0] > 0.0 { + let ratio = emis[0] / abso[0]; + // Planck 函数应该是正的 + assert!(ratio > 0.0); + } + } +} diff --git a/src/math/phtx.rs b/src/math/phtx.rs new file mode 100644 index 0000000..83ba2f7 --- /dev/null +++ b/src/math/phtx.rs @@ -0,0 +1,565 @@ +//! 详细光致电离不透明度计算（从 SIGAVS 读取的表格）。 +//! +//! 重构自 SYNSPEC `phtx.f` +//! +//! 计算两种类型的光致电离不透明度： +//! 1. 能级光致电离（NLTE 布居数） +//! 2. 氢/氦光致电离（LTE 布居数） + +use crate::state::constants::{HK, BN}; + +/// c3 常数：hc/k [cm K] +const C3: f64 = 1.4387886; + +/// 能级光致电离参数 +#[derive(Debug, Clone)] +pub struct LevelPhotoData { + /// 最大截面 CRMX(MLEVEL) + pub crmx: Vec, + /// 频率阈值索引 NFCR(MLEVEL) + pub nfcr: Vec, + /// 频率阈值 FRECR(MLEVEL, MCROSS) + pub frecr: Vec>, + /// 截面 CROSR(MLEVEL, MCROSS) + pub crosr: Vec>, +} + +impl Default for LevelPhotoData { + fn default() -> Self { + Self { + crmx: Vec::new(), + nfcr: Vec::new(), + frecr: Vec::new(), + crosr: Vec::new(), + } + } +} + +/// 氢/氦光致电离参数 +#[derive(Debug, Clone)] +pub struct HhePhotoData { + /// 最大截面 CRMY(MPHOT) + pub crmy: Vec, + /// 频率阈值索引 NFQHT(MPHOT) + pub nfqht: Vec, + /// 频率阈值 FRECQ(MPHOT, MCROSS) + pub frecq: Vec>, + /// 截面 QHOT(MPHOT, MCROSS) + pub qhot: Vec>, + /// 原子/离子信息 AQHT(MPHOT) + pub aqht: Vec, + /// 激发能 EQHT(MPHOT) [cm^-1] + pub eqht: Vec, + /// 统计权重 GQHT(MPHOT) + pub gqht: Vec, +} + +impl Default for HhePhotoData { + fn default() -> Self { + Self { + crmy: Vec::new(), + nfqht: Vec::new(), + frecq: Vec::new(), + qhot: Vec::new(), + aqht: Vec::new(), + eqht: Vec::new(), + gqht: Vec::new(), + } + } +} + +/// PHTX 输入参数 +pub struct PhtxParams<'a> { + /// 深度点温度 + pub temp: f64, + /// 深度索引（从 1 开始，用于判断是否在第一个深度点） + pub depth_id: usize, + /// 频率数组 + pub fre: &'a [f64], + /// 频率数 + pub nfre: usize, + /// 总频率数 + pub nfreq: usize, + /// 连续谱频率数 + pub nfreqc: usize, + /// 模式标志（1=连续谱模式） + pub icon: i32, + /// IASV 标志（能级光致电离） + pub iasv: i32, + /// NQHT 标志（氢/氦光致电离） + pub nqht: i32, + /// 能级布居数 POPUL(MLEVEL, ND) + pub popul: &'a [f64], + /// 原子数密度 RRR(ND, MION, MATOM) + pub rrr: &'a [f64], + /// 能级光致电离数据 + pub level_photo: &'a LevelPhotoData, + /// 氢/氦光致电离数据 + pub hhe_photo: &'a HhePhotoData, +} + +/// PHTX 内部状态（用于保存插值结果） +#[derive(Debug, Clone, Default)] +pub struct PhtxState { + /// 插值后的光致电离截面 PHOTI(MCROSS, MFREQ) + pub photi: Vec>, + /// 能级频率索引缓存 IJP(MLEVEL) + pub ijp: Vec, + /// 氢/氦频率索引缓存 IJQ(MPHOT) + pub ijq: Vec, + /// 是否已初始化 + pub initialized: bool, +} + +impl PhtxState { + /// 创建新的状态 + pub fn new(nlevel: usize, nphot: usize, nfreq: usize) -> Self { + Self { + photi: vec![vec![0.0; nfreq]; nlevel.max(nphot)], + ijp: vec![0; nlevel], + ijq: vec![0; nphot], + initialized: false, + } + } +} + +/// PHTX 输入/输出参数 +pub struct PhtxOutput<'a> { + /// 吸收系数 + pub abso: &'a mut [f64], + /// 发射系数 + pub emis: &'a mut [f64], +} + +/// 计算光致电离不透明度。 +/// +/// # 参数 +/// * `params` - 输入参数 +/// * `output` - 输出数组 +/// * `state` - 内部状态（保存插值结果） +/// +/// # 说明 +/// 如果 IASV=0 且 NQHT=0，函数直接返回。 +/// +/// 在第一个深度点（ID=1）时，进行截面插值计算。 +pub fn phtx(params: &PhtxParams, output: &mut PhtxOutput, state: &mut PhtxState) { + // 如果没有光致电离数据，直接返回 + if params.iasv == 0 && params.nqht == 0 { + return; + } + + let t = params.temp; + let nfre = params.nfre; + let nfreq = params.nfreq; + let nfreqc = params.nfreqc; + + // 确定频率范围 + let (ij0, actual_nfre) = if params.icon == 1 { + (0, nfreqc) + } else { + (2, nfre) + }; + + // 预计算 Planck 函数和受激因子 + let mut planf = vec![0.0; nfre]; + let mut stimu = vec![0.0; nfre]; + + for ij in 0..nfre { + let xx = params.fre[ij]; + let x15 = xx * 1.0e-15; + let bnu = BN * x15 * x15 * x15; + let hkf = HK * xx; + let exh = (hkf / t).exp(); + planf[ij] = bnu / (exh - 1.0); + stimu[ij] = 1.0 - 1.0 / exh; + } + + // 处理能级光致电离（IASV） + if params.iasv != 0 { + let level_photo = params.level_photo; + let nlevel = level_photo.crmx.len(); + + // 在第一个深度点，计算截面插值 + if params.depth_id == 1 { + for i in 0..nlevel { + if level_photo.crmx[i] == 0.0 { + continue; + } + + let nfcr_i = level_photo.nfcr[i] as usize; + let ik1_start = state.ijp[i].max(2) as usize; + + // 从索引 2 开始（对应 Fortran 的 IJ=3） + for ij in 2..nfre.min(nfreq) { + let mut ik2 = ik1_start; + + // 找到频率区间 + for ik in ik1_start..nfcr_i { + if level_photo.frecr[i][ik] < params.fre[ij] { + ik2 = ik; + break; + } + } + + let ik1 = ik2; + if ij == 2 { + state.ijp[i] = ik1 as i32; + } + + // 线性插值 + let dfr = if ik1 > 0 && ik1 < nfcr_i { + let denom = level_photo.frecr[i][ik1 - 1] - level_photo.frecr[i][ik1]; + if denom.abs() > 1e-30 { + (params.fre[ij] - level_photo.frecr[i][ik1]) / denom + } else { + 0.0 + } + } else { + 0.0 + }; + + state.photi[i][ij] = level_photo.crosr[i][ik1] + + dfr * (level_photo.crosr[i][ik1 - 1] - level_photo.crosr[i][ik1]); + } + + // 边界处理 + if nfre > 2 { + state.photi[i][0] = state.photi[i][2]; + state.photi[i][1] = state.photi[i][nfreq - 1]; + } + } + } + + // 累加吸收和发射系数 + for i in 0..nlevel { + if level_photo.crmx[i] == 0.0 { + continue; + } + + // 获取能级布居数 + let pop = if i < params.popul.len() { + params.popul[i] + } else { + continue; + }; + + for ij in ij0..actual_nfre.min(nfre) { + let ab = state.photi[i][ij] * pop * stimu[ij]; + output.abso[ij] += ab; + output.emis[ij] += ab * planf[ij]; + } + } + } + + // 处理氢/氦光致电离（NQHT） + if params.nqht == 0 { + return; + } + + let hhe_photo = params.hhe_photo; + let nqht = hhe_photo.crmy.len(); + + // 在第一个深度点，计算截面插值 + if params.depth_id == 1 { + for i in 0..nqht { + if hhe_photo.crmy[i] == 0.0 { + continue; + } + + let nfqht_i = hhe_photo.nfqht[i] as usize; + let ik1_start = state.ijq[i].max(2) as usize; + + for ij in 2..nfre.min(nfreq) { + let mut ik2 = ik1_start; + + for ik in ik1_start..nfqht_i { + if hhe_photo.frecq[i][ik] < params.fre[ij] { + ik2 = ik; + break; + } + } + + let ik1 = ik2; + if ij == 2 { + state.ijq[i] = ik1 as i32; + } + + let dfr = if ik1 > 0 && ik1 < nfqht_i { + let denom = hhe_photo.frecq[i][ik1 - 1] - hhe_photo.frecq[i][ik1]; + if denom.abs() > 1e-30 { + (params.fre[ij] - hhe_photo.frecq[i][ik1]) / denom + } else { + 0.0 + } + } else { + 0.0 + }; + + state.photi[i][ij] = + hhe_photo.qhot[i][ik1] + dfr * (hhe_photo.qhot[i][ik1 - 1] - hhe_photo.qhot[i][ik1]); + } + } + } + + // 累加吸收和发射系数（LTE 布居数） + for i in 0..nqht { + if hhe_photo.crmy[i] == 0.0 { + continue; + } + + // 计算 LTE 布居数 + let iat = hhe_photo.aqht[i] as i32; + let x = (hhe_photo.aqht[i] - iat as f64 + 1.0e-4) * 100.0; + let ion = x as i32 as usize + 1; + + // 简化的 RRR 索引 + let pop = if ion < params.rrr.len() { + params.rrr[ion] * hhe_photo.gqht[i] * (-hhe_photo.eqht[i] * C3 / t).exp() + } else { + continue; + }; + + for ij in 2..actual_nfre.min(nfre) { + let ab = state.photi[i][ij] * pop * stimu[ij]; + output.abso[ij] += ab; + output.emis[ij] += ab * planf[ij]; + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn create_test_level_photo() -> LevelPhotoData { + LevelPhotoData { + crmx: vec![1.0e-18, 0.0, 2.0e-18], // 能级 0 和 2 有光致电离 + nfcr: vec![10, 0, 10], + frecr: vec![ + vec![1.0e15; 10], + vec![], + vec![1.5e15; 10], + ], + crosr: vec![ + vec![1.0e-18; 10], + vec![], + vec![2.0e-18; 10], + ], + } + } + + fn create_test_hhe_photo() -> HhePhotoData { + HhePhotoData { + crmy: vec![1.0e-18], + nfqht: vec![10], + frecq: vec![vec![1.0e15; 10]], + qhot: vec![vec![1.0e-18; 10]], + aqht: vec![1.01], // H I + eqht: vec![0.0], + gqht: vec![2.0], + } + } + + #[test] + fn test_phtx_no_data() { + // 测试无数据时直接返回 + let level_photo = LevelPhotoData::default(); + let hhe_photo = HhePhotoData::default(); + let fre = vec![1.0e15, 2.0e15, 3.0e15]; + + let params = PhtxParams { + temp: 10000.0, + depth_id: 1, + fre: &fre, + nfre: 3, + nfreq: 3, + nfreqc: 3, + icon: 0, + iasv: 0, + nqht: 0, + popul: &[], + rrr: &[], + level_photo: &level_photo, + hhe_photo: &hhe_photo, + }; + + let mut state = PhtxState::default(); + let mut abso = vec![0.0; 3]; + let mut emis = vec![0.0; 3]; + let mut output = PhtxOutput { + abso: &mut abso, + emis: &mut emis, + }; + + phtx(¶ms, &mut output, &mut state); + + // 应该没有变化 + assert_eq!(abso, vec![0.0, 0.0, 0.0]); + assert_eq!(emis, vec![0.0, 0.0, 0.0]); + } + + #[test] + fn test_phtx_level_photo() { + // 测试能级光致电离 + let level_photo = create_test_level_photo(); + let hhe_photo = HhePhotoData::default(); + let fre: Vec = (1..=10).map(|i| i as f64 * 1.0e14).collect(); + + let popul = vec![1.0e8, 0.0, 2.0e8]; + + let params = PhtxParams { + temp: 10000.0, + depth_id: 1, + fre: &fre, + nfre: 10, + nfreq: 10, + nfreqc: 10, + icon: 0, + iasv: 1, + nqht: 0, + popul: &popul, + rrr: &[], + level_photo: &level_photo, + hhe_photo: &hhe_photo, + }; + + let mut state = PhtxState::new(3, 0, 10); + let mut abso = vec![0.0; 10]; + let mut emis = vec![0.0; 10]; + let mut output = PhtxOutput { + abso: &mut abso, + emis: &mut emis, + }; + + phtx(¶ms, &mut output, &mut state); + + // 检查有吸收系数被计算 + let total_abso: f64 = abso.iter().sum(); + assert!(total_abso > 0.0); + } + + #[test] + fn test_phtx_hhe_photo() { + // 测试氢/氦光致电离 + let level_photo = LevelPhotoData::default(); + let hhe_photo = create_test_hhe_photo(); + let fre: Vec = (1..=10).map(|i| i as f64 * 1.0e14).collect(); + + // APHT = 1.01 -> IAT=1, X=1.01, ION=2 -> rrr[2] + let rrr = vec![0.0, 0.0, 1.0e10]; + + let params = PhtxParams { + temp: 10000.0, + depth_id: 1, + fre: &fre, + nfre: 10, + nfreq: 10, + nfreqc: 10, + icon: 0, + iasv: 0, + nqht: 1, + popul: &[], + rrr: &rrr, + level_photo: &level_photo, + hhe_photo: &hhe_photo, + }; + + let mut state = PhtxState::new(0, 1, 10); + let mut abso = vec![0.0; 10]; + let mut emis = vec![0.0; 10]; + let mut output = PhtxOutput { + abso: &mut abso, + emis: &mut emis, + }; + + phtx(¶ms, &mut output, &mut state); + + // 检查有吸收系数被计算 + let total_abso: f64 = abso.iter().sum(); + assert!(total_abso > 0.0); + } + + #[test] + fn test_phtx_icon_mode() { + // 测试 icon=1 模式（连续谱） + let level_photo = create_test_level_photo(); + let hhe_photo = HhePhotoData::default(); + let fre: Vec = (1..=10).map(|i| i as f64 * 1.0e14).collect(); + + let popul = vec![1.0e8, 0.0, 2.0e8]; + + let params = PhtxParams { + temp: 10000.0, + depth_id: 1, + fre: &fre, + nfre: 10, + nfreq: 10, + nfreqc: 5, + icon: 1, // 连续谱模式 + iasv: 1, + nqht: 0, + popul: &popul, + rrr: &[], + level_photo: &level_photo, + hhe_photo: &hhe_photo, + }; + + let mut state = PhtxState::new(3, 0, 10); + let mut abso = vec![0.0; 10]; + let mut emis = vec![0.0; 10]; + let mut output = PhtxOutput { + abso: &mut abso, + emis: &mut emis, + }; + + phtx(¶ms, &mut output, &mut state); + + // 检查有吸收系数被计算 + let total_abso: f64 = abso.iter().sum(); + assert!(total_abso > 0.0); + } + + #[test] + fn test_phtx_second_depth() { + // 测试非第一个深度点（不进行插值） + let level_photo = create_test_level_photo(); + let hhe_photo = HhePhotoData::default(); + let fre: Vec = (1..=10).map(|i| i as f64 * 1.0e14).collect(); + + let popul = vec![1.0e8, 0.0, 2.0e8]; + + let params = PhtxParams { + temp: 10000.0, + depth_id: 2, // 非第一个深度点 + fre: &fre, + nfre: 10, + nfreq: 10, + nfreqc: 10, + icon: 0, + iasv: 1, + nqht: 0, + popul: &popul, + rrr: &[], + level_photo: &level_photo, + hhe_photo: &hhe_photo, + }; + + let mut state = PhtxState::new(3, 0, 10); + // 预填充状态 + state.photi[0][2] = 1.0e-18; + state.photi[2][2] = 2.0e-18; + + let mut abso = vec![0.0; 10]; + let mut emis = vec![0.0; 10]; + let mut output = PhtxOutput { + abso: &mut abso, + emis: &mut emis, + }; + + phtx(¶ms, &mut output, &mut state); + + // 应该使用预填充的值 + // 由于 ij0=2，只有 abso[2] 会被更新 + assert!(abso[2] > 0.0); + } +} diff --git a/src/math/sgmerg.rs b/src/math/sgmerg.rs new file mode 100644 index 0000000..4e7c970 --- /dev/null +++ b/src/math/sgmerg.rs @@ -0,0 +1,274 @@ +//! 合并能级光致电离截面计算。 +//! +//! 重构自 SYNSPEC `sgmerg.f` +//! +//! 计算合并能级的光致电离截面（Stark 加宽相关）。 +//! 注：此函数在 SYNSPEC 中实际上并未使用（formal routine）。 + +/// 氢 Rydberg 频率 [Hz] +const FRH: f64 = 3.28805e15; +/// 光致电离常数 [cm^2 Hz^3] +const PH2: f64 = 2.815e29 * 2.0; +/// 氢电离能 [cm^-1] +const EHB: f64 = 157802.77355; + +/// SGMERG 输入参数 +pub struct SgmergParams<'a> { + /// 能级索引 + pub ii: usize, + /// 深度索引 + pub id: usize, + /// 频率 [Hz] + pub fr: f64, + /// 元素索引 IEL(MLEVEL) + pub iel: &'a [i32], + /// 电荷 IZ(MATOM) + pub iz: &'a [f64], + /// 温度 TEMP(ND) + pub temp: &'a [f64], + /// 合并能级索引 IMRG(MLEVEL) + pub imrg: &'a [i32], + /// 统计权重 G(MLEVEL) + pub g: &'a mut [f64], + /// 合并能级贡献函数 GMER(NMRG, ND) + pub gmer: &'a [f64], + /// 主量子数边界 NQUANT(MLEVEL) + pub nquant: &'a [i32], + /// 氢 Stark 加宽积分 WNHINT(NLMX, ND) + pub wnhint: &'a [f64], + /// 最大主量子数 + pub nlmx: usize, +} + +/// 计算合并能级光致电离截面。 +/// +/// # 参数 +/// * `params` - 输入参数 +/// +/// # 返回 +/// 光致电离截面 [cm^2] +/// +/// # 说明 +/// 此函数计算高里德伯态（合并能级）对不透明度的贡献。 +/// 使用氢 Stark 加宽积分数据。 +pub fn sgmerg(params: &mut SgmergParams) -> f64 { + let ii = params.ii; + let id = params.id; + + // 获取元素和电荷 + let ie = params.iel[ii] as usize; + let ch = params.iz[ie] * params.iz[ie]; + + // 获取合并能级的统计权重 + if ii < params.imrg.len() && params.imrg[ii] > 0 { + let imrg_idx = (params.imrg[ii] - 1) as usize; // Fortran 1-indexed + if id > 0 && imrg_idx < params.gmer.len() { + params.g[ii] = params.gmer[imrg_idx]; + } + } + + let t = params.temp[id]; + let t1 = 1.0 / t; + let ex = EHB * ch * t1; + + // 获取主量子数范围 + let ii0 = if ii > 0 { + params.nquant[ii - 1] as usize + 1 + } else { + 1 + }; + + let mut sum = 0.0; + + // 遍历所有主量子数 + for i in ii0..=params.nlmx { + let x = i as f64; + let xi = 1.0 / (x * x); + let fredg = FRH * ch * xi; + + // 频率必须大于阈值 + if params.fr < fredg { + continue; + } + + let exi = (ex * xi).exp(); + + // 获取氢 Stark 积分 + let wnhint_val = if i < params.wnhint.len() { + params.wnhint[i] + } else { + 0.0 + }; + + let s = exi * wnhint_val * xi / x; + sum += s; + } + + // 计算最终截面 + let g_ii = if ii < params.g.len() { params.g[ii] } else { 1.0 }; + let sg0 = PH2 / (params.fr * params.fr * params.fr * g_ii) * ch * ch; + + sum * sg0 +} + +/// 纯计算版本（不修改状态） +/// +/// # 参数 +/// * `ii` - 能级索引 +/// * `id` - 深度索引 +/// * `fr` - 频率 +/// * `charge` - 电荷（Z^2） +/// * `temp` - 温度 +/// * `n0` - 起始主量子数 +/// * `nlmx` - 最大主量子数 +/// * `g` - 统计权重 +/// * `wnhint` - 氢 Stark 积分数组 +/// +/// # 返回 +/// 光致电离截面 +pub fn sgmerg_pure( + _ii: usize, + _id: usize, + fr: f64, + charge: f64, + temp: f64, + n0: usize, + nlmx: usize, + g: f64, + wnhint: &[f64], +) -> f64 { + let ch = charge; + let t1 = 1.0 / temp; + let ex = EHB * ch * t1; + + let mut sum = 0.0; + + for i in n0..=nlmx { + let x = i as f64; + let xi = 1.0 / (x * x); + let fredg = FRH * ch * xi; + + if fr < fredg { + continue; + } + + let exi = (ex * xi).exp(); + let wnhint_val = if i < wnhint.len() { wnhint[i] } else { 0.0 }; + + let s = exi * wnhint_val * xi / x; + sum += s; + } + + let sg0 = PH2 / (fr * fr * fr * g) * ch * ch; + sum * sg0 +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_constants() { + // 验证常量值 + assert!((FRH - 3.28805e15).abs() < 1e10); + assert!((PH2 - 5.63e29).abs() < 1e27); + assert!((EHB - 157802.77355).abs() < 1e-5); + } + + #[test] + fn test_sgmerg_pure_basic() { + // 基本测试 + let fr = 1.0e16; // 1e16 Hz + let charge = 1.0; // 氢 + let temp = 10000.0; + let n0 = 1; + let nlmx = 10; + let g = 2.0; + let wnhint = vec![0.0; 20]; // 全零积分 + + let result = sgmerg_pure(0, 0, fr, charge, temp, n0, nlmx, g, &wnhint); + + // 全零积分应该返回 0 + assert_eq!(result, 0.0); + } + + #[test] + fn test_sgmerg_pure_nonzero() { + // 非零测试 + let fr = 1.0e16; + let charge = 1.0; + let temp = 10000.0; + let n0 = 1; + let nlmx = 10; + let g = 2.0; + let mut wnhint = vec![0.0; 20]; + wnhint[5] = 1.0e10; // 非零积分 + + let result = sgmerg_pure(0, 0, fr, charge, temp, n0, nlmx, g, &wnhint); + + // 应该有非零结果 + assert!(result > 0.0); + } + + #[test] + fn test_sgmerg_pure_frequency_cutoff() { + // 频率截止测试 + let charge = 1.0; + let temp = 10000.0; + let n0 = 1; + let nlmx = 5; + let g = 2.0; + let mut wnhint = vec![1.0e10; 20]; // 全部非零 + + // 低频：应该只有少数能级贡献（n=1 阈值约 3.29e15 Hz） + let fr_low = 4.0e15; + let result_low = sgmerg_pure(0, 0, fr_low, charge, temp, n0, nlmx, g, &wnhint); + + // 高频：所有能级都应该贡献 + let fr_high = 1.0e16; + let result_high = sgmerg_pure(0, 0, fr_high, charge, temp, n0, nlmx, g, &wnhint); + + // 两者都应该为正 + assert!(result_low > 0.0); + assert!(result_high > 0.0); + + // 注意：高频在分母，所以截面可能更小 + // 但高频允许更多能级贡献，所以具体取决于参数 + // 这里只验证两者都为正 + } + + #[test] + fn test_sgmerg_pure_charge_scaling() { + // 电荷依赖测试 + let fr = 1.0e16; + let temp = 10000.0; + let n0 = 1; + let nlmx = 5; + let g = 2.0; + let wnhint = vec![1.0e10; 20]; + + let result_h = sgmerg_pure(0, 0, fr, 1.0, temp, n0, nlmx, g, &wnhint); + let result_he = sgmerg_pure(0, 0, fr, 4.0, temp, n0, nlmx, g, &wnhint); // Z=2, Z^2=4 + + // He 应该有更大的截面（电荷平方效应） + assert!(result_he > result_h); + } + + #[test] + fn test_sgmerg_pure_temperature_effect() { + // 温度效应测试 + let fr = 1.0e16; + let charge = 1.0; + let n0 = 1; + let nlmx = 5; + let g = 2.0; + let wnhint = vec![1.0e10; 20]; + + let result_low_t = sgmerg_pure(0, 0, fr, charge, 5000.0, n0, nlmx, g, &wnhint); + let result_high_t = sgmerg_pure(0, 0, fr, charge, 20000.0, n0, nlmx, g, &wnhint); + + // 两者都应该为正 + assert!(result_low_t > 0.0); + assert!(result_high_t > 0.0); + } +}