feat: 添加 lyahhe 函数 - Lyman alpha 氦展宽插值

- 从 synspec54.f:12768 翻译
- 使用 OnceLock 实现延迟初始化
- 支持二分查找和线性插值
- 3 个单元测试通过

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

src/synspec/math/lyahhe.rs

@@ -0,0 +1,186 @@
+//! Lyman alpha broadening by helium.
+//!
+//! Translated from SYNSPEC `lyahhe` subroutine (synspec54.f:12768).
+//!
+//! Calculates the Lyman alpha profile broadened by helium collisions,
+//! using cross-section data from N. Allard.
+
+use std::sync::OnceLock;
+
+/// Maximum number of cross-section data points
+const NXMAX: usize = 1000;
+
+/// Cross-section normalization factor for He-H broadening
+static STHE: OnceLock<f64> = OnceLock::new();
+
+/// Flag for wavelength unit conversion (1 = convert from Angstrom)
+static NUNHHE: OnceLock<i32> = OnceLock::new();
+
+/// Cached cross-section data
+struct HeHData {
+    xlhhe: Vec<f64>,
+    sighhe: Vec<f64>,
+    nxhhe: usize,
+}
+
+static HEH_DATA: OnceLock<HeHData> = OnceLock::new();
+
+/// Initialize the He-H broadening data from a file.
+///
+/// This should be called once before using `lyahhe`. The data is read from
+/// unit 67 (typically `siglyhhe_21_T14500.lam`).
+///
+/// # Arguments
+/// * `sthe_val` - Cross-section normalization factor
+/// * `nunhhe_val` - Flag for wavelength unit conversion
+/// * `data_lines` - Iterator of (wavelength, cross-section) pairs from the data file
+pub fn lyahhe_init<I>(sthe_val: f64, nunhhe_val: i32, data_lines: I)
+where
+    I: IntoIterator<Item = (f64, f64)>,
+{
+    let _ = STHE.set(sthe_val);
+    let _ = NUNHHE.set(nunhhe_val);
+
+    let nunhhe = nunhhe_val;
+    let mut xlhh0 = Vec::with_capacity(NXMAX);
+    let mut sighh0 = Vec::with_capacity(NXMAX);
+
+    for (xl, sig) in data_lines.into_iter().take(NXMAX) {
+        let xl_conv = if nunhhe == 1 {
+            1.0 / (1.0e-8 * xl + 1.0 / 1215.67)
+        } else {
+            xl
+        };
+        xlhh0.push(xl_conv);
+        sighh0.push(sig);
+    }
+
+    let nxhhe = xlhh0.len();
+
+    // Reverse the arrays (Fortran reads in reverse order)
+    let mut xlhhe: Vec<f64> = xlhh0.into_iter().rev().collect();
+    let mut sighhe: Vec<f64> = sighh0.into_iter().rev().collect();
+
+    // Ensure sorted for binary search
+    // After reversal, data should be in ascending wavelength order
+    // but let's verify and sort if needed
+    if nxhhe > 1 && xlhhe[0] > xlhhe[nxhhe - 1] {
+        xlhhe.reverse();
+        sighhe.reverse();
+    }
+
+    let _ = HEH_DATA.set(HeHData {
+        xlhhe,
+        sighhe,
+        nxhhe,
+    });
+}
+
+/// Calculate Lyman alpha profile broadened by helium.
+///
+/// Uses binary search and linear interpolation on pre-loaded cross-section data.
+///
+/// # Arguments
+/// * `xl` - Wavelength (Å)
+/// * `ahe` - He atom density
+///
+/// # Returns
+/// Profile value (0.0 if wavelength is outside data range)
+pub fn lyahhe(xl: f64, ahe: f64) -> f64 {
+    let data = match HEH_DATA.get() {
+        Some(d) => d,
+        None => return 0.0, // Data not initialized
+    };
+
+    let sthe = match STHE.get() {
+        Some(&s) => s,
+        None => return 0.0,
+    };
+
+    if data.nxhhe == 0 {
+        return 0.0;
+    }
+
+    let prof = 0.0_f64;
+
+    // Check bounds
+    if xl > data.xlhhe[data.nxhhe - 1] {
+        return prof;
+    }
+
+    // Binary search for the interval
+    let mut jl = 0usize;
+    let mut ju = data.nxhhe;
+
+    while ju - jl > 1 {
+        let jm = (ju + jl) / 2;
+        if (data.xlhhe[data.nxhhe - 1] > data.xlhhe[0]) == (xl > data.xlhhe[jm]) {
+            jl = jm;
+        } else {
+            ju = jm;
+        }
+    }
+
+    let mut j = jl;
+    if j == 0 {
+        j = 1;
+    }
+    if j >= data.nxhhe {
+        j = data.nxhhe - 1;
+    }
+
+    // Linear interpolation
+    let denom = data.xlhhe[j] - data.xlhhe[j - 1];
+    if denom.abs() < 1.0e-30 {
+        return 0.0;
+    }
+
+    let a1 = (xl - data.xlhhe[j - 1]) / denom;
+    let s1 = (1.0 - a1) * data.sighhe[j - 1] + a1 * data.sighhe[j];
+    let prof = s1 * ahe / sthe * 6.2831855;
+
+    prof
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_lyahhe_not_initialized() {
+        // Without initialization, should return 0.0
+        let result = lyahhe(1215.67, 1.0e20);
+        assert_eq!(result, 0.0);
+    }
+
+    #[test]
+    fn test_lyahhe_basic() {
+        // Initialize with synthetic data
+        let data = vec![
+            (1215.0, 1.0e-20),
+            (1215.5, 2.0e-20),
+            (1215.67, 5.0e-20),
+            (1216.0, 3.0e-20),
+            (1216.5, 1.0e-20),
+        ];
+        lyahhe_init(1.0e21, 0, data);
+
+        // Test interpolation at known point
+        let result = lyahhe(1215.67, 1.0e20);
+        assert!(result > 0.0);
+        assert!(result.is_finite());
+    }
+
+    #[test]
+    fn test_lyahhe_outside_range() {
+        let data = vec![
+            (1215.0, 1.0e-20),
+            (1216.0, 3.0e-20),
+        ];
+        lyahhe_init(1.0e21, 0, data);
+
+        // Wavelength outside range should return 0.0
+        let result = lyahhe(1220.0, 1.0e20);
+        assert_eq!(result, 0.0);
+    }
+}

src/synspec/math/mod.rs

@@ -2,58 +2,210 @@
 //!
 //! 重构自 SYNSPEC Fortran 代码。
 
+mod carbon;
+mod chckab;
 mod count_words;
 mod crosew;
 mod divhe2;
+mod densit;
+mod eldens;
+mod divstr;
+mod dwnfr0;
+mod dwnfr1;
 mod eps;
+mod expint;
 mod extprf;
 mod feautr;
 mod gamhe;
+mod gaunt;
+mod getwrd;
+mod gfree;
+mod ghydop;
 mod griem;
+mod gvdw;
+mod he2lin;
 mod he2sew;
+mod he2set;
+mod he1ini;
+mod he2ini;
+mod hephot;
+mod hidalg;
 mod heset;
+mod hydlin;
+mod hydini;
+mod hydliw;
+mod hydtab;
 mod hylset;
+mod hylsew;
 mod inibla;
+mod iniblh;
 mod iniblm;
+mod inimod;
+mod interp;
+mod inthe2;
+mod inthyd;
+mod intxen;
 mod intrp;
+mod linop;
+mod linopw;
+mod levsol;
+mod lineqs;
 mod ispec;
+mod locate;
+mod lyahhe;
+mod lymlin;
+mod matinv;
+mod molini;
 mod molop;
+mod molset;
+mod opac;
+mod opacon;
+mod ratmat;
+mod opadd;
+mod ougrid;
 mod partdv;
+mod phe1;
 mod phe2;
 mod phtion;
+mod resolv;
+mod rhonen;
+mod rte;
 mod phtx;
+mod radtem;
 mod pretab;
+mod reiman;
+mod sbfhe1;
+mod sbfhmi;
+mod sffhmi;
 mod sffhmi_old;
+mod sghe12;
 mod sgmerg;
+mod setwin;
+mod sigk;
+mod spsigk;
+mod stark0;
+mod starka;
 mod starkir;
+mod timing;
 mod tint;
+mod todens;
+mod tridag;
+mod voigte;
 mod voigtk;
+mod vopf;
+mod h2minus;
+mod h2opf;
+mod wn;
+mod wnstor;
 mod wtot;
+mod pfni;
+mod nlte;
+mod wgtjh1;
+mod xk2dop;
+mod xenini;
+mod yint;
+mod ylintp;
 
+pub use carbon::carbon;
+pub use chckab::{chckab, ChckabParams, ChckabResult};
 pub use count_words::count_words;
 pub use crosew::{croset, crosew, CrosetParams, CrosewParams};
 pub use divhe2::divhe2;
+pub use densit::{densit, DensitMode, DensitParams, DensitResult};
+pub use eldens::{eldens, EldensParams, EldensResult};
+pub use divstr::divstr;
+pub use dwnfr0::dwnfr0;
+pub use dwnfr1::{dwnfr1, Dwnfr1Params};
 pub use eps::eps;
+pub use expint::expint;
 pub use extprf::extprf;
 pub use feautr::{feautr, FeautrParams};
 pub use gamhe::{gamhe, GamheData, GamheParams};
+pub use gaunt::{gaunt, gntk};
+pub use getwrd::getwrd;
+pub use gfree::gfree;
+pub use ghydop::{ghydop, GhydopParams, GhydopResult};
 pub use griem::{griem, GriemParams};
+pub use gvdw::{gvdw, GvdwParams};
+pub use he2lin::{he2lin, he2liw, He2Common, He2linParams, He2linResult, He2liwParams, He2liwWindowParams};
 pub use he2sew::{he2sew, He2WindowParams};
+pub use he2set::{he2set, He2setParams, He2setOutput};
+pub use he1ini::{he1ini, He1ProfileData, He1Profile4471, He1ProfileOther};
+pub use he2ini::{he2ini, He2InitResult, He2ProfileTable};
+pub use hephot::hephot;
+pub use hidalg::hidalg;
+pub use hydlin::{hydlin, HydlinParams, HydlinResult};
+pub use hydliw::{hydliw, HydliwCommon, HydliwParams, HydliwResult, HydliwWindowParams};
+pub use hydini::{hydini, HydiniParams, HydInitResult, HydProfileTable, HydTableSource};
+pub use hydtab::{hydtab, HydtabParams};
 pub use heset::{heset, HesetInput, HesetOutput};
 pub use hylset::{hylset, HylsetParams, HylsetOutput};
+pub use hylsew::{hylsew, HylsewOutput};
 pub use inibla::{inibla, IniblaParams, IniblaOutput, compute_doppler_width, compute_vdw_width, compute_planck, CL, H, BOLK, BN, HK};
+pub use iniblh::{iniblh, IniblhParams, IniblhOutput, HydrogenLineInfo};
 pub use iniblm::{iniblm, IniblmParams, IniblmOutput, compute_molecular_doppler_width, compute_molecular_planck};
+pub use inimod::{inimod, InimodParams, InimodResult};
+pub use interp::interp;
+pub use inthe2::{inthe2, Inthe2Params, Inthe2Result};
+pub use inthyd::{inthyd, InthydParams, InthydResult};
+pub use intxen::{intxen, IntxenParams, IntxenResult};
 pub use intrp::{intrp, intrp_to_vec};
+pub use linop::{linop, LinopParams, LinopResult, LineData, Phe1Data, Phe2Common};
+pub use linopw::{linopw, LinopwParams, LinopwResult, LinopwLineData, Phe1DataW, Phe2CommonW};
+pub use levsol::levsol;
+pub use lineqs::lineqs;
 pub use ispec::{ispec, PROFILE_VOIGT, PROFILE_HYDROGEN, PROFILE_HEI_4471, PROFILE_HEI_4388, PROFILE_HEI_4026};
+pub use locate::locate;
+pub use lyahhe::{lyahhe, lyahhe_init};
+pub use lymlin::{lymlin, feautr_lyman, LymlinParams, LymlinResult};
+pub use matinv::matinv;
+pub use molini::{molini, MoliniParams, MoliniResult};
 pub use molop::{molop, MolopConfig, MolopModelState, MolopFreqParams, MolLineData, MolModelState, MolopOutput};
+pub use molset::{molset, MolsetParams, MolsetOutput};
+pub use opac::{opac, OpacParams, OpacResult};
+pub use opacon::{opacon, OpaconParams, OpaconResult};
+pub use opadd::{opadd, OpaddParams, OpaddResult};
+pub use ratmat::ratmat;
+pub use ougrid::{ougrid, OugridParams, OugridOutput};
 pub use partdv::{partdv, partdv_with_params, PartdvParams};
+pub use phe1::{phe1, Phe1Params};
 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 radtem::radtem;
+pub use resolv::{resolv, ResolvParams, ResolvResult};
+pub use rhonen::{rhonen, RhonenParams, RhonenResult};
+pub use rte::{rte, RteParams, RteResult};
 pub use pretab::{pretab, VoigtTables};
+pub use reiman::reiman;
+pub use sbfhe1::sbfhe1;
+pub use sbfhmi::sbfhmi;
+pub use sffhmi::sffhmi;
 pub use sffhmi_old::sffhmi_old;
+pub use sghe12::sghe12;
 pub use sgmerg::{sgmerg, sgmerg_pure, SgmergParams};
+pub use setwin::{setwin, SetWinParams};
+pub use sigk::{sigk, SigkParams};
+pub use spsigk::spsigk;
+pub use stark0::{stark0, Stark0Result};
+pub use starka::starka;
 pub use starkir::starkir;
+pub use timing::{timing, TimingResult, TIMING_TABLE, TIMING_FINAL};
 pub use tint::{tint, TintResult};
+pub use todens::{todens, TodensParams, TodensResult};
+pub use tridag::tridag;
+pub use voigte::voigte;
 pub use voigtk::{voigtk, MVOI};
+pub use vopf::vopf;
+pub use h2minus::h2minus;
+pub use h2opf::h2opf;
+pub use wn::wn;
+pub use wgtjh1::{wgtjh1, Wgtjh1Params, Wgtjh1Result};
+pub use wnstor::{wnstor, WnstorParams, WnstorOutput};
 pub use wtot::{wtot, LINE_4471, LINE_4387, LINE_4026, LINE_4922};
+pub use pfni::pfni;
+pub use nlte::{nlte, NlteParams, NlteOutput};
+pub use xk2dop::xk2dop;
+pub use xenini::{xenini, XenInitResult, XenProfileData};
+pub use yint::yint;
+pub use ylintp::ylintp;