SpectraRust/src/tlusty/math/bre.rs

//! 辐射平衡方程矩阵计算。
//!
//! 重构自 TLUSTY `bre.f`
//!
//! 计算辐射平衡方程对应的矩阵 A, B, C 部分（第 NFREQE+INRE 行）。
//! 包含积分方程部分和微分方程部分。

use crate::tlusty::state::constants::{HALF, SIG4P, UN};

// ============================================================================
// 常量
// ============================================================================

/// 最大线性化能级数
const MLVEXP: usize = 233;

// ============================================================================
// BRE 参数结构体
// ============================================================================

/// BRE 输入参数。
///
/// 包含所有来自 COMMON 块的必要数据。
#[derive(Debug)]
pub struct BreParams<'a> {
    // ==================== 基本参数 ====================
    /// 深度索引 ID (1-indexed)
    pub id: usize,

    // ==================== 维度参数 ====================
    /// 深度点数 ND
    pub nd: usize,
    /// 线性化频率数 NFREQE
    pub nfreqe: usize,
    /// 总频率数 NFREQ
    pub nfreq: usize,
    /// 线性化能级数 NLVEXP
    pub nlvexp: usize,

    // ==================== 索引参数 ====================
    /// 氦索引 INHE (0 表示无氦)
    pub inhe: usize,
    /// 辐射平衡索引 INRE
    pub inre: usize,
    /// 电子密度索引 INPC
    pub inpc: usize,
    /// 质量密度索引 INMP
    pub inmp: usize,

    // ==================== 控制参数 ====================
    /// Compton 散射标志
    pub icompt: i32,
    /// Compton 边界条件标志
    pub icombc: i32,
    /// Compton 密度导数标志
    pub icmdra: i32,
    /// 迭代次数
    pub iter: i32,
    /// 辐射平衡温度控制
    pub nretc: i32,
    /// 盘模型标志 (0=无盘, 1=有盘)
    pub idisk: i32,
    /// ALI 标志 (>5 启用完整 ALI)
    pub ifali: i32,

    // ==================== 频率映射 ====================
    /// 频率索引映射 [MFREQE]
    pub ijfr: &'a [usize],
    /// 频率起源索引 [MFREQ]
    pub ijorig: &'a [usize],
    /// 频率反转映射 [MFREQ]
    pub kij: &'a [usize],
    /// ALI 频率索引 [MFREQ]
    pub ijex: &'a [i32],

    // ==================== 模型状态 ====================
    /// 温度 [MDEPTH] (K)
    pub temp: &'a [f64],
    /// 电子密度 [MDEPTH]
    pub elec: &'a [f64],
    /// 柱质量密度 [MDEPTH]
    pub dm: &'a [f64],
    /// 密度倒数 [MDEPTH]
    pub dens1: &'a [f64],
    /// 分子权重 [MDEPTH]
    pub wmm: &'a [f64],
    /// 粘性加热 [MDEPTH]
    pub tvisc: &'a [f64],
    /// 辐射积分因子 [MDEPTH]
    pub reint: &'a [f64],
    /// THETAV 速度场 [MDEPTH]
    pub thetav: &'a [f64],

    // ==================== 辐射场 ====================
    /// 当前深度辐射 [MFREQE]
    pub rad0: &'a [f64],
    /// 前深度辐射 [MFREQE]
    pub radm: &'a [f64],
    /// FK 当前 [MFREQE]
    pub fk0: &'a [f64],
    /// FK 前深度 [MFREQE]
    pub fkm: &'a [f64],

    // ==================== 吸收/发射系数 ====================
    /// 吸收系数当前 [MFREQE]
    pub abso0: &'a [f64],
    /// 吸收系数前深度 [MFREQE]
    pub absom: &'a [f64],
    /// 发射系数当前 [MFREQE]
    pub emis0: &'a [f64],
    /// 发射系数前深度 [MFREQE]
    pub emism: &'a [f64],
    /// 散射系数当前 [MFREQE]
    pub scat0: &'a [f64],

    // ==================== 导数 ====================
    /// 吸收系数 T 导数当前 [MFREQE]
    pub dabt0: &'a [f64],
    /// 吸收系数 T 导数前深度 [MFREQE]
    pub dabtm: &'a [f64],
    /// 吸收系数 N 导数当前 [MFREQE]
    pub dabn0: &'a [f64],
    /// 吸收系数 N 导数前深度 [MFREQE]
    pub dabnm: &'a [f64],
    /// 吸收系数 M 导数当前 [MFREQE]
    pub dabm0: &'a [f64],
    /// 发射系数 T 导数当前 [MFREQE]
    pub demt0: &'a [f64],
    /// 发射系数 T 导数前深度 [MFREQE]
    pub demtm: &'a [f64],
    /// 发射系数 N 导数当前 [MFREQE]
    pub demn0: &'a [f64],
    /// 发射系数 N 导数前深度 [MFREQE]
    pub demnm: &'a [f64],
    /// 发射系数 M 导数当前 [MFREQE]
    pub demm0: &'a [f64],

    // ==================== 能级导数 ====================
    /// 能级导数当前 [MLVEXP × MFREQE]
    pub drch0: &'a [Vec<f64>],
    /// 能级导数前深度 [MLVEXP × MFREQE]
    pub drchm: &'a [Vec<f64>],
    /// 能级发射导数当前 [MLVEXP × MFREQE]
    pub dret0: &'a [Vec<f64>],

    // ==================== 深度权重 ====================
    /// 深度权重当前 [MFREQE]
    pub wdep0: &'a [f64],

    // ==================== 冷却率 ====================
    /// ALI 冷却率 [MDEPTH]
    pub fcool: &'a [f64],

    // ==================== ALI 辐射等效项 ====================
    /// REIT [MDEPTH]
    pub reit: &'a [f64],
    /// REIN [MDEPTH]
    pub rein: &'a [f64],
    /// REIM [MDEPTH]
    pub reim: &'a [f64],
    /// REIX [MDEPTH]
    pub reix: &'a [f64],
    /// REIP [MLVEXP × MDEPTH]
    pub reip: &'a [Vec<f64>],

    // ==================== ALI A 矩阵项 ====================
    /// AREIT [MDEPTH]
    pub areit: &'a [f64],
    /// AREIN [MDEPTH]
    pub arein: &'a [f64],
    /// AREIM [MDEPTH]
    pub areim: &'a [f64],
    /// AREIP [MLVEXP × MDEPTH]
    pub areip: &'a [Vec<f64>],

    // ==================== ALI C 矩阵项 ====================
    /// CREIT [MDEPTH]
    pub creit: &'a [f64],
    /// CREIN [MDEPTH]
    pub crein: &'a [f64],
    /// CREIM [MDEPTH]
    pub creim: &'a [f64],
    /// CREIX [MDEPTH]
    pub creix: &'a [f64],
    /// CREIP [MLVEXP × MDEPTH]
    pub creip: &'a [Vec<f64>],

    // ==================== RED 微分方程项 ====================
    /// REDIF [MDEPTH]
    pub redif: &'a [f64],
    /// REDT [MDEPTH]
    pub redt: &'a [f64],
    /// REDTP [MDEPTH]
    pub redtp: &'a [f64],
    /// REDX [MDEPTH]
    pub redx: &'a [f64],
    /// REDXP [MDEPTH]
    pub redxp: &'a [f64],
    /// REDN [MDEPTH]
    pub redn: &'a [f64],
    /// REDP [MLVEXP × MDEPTH]
    pub redp: &'a [Vec<f64>],

    // ==================== RED 前深度项 ====================
    /// REDTM [MDEPTH]
    pub redtm: &'a [f64],
    /// REDXM [MDEPTH]
    pub redxm: &'a [f64],
    /// REDNM [MDEPTH]
    pub rednm: &'a [f64],
    /// REDPM [MLVEXP × MDEPTH]
    pub redpm: &'a [Vec<f64>],

    // ==================== RED 后深度项 ====================
    /// REDNP [MDEPTH]
    pub rednp: &'a [f64],
    /// REDPP [MLVEXP × MDEPTH]
    pub redpp: &'a [Vec<f64>],

    // ==================== 盘模型项 ====================
    /// DTVIST [MDEPTH]
    pub dtvist: &'a [f64],
    /// DTVISR [MDEPTH]
    pub dtvisr: &'a [f64],
    /// DTVISN [MDEPTH]
    pub dtvisn: &'a [f64],

    // ==================== 边界条件 ====================
    /// FH [MFREQ]
    pub fh: &'a [f64],
    /// HEXTRD [MFREQ]
    pub hextrd: &'a [f64],

    // ==================== 有效温度 ====================
    pub teff: f64,
    /// 氢原子质量
    pub hmass: f64,

    // ==================== 电子散射截面 ====================
    /// SIGEC [MFREQ]
    pub sigec: &'a [f64],
    /// 汤姆逊散射截面
    pub sige: f64,
    /// CMD - Compton 密度导数
    pub cmd: f64,
}

// ============================================================================
// BRE 可变状态结构体
// ============================================================================

/// BRE 可变状态（矩阵和向量）。
#[derive(Debug)]
pub struct BreState<'a> {
    /// A 矩阵 [MTOT × MTOT]
    pub a: &'a mut [Vec<f64>],
    /// B 矩阵 [MTOT × MTOT]
    pub b: &'a mut [Vec<f64>],
    /// C 矩阵 [MTOT × MTOT]
    pub c: &'a mut [Vec<f64>],
    /// 左向量 VECL [MTOT]
    pub vecl: &'a mut [f64],
    /// REX 辅助数组 [MLEVEL]
    pub rex: &'a mut [f64],
}

// ============================================================================
// Compton 辅助计算
// ============================================================================

/// 计算 Compton 散射辅助量（简化版）。
fn compt0_bre(
    _ij: usize,
    _id: usize,
    ab: f64,
    nfreq: usize,
    kij: &[usize],
    elec: &[f64],
    sige: f64,
    ij_idx: usize,
    id_idx: usize,
) -> (f64, f64, f64, f64, f64, f64) {
    // IJI = NFREQ - KIJ(IJ) + 1
    let iji = nfreq - kij[ij_idx] + 1;

    if iji == 1 {
        return (0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
    }

    // 简化计算 - 完整实现需要调用 compt0 函数
    let ss0 = elec[id_idx] * sige / ab;

    // 返回 (CMA, CMB, CMC, CME, CMS, CMD)
    (0.0, 0.0, 0.0, 0.0, ss0, 0.0)
}

// ============================================================================
// BRE 主函数
// ============================================================================

/// 计算辐射平衡方程的矩阵 A, B, C 部分。
///
/// # 参数
///
/// * `params` - 输入参数
/// * `state` - 可变状态（矩阵 A, B, C, VECL）
///
/// # 说明
///
/// 此函数修改矩阵 A, B, C 的第 (NFREQE+INRE) 行，
/// 对应辐射平衡方程的积分部分和微分部分。
pub fn bre(params: &BreParams, state: &mut BreState) {
    let id = params.id;
    let id_idx = id - 1; // 0-indexed

    // 计算矩阵列索引
    let nhe = params.nfreqe + params.inhe;
    let nre = params.nfreqe + params.inre;
    let npc = params.nfreqe + params.inpc;
    let nmp = params.nfreqe + params.inmp;
    let nse = params.nfreqe + params.inse() - 1;

    // IJ1 = 1 或 2 (Compton 散射时从 2 开始)
    let mut ij1 = 1;
    if params.icompt > 0 && params.icombc > 0 && !params.ijex.is_empty() && params.ijex[0] > 0 {
        ij1 = 2;
    }

    // 检查温度控制
    let ittc = (params.nretc.abs() / 100) as i32;
    if params.iter > ittc {
        let mod_val = (params.nretc.abs() % 100) as usize;
        if id <= mod_val {
            state.b[nre - 1][nre - 1] = 1.0;
            if params.nretc < 0 {
                state.c[nre - 1][nre - 1] = -1.0;
                if id_idx + 1 < params.temp.len() {
                    state.vecl[nre - 1] = params.temp[id_idx + 1] - params.temp[id_idx];
                }
            }
            return;
        }
    }

    // RHS 向量初始化（ALI 冷却率）
    state.vecl[nre - 1] = params.fcool[id_idx];
    if params.idisk == 1 {
        state.vecl[nre - 1] -= params.reint[id_idx] * params.tvisc[id_idx];
    }

    if params.reint[id_idx] <= 0.0 {
        // 跳转到微分方程部分
        bre_differential(params, state, id, nre, nhe, npc, nmp, nse);
        return;
    }

    // ==================== 积分方程部分 ====================

    let mut brepc = 0.0;
    let mut bremp = 0.0;

    // 初始化 REX
    for i in 0..params.nlvexp.min(state.rex.len()) {
        state.rex[i] = 0.0;
    }

    if params.nfreqe > 0 {
        for ij in ij1..=params.nfreqe {
            let ij_idx = ij - 1;
            let ijt = params.ijfr[ij_idx];

            // 累积积分项
            let sigec_val = if ijt > 0 && ijt <= params.sigec.len() {
                params.sigec[ijt - 1]
            } else {
                0.0
            };

            brepc += ((params.dabn0[ij_idx] - sigec_val) * params.rad0[ij_idx]
                - params.demn0[ij_idx])
                * params.wdep0[ij_idx];
            bremp += (params.dabm0[ij_idx] * params.rad0[ij_idx] - params.demm0[ij_idx])
                * params.wdep0[ij_idx];

            for i in 0..params.nlvexp.min(state.rex.len()) {
                state.rex[i] += (params.drch0[i][ij_idx] * params.rad0[ij_idx]
                    - params.dret0[i][ij_idx])
                    * params.wdep0[ij_idx];
            }

            // B 矩阵对角项
            state.b[nre - 1][nre - 1] += (params.dabt0[ij_idx] * params.rad0[ij_idx]
                - params.demt0[ij_idx])
                * params.wdep0[ij_idx]
                * params.reint[id_idx];

            // 加热项
            let heat = params.abso0[ij_idx] - params.scat0[ij_idx];
            state.b[nre - 1][ij - 1] = params.wdep0[ij_idx] * heat * params.reint[id_idx];

            // RHS 向量
            state.vecl[nre - 1] -= (heat * params.rad0[ij_idx] - params.emis0[ij_idx])
                * params.wdep0[ij_idx]
                * params.reint[id_idx];

            // Compton 散射项
            if params.icompt > 5 {
                let (_cma, cmb, _cmc, cme, cms, _cmd) = compt0_bre(
                    ijt,
                    id,
                    params.abso0[ij_idx],
                    params.nfreq,
                    params.kij,
                    params.elec,
                    params.sige,
                    ij_idx,
                    id_idx,
                );

                state.vecl[nre - 1] +=
                    params.abso0[ij_idx] * cms * params.wdep0[ij_idx] * params.reint[id_idx];

                if params.icompt > 6 {
                    if params.icmdra > 0 {
                        state.b[nre - 1][ij - 1] -=
                            params.abso0[ij_idx] * (cmb + cme) * params.wdep0[ij_idx]
                                * params.reint[id_idx];
                    } else {
                        state.b[nre - 1][ij - 1] -=
                            params.abso0[ij_idx] * (cmb + cme) * params.reint[id_idx];
                    }
                    // 注：完整的 Compton 处理需要更多代码
                }
            }
        }
    }

    // ALI 修正项
    state.b[nre - 1][nre - 1] += params.reit[id_idx] * params.reint[id_idx];
    if params.inpc > 0 {
        state.b[nre - 1][npc - 1] += (brepc + params.rein[id_idx]) * params.reint[id_idx];
    }
    if params.inmp > 0 {
        state.b[nre - 1][nmp - 1] += (bremp + params.reim[id_idx]) * params.reint[id_idx];
    }
    if params.inhe > 0 {
        state.b[nre - 1][nhe - 1] = params.reix[id_idx] * params.reint[id_idx];
    }

    // IFALI > 5 时的 A 和 C 矩阵项
    if params.ifali > 5 {
        state.a[nre - 1][nre - 1] = params.areit[id_idx] * params.reint[id_idx];
        if params.inpc > 0 {
            state.a[nre - 1][npc - 1] = params.arein[id_idx] * params.reint[id_idx];
        }
        if params.inmp > 0 {
            state.a[nre - 1][nmp - 1] = params.areim[id_idx] * params.reint[id_idx];
        }
        state.c[nre - 1][nre - 1] = params.creit[id_idx] * params.reint[id_idx];
        if params.inpc > 0 {
            state.c[nre - 1][npc - 1] = params.crein[id_idx] * params.reint[id_idx];
        }
        if params.inmp > 0 {
            state.c[nre - 1][nmp - 1] = params.creim[id_idx] * params.reint[id_idx];
        }
        if params.inhe > 0 {
            state.c[nre - 1][nhe - 1] = params.creix[id_idx] * params.reint[id_idx];
        }
    }

    // 盘模型项
    if params.idisk == 1 {
        state.b[nre - 1][nre - 1] += params.dtvist[id_idx] * params.reint[id_idx];
        if params.inpc > 0 {
            state.b[nre - 1][npc - 1] -= params.dtvisr[id_idx] * params.reint[id_idx];
        }
        if params.inhe > 0 {
            state.b[nre - 1][nhe - 1] =
                (params.dtvisr[id_idx] + params.dtvisn[id_idx]) * params.reint[id_idx];
        }
        if params.inmp > 0 {
            state.b[nre - 1][nmp - 1] = params.dtvisr[id_idx] * params.hmass
                / params.wmm[id_idx]
                * params.reint[id_idx];
        }
    }

    // 能级相关项
    for ii in 0..params.nlvexp {
        if nse + ii < state.b[nre - 1].len() {
            state.b[nre - 1][nse + ii] +=
                (state.rex[ii] + params.reip[ii][id_idx]) * params.reint[id_idx];
        }
    }

    if params.ifali > 5 && id > 1 {
        for ii in 0..params.nlvexp {
            if nse + ii < state.a[nre - 1].len() {
                state.a[nre - 1][nse + ii] += params.areip[ii][id_idx] * params.reint[id_idx];
            }
        }
    }

    if params.ifali > 5 && id < params.nd {
        for ii in 0..params.nlvexp {
            if nse + ii < state.c[nre - 1].len() {
                state.c[nre - 1][nse + ii] += params.creip[ii][id_idx] * params.reint[id_idx];
            }
        }
    }

    // ==================== 微分方程部分 ====================
    bre_differential(params, state, id, nre, nhe, npc, nmp, nse);
}

/// 计算微分方程部分。
fn bre_differential(
    params: &BreParams,
    state: &mut BreState,
    id: usize,
    nre: usize,
    nhe: usize,
    npc: usize,
    nmp: usize,
    nse: usize,
) {
    let id_idx = id - 1;

    if params.redif[id_idx] == 0.0 {
        return;
    }

    // TEFF^4 项
    let mut teffd = params.teff.powi(4);
    if params.idisk == 1 {
        teffd *= UN - params.thetav[id_idx];
    }
    state.vecl[nre - 1] += SIG4P * teffd * params.redif[id_idx];

    if id == 1 {
        // 上边界条件
        bre_upper_boundary(params, state, nre, nhe, npc, nse);
        return;
    }

    // ==================== 内部深度点 ====================

    let ddm = (params.dm[id_idx] - params.dm[id_idx - 1]) * HALF;
    let mut aren = 0.0;
    let mut bren = 0.0;
    let mut arepc = 0.0;
    let mut brepc = 0.0;

    // GP, GN 几何因子
    let (gp, gn) = if params.inmp > 0 { (UN, 0.0) } else { (0.0, UN) };

    // 初始化辅助数组
    let mut rexa = vec![0.0; params.nlvexp];
    let mut rexb = vec![0.0; params.nlvexp];

    if params.nfreqe > 0 {
        for ij in 1..=params.nfreqe {
            let ij_idx = ij - 1;

            let omeg0 = params.abso0[ij_idx] * params.dens1[id_idx];
            let omegm = params.absom[ij_idx] * params.dens1[id_idx - 1];
            let dtaum = (omeg0 + omegm) * ddm;

            if dtaum.abs() < 1e-30 {
                continue;
            }

            let frd = params.fk0[ij_idx] * params.rad0[ij_idx]
                - params.fkm[ij_idx] * params.radm[ij_idx];
            let gamr = frd / dtaum;
            let a1 = gamr / (omeg0 + omegm);
            let a3r = a1 * params.dens1[id_idx - 1] * params.wdep0[ij_idx];
            let b3r = a1 * params.dens1[id_idx] * params.wdep0[ij_idx];

            // A 矩阵项
            state.a[nre - 1][ij - 1] =
                -params.wdep0[ij_idx] * params.fkm[ij_idx] / dtaum * params.redif[id_idx];
            let rtr = omegm * params.wmm[id_idx - 1] * a3r;
            aren += rtr * gn;
            arepc -= a3r * params.dabnm[ij_idx] + rtr * gn;

            if params.inmp != 0 {
                state.a[nre - 1][nmp - 1] += rtr * gp * params.redif[id_idx];
            }
            state.a[nre - 1][nre - 1] -= a3r * params.dabtm[ij_idx] * params.redif[id_idx];

            // B 矩阵项
            state.b[nre - 1][ij - 1] +=
                params.wdep0[ij_idx] * params.fk0[ij_idx] / dtaum * params.redif[id_idx];
            let rtr = omeg0 * params.wmm[id_idx] * b3r;
            bren += rtr * gn;
            brepc -= b3r * params.dabn0[ij_idx] - rtr * gn;

            if params.inmp != 0 {
                state.b[nre - 1][nmp - 1] +=
                    (rtr + params.redx[id_idx]) * gp * params.redif[id_idx];
            }

            // 温度列
            state.b[nre - 1][nre - 1] -= b3r * params.dabt0[ij_idx] * params.redif[id_idx];

            // 能级导数
            for i in 0..params.nlvexp {
                rexa[i] -= a3r * params.drchm[i][ij_idx];
                rexb[i] -= b3r * params.drch0[i][ij_idx];
            }

            // RHS
            state.vecl[nre - 1] -= params.wdep0[ij_idx] * gamr * params.redif[id_idx];
        }
    }

    // N 列（氦）
    if params.inhe != 0 {
        state.a[nre - 1][nhe - 1] = (aren + params.redxm[id_idx]) * params.redif[id_idx];
        state.b[nre - 1][nhe - 1] += (bren + params.redx[id_idx]) * params.redif[id_idx];
    }

    // 温度列
    state.a[nre - 1][nre - 1] += params.redtm[id_idx] * params.redif[id_idx];
    state.b[nre - 1][nre - 1] += params.redt[id_idx] * params.redif[id_idx];
    state.c[nre - 1][nre - 1] += params.redtp[id_idx] * params.redif[id_idx];

    // 电子密度列
    if params.inpc != 0 {
        state.a[nre - 1][npc - 1] +=
            (arepc + params.rednm[id_idx] - params.redxm[id_idx]) * params.redif[id_idx];
        state.b[nre - 1][npc - 1] +=
            (brepc + params.redn[id_idx] - params.redx[id_idx]) * params.redif[id_idx];
        state.c[nre - 1][npc - 1] += params.rednp[id_idx] * params.redif[id_idx];
    }

    // 能级列
    for ii in 0..params.nlvexp {
        if nse + ii < state.a[nre - 1].len() {
            state.a[nre - 1][nse + ii] +=
                (rexa[ii] + params.redpm[ii][id_idx]) * params.redif[id_idx];
        }
        if nse + ii < state.b[nre - 1].len() {
            state.b[nre - 1][nse + ii] +=
                (rexb[ii] + params.redp[ii][id_idx]) * params.redif[id_idx];
        }
        if nse + ii < state.c[nre - 1].len() {
            state.c[nre - 1][nse + ii] += params.redpp[ii][id_idx] * params.redif[id_idx];
        }
    }
}

/// 上边界条件（ID = 1）。
fn bre_upper_boundary(
    params: &BreParams,
    state: &mut BreState,
    nre: usize,
    nhe: usize,
    npc: usize,
    nse: usize,
) {
    let id_idx = 0; // ID = 1

    if params.nfreqe > 0 {
        for ij in 1..=params.nfreqe {
            let ij_idx = ij - 1;
            let ijt = params.ijfr[ij_idx];

            let fh_val = if ijt > 0 && ijt <= params.fh.len() {
                params.fh[ijt - 1]
            } else {
                0.0
            };
            let hextrd_val = if ijt > 0 && ijt <= params.hextrd.len() {
                params.hextrd[ijt - 1]
            } else {
                0.0
            };

            let wf = params.wdep0[ij_idx] * fh_val * params.redif[id_idx];
            state.b[nre - 1][ij - 1] += wf;
            state.vecl[nre - 1] -=
                wf * params.rad0[ij_idx] + params.wdep0[ij_idx] * hextrd_val * params.redif[id_idx];
        }
    }

    // 温度列
    state.b[nre - 1][nre - 1] += params.redt[id_idx] * params.redif[id_idx];
    state.c[nre - 1][nre - 1] += params.redtp[id_idx] * params.redif[id_idx];

    // N 和电子密度列
    if params.inhe != 0 {
        state.b[nre - 1][nhe - 1] += params.redx[id_idx] * params.redif[id_idx];
    }
    if params.inpc != 0 {
        state.b[nre - 1][npc - 1] += params.redn[id_idx] * params.redif[id_idx];
    }
    if params.inhe != 0 {
        state.c[nre - 1][nhe - 1] += params.redxp[id_idx] * params.redif[id_idx];
    }
    if params.inpc != 0 {
        state.c[nre - 1][npc - 1] += params.rednp[id_idx] * params.redif[id_idx];
    }

    // 能级列
    for ii in 0..params.nlvexp {
        if nse + ii < state.b[nre - 1].len() {
            state.b[nre - 1][nse + ii] += params.redp[ii][id_idx] * params.redif[id_idx];
        }
        if nse + ii < state.c[nre - 1].len() {
            state.c[nre - 1][nse + ii] += params.redpp[ii][id_idx] * params.redif[id_idx];
        }
    }
}

// ============================================================================
// 辅助方法
// ============================================================================

impl<'a> BreParams<'a> {
    /// 计算 INSE（谱线起始索引）
    fn inse(&self) -> usize {
        1
    }
}

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

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

    #[test]
    fn test_bre_compile() {
        // 编译测试 - 验证类型签名正确
        assert!(true);
    }

    #[test]
    fn test_compt0_bre_iji_1() {
        // 当 IJI = 1 时，所有输出应为 0
        let kij = vec![100]; // NFREQ - 100 + 1 = 1
        let elec = vec![1e12];
        let result = compt0_bre(1, 1, 1e-8, 100, &kij, &elec, 6.6516e-25, 0, 0);
        assert!((result.0).abs() < 1e-15);
        assert!((result.1).abs() < 1e-15);
        assert!((result.2).abs() < 1e-15);
    }
}