SpectraRust/src/synspec/math/inibla.rs

//! 线列表初始化驱动程序。
//!
//! 重构自 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<f64>,

    /// Planck 函数值 - PLAN
    pub plan: Vec<f64>,

    /// 受激因子 (1 - e^(-hν/kT)) - STIM
    pub stim: Vec<f64>,

    /// Van der Waals 宽度 - VDWC
    pub vdwc: Vec<f64>,

    /// Doppler 宽度数组 (原子数 × 深度点数) - DOPA1
    pub dopa1: Vec<f64>,
}

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

        // 验证输出数组长度
        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);
    }
}