From b71930cf8e2943f8ea8c4f2ecd5246a34e6136ed Mon Sep 17 00:00:00 2001 From: Asfmq <2696428814@qq.com> Date: Wed, 25 Mar 2026 12:33:47 +0800 Subject: [PATCH] =?UTF-8?q?=E6=A8=A1=E5=9D=97=E9=87=8D=E6=9E=84=E5=A4=A7?= =?UTF-8?q?=E9=83=A8=E5=88=86=E5=B7=B2=E5=AE=8C=E6=88=90?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/math/mod.rs | 4 +- src/math/pretab.rs | 199 +++++++++++++++++++++++++++++++++++++++++ src/state/constants.rs | 4 + 3 files changed, 206 insertions(+), 1 deletion(-) create mode 100644 src/math/pretab.rs diff --git a/src/math/mod.rs b/src/math/mod.rs index 985f520..bb90ca0 100644 --- a/src/math/mod.rs +++ b/src/math/mod.rs @@ -187,6 +187,7 @@ mod princ; mod pzeval; mod prnt; mod prsent; +mod pretab; mod profil; mod profsp; mod quartc; @@ -397,7 +398,7 @@ pub use gamhe::{gamhe, GamheData, GamheParams}; pub use gami::gami; pub use getlal::{getlal, GetlalParams, GetlalResult}; pub use gamsp::gamsp; -pub use gfree::{gfree0, gfreed}; +pub use gfree::{gfree0, gfree1, gfreed}; pub use ghydop::{ghydop, GhydopParams, GhydopResult}; pub use gaunt::gaunt; pub use gntk::gntk; @@ -507,6 +508,7 @@ pub use prchan::{prchan, PrchanParams, PrchanOutput, format_change_report}; pub use princ::{princ_pure, PrincParams, PrincOutput, PrincTransResult, PrincDepthResult, format_princ_header, format_princ_line}; pub use prsent::{prsent, PrsentParams, PrsentOutput, ThermTables}; +pub use pretab::{pretab, VoigtTables}; pub use profil::{profil, ProfilParams}; pub use profsp::{profsp, ProfspParams}; pub use pfni::pfni; diff --git a/src/math/pretab.rs b/src/math/pretab.rs new file mode 100644 index 0000000..d64bbdf --- /dev/null +++ b/src/math/pretab.rs @@ -0,0 +1,199 @@ +//! Voigt 函数展开系数预计算表。 +//! +//! 重构自 SYNSPEC `pretab.f`。 +//! +//! 预计算 Voigt 函数的展开系数表： +//! - 200 步/Doppler 宽度 +//! - 最多 10 个 Doppler 宽度 +//! - 共 MVOI=2001 个点 + +use crate::math::interp::interp; + +/// Voigt 表步长 (每 Doppler 宽度的步数) +const VSTEPS: f64 = 200.0; + +/// TABVI 数据: 速度值 (0 到 12, 81 点) +static TABVI: [f64; 81] = [ + 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, + 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, + 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, + 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, + 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, 5.8, + 6.0, 6.2, 6.4, 6.6, 6.8, 7.0, 7.2, 7.4, 7.6, 7.8, + 8.0, 8.2, 8.4, 8.6, 8.8, 9.0, 9.2, 9.4, 9.6, 9.8, + 10.0, 10.2, 10.4, 10.6, 10.8, 11.0, 11.2, 11.4, 11.6, 11.8, + 12.0, +]; + +/// TABH1 数据: H1 展开系数 (81 点) +static TABH1: [f64; 81] = [ + -1.12838, -1.10596, -1.04048, -0.93703, -0.80346, -0.64945, + -0.48552, -0.32192, -0.16772, -0.03012, 0.08594, 0.17789, + 0.24537, 0.28981, 0.31394, 0.32130, 0.31573, 0.30094, 0.28027, + 0.25648, 0.231726, 0.207528, 0.184882, 0.164341, 0.146128, + 0.130236, 0.116515, 0.104739, 0.094653, 0.086005, 0.078565, + 0.072129, 0.066526, 0.061615, 0.057281, 0.053430, 0.049988, + 0.046894, 0.044098, 0.041561, 0.039250, 0.035195, 0.031762, + 0.028824, 0.026288, 0.024081, 0.022146, 0.020441, 0.018929, + 0.017582, 0.016375, 0.015291, 0.014312, 0.013426, 0.012620, + 0.0118860, 0.0112145, 0.0105990, 0.0100332, 0.0095119, + 0.0090306, 0.0085852, 0.0081722, 0.0077885, 0.0074314, + 0.0070985, 0.0067875, 0.0064967, 0.0062243, 0.0059688, + 0.0057287, 0.0055030, 0.0052903, 0.0050898, 0.0049006, + 0.0047217, 0.0045526, 0.0043924, 0.0042405, 0.0040964, + 0.0039595, +]; + +/// Voigt 函数预计算表。 +/// +/// 包含 H0, H1, H2 三个展开系数数组。 +#[derive(Debug, Clone)] +pub struct VoigtTables { + /// H0 展开系数 (Gauss 轮廓) + pub h0tab: Vec, + /// H1 展开系数 (Voigt 轮廓一阶修正) + pub h1tab: Vec, + /// H2 展开系数 (Voigt 轮廓二阶修正) + pub h2tab: Vec, +} + +impl Default for VoigtTables { + fn default() -> Self { + Self::new() + } +} + +impl VoigtTables { + /// 创建新的 Voigt 表 (未初始化)。 + pub fn new() -> Self { + Self { + h0tab: vec![0.0; crate::state::MVOI], + h1tab: vec![0.0; crate::state::MVOI], + h2tab: vec![0.0; crate::state::MVOI], + } + } + + /// 预计算 Voigt 展开系数表。 + /// + /// 重构自 SYNSPEC `PRETAB` 子程序。 + /// + /// # 算法 + /// 1. 设置 H0TAB(i) = (i-1)/VSTEPS 作为速度值 + /// 2. 使用二次插值从 TABH1 计算 H1TAB + /// 3. 计算 H0TAB(i) = exp(-v²) (Gauss 轮廓) + /// 4. 计算 H2TAB(i) = H0TAB(i) * (1 - 2v²) (二阶导数) + pub fn compute(&mut self) { + let n = crate::state::MVOI; + + // 准备插值输入 + let mut x = TABVI.to_vec(); + let mut y = TABH1.to_vec(); + let mut xx = vec![0.0; n]; + let mut yy = vec![0.0; n]; + + // 步骤 1: 设置速度值 + for i in 0..n { + xx[i] = (i as f64) / VSTEPS; + } + + // 步骤 2: 插值计算 H1TAB + // INTERP(X, Y, XX, YY, NX, NXX, NPOL, ILOGX, ILOGY) + // NPOL=2 表示线性插值 (NPOL-1 = 1 阶) + interp(&mut x, &mut y, &mut xx, &mut yy, 81, n, 2, 0, 0); + + // 复制插值结果到 H1TAB + self.h1tab.copy_from_slice(&yy); + + // 步骤 3 & 4: 计算 H0TAB 和 H2TAB + for i in 0..n { + let vv = ((i as f64) / VSTEPS) * ((i as f64) / VSTEPS); + self.h0tab[i] = (-vv).exp(); + self.h2tab[i] = self.h0tab[i] - (vv + vv) * self.h0tab[i]; + } + } +} + +/// 预计算 Voigt 展开系数表 (便捷函数)。 +/// +/// # 返回值 +/// 包含 H0TAB, H1TAB, H2TAB 三个数组的结构体 +pub fn pretab() -> VoigtTables { + let mut tables = VoigtTables::new(); + tables.compute(); + tables +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn test_pretab_basic() { + let tables = pretab(); + + // 检查数组大小 + assert_eq!(tables.h0tab.len(), crate::state::MVOI); + assert_eq!(tables.h1tab.len(), crate::state::MVOI); + assert_eq!(tables.h2tab.len(), crate::state::MVOI); + + // H0TAB 在 v=0 时应为 1 (exp(0) = 1) + assert_relative_eq!(tables.h0tab[0], 1.0, epsilon = 1e-15); + + // H2TAB 在 v=0 时应为 1 (1 - 0*2 = 1) + assert_relative_eq!(tables.h2tab[0], 1.0, epsilon = 1e-15); + + // H0TAB 应随 v 增大而减小 + assert!(tables.h0tab[100] < tables.h0tab[50]); + assert!(tables.h0tab[200] < tables.h0tab[100]); + } + + #[test] + fn test_h0tab_values() { + let tables = pretab(); + + // 检查一些已知值 + // v = 0.5 (i=100), H0 = exp(-0.25) ≈ 0.7788 + let v_100 = 100.0_f64 / VSTEPS; + let expected_h0_100 = (-(v_100 * v_100)).exp(); + assert_relative_eq!(tables.h0tab[100], expected_h0_100, epsilon = 1e-10); + + // v = 1.0 (i=200), H0 = exp(-1) ≈ 0.3679 + let v_200 = 200.0_f64 / VSTEPS; + let expected_h0_200 = (-(v_200 * v_200)).exp(); + assert_relative_eq!(tables.h0tab[200], expected_h0_200, epsilon = 1e-10); + } + + #[test] + fn test_h2tab_values() { + let tables = pretab(); + + // H2 = H0 * (1 - 2*v²) + // 在 v = 0.5 时, H2 = H0 * (1 - 0.5) = 0.5 * H0 + let v_100 = 100.0_f64 / VSTEPS; + let vv = v_100 * v_100; + let expected_h2 = tables.h0tab[100] * (1.0 - 2.0 * vv); + assert_relative_eq!(tables.h2tab[100], expected_h2, epsilon = 1e-10); + } + + #[test] + fn test_h1tab_interpolation() { + let tables = pretab(); + + // H1TAB 应该是从 TABH1 插值得到的 + // 在 v = 0 时, 应该接近 TABH1[0] = -1.12838 + assert_relative_eq!(tables.h1tab[0], TABH1[0], epsilon = 1e-5); + + // 在 v = 1.0 时 (i=200), 应该在 TABH1[10] 和 TABH1[11] 之间 + // TABH1[10] = 0.08594 (v=1.0) + assert_relative_eq!(tables.h1tab[200], TABH1[10], epsilon = 1e-5); + } + + #[test] + fn test_voigt_tables_default() { + let tables = VoigtTables::default(); + assert_eq!(tables.h0tab.len(), crate::state::MVOI); + assert_eq!(tables.h1tab.len(), crate::state::MVOI); + assert_eq!(tables.h2tab.len(), crate::state::MVOI); + } +} diff --git a/src/state/constants.rs b/src/state/constants.rs index 5883126..4e93cac 100644 --- a/src/state/constants.rs +++ b/src/state/constants.rs @@ -103,6 +103,10 @@ pub const MCFE: usize = 7824000; /// ODF 频率范围 pub const MFRO: usize = MFREQL; +// SYNSPEC 相关 +/// Voigt 表大小 (200 steps per doppler width, up to 10 Doppler widths) +pub const MVOI: usize = 2001; + // ============================================================================ // 物理常数 (CGS 单位) // ============================================================================