//! 辐射转移方程矩阵（几何深度版本）。
//!
//! 重构自 TLUSTY `brtez.f`。
//!
//! 使用 ZD（几何深度）而不是 DM（柱质量密度）来计算矩阵元素。

const UN: f64 = 1.0;
const HALF: f64 = 0.5;
const THIRD: f64 = 1.0 / 3.0;
const SIXTH: f64 = 1.0 / 6.0;
const XCON: f64 = 8.0935e-21;
const YCON: f64 = 1.68638e-10;
const SIGE: f64 = 6.6516e-25;  // Thomson 散射截面

/// 简化版 Compton 计算（BRTEZ 使用）。
///
/// 返回 (cma, cmb, cmc, cme, cms, cmd)
fn compt0_brtez(
    ij: usize,
    id: usize,
    ab: f64,
    nfreq: usize,
    kij: &[usize],
    elec: &[f64],
) -> (f64, f64, f64, f64, f64, f64) {
    let iji = nfreq - kij[ij - 1] + 1;
    if iji == 1 {
        return (0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
    }
    let ss0 = elec[id - 1] * SIGE / ab;
    (0.0, 0.0, 0.0, 0.0, ss0, 0.0)
}

/// BRTEZ 参数结构体
#[allow(non_snake_case)]
pub struct BrtezParams<'a> {
    // 标量参数
    pub id: usize,
    pub nfreqe: usize,
    pub nd: usize,
    pub ndre: usize,
    pub inhe: usize,
    pub inre: usize,
    pub inpc: usize,
    pub inse: usize,
    pub inmp: usize,
    pub nlvexp: usize,
    pub nn0: usize,
    pub isplin: i32,
    pub icompt: i32,
    pub icombc: i32,
    pub ichcoo: i32,
    pub idisk: i32,
    pub ifz0: i32,
    pub ibc: i32,
    pub iwinbl: i32,
    pub radzer: f64,
    pub nfreq: usize,
    pub dst: f64,
    pub dsn: f64,
    pub bn: f64,
    pub hk: f64,

    // 数组参数
    pub zd: &'a [f64],
    pub temp: &'a [f64],
    pub dens: &'a [f64],
    pub elec: &'a [f64],
    pub freq: &'a [f64],
    pub dlnfr: &'a [f64],
    pub delj: &'a [f64],
    pub kij: &'a [usize],
    pub ijex: &'a [usize],
    pub ijorig: &'a [usize],
    pub ijfr: &'a [usize],
    pub sigec: &'a [f64],
    pub fhd: &'a [f64],
    pub hextrd: &'a [f64],

    // 辐射场和相关量
    pub rad: &'a [f64],       // nfreq x nd
    pub radex: &'a [f64],     // nfreqe x nd
    pub fk0: &'a [f64],
    pub fkp: &'a [f64],
    pub fkm: &'a [f64],
    pub abso0: &'a [f64],
    pub absop: &'a [f64],
    pub absom: &'a [f64],
    pub scat0: &'a [f64],
    pub scatp: &'a [f64],
    pub scatm: &'a [f64],
    pub emis0: &'a [f64],
    pub emisp: &'a [f64],
    pub emism: &'a [f64],
    pub q0: &'a [f64],
    pub uu0: &'a [f64],
    pub wmm: &'a [f64],

    // 矩阵导数项
    pub dabt0: &'a [f64],
    pub dabtp: &'a [f64],
    pub dabtm: &'a [f64],
    pub dabn0: &'a [f64],
    pub dabnp: &'a [f64],
    pub dabnm: &'a [f64],
    pub dabm0: &'a [f64],
    pub dabmp: &'a [f64],
    pub dabmm: &'a [f64],
    pub demt0: &'a [f64],
    pub demtp: &'a [f64],
    pub demtm: &'a [f64],
    pub demn0: &'a [f64],
    pub demnp: &'a [f64],
    pub demnm: &'a [f64],
    pub demm0: &'a [f64],
    pub demmp: &'a [f64],
    pub demmm: &'a [f64],
    pub drch0: &'a [f64],  // nlvexp x nfreqe
    pub drchp: &'a [f64],
    pub drchm: &'a [f64],
    pub dret0: &'a [f64],
    pub dretp: &'a [f64],
    pub dretm: &'a [f64],

    // 辐射场指针
    pub rad0: &'a [f64],
    pub radp: &'a [f64],
    pub radm: &'a [f64],
}

/// BRTEZ 状态结构体（可变矩阵）
#[allow(non_snake_case)]
pub struct BrtezState<'a> {
    pub a: &'a mut [&'a mut [f64]],
    pub b: &'a mut [&'a mut [f64]],
    pub c: &'a mut [&'a mut [f64]],
    pub vecl: &'a mut [f64],
}

/// 辐射转移方程矩阵（几何深度版本）。
///
/// # 参数
/// * `params` - 输入参数
/// * `state` - 可变矩阵状态
#[allow(non_snake_case)]
pub fn brtez(params: &BrtezParams, state: &mut BrtezState) {
    if params.nfreqe == 0 {
        return;
    }

    let id = params.id;
    let id_idx = id - 1;  // 0-indexed

    let mut ispl = params.isplin;
    if ispl >= 5 {
        ispl -= 5;
    }

    let nhe = params.nfreqe + params.inhe;
    let nre = params.nfreqe + params.inre;
    let npc = params.nfreqe + params.inpc;
    let nse = params.nfreqe + params.inse - 1;
    let nmp = params.nfreqe + params.inmp;

    let mut gp = 0.0;
    let mut gn = UN;
    if params.inmp > 0 {
        gp = UN;
        gn = 0.0;
    }

    // Compton 散射边界条件
    let mut ij1 = 1;
    if params.icompt > 0 && params.icombc > 0 && params.ijex[0] > 0 {
        ij1 = 2;
        let ij = 1;
        let iji = params.nfreq;
        let zj1 = (-params.hk * params.freq[ij - 1] / params.temp[id_idx]).exp();
        let zj2 = (-params.hk * params.freq[ij] / params.temp[id_idx]).exp();
        let dlt = params.delj[iji - 2 + id_idx * (params.nfreq - 1)];

        let (combid, comaid) = if params.ichcoo == 0 {
            let zj0 = UN / (params.hk * (params.freq[ij - 1] * params.freq[ij]).sqrt() / params.temp[id_idx]);
            let zxx = UN - 3.0 * zj0 + (UN - dlt) * zj1 + dlt * zj2;
            let combid = zj0 / params.dlnfr[iji - 2] + (UN - dlt) * zxx;
            let comaid = -zj0 / params.dlnfr[iji - 2] + dlt * zxx;
            (combid, comaid)
        } else {
            let e2 = YCON * params.temp[id_idx];
            let zxx0 = XCON * params.freq[ij - 1] * (UN + zj1) - 3.0 * e2;
            let zxxm = XCON * params.freq[ij] * (UN + zj2) - 3.0 * e2;
            let zxx = (UN - dlt) * zxx0 + dlt * zxxm;
            let combid = e2 / params.dlnfr[iji - 2] + (UN - dlt) * zxx;
            let comaid = -e2 / params.dlnfr[iji - 2] + dlt * zxx;
            (combid, comaid)
        };

        state.b[ij - 1][ij - 1] = combid;
        state.b[ij - 1][ij] = comaid;
        // RAD 数组索引: rad(iji, id) -> rad[(iji-1) + id_idx * nfreq]
        let rad_idx = (iji - 1) + id_idx * params.nfreq;
        let rad_idx_prev = (iji - 2) + id_idx * params.nfreq;
        state.vecl[ij - 1] = -state.b[ij - 1][ij - 1] * params.rad[rad_idx]
                           - state.b[ij - 1][ij] * params.rad[rad_idx_prev];
    }

    // ID = 1: 上边界条件
    if id > 1 {
        brtez_internal(params, state, id, ij1, nhe, nre, npc, nse, nmp, gn, gp, ispl);
        return;
    }

    brtez_upper_boundary(params, state, id_idx, ij1, nhe, nre, npc, nse, nmp, gn, gp, ispl);
}

/// 上边界条件 (ID = 1)
#[allow(non_snake_case)]
fn brtez_upper_boundary(
    params: &BrtezParams,
    state: &mut BrtezState,
    id_idx: usize,
    ij1: usize,
    nhe: usize,
    nre: usize,
    npc: usize,
    nse: usize,
    nmp: usize,
    gn: f64,
    gp: f64,
    ispl: i32,
) {
    let ddp = (params.zd[0] - params.zd[1]) * HALF;

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

        let omeg0 = params.abso0[ij_idx];
        let omegp = params.absop[ij_idx];
        let dzp = omeg0 + omegp;
        let dtaup = dzp * ddp;

        let alf1 = (params.fk0[ij_idx] * params.rad0[ij_idx]
                  - params.fkp[ij_idx] * params.radp[ij_idx]) / dtaup;

        let chiel0 = params.scat0[ij_idx];
        let chielp = params.scatp[ij_idx];
        let s0 = (params.emis0[ij_idx] + chiel0 * params.rad0[ij_idx]) / params.abso0[ij_idx];

        let mut bs = HALF * dtaup;
        let mut cs = 0.0;
        let mut c2 = 0.0;
        let mut gam2 = 0.0;
        let mut sp = 0.0;

        // Compton 散射
        let (cma, cmb, cmc, cme, cms, cmd) = if params.icompt > 0 {
            compt0_brtez(ijt + 1, id_idx + 1, params.abso0[ij_idx], params.nfreq, params.kij, params.elec)
        } else {
            (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
        };

        let s0 = if params.icompt > 0 { s0 + cms } else { s0 };

        if ispl % 3 > 0 {
            // Spline 配点法或 Hermite 方法
            bs = dtaup * THIRD;
            cs = HALF * bs;
            sp = (params.emisp[ij_idx] + chielp * params.radp[ij_idx]) / params.absop[ij_idx];
            c2 = cs / params.absop[ij_idx];
            gam2 = cs * (params.radp[ij_idx] - sp);
        }

        // 辅助量
        let alf2 = bs * (params.rad0[ij_idx] - s0);
        let bet2 = alf2 + gam2;
        let x1 = (alf1 - bet2) / dzp;
        let b2 = (bs + params.q0[ij_idx]) / params.abso0[ij_idx];
        let mut b1 = x1;
        b1 = b1 + params.uu0[ij_idx] * s0 * params.zd[0] / params.dens[0];  // 使用 ZD(1) 替代 DM(1)
        let mut c1 = x1;
        b1 = b1 - b2 * s0;
        c1 = c1 - c2 * sp;

        // 矩阵元素
        let rtn = omeg0 * params.wmm[id_idx] * b1;
        let rtnc = omegp * params.wmm[id_idx + 1] * c1;

        state.b[ij_idx][nre - 1] = b1 * params.dabt0[ij_idx]
                                 + b2 * (params.demt0[ij_idx] + params.dst * params.rad0[ij_idx]);
        state.c[ij_idx][nre - 1] = c1 * params.dabtp[ij_idx]
                                 + c2 * (params.demtp[ij_idx] + params.dst * params.radp[ij_idx]);

        let sigec_val = if ijt < params.sigec.len() { params.sigec[ijt] } else { 0.0 };

        state.b[ij_idx][npc - 1] = b1 * params.dabn0[ij_idx]
                                 + b2 * (params.demn0[ij_idx] + (params.dsn + sigec_val) * params.rad0[ij_idx]);
        state.c[ij_idx][npc - 1] = c1 * params.dabnp[ij_idx]
                                 + c2 * (params.demnp[ij_idx] + (params.dsn + sigec_val) * params.radp[ij_idx]);

        state.b[ij_idx][nmp - 1] = b1 * params.dabm0[ij_idx] + b2 * params.demm0[ij_idx] - gp * rtn;
        state.c[ij_idx][nmp - 1] = c1 * params.dabmp[ij_idx] + c2 * params.demmp[ij_idx] - gp * rtnc;

        // 能级列
        for ii in 0..params.nlvexp {
            let drch0_val = params.drch0[ii * params.nfreqe + ij_idx];
            let drchp_val = params.drchp[ii * params.nfreqe + ij_idx];
            let dret0_val = params.dret0[ii * params.nfreqe + ij_idx];
            let dretp_val = params.dretp[ii * params.nfreqe + ij_idx];
            state.b[ij_idx][nse + ii] = state.b[ij_idx][nse + ii] + b1 * drch0_val + b2 * dret0_val;
            state.c[ij_idx][nse + ii] = state.c[ij_idx][nse + ii] + c1 * drchp_val + c2 * dretp_val;
        }

        state.b[ij_idx][params.nfreqe - 1] = 0.0;
        state.b[ij_idx][ij_idx] = -params.fk0[ij_idx] / dtaup - params.fhd[ijt]
                                 - bs * (UN - chiel0 / params.abso0[ij_idx])
                                 + params.q0[ij_idx] * chiel0 / params.abso0[ij_idx];
        state.c[ij_idx][params.nfreqe - 1] = 0.0;
        state.c[ij_idx][ij_idx] = params.fkp[ij_idx] / dtaup - cs * (UN - chielp / params.absop[ij_idx]);

        // RHS 向量
        state.vecl[ij_idx] = alf1 + bet2 + params.fhd[ijt] * params.rad0[ij_idx] - s0 * params.q0[ij_idx];
        if params.iwinbl < 0 {
            state.vecl[ij_idx] -= params.hextrd[ijt];
        }

        // Compton 散射附加项
        if params.icompt > 4 {
            let iji = params.nfreq - params.kij[ijt] + 1;
            state.b[ij_idx][ij_idx] += bs * (cmb + cme);
            if iji > 1 {
                let ijm = params.ijex[params.ijorig[iji - 2] - 1];
                if ijm > 0 {
                    state.b[ij_idx][ijm - 1] += bs * cma;
                }
            }
            if iji < params.nfreq {
                let ijp_idx = params.ijex[params.ijorig[iji] - 1];
                if ijp_idx > 0 {
                    state.b[ij_idx][ijp_idx - 1] += bs * cmc;
                }
            }
            if params.inre > 0 {
                state.b[ij_idx][nre - 1] += cmd * bs;
            }
            if params.inpc > 0 {
                state.b[ij_idx][npc - 1] += cms * bs / params.elec[id_idx];
            }
        }
    }
}

/// 内部深度点 (1 < ID < ND)
#[allow(non_snake_case)]
fn brtez_internal(
    params: &BrtezParams,
    state: &mut BrtezState,
    id: usize,
    ij1: usize,
    nhe: usize,
    nre: usize,
    npc: usize,
    nse: usize,
    nmp: usize,
    gn: f64,
    gp: f64,
    ispl: i32,
) {
    let id_idx = id - 1;
    let ddm = (params.zd[id_idx - 1] - params.zd[id_idx]) * HALF;

    if id == params.nd {
        brtez_lower_boundary(params, state, id_idx, ij1, nhe, nre, npc, nse, nmp, gn, gp, ddm, ispl);
        return;
    }

    let ddp = (params.zd[id_idx] - params.zd[id_idx + 1]) * HALF;

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

        let omeg0 = params.abso0[ij_idx];
        let omegp = params.absop[ij_idx];
        let omegm = params.absom[ij_idx];
        let dzp = omeg0 + omegp;
        let dzm = omeg0 + omegm;
        let dtaup = dzp * ddp;
        let dtaum = dzm * ddm;
        let dtau0 = HALF * (dtaup + dtaum);

        let frd = params.fk0[ij_idx] * params.rad0[ij_idx];
        let alf1 = (frd - params.fkp[ij_idx] * params.radp[ij_idx]) / dtaup / dtau0;
        let gam1 = (frd - params.fkm[ij_idx] * params.radm[ij_idx]) / dtaum / dtau0;
        let bet1 = alf1 + gam1;
        let x1 = HALF * bet1 / dtau0;

        let mut a1 = (gam1 + x1 * dtaum) / dzm;
        let mut c1 = (alf1 + x1 * dtaup) / dzp;
        let mut b1 = a1 + c1;

        let mut bs = UN;
        let chielm = params.scatm[ij_idx];
        let chiel0 = params.scat0[ij_idx];
        let chielp = params.scatp[ij_idx];
        let s0 = (params.emis0[ij_idx] + chiel0 * params.rad0[ij_idx]) / params.abso0[ij_idx];

        let mut as_val = 0.0;
        let mut cs = 0.0;
        let mut a2 = 0.0;
        let mut c2 = 0.0;
        let mut a3 = 0.0;
        let mut c3 = 0.0;
        let mut bet2 = 0.0;
        let mut sm = 0.0;
        let mut sp = 0.0;

        // Compton 散射
        let (cma, cmb, cmc, cme, cms, cmd) = if params.icompt > 0 {
            compt0_brtez(ijt + 1, id, params.abso0[ij_idx], params.nfreq, params.kij, params.elec)
        } else {
            (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
        };

        let s0 = if params.icompt > 0 { s0 + cms } else { s0 };

        if ispl % 3 != 0 {
            sm = (params.emism[ij_idx] + params.radm[ij_idx] * chielm) / params.absom[ij_idx];
            sp = (params.emisp[ij_idx] + params.radp[ij_idx] * chielp) / params.absop[ij_idx];

            if ispl == 1 {
                // Spline 配点法
                as_val = dtaum / dtau0 * SIXTH;
                cs = dtaup / dtau0 * SIXTH;
                bs = 0.666666666666667_f64;
                let alf2 = as_val * (params.radm[ij_idx] - sm);
                let gam2 = cs * (params.radp[ij_idx] - sp);
                bet2 = alf2 + gam2;
                let x = HALF * bet2 / dtau0;
                a2 = (gam2 - x * dtaum) / dzm;
                c2 = (alf2 - x * dtaup) / dzp;
            } else {
                // Hermite 方法
                as_val = dtaup * dtaup / dtaum / dtau0;
                cs = dtaum * dtaum / dtaup / dtau0;
                let al3 = (params.radp[ij_idx] - sp - params.rad0[ij_idx] + s0) * SIXTH;
                let ga3 = (params.radm[ij_idx] - sm - params.rad0[ij_idx] + s0) * SIXTH;
                let av = al3 * cs;
                let cv = ga3 * as_val;
                as_val = (UN - HALF * as_val) * SIXTH;
                cs = (UN - HALF * cs) * SIXTH;
                bs = UN - as_val - cs;
                let x = (av + cv) / dtau0 / 4.0;
                a2 = (x * dtaum + HALF * cv - av) / dzm;
                c2 = (x * dtaup + HALF * av - cv) / dzp;
                bet2 = as_val * (params.radm[ij_idx] - sm) + cs * (params.radp[ij_idx] - sp);
            }
        }

        // 辅助量
        b1 = b1 - (a2 + c2);
        a1 = a1 - a2;
        c1 = c1 - c2;
        a2 = as_val / params.absom[ij_idx];
        c2 = cs / params.absop[ij_idx];
        a3 = a2 * sm;
        c3 = c2 * sp;

        let b2 = bs / params.abso0[ij_idx];
        let b3 = b2 * s0;

        a1 = a1 - a3;
        b1 = b1 - b3;
        c1 = c1 - c3;

        // 矩阵元素
        let rtna = omegm * params.wmm[id_idx - 1] * a1;
        let rtn = omeg0 * params.wmm[id_idx] * b1;
        let rtnc = omegp * params.wmm[id_idx + 1] * c1;

        state.a[ij_idx][nhe - 1] = -gn * rtna;
        state.b[ij_idx][nhe - 1] = -gn * rtn;
        state.c[ij_idx][nhe - 1] = -gn * rtnc;

        let a1_adj = a1 - a3;
        let b1_adj = b1 - b3;
        let c1_adj = c1 - c3;

        state.a[ij_idx][nre - 1] = a1_adj * params.dabtm[ij_idx]
                                 + a2 * (params.demtm[ij_idx] + params.dst * params.radm[ij_idx]);
        state.b[ij_idx][nre - 1] = b1_adj * params.dabt0[ij_idx]
                                 + b2 * (params.demt0[ij_idx] + params.dst * params.rad0[ij_idx]);
        state.c[ij_idx][nre - 1] = c1_adj * params.dabtp[ij_idx]
                                 + c2 * (params.demtp[ij_idx] + params.dst * params.radp[ij_idx]);

        let sigec_val = if ijt < params.sigec.len() { params.sigec[ijt] } else { 0.0 };

        state.a[ij_idx][npc - 1] = a1_adj * params.dabnm[ij_idx]
                                 + a2 * (params.demnm[ij_idx] + (params.dsn + sigec_val) * params.radm[ij_idx]);
        state.b[ij_idx][npc - 1] = b1_adj * params.dabn0[ij_idx]
                                 + b2 * (params.demn0[ij_idx] + (params.dsn + sigec_val) * params.rad0[ij_idx]);
        state.c[ij_idx][npc - 1] = c1_adj * params.dabnp[ij_idx]
                                 + c2 * (params.demnp[ij_idx] + (params.dsn + sigec_val) * params.radp[ij_idx]);

        state.a[ij_idx][nmp - 1] = a1_adj * params.dabmm[ij_idx] + a2 * params.demmm[ij_idx] - gp * rtna;
        state.b[ij_idx][nmp - 1] = b1_adj * params.dabm0[ij_idx] + b2 * params.demm0[ij_idx] - gp * rtn;
        state.c[ij_idx][nmp - 1] = c1_adj * params.dabmp[ij_idx] + c2 * params.demmp[ij_idx] - gp * rtnc;

        // 能级列
        for ii in 0..params.nlvexp {
            let drchm_val = params.drchm[ii * params.nfreqe + ij_idx];
            let drch0_val = params.drch0[ii * params.nfreqe + ij_idx];
            let drchp_val = params.drchp[ii * params.nfreqe + ij_idx];
            let dretm_val = params.dretm[ii * params.nfreqe + ij_idx];
            let dret0_val = params.dret0[ii * params.nfreqe + ij_idx];
            let dretp_val = params.dretp[ii * params.nfreqe + ij_idx];

            state.a[ij_idx][nse + ii] = state.a[ij_idx][nse + ii] + a1_adj * drchm_val + a2 * dretm_val;
            state.b[ij_idx][nse + ii] = state.b[ij_idx][nse + ii] + b1_adj * drch0_val + b2 * dret0_val;
            state.c[ij_idx][nse + ii] = state.c[ij_idx][nse + ii] + c1_adj * drchp_val + c2 * dretp_val;
        }

        state.a[ij_idx][params.nfreqe - 1] = 0.0;
        state.a[ij_idx][ij_idx] = params.fkm[ij_idx] / dtaum / dtau0 - as_val * (UN - chielm / params.absom[ij_idx]);
        state.b[ij_idx][params.nfreqe - 1] = 0.0;
        state.b[ij_idx][ij_idx] = -params.fk0[ij_idx] / dtau0 * (UN / dtaup + UN / dtaum)
                                 - bs * (UN - chiel0 / params.abso0[ij_idx]);
        state.c[ij_idx][params.nfreqe - 1] = 0.0;
        state.c[ij_idx][ij_idx] = params.fkp[ij_idx] / dtaup / dtau0 - cs * (UN - chielp / params.absop[ij_idx]);

        // RHS 向量
        state.vecl[ij_idx] = bet1 + bet2 + bs * (params.rad0[ij_idx] - s0);

        // Compton 散射附加项
        if params.icompt > 4 {
            let iji = params.nfreq - params.kij[ijt] + 1;
            state.b[ij_idx][ij_idx] += bs * (cmb + cme);
            if iji > 1 {
                let ijm = params.ijex[params.ijorig[iji - 2] - 1];
                if ijm > 0 {
                    state.b[ij_idx][ijm - 1] += bs * cma;
                }
            }
            if iji < params.nfreq {
                let ijp_idx = params.ijex[params.ijorig[iji] - 1];
                if ijp_idx > 0 {
                    state.b[ij_idx][ijp_idx - 1] += bs * cmc;
                }
            }
            if params.inre > 0 {
                state.b[ij_idx][nre - 1] += cmd * bs;
            }
            if params.inpc > 0 {
                state.b[ij_idx][npc - 1] += cms * bs / params.elec[id_idx];
            }
        }
    }
}

/// 下边界条件 (ID = ND)
#[allow(non_snake_case)]
fn brtez_lower_boundary(
    params: &BrtezParams,
    state: &mut BrtezState,
    id_idx: usize,
    ij1: usize,
    _nhe: usize,
    nre: usize,
    npc: usize,
    nse: usize,
    nmp: usize,
    gn: f64,
    gp: f64,
    ddm: f64,
    ispl: i32,
) {
    if params.idisk == 0 || params.ifz0 < 0 {
        // 非盘模型
        brtez_lower_nodisk(params, state, id_idx, ij1, nre, npc, nse, nmp, gn, gp, ddm, ispl);
    } else {
        // 盘模型
        brtez_lower_disk(params, state, id_idx, ij1, nre, npc, nse, nmp, gn, gp, ddm);
    }
}

/// 下边界条件（非盘模型）
#[allow(non_snake_case)]
fn brtez_lower_nodisk(
    params: &BrtezParams,
    state: &mut BrtezState,
    id_idx: usize,
    ij1: usize,
    nre: usize,
    npc: usize,
    nse: usize,
    nmp: usize,
    gn: f64,
    gp: f64,
    ddm: f64,
    _ispl: i32,
) {
    let t = params.temp[id_idx];
    let tm = params.temp[id_idx - 1];
    let hkt = params.hk / t;
    let hktm = params.hk / tm;

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

        let chielm = params.scatm[ij_idx];
        let chiel0 = params.scat0[ij_idx];
        let omegm = params.absom[ij_idx];
        let omeg0 = params.abso0[ij_idx];
        let dzm = omeg0 + omegm;
        let dtaum = dzm * ddm;

        let frd = params.fk0[ij_idx] * params.rad0[ij_idx] - params.fkm[ij_idx] * params.radm[ij_idx];
        let gam1 = frd / dtaum;
        let mut a1 = gam1 / dzm;

        let mut as_val = 0.0;
        let mut bs = 0.0;
        let mut a2 = 0.0;
        let mut b2 = 0.0;
        let mut a3 = 0.0;
        let mut b3 = 0.0;
        let mut alf2 = 0.0;
        let mut bet2 = 0.0;
        let mut gam2 = 0.0;

        // 二阶边界条件
        if params.ibc > 0 && params.ibc < 4 {
            bs = dtaum * HALF;
            let s0 = (params.emis0[ij_idx] + chiel0 * params.rad0[ij_idx]) / params.abso0[ij_idx];

            // Compton 散射
            let s0 = if params.icompt > 0 {
                let (_, _, _, _, cms, _) = compt0_brtez(ijt + 1, id_idx + 1, params.abso0[ij_idx], params.nfreq, params.kij, params.elec);
                s0 + cms
            } else {
                s0
            };

            gam2 = bs * (params.rad0[ij_idx] - s0);
            bet2 = gam2;
            let x1 = bet2 / dzm;
            a1 = a1 - x1;
            b2 = bs / params.abso0[ij_idx];
            b3 = b2 * s0;
        }

        // 辅助参数
        let fr = params.freq[ijt];
        let fr15 = fr * 1e-15;
        let x = hkt * fr;
        let ex = x.exp();
        let xm = hktm * fr;
        let exm = xm.exp();
        let plan = params.bn * fr15 * fr15 * fr15 / (ex - UN);

        let mut planm = params.bn * fr15 * fr15 * fr15 / (exm - UN);

        if params.inre == 0 || params.id >= params.ndre {
            planm = params.bn * fr15 * fr15 * fr15 / (exm - UN);
            let gam3 = (plan - planm) / dtaum * THIRD;
            a1 = a1 - gam3 / dzm;
            // gam1 -= gam3; // 在原代码中 gam1 未使用
        }

        let c1 = a1;
        let b1 = c1;
        let a1_adj = a1 - a3;
        let b1_adj = b1 - b3;

        // 矩阵元素
        if params.inre == 0 || params.id >= params.ndre {
            let dplanm = planm * xm / tm / (UN - UN / exm);
            state.a[ij_idx][nre - 1] = a1_adj * params.dabtm[ij_idx]
                                     + a2 * (params.demtm[ij_idx] + params.dst * params.radm[ij_idx])
                                     - dplanm / dtaum * THIRD;

            let sigec_val = if ijt < params.sigec.len() { params.sigec[ijt] } else { 0.0 };
            state.a[ij_idx][npc - 1] = a1_adj * params.dabnm[ij_idx]
                                     + a2 * (params.demnm[ij_idx] + (params.dsn + sigec_val) * params.radm[ij_idx]);

            let bb = HALF + THIRD / dtaum;
            let dplan = plan * x / t / (UN - UN / ex);
            state.b[ij_idx][nre - 1] = b1_adj * params.dabt0[ij_idx]
                                     + b2 * (params.demt0[ij_idx] + params.dst * params.rad0[ij_idx])
                                     + bb * dplan;
            state.b[ij_idx][npc - 1] = b1_adj * params.dabn0[ij_idx]
                                     + b2 * (params.demn0[ij_idx] + (params.dsn + sigec_val) * params.rad0[ij_idx]);

            let rtna = omegm * params.wmm[id_idx - 1] * a1_adj;
            let rtn = omeg0 * params.wmm[id_idx] * b1_adj;

            state.a[ij_idx][nmp - 1] = a1_adj * params.dabmm[ij_idx] + a2 * params.demmm[ij_idx] - gp * rtna;
            state.b[ij_idx][nmp - 1] = b1_adj * params.dabm0[ij_idx] + b2 * params.demm0[ij_idx] - gp * rtn;

            // 能级列
            for ii in 0..params.nlvexp {
                let drchm_val = params.drchm[ii * params.nfreqe + ij_idx];
                let drch0_val = params.drch0[ii * params.nfreqe + ij_idx];
                let dretm_val = params.dretm[ii * params.nfreqe + ij_idx];
                let dret0_val = params.dret0[ii * params.nfreqe + ij_idx];
                state.a[ij_idx][nse + ii] = state.a[ij_idx][nse + ii] + a1_adj * drchm_val + a2 * dretm_val;
                state.b[ij_idx][nse + ii] = state.b[ij_idx][nse + ii] + b1_adj * drch0_val + b2 * dret0_val;
            }

            state.a[ij_idx][params.nfreqe - 1] = 0.0;
            state.a[ij_idx][ij_idx] = params.fkm[ij_idx] / dtaum - as_val * (UN - chielm / params.absom[ij_idx]);
            state.b[ij_idx][params.nfreqe - 1] = 0.0;
        }

        // RHS 向量
        if params.ibc == 0 || params.ibc == 4 {
            state.b[ij_idx][ij_idx] = state.b[ij_idx][ij_idx] - params.fk0[ij_idx] / dtaum
                                    - bs * (UN - chiel0 / params.abso0[ij_idx]) - HALF;
            state.vecl[ij_idx] = gam1 + bet2 - HALF * (plan - params.rad0[ij_idx]);
        } else {
            state.b[ij_idx][ij_idx] = state.b[ij_idx][ij_idx] - params.fk0[ij_idx] / dtaum
                                    - bs * (UN - chiel0 / params.abso0[ij_idx]) - params.fhd[ijt];
            state.vecl[ij_idx] = gam1 + bet2 - HALF * plan + params.fhd[ijt] * params.rad0[ij_idx];
        }

        // Compton 散射附加项
        if params.icompt > 4 {
            let (cma, cmb, cmc, cme, cms, cmd) = compt0_brtez(ijt + 1, id_idx + 1, params.abso0[ij_idx], params.nfreq, params.kij, params.elec);
            let iji = params.nfreq - params.kij[ijt] + 1;
            state.b[ij_idx][ij_idx] += bs * (cmb + cme);
            if iji > 1 {
                let ijm = params.ijex[params.ijorig[iji - 2] - 1];
                if ijm > 0 {
                    state.b[ij_idx][ijm - 1] += bs * cma;
                }
            }
            if iji < params.nfreq {
                let ijp_idx = params.ijex[params.ijorig[iji] - 1];
                if ijp_idx > 0 {
                    state.b[ij_idx][ijp_idx - 1] += bs * cmc;
                }
            }
            if params.inre > 0 {
                state.b[ij_idx][nre - 1] += cmd * bs;
            }
            if params.inpc > 0 {
                state.b[ij_idx][npc - 1] += cms * bs / params.elec[id_idx];
            }
        }
    }
}

/// 下边界条件（盘模型）
#[allow(non_snake_case)]
fn brtez_lower_disk(
    params: &BrtezParams,
    state: &mut BrtezState,
    id_idx: usize,
    ij1: usize,
    nre: usize,
    npc: usize,
    nse: usize,
    nmp: usize,
    gn: f64,
    gp: f64,
    ddm: f64,
) {
    for ij in ij1..=params.nfreqe {
        let ij_idx = ij - 1;
        let ijt = params.ijfr[ij_idx] - 1;

        let chielm = params.scatm[ij_idx];
        let chiel0 = params.scat0[ij_idx];
        let omegm = params.absom[ij_idx];
        let omeg0 = params.abso0[ij_idx];
        let dzm = omeg0 + omegm;
        let dtaum = dzm * ddm;

        let frd = params.fk0[ij_idx] * params.rad0[ij_idx] - params.fkm[ij_idx] * params.radm[ij_idx];
        let gam1 = frd / dtaum;
        let mut a1 = gam1 / dzm;

        let mut as_val = 0.0;
        let mut a2 = 0.0;
        let mut a3 = 0.0;
        let mut alf2 = 0.0;
        let mut gam2 = 0.0;

        let bs = dtaum * HALF;
        let s0 = (params.emis0[ij_idx] + chiel0 * params.rad0[ij_idx]) / params.abso0[ij_idx];

        // Compton 散射
        let s0 = if params.icompt > 0 {
            let (_, _, _, _, cms, _) = compt0_brtez(ijt + 1, id_idx + 1, params.abso0[ij_idx], params.nfreq, params.kij, params.elec);
            s0 + cms
        } else {
            s0
        };

        gam2 = bs * (params.rad0[ij_idx] - s0);
        let bet2 = alf2 + gam2;
        let x1 = bet2 / dzm;
        a1 = a1 - x1;
        let b2 = bs / params.abso0[ij_idx];
        let b3 = b2 * s0;
        let b1 = a1;
        let a1_adj = a1 - a3;
        let b1_adj = b1 - b3;

        // 矩阵 A 元素
        state.a[ij_idx][nre - 1] = a1_adj * params.dabtm[ij_idx]
                                 + a2 * (params.demtm[ij_idx] + params.dst * params.radm[ij_idx]);

        let sigec_val = if ijt < params.sigec.len() { params.sigec[ijt] } else { 0.0 };
        state.a[ij_idx][npc - 1] = a1_adj * params.dabnm[ij_idx]
                                 + a2 * (params.demnm[ij_idx] + (params.dsn + sigec_val) * params.radm[ij_idx]);

        let rtna = omegm * params.wmm[id_idx - 1] * a1_adj;
        state.a[ij_idx][nmp - 1] = a1_adj * params.dabmm[ij_idx] + a2 * params.demmm[ij_idx] - gp * rtna;

        // 能级列
        for ii in 0..params.nlvexp {
            let drchm_val = params.drchm[ii * params.nfreqe + ij_idx];
            let dretm_val = params.dretm[ii * params.nfreqe + ij_idx];
            state.a[ij_idx][nse + ii] = a1_adj * drchm_val + a2 * dretm_val;
        }

        state.a[ij_idx][params.nfreqe - 1] = 0.0;
        state.a[ij_idx][ij_idx] = params.fkm[ij_idx] / dtaum - as_val * (UN - chielm / params.absom[ij_idx]);

        // 矩阵 B 元素
        state.b[ij_idx][nre - 1] = b1_adj * params.dabt0[ij_idx]
                                 + b2 * (params.demt0[ij_idx] + params.dst * params.rad0[ij_idx]);
        state.b[ij_idx][npc - 1] = b1_adj * params.dabn0[ij_idx]
                                 + b2 * (params.demn0[ij_idx] + (params.dsn + sigec_val) * params.rad0[ij_idx]);

        let rtn = omeg0 * params.wmm[id_idx] * b1_adj;
        state.b[ij_idx][nmp - 1] = b1_adj * params.dabm0[ij_idx] + b2 * params.demm0[ij_idx] - gp * rtn;

        // 能级列
        for ii in 0..params.nlvexp {
            let drch0_val = params.drch0[ii * params.nfreqe + ij_idx];
            let dret0_val = params.dret0[ii * params.nfreqe + ij_idx];
            state.b[ij_idx][nse + ii] = b1_adj * drch0_val + b2 * dret0_val;
        }

        state.b[ij_idx][params.nfreqe - 1] = 0.0;
        state.b[ij_idx][ij_idx] = -params.fk0[ij_idx] / dtaum - bs * (UN - chiel0 / params.abso0[ij_idx]);

        // RHS 向量
        state.vecl[ij_idx] = gam1 + bet2;

        // Compton 散射附加项
        if params.icompt > 4 {
            let (cma, cmb, cmc, cme, cms, cmd) = compt0_brtez(ijt + 1, id_idx + 1, params.abso0[ij_idx], params.nfreq, params.kij, params.elec);
            let iji = params.nfreq - params.kij[ijt] + 1;
            state.b[ij_idx][ij_idx] += bs * (cmb + cme);
            if iji > 1 {
                let ijm = params.ijex[params.ijorig[iji - 2] - 1];
                if ijm > 0 {
                    state.b[ij_idx][ijm - 1] += bs * cma;
                }
            }
            if iji < params.nfreq {
                let ijp_idx = params.ijex[params.ijorig[iji] - 1];
                if ijp_idx > 0 {
                    state.b[ij_idx][ijp_idx - 1] += bs * cmc;
                }
            }
            if params.inre > 0 {
                state.b[ij_idx][nre - 1] += cmd * bs;
            }
            if params.inpc > 0 {
                state.b[ij_idx][npc - 1] += cms * bs / params.elec[id_idx];
            }
        }
    }
}

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

    #[test]
    fn test_brtez_constants() {
        assert_relative_eq!(UN, 1.0, epsilon = 1e-15);
        assert_relative_eq!(HALF, 0.5, epsilon = 1e-15);
        assert_relative_eq!(THIRD, 1.0 / 3.0, epsilon = 1e-15);
        assert_relative_eq!(SIXTH, 1.0 / 6.0, epsilon = 1e-15);
        assert_relative_eq!(XCON, 8.0935e-21, epsilon = 1e-25);
        assert_relative_eq!(YCON, 1.68638e-10, epsilon = 1e-15);
    }

    #[test]
    fn test_brtez_early_return() {
        // 当 nfreqe = 0 时应直接返回
        let params = BrtezParams {
            id: 1,
            nfreqe: 0,
            nd: 10,
            ndre: 5,
            inhe: 1,
            inre: 1,
            inpc: 1,
            inse: 1,
            inmp: 0,
            nlvexp: 0,
            nn0: 10,
            isplin: 0,
            icompt: 0,
            icombc: 0,
            ichcoo: 0,
            idisk: 0,
            ifz0: 0,
            ibc: 0,
            iwinbl: 0,
            radzer: 0.0,
            nfreq: 100,
            dst: 0.0,
            dsn: 0.0,
            bn: 0.0,
            hk: 0.0,
            zd: &[],
            temp: &[],
            dens: &[],
            elec: &[],
            freq: &[],
            dlnfr: &[],
            delj: &[],
            kij: &[],
            ijex: &[],
            ijorig: &[],
            ijfr: &[],
            sigec: &[],
            fhd: &[],
            hextrd: &[],
            rad: &[],
            radex: &[],
            fk0: &[],
            fkp: &[],
            fkm: &[],
            abso0: &[],
            absop: &[],
            absom: &[],
            scat0: &[],
            scatp: &[],
            scatm: &[],
            emis0: &[],
            emisp: &[],
            emism: &[],
            q0: &[],
            uu0: &[],
            wmm: &[],
            dabt0: &[],
            dabtp: &[],
            dabtm: &[],
            dabn0: &[],
            dabnp: &[],
            dabnm: &[],
            dabm0: &[],
            dabmp: &[],
            dabmm: &[],
            demt0: &[],
            demtp: &[],
            demtm: &[],
            demn0: &[],
            demnp: &[],
            demnm: &[],
            demm0: &[],
            demmp: &[],
            demmm: &[],
            drch0: &[],
            drchp: &[],
            drchm: &[],
            dret0: &[],
            dretp: &[],
            dretm: &[],
            rad0: &[],
            radp: &[],
            radm: &[],
        };

        // 空矩阵
        let mut a_data: Vec<Vec<f64>> = vec![];
        let mut b_data: Vec<Vec<f64>> = vec![];
        let mut c_data: Vec<Vec<f64>> = vec![];
        let mut vecl_data: Vec<f64> = vec![];

        let mut a_rows: Vec<&mut [f64]> = vec![];
        let mut b_rows: Vec<&mut [f64]> = vec![];
        let mut c_rows: Vec<&mut [f64]> = vec![];

        let mut state = BrtezState {
            a: &mut a_rows,
            b: &mut b_rows,
            c: &mut c_rows,
            vecl: &mut vecl_data,
        };

        brtez(&params, &mut state);
        // 应该不崩溃
    }
}