SpectraRust/src/tlusty/math/hydrogen/bhe.rs

//! 流体静力学平衡方程矩阵填充。
//!
//! 重构自 TLUSTY 的 BHE, BHED, BHEZ 子程序。
//!
//! 这些函数填充矩阵 A, B, C 的流体静力学平衡行
//! (NFREQE+INHE)-th row。

use crate::tlusty::state::constants::{BOLK, HALF, MDEPTH, MFREQ, MLEVEL, MTOT, TWO, UN};

// ============================================================================
// 控制参数
// ============================================================================

/// 矩阵索引控制参数。
/// 对应 COMMON /MATKEY/
#[derive(Debug, Clone, Default)]
pub struct MatKey {
    /// 方程总数
    pub nn: usize,
    /// NN0
    pub nn0: usize,
    /// 流体静力学平衡方程索引偏移
    pub inhe: isize,
    /// 辐射平衡方程索引偏移
    pub inre: isize,
    /// 粒子守恒方程索引偏移
    pub inpc: isize,
    /// 统计平衡方程索引偏移
    pub inse: isize,
    /// z-距离方程索引偏移
    pub inzd: isize,
    /// 虚拟大质量粒子密度索引偏移
    pub inmp: isize,
    /// 辐射平衡深度点数
    pub ndre: usize,
}

/// 输入参数控制。
/// 对应部分 COMMON /INPPAR/ 和 /FIXDEN/
#[derive(Debug, Clone, Default)]
pub struct InputControl {
    /// 边界条件类型 (0: 恒星大气, 1: 盘新版, 2: 盘旧版, 3: 简单气压边界)
    pub ibche: i32,
    /// 固定密度标志
    pub ifixde: i32,
    /// 辐射压力标志
    pub ifprad: i32,
    /// 表面气压 (用于 ibche=3)
    pub pgas0: f64,
    /// 重力加速度缩放因子
    pub qgrav: f64,
}

/// 模型维度参数。
#[derive(Debug, Clone)]
pub struct ModelDims {
    /// 显式频率数
    pub nfreqe: usize,
    /// 深度点数
    pub nd: usize,
    /// 线性化能级数
    pub nlvexp: usize,
}

impl Default for ModelDims {
    fn default() -> Self {
        Self {
            nfreqe: 0,
            nd: 1,
            nlvexp: 0,
        }
    }
}

// ============================================================================
// BHE 参数结构体
// ============================================================================

/// BHE 输入参数。
#[derive(Debug, Clone)]
pub struct BheParams {
    /// 深度索引 (1-based in Fortran)
    pub id: usize,

    /// 模型维度
    pub dims: ModelDims,

    /// 矩阵索引控制
    pub matkey: MatKey,

    /// 输入控制
    pub ctrl: InputControl,

    // 模型状态 (深度相关)
    /// 温度 (K) [MDEPTH]
    pub temp: Vec<f64>,
    /// 总粒子密度 (cm⁻³) [MDEPTH]
    pub dens: Vec<f64>,
    /// 总粒子数 [MDEPTH]
    pub totn: Vec<f64>,
    /// 电子密度 (cm⁻³) [MDEPTH]
    pub elec: Vec<f64>,
    /// 平均分子量倒数 [MDEPTH]
    pub wmm: Vec<f64>,
    /// 柱质量密度 (g/cm²) [MDEPTH]
    pub dm: Vec<f64>,
    /// 湍流速度 [MDEPTH]
    pub vturb: Vec<f64>,
    /// z-距离 [MDEPTH]
    pub zd: Vec<f64>,

    // 频率相关
    /// 频率权重 [MFREQ]
    pub w: Vec<f64>,
    /// 频率索引映射 [MFREQ]
    pub ijfr: Vec<usize>,
    /// 跳过标志 [MDEPTH][MFREQ]
    pub lskip: Vec<Vec<bool>>,
    /// 表面通量权重 [MFREQ]
    pub fh: Vec<f64>,
    /// 外辐射 (氦) [MFREQ]
    pub hextrd: Vec<f64>,

    // 辐射相关 (显式频率)
    /// 辐射场 (当前) [MFREQ]
    pub rad0: Vec<f64>,
    /// 辐射场 (前) [MFREQ]
    pub radm: Vec<f64>,
    /// FK 算子 (当前) [MFREQ]
    pub fk0: Vec<f64>,
    /// FK 算子 (前) [MFREQ]
    pub fkm: Vec<f64>,
    /// 吸收系数 (当前) [MFREQ]
    pub abso0: Vec<f64>,
    /// 吸收系数温度导数 (当前) [MFREQ]
    pub dabt0: Vec<f64>,
    /// 吸收系数密度导数 (当前) [MFREQ]
    pub dabn0: Vec<f64>,
    /// 发射系数温度导数 (当前) [MFREQ]
    pub demt0: Vec<f64>,
    /// 深度权重 (当前) [MFREQ]
    pub wdep0: Vec<f64>,
    /// 能级导数 [MLEVEL][MFREQ]
    pub drch0: Vec<Vec<f64>>,

    // HE 相关数组 (深度相关)
    /// HEIT - 氦温度相关 [MDEPTH]
    pub heit: Vec<f64>,
    /// HEIN - 氦密度相关 [MDEPTH]
    pub hein: Vec<f64>,
    /// HEITM - 氦温度相关 (前) [MDEPTH]
    pub heitm: Vec<f64>,
    /// HEINM - 氦密度相关 (前) [MDEPTH]
    pub heinm: Vec<f64>,
    /// HEITP - 氦温度相关导数 [MDEPTH]
    pub heitp: Vec<f64>,
    /// HEINP - 氦密度相关导数 [MDEPTH]
    pub heinp: Vec<f64>,

    // HE 能级相关 [nlvexp][MDEPTH]
    pub heip: Vec<Vec<f64>>,
    pub heipm: Vec<Vec<f64>>,
    pub heipp: Vec<Vec<f64>>,

    // 辐射压力
    /// 固定频率辐射压力 [MDEPTH]
    pub fprd: Vec<f64>,

    // Psi 向量
    /// Ψ (当前) [MTOT]
    pub psi0: Vec<f64>,
    /// Ψ (前) [MTOT]
    pub psim: Vec<f64>,
    /// Ψ (后) [MTOT]
    pub psip: Vec<f64>,
}

impl Default for BheParams {
    fn default() -> Self {
        Self {
            id: 1,
            dims: ModelDims::default(),
            matkey: MatKey::default(),
            ctrl: InputControl::default(),
            temp: vec![0.0; MDEPTH],
            dens: vec![0.0; MDEPTH],
            totn: vec![0.0; MDEPTH],
            elec: vec![0.0; MDEPTH],
            wmm: vec![0.0; MDEPTH],
            dm: vec![0.0; MDEPTH],
            vturb: vec![0.0; MDEPTH],
            zd: vec![0.0; MDEPTH],
            w: vec![0.0; MFREQ],
            ijfr: vec![0; MFREQ],
            lskip: vec![vec![false; MFREQ]; MDEPTH],
            fh: vec![0.0; MFREQ],
            hextrd: vec![0.0; MFREQ],
            rad0: vec![0.0; MFREQ],
            radm: vec![0.0; MFREQ],
            fk0: vec![0.0; MFREQ],
            fkm: vec![0.0; MFREQ],
            abso0: vec![0.0; MFREQ],
            dabt0: vec![0.0; MFREQ],
            dabn0: vec![0.0; MFREQ],
            demt0: vec![0.0; MFREQ],
            wdep0: vec![0.0; MFREQ],
            drch0: vec![vec![0.0; MFREQ]; MLEVEL],
            heit: vec![0.0; MDEPTH],
            hein: vec![0.0; MDEPTH],
            heitm: vec![0.0; MDEPTH],
            heinm: vec![0.0; MDEPTH],
            heitp: vec![0.0; MDEPTH],
            heinp: vec![0.0; MDEPTH],
            heip: vec![vec![0.0; MDEPTH]; MLEVEL],
            heipm: vec![vec![0.0; MDEPTH]; MLEVEL],
            heipp: vec![vec![0.0; MDEPTH]; MLEVEL],
            fprd: vec![0.0; MDEPTH],
            psi0: vec![0.0; MTOT],
            psim: vec![0.0; MTOT],
            psip: vec![0.0; MTOT],
        }
    }
}

// ============================================================================
// BHE 可变状态
// ============================================================================

/// BHE 可变状态 (矩阵和向量)。
#[derive(Debug, Clone, Default)]
pub struct BheState {
    /// A 矩阵 [MTOT][MTOT]
    pub a: Vec<Vec<f64>>,
    /// B 矩阵 [MTOT][MTOT]
    pub b: Vec<Vec<f64>>,
    /// C 矩阵 [MTOT][MTOT]
    pub c: Vec<Vec<f64>>,
    /// 右端向量 [MTOT]
    pub vecl: Vec<f64>,

    // CMATZD - z-m 关系的 C 矩阵元素 (BHED 专用)
    pub czz: f64,
    pub czn: f64,
    pub cze: f64,
    pub czm: f64,
}

impl BheState {
    pub fn new() -> Self {
        Self {
            a: vec![vec![0.0; MTOT]; MTOT],
            b: vec![vec![0.0; MTOT]; MTOT],
            c: vec![vec![0.0; MTOT]; MTOT],
            vecl: vec![0.0; MTOT],
            czz: 0.0,
            czn: 0.0,
            cze: 0.0,
            czm: 0.0,
        }
    }
}

// ============================================================================
// 物理常量
// ============================================================================

/// GN = 重力加速度相关常数
const GN: f64 = 1.0;
/// GP = 重力加速度相关常数 (大质量粒子)
const GP: f64 = 1.0;
/// PCK = 辐射压力常数
const PCK: f64 = 4.0 * std::f64::consts::PI / 3.0;

// ============================================================================
// 辅助函数
// ============================================================================

/// 计算 erfcx(x) = exp(x²) * erfc(x)
/// 使用 Fortran 代码中的系数
fn erfcx(x: f64) -> f64 {
    if x < 3.0 {
        // 使用近似公式
        8.86226925e-1 * (x * x).exp() * erfc_approx(x)
    } else {
        // 渐近展开
        HALF * (UN - HALF / (x * x)) / x
    }
}

/// erfc 近似
fn erfc_approx(x: f64) -> f64 {
    // 简单近似，精度足够
    if x < 0.0 {
        2.0 - erfc_approx(-x)
    } else if x < 6.0 {
        let t = 1.0 / (1.0 + 0.5 * x);
        let poly = t * (0.17087277 + t
            * (-0.82215223 + t
                * (1.48851587 + t
                    * (-1.13520398 + t
                        * (0.27886807 + t * (-0.18628806 + t * 0.4206475))))));
        t * (-x * x).exp() * poly
    } else {
        0.0
    }
}

// ============================================================================
// BHE 主函数
// ============================================================================

/// 流体静力学平衡方程矩阵填充。
///
/// 填充矩阵 A 和 B 的 (NFREQE+INHE)-th 行，
/// 对应流体静力学平衡方程。
///
/// # 参数
/// - `params`: 输入参数
/// - `state`: 可变状态 (矩阵和向量)
pub fn bhe(params: &BheParams, state: &mut BheState) {
    let id = params.id;
    let dims = &params.dims;
    let matkey = &params.matkey;

    // 计算行索引 (0-based)
    let nhe = dims.nfreqe + matkey.inhe as usize;
    let nre = dims.nfreqe + matkey.inre as usize;
    let npc = dims.nfreqe + matkey.inpc as usize;
    let nse = dims.nfreqe + matkey.inse as usize - 1; // Fortran: NSE = NFREQE+INSE-1

    // 固定密度情况
    if params.ctrl.ifixde > 0 {
        state.b[nhe][nhe] = UN;
        state.b[nhe][npc] = -UN;
        state.vecl[nhe] = params.dens[id - 1] * params.wmm[id - 1]
            + params.elec[id - 1]
            - params.totn[id - 1];
        return;
    }

    // 虚拟大质量粒子密度方程
    if matkey.inmp > 0 {
        let nmp = dims.nfreqe + matkey.inmp as usize;
        state.b[nmp][nmp] = -UN;
        state.b[nmp][nhe] = UN;
        if matkey.inpc > 0 {
            state.b[nmp][npc] = -UN;
        }
    }

    // 初始化累积变量
    let mut hext: f64 = 0.0;
    let mut hexn: f64 = 0.0;
    let mut grd: f64 = 0.0;
    let mut hex = vec![0.0; dims.nlvexp];

    // 上边界条件 (id = 1)
    if id > 1 {
        // 正常深度点 - 跳转到标号 50
        fill_normal_depth(params, state, id, nhe, nre, npc, nse, &mut grd);
        return;
    }

    // === 上边界条件 (ID = 1) ===
    // 基于 Mihalas (1978) 的线性化方程 (7-10)

    let mut x1 = 0.0;
    if dims.nfreqe > 0 && params.ctrl.ifprad > 0 {
        x1 = PCK / params.dens[id - 1];

        for ij in 1..=dims.nfreqe {
            let ijt = params.ijfr[ij - 1];
            if !params.lskip[id - 1][ijt] {
                let fluxw = params.w[ijt]
                    * (params.fh[ijt] * params.rad0[ij - 1] - params.hextrd[ijt]);
                grd += fluxw * params.abso0[ij - 1];
                hexn += fluxw * params.dabn0[ij - 1];
                hext += fluxw * params.dabt0[ij - 1];

                for i in 0..dims.nlvexp {
                    hex[i] += fluxw * params.drch0[i][ij - 1];
                }

                // 平均强度列
                state.b[nhe][ij - 1] = x1 * params.w[ijt] * params.fh[ijt] * params.abso0[ij - 1];
            }
        }
    }

    let rtn = x1 * params.wmm[id - 1] / params.dens[id - 1] * (grd + params.fprd[id - 1]);
    let vt0 = HALF * params.vturb[id - 1] * params.vturb[id - 1] / params.dm[id - 1]
        * params.wmm[id - 1];

    // 总粒子密度、虚拟大质量粒子密度、温度、电子密度列
    state.b[nhe][nhe] = BOLK * params.temp[id - 1] / params.dm[id - 1] - GN * (rtn - vt0);

    if matkey.inmp > 0 {
        let nmp = dims.nfreqe + matkey.inmp as usize;
        state.b[nhe][nmp] = GP * (vt0 - rtn);
    }

    if matkey.inre > 0 {
        state.b[nhe][nre] = BOLK * params.totn[id - 1] / params.dm[0]
            + x1 * (hext + params.heit[id - 1]);
        state.c[nhe][nre] = x1 * params.heitp[id - 1];
    }

    if matkey.inpc > 0 {
        state.b[nhe][npc] =
            x1 * (hexn + params.hein[id - 1]) + GN * (rtn - vt0);
        state.c[nhe][npc] = x1 * params.heinp[id - 1];
    }

    // 能级占据数列
    for ii in 0..dims.nlvexp {
        state.b[nhe][nse + ii] += x1 * (hex[ii] + params.heip[ii][id - 1]);
        state.c[nhe][nse + ii] += x1 * params.heipp[ii][id - 1];
    }

    // 右端向量
    let grav = params.ctrl.qgrav * params.zd[id - 1];
    state.vecl[nhe] = grav
        - BOLK * params.temp[id - 1] * params.totn[id - 1] / params.dm[id - 1]
        - x1 * (grd + params.fprd[id - 1])
        - vt0 / params.wmm[id - 1] * params.dens[id - 1];
}

/// 填充正常深度点 (ID > 1) 的矩阵
fn fill_normal_depth(
    params: &BheParams,
    state: &mut BheState,
    id: usize,
    nhe: usize,
    nre: usize,
    npc: usize,
    nse: usize,
    grd: &mut f64,
) {
    let dims = &params.dims;
    let matkey = &params.matkey;

    // 平均强度列
    if dims.nfreqe > 0 && params.ctrl.ifprad > 0 {
        for ij in 1..=dims.nfreqe {
            if !params.lskip[id - 1][params.ijfr[ij - 1]] {
                *grd += (params.fk0[ij - 1] * params.rad0[ij - 1]
                    - params.fkm[ij - 1] * params.radm[ij - 1])
                    * params.w[params.ijfr[ij - 1]];
                state.a[nhe][ij - 1] = -PCK * params.w[params.ijfr[ij - 1]] * params.fkm[ij - 1];
                state.b[nhe][ij - 1] = PCK * params.w[params.ijfr[ij - 1]] * params.fk0[ij - 1];
            }
        }
    }

    let vt0 = HALF * params.vturb[id - 1] * params.vturb[id - 1] * params.wmm[id - 1];
    let vtm = HALF * params.vturb[id - 2] * params.vturb[id - 2] * params.wmm[id - 2];

    // 总粒子密度列
    state.a[nhe][nhe] = -BOLK * params.temp[id - 2] - GN * vtm;
    state.b[nhe][nhe] = BOLK * params.temp[id - 1] + GN * vt0;

    // 温度列
    if matkey.inre > 0 {
        state.a[nhe][nre] = -BOLK * params.totn[id - 2] + PCK * params.heitm[id - 1];
        state.b[nhe][nre] = BOLK * params.totn[id - 1] + PCK * params.heit[id - 1];
        state.c[nhe][nre] = PCK * params.heitp[id - 1];
    }

    // 电子密度列
    if matkey.inpc > 0 {
        state.a[nhe][npc] = GN * vtm + PCK * params.heinm[id - 1];
        state.b[nhe][npc] = -GN * vt0 + PCK * params.hein[id - 1];
        state.c[nhe][npc] = PCK * params.heinp[id - 1];
    }

    // 虚拟大质量粒子密度列
    if matkey.inmp > 0 {
        let nmp = dims.nfreqe + matkey.inmp as usize;
        state.a[nhe][nmp] = -GP * vtm;
        state.b[nhe][nmp] = GP * vt0;
    }

    // 能级占据数列
    for ii in 0..dims.nlvexp {
        state.a[nhe][nse + ii] += PCK * params.heipm[ii][id - 1];
        state.b[nhe][nse + ii] += PCK * params.heip[ii][id - 1];
        state.c[nhe][nse + ii] += PCK * params.heipp[ii][id - 1];
    }

    // 右端向量
    let grav = params.ctrl.qgrav * (params.dm[id - 1] - params.dm[id - 2]);
    state.vecl[nhe] = grav
        - BOLK * (params.temp[id - 1] * params.totn[id - 1]
            - params.temp[id - 2] * params.totn[id - 2])
        - PCK * (*grd + params.fprd[id - 1])
        - vt0 / params.wmm[id - 1] * params.dens[id - 1]
        + vtm / params.wmm[id - 2] * params.dens[id - 2];
}

// ============================================================================
// BHED 函数
// ============================================================================

/// 流体静力学平衡方程矩阵填充 (扩展版，含 z-m 关系)。
///
/// 与 BHE 类似，但增加:
/// - z-距离与质量深度关系 (NFREQE+INZD)-th 行
/// - 支持盘模型的边界条件
pub fn bhed(params: &BheParams, state: &mut BheState) {
    let id = params.id;
    let dims = &params.dims;
    let matkey = &params.matkey;

    let nhe = dims.nfreqe + matkey.inhe as usize;
    let nre = dims.nfreqe + matkey.inre as usize;
    let npc = dims.nfreqe + matkey.inpc as usize;
    let nse = dims.nfreqe + matkey.inse as usize - 1;

    if matkey.inhe <= 0 {
        // 跳过流体静力学平衡，只处理 z-m 关系
        fill_zm_relation(params, state, id, dims, matkey);
        return;
    }

    // 虚拟大质量粒子密度方程
    if matkey.inmp > 0 {
        let nmp = dims.nfreqe + matkey.inmp as usize;
        state.b[nmp][nmp] = -UN;
        state.b[nmp][nhe] = UN;
        if matkey.inpc > 0 {
            state.b[nmp][npc] = -UN;
        }
    }

    // 初始化
    let mut hext: f64 = 0.0;
    let mut hexn: f64 = 0.0;
    let mut grd: f64 = 0.0;
    let mut hex = vec![0.0; dims.nlvexp];

    if id > 1 {
        // 正常深度点
        fill_normal_depth_bhed(params, state, id, nhe, nre, npc, nse, &mut grd);
        fill_zm_relation(params, state, id, dims, matkey);
        return;
    }

    // === 上边界条件 (ID = 1) ===
    let ibche = params.ctrl.ibche;

    if ibche <= 0 {
        // 1. 恒星大气边界条件 (Mihalas 1978)
        fill_upper_boundary_stellar(params, state, nhe, nre, npc, nse, &mut hext, &mut hexn, &mut grd, &mut hex);
    } else if ibche == 1 {
        // 2. 盘模型边界条件 (Hubeny 1990 新版)
        fill_upper_boundary_disk_new(params, state, nhe, nre, npc, nse, &mut hext, &mut hexn, &mut grd, &mut hex);
    } else if ibche == 2 {
        // 3. 盘模型边界条件 (Hubeny 1990 旧版)
        fill_upper_boundary_disk_old(params, state, nhe, nre, &mut grd);
    } else {
        // 4. 简单气压边界条件 P_gas(ID=1) = PGAS0
        fill_upper_boundary_simple(params, state, nhe, nre);
    }

    // z-m 关系
    fill_zm_relation(params, state, id, dims, matkey);
}

/// 填充恒星大气上边界条件
fn fill_upper_boundary_stellar(
    params: &BheParams,
    state: &mut BheState,
    nhe: usize,
    nre: usize,
    npc: usize,
    nse: usize,
    hext: &mut f64,
    hexn: &mut f64,
    grd: &mut f64,
    hex: &mut [f64],
) {
    let dims = &params.dims;
    let matkey = &params.matkey;

    let x1 = PCK / params.dens[0];

    if dims.nfreqe > 0 {
        for ij in 1..=dims.nfreqe {
            let ijt = params.ijfr[ij - 1];
            if !params.lskip[0][ijt] {
                let fluxw = params.w[ijt]
                    * (params.fh[ijt] * params.rad0[ij - 1] - params.hextrd[ijt]);
                *grd += fluxw * params.abso0[ij - 1];
                *hexn += fluxw * params.dabn0[ij - 1];
                *hext += fluxw * params.dabt0[ij - 1];
                for i in 0..dims.nlvexp {
                    hex[i] += fluxw * params.drch0[i][ij - 1];
                }
                state.b[nhe][ij - 1] = x1 * params.wdep0[ij - 1] * params.fh[ijt] * params.abso0[ij - 1];
            }
        }
    }

    let rtn = x1 * params.wmm[0] / params.dens[0] * (*grd + params.fprd[0]);
    let vt0 = HALF * params.vturb[0] * params.vturb[0] / params.dm[0] * params.wmm[0];

    state.b[nhe][nhe] = BOLK * params.temp[0] / params.dm[0] - GN * (rtn - vt0);

    if matkey.inmp > 0 {
        let nmp = dims.nfreqe + matkey.inmp as usize;
        state.b[nhe][nmp] = GP * (vt0 - rtn);
    }

    if matkey.inre > 0 {
        state.b[nhe][nre] = BOLK * params.psi0[nhe] / params.dm[0]
            + x1 * (*hext + params.heit[0]);
        state.c[nhe][nre] = x1 * params.heitp[0];
    }

    if matkey.inpc > 0 {
        state.b[nhe][npc] = x1 * (*hexn + params.hein[0]) + GN * (rtn - vt0);
        state.c[nhe][npc] = x1 * params.heinp[0];
    }

    for ii in 0..dims.nlvexp {
        state.b[nhe][nse + ii] += x1 * (hex[ii] + params.heip[ii][0]);
        state.c[nhe][nse + ii] += x1 * params.heipp[ii][0];
    }

    let grav = params.ctrl.qgrav * params.zd[0];
    state.vecl[nhe] = grav
        - BOLK * params.temp[0] * params.psi0[nhe] / params.dm[0]
        - x1 * (*grd + params.fprd[0])
        - vt0 / params.wmm[0] * params.dens[0];
}

/// 填充盘模型上边界条件 (新版)
fn fill_upper_boundary_disk_new(
    params: &BheParams,
    state: &mut BheState,
    nhe: usize,
    nre: usize,
    npc: usize,
    nse: usize,
    hext: &mut f64,
    hexn: &mut f64,
    grd: &mut f64,
    hex: &mut [f64],
) {
    let dims = &params.dims;
    let matkey = &params.matkey;

    if dims.nfreqe > 0 {
        for ij in 1..=dims.nfreqe {
            let ijt = params.ijfr[ij - 1];
            if !params.lskip[0][ijt] {
                let fluxw = params.w[ijt]
                    * (params.fh[ijt] * params.rad0[ij - 1] - params.hextrd[ijt]);
                *grd += fluxw * params.abso0[ij - 1];
                *hexn += fluxw * params.dabn0[ij - 1];
                *hext += fluxw * params.dabt0[ij - 1];
                for i in 0..dims.nlvexp {
                    hex[i] += fluxw * params.drch0[i][ij - 1];
                }
            }
        }
    }

    let ccc = PCK / params.ctrl.qgrav;
    let hr1 = ccc * (*grd + params.fprd[0]) / params.dens[0];
    let pg1 = BOLK * params.psi0[nhe] * params.temp[0];
    let hg1 = (TWO * pg1 / params.dens[0] / params.ctrl.qgrav).sqrt();
    let mut x = (params.zd[0] - hr1) / hg1;

    let f1 = if x < 3.0 {
        if x < 0.0 {
            x = 0.0;
        }
        erfcx(x)
    } else {
        HALF * (UN - HALF / (x * x)) / x
    };

    let x1 = x * 1.01;
    let f1d = if x1 < 3.0 {
        erfcx(x1)
    } else {
        HALF * (UN - HALF / (x1 * x1)) / x1
    };

    let f1d = if x > 0.0 { (f1d - f1) * 100.0 / x } else { 0.0 };
    let ggg = params.dens[0] * hg1 * f1;
    let rf1 = params.dens[0] * f1d;
    let ccd = ccc * f1d;

    for ij in 1..=dims.nfreqe {
        state.b[nhe][ij - 1] = -ccd * params.wdep0[ij - 1] * params.fh[params.ijfr[ij - 1]] * params.abso0[ij - 1];
    }

    state.b[nhe][nhe] += (ggg + hr1 * rf1) / params.psi0[nhe];

    if matkey.inre > 0 {
        state.b[nhe][nre] = (ggg - rf1 * params.zd[0] + rf1 * hr1) * HALF / params.temp[0]
            - ccd * (*hext + params.heit[0]);
    }

    if matkey.inzd > 0 {
        let nzd = dims.nfreqe + matkey.inzd as usize;
        state.b[nhe][nzd] = rf1;
    }

    if matkey.inpc > 0 {
        state.b[nhe][npc] = -ccd * (*hexn + params.hein[0]);
    }

    for ii in 0..dims.nlvexp {
        state.b[nhe][nse + ii] = -ccd * (hex[ii] + params.heip[ii][0]);
    }

    state.vecl[nhe] = params.dm[0] - ggg;
}

/// 填充盘模型上边界条件 (旧版)
fn fill_upper_boundary_disk_old(
    params: &BheParams,
    state: &mut BheState,
    nhe: usize,
    nre: usize,
    grd: &mut f64,
) {
    let dims = &params.dims;
    let matkey = &params.matkey;

    if dims.nfreqe > 0 {
        for ij in 1..=dims.nfreqe {
            let ijt = params.ijfr[ij - 1];
            if !params.lskip[0][ijt] {
                let fluxw = params.w[ijt]
                    * (params.fh[ijt] * params.rad0[ij - 1] - params.hextrd[ijt]);
                *grd += fluxw * params.abso0[ij - 1];
            }
        }
    }

    let ccc = PCK / params.ctrl.qgrav;
    let pr1 = ccc * (*grd + params.fprd[0]) / params.dens[0];
    let pg1 = BOLK * params.psi0[nhe] * params.temp[0];
    let hg1 = (TWO * pg1 / params.dens[0] / params.ctrl.qgrav).sqrt();
    let mut x = (params.zd[0] - pr1) / hg1;

    let f1 = if x < 3.0 {
        if x < 0.0 {
            x = 0.0;
        }
        erfcx(x)
    } else {
        HALF * (UN - HALF / (x * x)) / x
    };

    let ggg = hg1 * params.ctrl.qgrav * HALF / f1;

    state.b[nhe][nhe] = BOLK * params.temp[0];
    if matkey.inre > 0 {
        state.b[nhe][dims.nfreqe + matkey.inre as usize] = pg1 / params.temp[0];
    }

    state.vecl[nhe] = params.dm[0] * ggg - pg1;
}

/// 填充简单气压边界条件
fn fill_upper_boundary_simple(
    params: &BheParams,
    state: &mut BheState,
    nhe: usize,
    nre: usize,
) {
    let matkey = &params.matkey;

    state.b[nhe][nhe] = BOLK * params.temp[0];
    if matkey.inre > 0 {
        state.b[nhe][nre] = BOLK * params.psi0[nhe];
    }
    state.vecl[nhe] = params.ctrl.pgas0 - BOLK * params.temp[0] * params.psi0[nhe];
}

/// 填充 BHED 正常深度点
fn fill_normal_depth_bhed(
    params: &BheParams,
    state: &mut BheState,
    id: usize,
    nhe: usize,
    nre: usize,
    npc: usize,
    nse: usize,
    grd: &mut f64,
) {
    let dims = &params.dims;
    let matkey = &params.matkey;

    if dims.nfreqe > 0 {
        for ij in 1..=dims.nfreqe {
            if !params.lskip[id - 1][params.ijfr[ij - 1]] {
                *grd += (params.fk0[ij - 1] * params.rad0[ij - 1]
                    - params.fkm[ij - 1] * params.radm[ij - 1])
                    * params.w[params.ijfr[ij - 1]];
                state.a[nhe][ij - 1] = -PCK * params.w[params.ijfr[ij - 1]] * params.fkm[ij - 1];
                state.b[nhe][ij - 1] = PCK * params.w[params.ijfr[ij - 1]] * params.fk0[ij - 1];
            }
        }
    }

    let vt0 = HALF * params.vturb[id - 1] * params.vturb[id - 1] * params.wmm[id - 1];
    let vtm = HALF * params.vturb[id - 2] * params.vturb[id - 2] * params.wmm[id - 2];

    state.a[nhe][nhe] = -BOLK * params.temp[id - 2] - GN * vtm;
    state.b[nhe][nhe] = BOLK * params.temp[id - 1] + GN * vt0;

    if matkey.inre > 0 {
        state.a[nhe][nre] = -BOLK * params.psim[nhe] + PCK * params.heitm[id - 1];
        state.b[nhe][nre] = BOLK * params.psi0[nhe] + PCK * params.heit[id - 1];
        state.c[nhe][nre] = PCK * params.heitp[id - 1];
    }

    if matkey.inpc > 0 {
        state.a[nhe][npc] = GN * vtm + PCK * params.heinm[id - 1];
        state.b[nhe][npc] = -GN * vt0 + PCK * params.hein[id - 1];
        state.c[nhe][npc] = PCK * params.heinp[id - 1];
    }

    if matkey.inmp > 0 {
        let nmp = dims.nfreqe + matkey.inmp as usize;
        state.a[nhe][nmp] = -GP * vtm;
        state.b[nhe][nmp] = GP * vt0;
    }

    if matkey.inzd > 0 {
        let nzd = dims.nfreqe + matkey.inzd as usize;
        let dgrav = -params.ctrl.qgrav * (params.dm[id - 1] - params.dm[id - 2]) * HALF;
        state.a[nhe][nzd] = dgrav;
        state.b[nhe][nzd] = dgrav;
    }

    for ii in 0..dims.nlvexp {
        state.a[nhe][nse + ii] += PCK * params.heipm[ii][id - 1];
        state.b[nhe][nse + ii] += PCK * params.heip[ii][id - 1];
        state.c[nhe][nse + ii] += PCK * params.heipp[ii][id - 1];
    }

    let grav = params.ctrl.qgrav * (params.zd[id - 1] + params.zd[id - 2]) * HALF;
    state.vecl[nhe] = grav * (params.dm[id - 1] - params.dm[id - 2])
        - BOLK * (params.temp[id - 1] * params.psi0[nhe]
            - params.temp[id - 2] * params.psim[nhe])
        - PCK * (*grd + params.fprd[id - 1])
        - vt0 / params.wmm[id - 1] * params.dens[id - 1]
        + vtm / params.wmm[id - 2] * params.dens[id - 2];
}

/// 填充 z-m 关系方程
///
/// # Fortran 索引说明
///
/// Fortran 使用 1-indexed，此函数中的 `id` 参数也是 1-indexed。
///
/// Fortran 代码使用向前差分:
/// - `DDP = (DM(ID+1) - DM(ID)) * HALF`
/// - `VECL(NZD) = ZD(ID+1) - ZD(ID) + ...`
///
/// 转换为 Rust 0-indexed:
/// - Fortran `DM(ID)` → Rust `dm[id-1]`
/// - Fortran `DM(ID+1)` → Rust `dm[id]`
fn fill_zm_relation(
    params: &BheParams,
    state: &mut BheState,
    id: usize,
    dims: &ModelDims,
    matkey: &MatKey,
) {
    if matkey.inzd <= 0 {
        return;
    }

    let nzd = dims.nfreqe + matkey.inzd as usize;

    // 下边界条件 z(ND) = 0
    state.b[nzd][nzd] = UN;
    if id == dims.nd {
        return;
    }

    // 正常深度点 - 使用向前差分
    // Fortran: DDP = (DM(ID+1) - DM(ID)) * HALF
    // Rust: ddp = (dm[id] - dm[id-1]) * HALF  (因为 Fortran ID 是 1-indexed)
    let ddp = (params.dm[id] - params.dm[id - 1]) * HALF;

    state.b[nzd][nzd] = UN;
    state.czz = -UN;

    // Fortran: X1 = GN * WMM(ID) * DDP
    // Rust: wmm[id-1] (0-indexed)
    let x1 = GN * params.wmm[id - 1] * ddp;

    if matkey.inhe > 0 {
        let nhe = dims.nfreqe + matkey.inhe as usize;
        // Fortran: B(NZD,NHE) = X1 / DENS(ID) / DENS(ID)
        // Fortran: CZN = X1 / DENS(ID+1) / DENS(ID+1)
        state.b[nzd][nhe] = x1 / params.dens[id - 1] / params.dens[id - 1];
        state.czn = x1 / params.dens[id] / params.dens[id];
    }

    if matkey.inpc > 0 {
        let npc = dims.nfreqe + matkey.inpc as usize;
        state.b[nzd][npc] = -x1 / params.dens[id - 1] / params.dens[id - 1];
        state.cze = -x1 / params.dens[id] / params.dens[id];
    }

    if matkey.inmp > 0 {
        let nmp = dims.nfreqe + matkey.inmp as usize;
        // Fortran: PSI0(NFREQE+INMP) 和 PSIP(NFREQE+INMP)
        state.b[nzd][nmp] = ddp / params.dens[id - 1] / params.psi0[nmp];
        state.czm = ddp / params.dens[id] / params.psip[nmp];
    }

    // Fortran: VECL(NZD) = ZD(ID+1) - ZD(ID) + DDP/DENS(ID) + DDP/DENS(ID+1)
    state.vecl[nzd] = params.zd[id] - params.zd[id - 1]
        + ddp / params.dens[id - 1]
        + ddp / params.dens[id];
}

// ============================================================================
// BHEZ 函数
// ============================================================================

/// 流体静力学平衡方程矩阵填充 (z 坐标版本)。
///
/// 专门用于 z 坐标情况，
/// 边界条件处理略有不同。
pub fn bhez(params: &BheParams, state: &mut BheState) {
    let id = params.id;
    let dims = &params.dims;
    let matkey = &params.matkey;

    let nhe = dims.nfreqe + matkey.inhe as usize;
    let nre = dims.nfreqe + matkey.inre as usize;
    let npc = dims.nfreqe + matkey.inpc as usize;
    let nzd = dims.nfreqe + matkey.inzd as usize;
    let nse = dims.nfreqe + matkey.inse as usize - 1;

    if matkey.inhe <= 0 {
        return;
    }

    // 虚拟大质量粒子密度方程
    if matkey.inmp > 0 {
        let nmp = dims.nfreqe + matkey.inmp as usize;
        state.b[nmp][nmp] = -UN;
        state.b[nmp][nhe] = UN;
        if matkey.inpc > 0 {
            state.b[nmp][npc] = -UN;
        }
    }

    // 初始化
    let mut hext: f64 = 0.0;
    let mut hexn: f64 = 0.0;
    let mut grd: f64 = 0.0;
    let mut hex = vec![0.0; dims.nlvexp];

    if id > 1 {
        // 正常深度点
        fill_normal_depth_bhez(params, state, id, nhe, nre, npc, nse, nzd, &mut grd);
        return;
    }

    // === 上边界条件 (ID = 1) ===
    let ibche = params.ctrl.ibche;

    if ibche == 0 {
        // 恒星大气边界条件
        fill_upper_boundary_stellar_bhez(params, state, nhe, nre, npc, nse, nzd, &mut hext, &mut hexn, &mut grd, &mut hex);
    } else if ibche == 1 {
        // 盘模型边界条件 (新版)
        fill_upper_boundary_disk_new_bhez(params, state, nhe, nre, npc, nse, nzd, &mut hext, &mut hexn, &mut grd, &mut hex);
    } else if ibche == 2 {
        // 盘模型边界条件 (旧版)
        fill_upper_boundary_disk_old_bhez(params, state, nhe, nre, nzd, &mut grd);
    } else {
        // 简单气压边界条件
        fill_upper_boundary_simple(params, state, nhe, nre);
    }
}

/// BHEZ 恒星大气上边界
fn fill_upper_boundary_stellar_bhez(
    params: &BheParams,
    state: &mut BheState,
    nhe: usize,
    nre: usize,
    npc: usize,
    nse: usize,
    nzd: usize,
    hext: &mut f64,
    hexn: &mut f64,
    grd: &mut f64,
    hex: &mut [f64],
) {
    // 与 BHED 的恒星大气边界相同
    fill_upper_boundary_stellar(params, state, nhe, nre, npc, nse, hext, hexn, grd, hex);
}

/// BHEZ 盘模型上边界 (新版)
fn fill_upper_boundary_disk_new_bhez(
    params: &BheParams,
    state: &mut BheState,
    nhe: usize,
    nre: usize,
    npc: usize,
    nse: usize,
    nzd: usize,
    hext: &mut f64,
    hexn: &mut f64,
    grd: &mut f64,
    hex: &mut [f64],
) {
    let dims = &params.dims;
    let matkey = &params.matkey;

    if dims.nfreqe > 0 {
        for ij in 1..=dims.nfreqe {
            let ijt = params.ijfr[ij - 1];
            if !params.lskip[0][ijt] {
                let fluxw = params.w[ijt]
                    * (params.fh[ijt] * params.rad0[ij - 1] - params.hextrd[ijt]);
                *grd += fluxw * params.abso0[ij - 1];
                *hexn += fluxw * params.dabn0[ij - 1];
                *hext += fluxw * params.dabt0[ij - 1];
                for i in 0..dims.nlvexp {
                    hex[i] += fluxw * params.drch0[i][ij - 1];
                }
            }
        }
    }

    let ccc = PCK / params.ctrl.qgrav;
    let hr1 = ccc * (*grd + params.fprd[0]) / params.dens[0];
    let pg1 = BOLK * params.psi0[nhe] * params.temp[0];
    let hg1 = (TWO * pg1 / params.dens[0] / params.ctrl.qgrav).sqrt();
    let mut x = (params.zd[0] - hr1) / hg1;

    let f1 = if x < 3.0 {
        if x < 0.0 {
            x = 0.0;
        }
        erfcx(x)
    } else {
        HALF * (UN - HALF / (x * x)) / x
    };

    let x1 = x * 1.01;
    let f1d = if x1 < 3.0 {
        erfcx(x1)
    } else {
        HALF * (UN - HALF / (x1 * x1)) / x1
    };

    let f1d = if x > 0.0 { (f1d - f1) * 100.0 / x } else { 0.0 };
    let ggg = params.dens[0] * hg1 * f1;
    let rf1 = params.dens[0] * f1d;
    let ccd = ccc * f1d;

    for ij in 1..=dims.nfreqe {
        state.b[nhe][ij - 1] = -ccd * params.wdep0[ij - 1] * params.fh[params.ijfr[ij - 1]] * params.abso0[ij - 1];
    }

    state.b[nhe][nhe] += (ggg + hr1 * rf1) / params.psi0[nhe];

    if matkey.inre > 0 {
        state.b[nhe][nre] = (ggg - rf1 * params.zd[0] + rf1 * hr1) * HALF / params.temp[0]
            - ccd * (*hext + params.heit[0]);
    }

    // Fortran: IF(INZD.GT.0) B(NHE,NZD)=RF1
    if matkey.inzd > 0 {
        state.b[nhe][nzd] = rf1;
    }

    if matkey.inpc > 0 {
        state.b[nhe][npc] = -ccd * (*hexn + params.hein[0]);
    }

    for ii in 0..dims.nlvexp {
        state.b[nhe][nse + ii] = -ccd * (hex[ii] + params.heip[ii][0]);
    }

    state.vecl[nhe] = params.dm[0] - ggg;
}

/// BHEZ 盘模型上边界 (旧版)
fn fill_upper_boundary_disk_old_bhez(
    params: &BheParams,
    state: &mut BheState,
    nhe: usize,
    nre: usize,
    nzd: usize,
    grd: &mut f64,
) {
    let dims = &params.dims;
    let matkey = &params.matkey;

    if dims.nfreqe > 0 {
        for ij in 1..=dims.nfreqe {
            let ijt = params.ijfr[ij - 1];
            if !params.lskip[0][ijt] {
                let fluxw = params.w[ijt]
                    * (params.fh[ijt] * params.rad0[ij - 1] - params.hextrd[ijt]);
                *grd += fluxw * params.abso0[ij - 1];
            }
        }
    }

    let ccc = PCK / params.ctrl.qgrav;
    let pr1 = ccc * (*grd + params.fprd[0]) / params.dens[0];
    let pg1 = BOLK * params.psi0[nhe] * params.temp[0];
    let hg1 = (TWO * pg1 / params.dens[0] / params.ctrl.qgrav).sqrt();
    let mut x = (params.zd[0] - pr1) / hg1;

    let f1 = if x < 3.0 {
        if x < 0.0 {
            x = 0.0;
        }
        erfcx(x)
    } else {
        HALF * (UN - HALF / (x * x)) / x
    };

    let ggg = hg1 * params.ctrl.qgrav * HALF / f1;

    state.b[nhe][nhe] = BOLK * params.temp[0];
    if matkey.inre > 0 {
        state.b[nhe][dims.nfreqe + matkey.inre as usize] = pg1 / params.temp[0];
    }

    state.vecl[nhe] = params.dm[0] * ggg - pg1;
}

/// BHEZ 正常深度点
fn fill_normal_depth_bhez(
    params: &BheParams,
    state: &mut BheState,
    id: usize,
    nhe: usize,
    nre: usize,
    npc: usize,
    nse: usize,
    nzd: usize,
    grd: &mut f64,
) {
    let dims = &params.dims;
    let matkey = &params.matkey;

    let grav = params.ctrl.qgrav * (params.zd[id - 1] + params.zd[id - 2]) * HALF;
    let gravz = grav * (params.zd[id - 1] - params.zd[id - 2]);
    let dgrv = gravz * HALF * params.wmm[id - 1];

    if dims.nfreqe > 0 {
        for ij in 1..=dims.nfreqe {
            if !params.lskip[id - 1][params.ijfr[ij - 1]] {
                *grd += (params.fk0[ij - 1] * params.rad0[ij - 1]
                    - params.fkm[ij - 1] * params.radm[ij - 1])
                    * params.w[params.ijfr[ij - 1]];
                state.a[nhe][ij - 1] = -PCK * params.w[params.ijfr[ij - 1]] * params.fkm[ij - 1];
                state.b[nhe][ij - 1] = PCK * params.w[params.ijfr[ij - 1]] * params.fk0[ij - 1];
            }
        }
    }

    let vt0 = HALF * params.vturb[id - 1] * params.vturb[id - 1] * params.wmm[id - 1];
    // Fortran: VTM=HALF*VTURB(ID-1)*VTURB(ID-1)*WMM(ID) — WMM 用 ID 不是 ID-1
    let vtm = HALF * params.vturb[id - 2] * params.vturb[id - 2] * params.wmm[id - 1];

    state.a[nhe][nhe] = -BOLK * params.temp[id - 2] - GN * (vtm + dgrv);
    state.b[nhe][nhe] = BOLK * params.temp[id - 1] + GN * (vt0 + dgrv);

    if matkey.inre > 0 {
        state.a[nhe][nre] = -BOLK * params.psim[nhe] + PCK * params.heitm[id - 1];
        state.b[nhe][nre] = BOLK * params.psi0[nhe] + PCK * params.heit[id - 1];
        state.c[nhe][nre] = PCK * params.heitp[id - 1];
    }

    if matkey.inpc > 0 {
        state.a[nhe][npc] = GN * (vtm + dgrv) + PCK * params.heinm[id - 1];
        state.b[nhe][npc] = -GN * (vt0 + dgrv) + PCK * params.hein[id - 1];
        state.c[nhe][npc] = PCK * params.heinp[id - 1];
    }

    if matkey.inmp > 0 {
        let nmp = dims.nfreqe + matkey.inmp as usize;
        state.a[nhe][nmp] = -GP * (vtm + dgrv);
        state.b[nhe][nmp] = GP * (vt0 + dgrv);
    }

    for ii in 0..dims.nlvexp {
        state.a[nhe][nse + ii] += PCK * params.heipm[ii][id - 1];
        state.b[nhe][nse + ii] += PCK * params.heip[ii][id - 1];
        state.c[nhe][nse + ii] += PCK * params.heipp[ii][id - 1];
    }

    state.vecl[nhe] = -gravz * (params.dens[id - 1] + params.dens[id - 2]) * HALF
        - BOLK * (params.temp[id - 1] * params.psi0[nhe]
            - params.temp[id - 2] * params.psim[nhe])
        - PCK * (*grd + params.fprd[id - 1])
        - vt0 / params.wmm[id - 1] * params.dens[id - 1]
        + vtm / params.wmm[id - 2] * params.dens[id - 2];
}

// ============================================================================
// 测试
// ============================================================================

#[cfg(test)]
mod tests {
    use super::*;

    fn create_test_params() -> BheParams {
        let mut params = BheParams::default();
        params.id = 1;
        params.dims = ModelDims {
            nfreqe: 10,
            nd: 5,
            nlvexp: 3,
        };
        params.matkey = MatKey {
            nn: 20,
            nn0: 20,
            inhe: 10,
            inre: 11,
            inpc: 12,
            inse: 14,
            inzd: 0,
            inmp: 0,
            ndre: 5,
        };
        params.ctrl = InputControl {
            ibche: 0,
            ifixde: 0,
            ifprad: 1,
            pgas0: 0.0,
            qgrav: 1.0,
        };

        // 设置简单的测试值
        params.temp[0] = 10000.0;
        params.dens[0] = 1e12;
        params.totn[0] = 1e12;
        params.elec[0] = 1e10;
        params.wmm[0] = 1.0;
        params.dm[0] = 1e-4;
        params.vturb[0] = 5e5;
        params.zd[0] = 1e8;

        for i in 0..10 {
            params.w[i] = 1.0;
            params.ijfr[i] = i;
            params.fh[i] = 0.5;
            params.hextrd[i] = 0.0;
            params.rad0[i] = 1e10;
            params.abso0[i] = 1e-8;
            params.dabt0[i] = 1e-10;
            params.dabn0[i] = 1e-10;
            params.wdep0[i] = 1.0;
        }

        for i in 0..3 {
            for j in 0..10 {
                params.drch0[i][j] = 1e-12;
            }
            params.heip[i][0] = 0.0;
            params.heipp[i][0] = 0.0;
        }

        params.heit[0] = 0.0;
        params.hein[0] = 0.0;
        params.heitp[0] = 0.0;
        params.heinp[0] = 0.0;
        params.fprd[0] = 0.0;
        params.psi0[20] = 1e12; // nhe index

        params
    }

    #[test]
    fn test_bhe_fixed_density() {
        let mut params = create_test_params();
        params.ctrl.ifixde = 1;
        params.ctrl.ifprad = 0;

        let mut state = BheState::new();
        bhe(&params, &mut state);

        let nhe = params.dims.nfreqe + params.matkey.inhe as usize;
        let npc = params.dims.nfreqe + params.matkey.inpc as usize;

        assert!((state.b[nhe][nhe] - UN).abs() < 1e-10);
        assert!((state.b[nhe][npc] - (-UN)).abs() < 1e-10);
    }

    #[test]
    fn test_bhe_upper_boundary() {
        let params = create_test_params();
        let mut state = BheState::new();
        bhe(&params, &mut state);

        let nhe = params.dims.nfreqe + params.matkey.inhe as usize;

        // 检查矩阵元素被填充
        assert!(state.b[nhe][nhe].abs() > 0.0 || state.vecl[nhe].abs() > 0.0);
    }

    #[test]
    fn test_bhed_simple() {
        let params = create_test_params();
        let mut state = BheState::new();
        bhed(&params, &mut state);

        let nhe = params.dims.nfreqe + params.matkey.inhe as usize;

        // 检查矩阵元素被填充
        assert!(state.b[nhe][nhe].abs() > 0.0 || state.vecl[nhe].abs() > 0.0);
    }

    #[test]
    fn test_bhez_simple() {
        let params = create_test_params();
        let mut state = BheState::new();
        bhez(&params, &mut state);

        let nhe = params.dims.nfreqe + params.matkey.inhe as usize;

        // 检查矩阵元素被填充
        assert!(state.b[nhe][nhe].abs() > 0.0 || state.vecl[nhe].abs() > 0.0);
    }

    #[test]
    fn test_erfcx() {
        // 测试 erfcx 函数
        let x = 1.0;
        let result = erfcx(x);
        assert!(result > 0.0);
        assert!(result < 2.0);
    }
}