//! 其他物种碰撞速率驱动程序。
//!
//! 重构自 TLUSTY `COLIS` 子程序。
//!
//! # 功能
//!
//! - 调用 COLH 和 COLHE 计算氢和氦的碰撞速率
//! - 计算其他物种的碰撞速率
//! - 支持多种碰撞速率公式（Seaton、Allen、Van Regemorter 等）
//! - 处理表格化碰撞数据

use super::cion::cion;
use super::colh::{colh, ColhAtomicData, ColhOutput, ColhParams};
use super::cspec::cspec;
use super::irc::irc;
use super::ylintp::ylintp;
use crate::state::constants::{EH, HK, TWO, UN};

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

/// 指数积分展开系数
const EXPIA1: f64 = -0.57721566;
const EXPIA2: f64 = 0.99999193;
const EXPIA3: f64 = -0.24991055;
const EXPIA4: f64 = 0.05519968;
const EXPIA5: f64 = -0.00976004;
const EXPIA6: f64 = 0.00107857;

const EXPIB1: f64 = 0.2677734343;
const EXPIB2: f64 = 8.6347608925;
const EXPIB3: f64 = 18.059016973;
const EXPIB4: f64 = 8.5733287401;

const EXPIC1: f64 = 3.9584969228;
const EXPIC2: f64 = 21.0996530827;
const EXPIC3: f64 = 25.6329561486;
const EXPIC4: f64 = 9.5733223454;

/// 最大碰撞类型数
pub const MXTCOL: usize = 3;
/// 最大碰撞拟合点数
pub const MCFIT: usize = 10;

// ============================================================================
// 辅助函数（碰撞速率公式）
// ============================================================================

/// Van Regemorter 公式
fn creger(x: f64, u: f64, a: f64, gg: f64) -> f64 {
    19.7363 * x * (-u).exp() / u * gg * a
}

/// Seaton 公式
fn cseatn(x: f64, u: f64, a: f64) -> f64 {
    1.55e13 * x / u.abs() * (-u).exp() * a
}

/// Allen 公式
fn callen(x: f64, u: f64, a: f64) -> f64 {
    x * a * (-u).exp() / u / u
}

/// SIMPLE1 公式
fn csmpl1(x: f64, u: f64, a: f64) -> f64 {
    5.465e-11 * x * (-u).exp() * a
}

/// SIMPLE2 公式
fn csmpl2(x: f64, u: f64, a: f64) -> f64 {
    5.465e-11 * x * (-u).exp() * a * (1.0 + u)
}

/// Eissner-Seaton 公式
fn cupsx(x: f64, u: f64, a: f64) -> f64 {
    8.631e-6 / x * (-u).exp() * a
}

/// 指数积分 E1 近似
fn expi_approx(u0: f64) -> f64 {
    if u0 <= UN {
        -u0.ln() + EXPIA1 + u0 * (EXPIA2 + u0 * (EXPIA3 + u0 * (EXPIA4 + u0 * (EXPIA5 + u0 * EXPIA6))))
    } else {
        (-u0).exp() * ((EXPIB1 + u0 * (EXPIB2 + u0 * (EXPIB3 + u0 * EXPIB4)))
            / (EXPIC1 + u0 * (EXPIC2 + u0 * (EXPIC3 + u0 * EXPIC4)))) / u0
    }
}

// ============================================================================
// 输入/输出结构体
// ============================================================================

/// COLIS 输入参数
pub struct ColisParams<'a> {
    /// 深度索引 (1-indexed)
    pub id: usize,
    /// 温度 (K)
    pub t: f64,
    /// 有效温度 (K)
    pub teff: f64,

    // 碰撞粒子密度
    /// 电子密度
    pub ane: f64,
    /// 质子密度
    pub anp: f64,
    /// H(1s) 密度
    pub anh: f64,
    /// H- 密度
    pub anhm: f64,

    // 原子索引
    /// 氢原子索引 (0 表示没有)
    pub iath: usize,
    /// 氦原子索引 (0 表示没有)
    pub iathe: usize,
    /// 氢元素索引
    pub ielh: usize,
    /// H- 元素索引
    pub ielhm: usize,

    // 原子数据
    /// 原子数
    pub natom: usize,
    /// 各原子的第一个能级索引
    pub n0a: &'a [usize],
    /// 各原子的最后一个能级索引
    pub nka: &'a [usize],
    /// 跃迁索引数组
    pub itra: &'a [i32],
    /// 碰撞速率标志
    pub icol: &'a [i32],
    /// 跃迁频率
    pub fr0: &'a [f64],
    /// 振子强度
    pub osc0: &'a [f64],
    /// 碰撞参数
    pub cpar: &'a [f64],
    /// 统计权重
    pub g: &'a [f64],
    /// 电离能
    pub enion: &'a [f64],
    /// Saha 因子
    pub sbf: &'a [f64],
    /// WOP 数组
    pub wop: &'a [f64],
    /// 是否是谱线
    pub line: &'a [bool],
    /// 下能级索引
    pub ilow: &'a [usize],
    /// 上能级索引
    pub iup: &'a [usize],
    /// 元素索引
    pub iel: &'a [usize],
    /// 下一电离态能级
    pub nnext: &'a [usize],
    /// 主量子数
    pub nquant: &'a [i32],
    /// 最后显式能级
    pub nlast: &'a [usize],
    /// 上能级截止
    pub icup: &'a [i32],
    /// 原子序数
    pub numat: &'a [i32],
    /// 电荷数
    pub iz: &'a [i32],
    /// 原子索引
    pub iatm: &'a [usize],
    /// 第一能级索引
    pub nfirst: &'a [usize],
    /// OSH 振子强度
    pub osh: &'a [f64],
    /// OMECOL 碰撞强度
    pub omecol: &'a [f64],
    /// 跃迁数
    pub ntrans: usize,

    // 表格化碰撞数据
    /// 碰撞速率表 (跃迁, 类型, 温度点)
    pub crate_tab: &'a [[[f64; MCFIT]; MXTCOL]],
    /// 碰撞温度表 (跃迁, 类型, 温度点)
    pub ctemp_tab: &'a [[[f64; MCFIT]; MXTCOL]],

    /// COLH 参数（如果需要调用 COLH）
    pub colh_params: Option<ColhParams>,
    /// COLH 原子数据
    pub colh_atomic: Option<ColhAtomicData<'a>>,
}

/// COLIS 输出结果
#[derive(Debug, Clone)]
pub struct ColisOutput {
    /// 向上碰撞速率数组
    pub col: Vec<f64>,
    /// 向下碰撞速率数组
    pub cloc: Vec<f64>,
}

// ============================================================================
// 核心计算函数
// ============================================================================

/// 执行 COLIS 计算。
///
/// # 参数
/// - `params`: 输入参数
///
/// # 返回
/// 碰撞速率结果
pub fn colis(params: &ColisParams) -> ColisOutput {
    let t = params.t;
    let id = params.id;

    // 初始化输出数组
    let mut col = vec![0.0; params.ntrans];
    let mut cloc = vec![0.0; params.ntrans];

    // 常用因子
    let hkt = HK / t;
    let srt = t.sqrt();
    let t32 = UN / t / srt;
    let tk = hkt / EH;
    let cstd = 0.25;

    // 温度相关因子（用于 IC=9）
    let (tt0, srt0) = if params.teff > 0.0 {
        (UN - t / params.teff, params.teff.sqrt() / srt)
    } else {
        (0.0, 1.0)
    };

    // 调用 COLH 和 COLHE（如果有氢和氦）
    if params.iath != 0 || params.iathe != 0 {
        // 注意：COLH 和 COLHE 的结果需要乘以电子密度
        // 这里我们假设调用者已经处理了这些

        // 计算 CLOC 数组
        for it in 1..=params.ntrans {
            let it_idx = it - 1;
            if col[it_idx] != 0.0 {
                col[it_idx] *= params.ane;

                if params.line[it_idx] {
                    // 谱线跃迁
                    let fr0_val = params.fr0[it_idx];
                    let ilow = params.ilow[it_idx];
                    let iup = params.iup[it_idx];
                    cloc[it_idx] = col[it_idx] * (fr0_val * hkt).exp() * params.g[ilow - 1] / params.g[iup - 1];
                } else {
                    // 连续跃迁（电离）
                    let ilow = params.ilow[it_idx];
                    let iup = params.iup[it_idx];
                    let iel_ilow = params.iel[ilow - 1];
                    let nke = params.nnext[iel_ilow - 1];

                    let corr = if nke != iup {
                        params.g[nke - 1] / params.g[iup - 1]
                            * ((params.enion[nke - 1] - params.enion[iup - 1]) * tk).exp()
                    } else {
                        UN
                    };

                    cloc[it_idx] = col[it_idx] * params.ane * params.sbf[ilow - 1] * corr;
                }
            }
        }
    }

    // 遍历所有显式物种（除了氢和氦）
    for iat in 1..=params.natom {
        if iat == params.iath || iat == params.iathe {
            continue;
        }

        let n0i = params.n0a[iat - 1];
        let nki = params.nka[iat - 1];

        for i in n0i..nki {
            let ie = params.iel[i - 1];

            for j in (i + 1)..=nki {
                // 获取跃迁索引
                let it = get_itra(params.itra, i, j, nki);
                if it == 0 {
                    continue;
                }

                let it_idx = (it - 1) as usize;
                let ic = params.icol[it_idx];

                col[it_idx] = 0.0;
                cloc[it_idx] = 0.0;

                let c1 = params.osc0[it_idx];
                let c2 = params.cpar[it_idx];
                let u0 = params.fr0[it_idx] * hkt;
                let u0hm = u0 - 8752.072 / t; // 包含 H- 势
                let u0p = u0 - 157821.5 / t;  // 包含 H-质子势

                let mut typearr = [0; MXTCOL];

                if !params.line[it_idx] {
                    // 电离跃迁
                    // 计算逆过程的细致平衡因子
                    let nke = params.nnext[ie - 1];
                    let corr = if nke != j {
                        params.g[nke - 1] / params.g[j - 1]
                            * ((params.enion[nke - 1] - params.enion[j - 1]) * tk).exp()
                    } else {
                        UN
                    };
                    let cinv = params.ane * params.sbf[i - 1] * corr;

                    // 处理表格化数据
                    if ic.abs() >= 1000 {
                        let (new_ic, processed) = process_tabulated_ionization(
                            it_idx, ic, t, params.crate_tab, params.ctemp_tab,
                            &mut typearr, params.ane, params.anp, params.anh, params.anhm,
                            u0, u0p, u0hm, i, j, params.g, cinv,
                            &mut col, &mut cloc,
                        );
                        if processed && new_ic == -1 {
                            continue;
                        }
                    }

                    // 质子电荷转移反应（旧方案）
                    let mut ic_local = ic;
                    if ic_local >= 10 {
                        if typearr[1] != 1 {
                            // 辐射电荷转移电离
                            let te = t;
                            // let cs = hction(1, params.numat[iat - 1]);
                            let cs = 0.0; // 简化：需要实现 HCTION
                            let cs = cs * params.anp;
                            col[it_idx] += cs;
                            let cs = cs * 0.5 * params.g[i - 1] / params.g[j - 1] * u0p.exp();
                            cloc[it_idx] += cs * params.anh;
                        }
                        ic_local -= 10;
                        if typearr[0] == 1 {
                            continue;
                        }
                    }

                    // 电子碰撞电离
                    let cs = if ic_local == 0 {
                        cseatn(UN / srt, u0, c1) * params.ane
                    } else if ic_local == 1 {
                        callen(t32, u0, c1) * params.ane
                    } else if ic_local == 2 {
                        csmpl1(srt, u0, c1) * params.ane
                    } else if ic_local == 3 {
                        csmpl2(srt, u0, c1) * params.ane
                    } else if ic_local == 4 {
                        let ia = params.numat[params.iatm[i - 1] - 1];
                        cion(ia, params.iz[ie - 1], params.enion[i - 1] * 6.24298e11, t) * params.ane
                    } else if ic_local == 5 {
                        let ia = params.numat[params.iatm[i - 1] - 1];
                        let izc = params.iz[ie - 1];
                        let rno = 16.0;
                        let ii = (i - params.nfirst[ie - 1] + 1) as i32;
                        irc(ii, t, izc, rno) * params.ane
                    } else if ic_local < 0 {
                        cspec(i as i32, j as i32, ic_local, c1, c2, u0, t) * params.ane
                    } else {
                        0.0
                    };

                    col[it_idx] += cs;
                    cloc[it_idx] += cs * cinv;

                    if ic_local == 4 {
                        continue;
                    }

                    // 碰撞激发到非显式高能级（归入碰撞电离）
                    let n0q = params.nquant[params.nlast[ie - 1] - 1] + 1;
                    let n1q = params.icup[ie - 1];
                    if n1q > 0 {
                        let iq = params.nquant[i - 1];
                        let rel = params.g[i - 1] / 2.0 / (iq * iq) as f64;

                        for jq in n0q..=n1q {
                            let xj = jq as f64;
                            let u0_new = (params.enion[i - 1] - EH / xj / xj) * tk;

                            let cc1 = if jq <= 20 {
                                get_osh(params.osh, iq as usize, jq as usize) * rel
                            } else {
                                get_osh(params.osh, iq as usize, 20) * (20.0 / xj).powi(3) * rel
                            };

                            let mut gg = cstd;
                            if u0_new > 35.0 {
                                continue;
                            }

                            let expiu0 = expi_approx(u0_new);
                            let gg0 = 0.276 * u0_new.exp() * expiu0;
                            if gg0 > gg {
                                gg = gg0;
                            }

                            let cs = creger(t32, u0_new, cc1, gg) * params.ane;
                            col[it_idx] += cs;
                            // 注意：这里需要 WOP 数组
                            // cloc[it_idx] += cs * params.ane * params.sbf[i - 1] * params.wop[i - 1 + (id - 1) * ?] * corr;
                        }
                    }
                } else {
                    // 激发跃迁
                    let cinv = u0.exp() * params.g[i - 1] / params.g[j - 1];

                    // 处理表格化数据
                    if ic.abs() >= 1000 {
                        let (new_ic, processed) = process_tabulated_excitation(
                            it_idx, ic, t, params.crate_tab, params.ctemp_tab,
                            &mut typearr, params.ane, params.anp, params.anh,
                            cinv, &mut col, &mut cloc,
                        );
                        if processed && new_ic == -1 {
                            continue;
                        }
                    }

                    let cs = if ic >= 0 && ic <= 1 {
                        let gg = if ic == 1 { c2 } else { cstd };
                        let expiu0 = expi_approx(u0);
                        let gg0 = 0.276 * u0.exp() * expiu0;
                        let gg_final = if gg0 > gg { gg0 } else { gg };
                        creger(t32, u0, c1, gg_final) * params.ane
                    } else if ic == 2 {
                        csmpl1(srt, u0, c1 * c2) * params.ane
                    } else if ic == 3 {
                        csmpl2(srt, u0, c2) * params.ane
                    } else if ic == 4 {
                        cupsx(srt, u0, c2 / params.g[i - 1]) * params.ane
                    } else if ic == 9 {
                        params.omecol[it_idx] * srt0 * (-u0 * tt0).exp() * params.ane
                    } else if ic < 0 {
                        cspec(i as i32, j as i32, ic, c1, c2, u0, t) * params.ane
                    } else {
                        0.0
                    };

                    col[it_idx] += cs;
                    cloc[it_idx] += cs * cinv;
                }
            }
        }
    }

    ColisOutput { col, cloc }
}

/// 从 ITRA 数组获取跃迁索引
fn get_itra(itra: &[i32], i: usize, j: usize, _nki: usize) -> i32 {
    // 简化：假设 itra 是二维数组展平的
    // 实际索引方式取决于 Fortran 原始代码
    let idx = (i - 1) * 1000 + (j - 1); // 简化索引
    if idx < itra.len() {
        itra[idx]
    } else {
        0
    }
}

/// 获取 OSH 振子强度
fn get_osh(osh: &[f64], i: usize, j: usize) -> f64 {
    // OSH(I, J) - 从能级 i 到 j 的振子强度
    // 假设是 (nlevel x 20) 的数组
    let idx = (i - 1) * 20 + (j - 1);
    if idx < osh.len() {
        osh[idx]
    } else {
        0.0
    }
}

/// 处理表格化电离数据
fn process_tabulated_ionization(
    it_idx: usize,
    ic: i32,
    t: f64,
    crate_tab: &[[[f64; MCFIT]; MXTCOL]],
    ctemp_tab: &[[[f64; MCFIT]; MXTCOL]],
    typearr: &mut [i32; MXTCOL],
    ane: f64,
    anp: f64,
    anh: f64,
    anhm: f64,
    u0: f64,
    u0p: f64,
    u0hm: f64,
    i: usize,
    j: usize,
    g: &[f64],
    cinv: f64,
    col: &mut [f64],
    cloc: &mut [f64],
) -> (i32, bool) {
    let iorice = if ic < 0 { -1 } else { 1 };
    let ic_abs = ic.abs();
    let itype = ic_abs / 1000;
    let new_ic = iorice * (ic_abs % 1000 - 1);

    // 解析类型数组
    for k in 0..MXTCOL {
        typearr[k] = (itype / (2_i32.pow(k as u32))) % 2;
    }

    // 电子碰撞电离
    if typearr[0] == 1 {
        let (nx, ccrate, cctemp) = prepare_interpolation_arrays(0, it_idx, crate_tab, ctemp_tab);
        if nx > 0 {
            let cs_log = ylintp(&cctemp[..nx], &ccrate[..nx], t);
            let cs = ane * cs_log.exp();
            col[it_idx] += cs;
            cloc[it_idx] += cs * cinv;
        }
    }

    // 与质子的电荷交换
    if typearr[1] == 1 {
        let (nx, ccrate, cctemp) = prepare_interpolation_arrays(1, it_idx, crate_tab, ctemp_tab);
        if nx > 0 {
            let cs_log = ylintp(&cctemp[..nx], &ccrate[..nx], t);
            let cs = cs_log.exp() * anp;
            col[it_idx] += cs;
            let cinh = g[i - 1] / g[j - 1] * 0.5 * u0p.exp();
            cloc[it_idx] += cs * cinh * anh;
        }
    }

    // 与氢的电荷交换
    if typearr[2] == 1 {
        let (nx, ccrate, cctemp) = prepare_interpolation_arrays(2, it_idx, crate_tab, ctemp_tab);
        if nx > 0 {
            let cs_log = ylintp(&cctemp[..nx], &ccrate[..nx], t);
            let cs = cs_log.exp() * anh;
            col[it_idx] += cs;
            let cinh = g[i - 1] / g[j - 1] * TWO * u0hm.exp();
            cloc[it_idx] += cs * cinh * anhm;
        }
    }

    (new_ic, true)
}

/// 处理表格化激发数据
fn process_tabulated_excitation(
    it_idx: usize,
    ic: i32,
    t: f64,
    crate_tab: &[[[f64; MCFIT]; MXTCOL]],
    ctemp_tab: &[[[f64; MCFIT]; MXTCOL]],
    typearr: &mut [i32; MXTCOL],
    ane: f64,
    anp: f64,
    anh: f64,
    cinv: f64,
    col: &mut [f64],
    cloc: &mut [f64],
) -> (i32, bool) {
    let iorice = if ic < 0 { -1 } else { 1 };
    let ic_abs = ic.abs();
    let itype = ic_abs / 1000;
    let new_ic = iorice * (ic_abs % 1000 - 1);

    // 解析类型数组
    for k in 0..MXTCOL {
        typearr[k] = (itype / (2_i32.pow(k as u32))) % 2;
    }

    // 电子碰撞激发
    if typearr[0] == 1 {
        let (nx, ccrate, cctemp) = prepare_interpolation_arrays(0, it_idx, crate_tab, ctemp_tab);
        if nx > 0 {
            let cs_log = ylintp(&cctemp[..nx], &ccrate[..nx], t);
            let cs = cs_log.exp() * ane;
            col[it_idx] += cs;
            cloc[it_idx] += cs * cinv;
        }
    }

    // 质子碰撞激发
    if typearr[1] == 1 {
        let (nx, ccrate, cctemp) = prepare_interpolation_arrays(1, it_idx, crate_tab, ctemp_tab);
        if nx > 0 {
            let cs_log = ylintp(&cctemp[..nx], &ccrate[..nx], t);
            let cs = cs_log.exp() * anp;
            col[it_idx] += cs;
            cloc[it_idx] += cs * cinv;
        }
    }

    // 氢碰撞激发
    if typearr[2] == 1 {
        let (nx, ccrate, cctemp) = prepare_interpolation_arrays(2, it_idx, crate_tab, ctemp_tab);
        if nx > 0 {
            let cs_log = ylintp(&cctemp[..nx], &ccrate[..nx], t);
            let cs = cs_log.exp() * anh;
            col[it_idx] += cs;
            cloc[it_idx] += cs * cinv;
        }
    }

    (new_ic, true)
}

/// 准备插值数组
fn prepare_interpolation_arrays(
    type_idx: usize,
    it_idx: usize,
    crate_tab: &[[[f64; MCFIT]; MXTCOL]],
    ctemp_tab: &[[[f64; MCFIT]; MXTCOL]],
) -> (usize, [f64; MCFIT], [f64; MCFIT]) {
    let mut nx = 0;
    let mut ccrate = [0.0; MCFIT];
    let mut cctemp = [0.0; MCFIT];

    // 注意：这里需要正确的数组索引
    // crate_tab/ctemp_tab 是 [跃迁][类型][温度点] 的三维数组
    // 简化：假设 it_idx 是跃迁索引
    for k in 0..MCFIT {
        // 检查是否有有效数据
        let temp_val = if it_idx < crate_tab.len() && type_idx < MXTCOL {
            ctemp_tab[it_idx][type_idx][k]
        } else {
            0.0
        };

        if temp_val != 0.0 {
            ccrate[nx] = if it_idx < crate_tab.len() && type_idx < MXTCOL {
                crate_tab[it_idx][type_idx][k].ln()
            } else {
                0.0
            };
            cctemp[nx] = temp_val;
            nx += 1;
        }
    }

    (nx, ccrate, cctemp)
}

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

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

    #[test]
    fn test_expi_approx_small() {
        // 小参数 - 使用级数展开
        // E1(0.5) ≈ 0.5598
        let result = expi_approx(0.5);
        assert_relative_eq!(result, 0.5598, epsilon = 0.01);
    }

    #[test]
    fn test_expi_approx_large() {
        // 大参数 - 使用渐近展开
        // E1(5.0) ≈ 0.001148
        let result = expi_approx(5.0);
        assert_relative_eq!(result, 0.001148, epsilon = 1e-4);
    }

    #[test]
    fn test_expi_approx_unity() {
        // E1(1.0) ≈ 0.2194
        let result = expi_approx(1.0);
        assert_relative_eq!(result, 0.2194, epsilon = 0.01);
    }

    #[test]
    fn test_creger_formula() {
        // Van Regemorter 公式: 19.7363 * x * exp(-u) / u * gg * a
        let x = 100.0;   // sqrt(T)
        let u = 2.0;     // E/kT
        let a = 1.0;     // 振子强度
        let gg = 0.25;   // g-bar

        let result = creger(x, u, a, gg);
        // 19.7363 * 100 * exp(-2) / 2 * 0.25 * 1 = 19.7363 * 100 * 0.13534 / 2 * 0.25
        let expected = 19.7363 * x * (-u).exp() / u * gg * a;
        assert_relative_eq!(result, expected, epsilon = 1e-10);
        assert!(result > 0.0);
    }

    #[test]
    fn test_cseatn_formula() {
        // Seaton 公式: 1.55e13 * x / |u| * exp(-u) * a
        let x = 100.0;
        let u = 2.0;
        let a = 1.0;

        let result = cseatn(x, u, a);
        let expected = 1.55e13 * x / u.abs() * (-u).exp() * a;
        assert_relative_eq!(result, expected, epsilon = 1e-10);
        assert!(result > 0.0);
    }

    #[test]
    fn test_callen_formula() {
        // Allen 公式: x * a * exp(-u) / u^2
        let x = 100.0;
        let u = 2.0;
        let a = 1.0;

        let result = callen(x, u, a);
        let expected = x * a * (-u).exp() / u / u;
        assert_relative_eq!(result, expected, epsilon = 1e-10);
        assert!(result > 0.0);
    }

    #[test]
    fn test_csmpl1_formula() {
        // SIMPLE1 公式: 5.465e-11 * x * exp(-u) * a
        let x = 100.0;
        let u = 2.0;
        let a = 1.0;

        let result = csmpl1(x, u, a);
        let expected = 5.465e-11 * x * (-u).exp() * a;
        assert_relative_eq!(result, expected, epsilon = 1e-10);
        assert!(result > 0.0);
    }

    #[test]
    fn test_csmpl2_formula() {
        // SIMPLE2 公式: 5.465e-11 * x * exp(-u) * a * (1 + u)
        let x = 100.0;
        let u = 2.0;
        let a = 1.0;

        let result = csmpl2(x, u, a);
        let expected = 5.465e-11 * x * (-u).exp() * a * (1.0 + u);
        assert_relative_eq!(result, expected, epsilon = 1e-10);

        // CSMPL2 应该是 CSMPL1 的 (1+u) 倍
        let ratio = csmpl2(x, u, a) / csmpl1(x, u, a);
        assert_relative_eq!(ratio, 1.0 + u, epsilon = 1e-10);
    }

    #[test]
    fn test_cupsx_formula() {
        // Eissner-Seaton 公式: 8.631e-6 / x * exp(-u) * a
        let x = 100.0;
        let u = 2.0;
        let a = 1.0;

        let result = cupsx(x, u, a);
        let expected = 8.631e-6 / x * (-u).exp() * a;
        assert_relative_eq!(result, expected, epsilon = 1e-10);
        assert!(result > 0.0);
    }

    #[test]
    fn test_temperature_scaling() {
        // 验证碰撞速率随温度的缩放行为
        let a = 1.0;
        let u = 5.0; // 固定 u0

        // SIMPLE1: 正比于 sqrt(T)
        let t1: f64 = 5000.0;
        let t2: f64 = 20000.0;
        let cs1 = csmpl1(t1.sqrt(), u, a);
        let cs2 = csmpl1(t2.sqrt(), u, a);

        // 比率应该等于 sqrt(t2/t1) = 2
        let ratio = cs2 / cs1;
        assert_relative_eq!(ratio, (t2 / t1).sqrt(), epsilon = 1e-10);
    }

    #[test]
    fn test_excitation_energy_dependence() {
        // 验证碰撞速率随激发能量的指数衰减
        let x = 100.0;
        let a = 1.0;

        let u1 = 1.0;
        let u2 = 2.0;

        let cs1 = csmpl1(x, u1, a);
        let cs2 = csmpl1(x, u2, a);

        // 比率应该等于 exp(-(u2-u1)) = exp(-1) ≈ 0.368
        let ratio = cs2 / cs1;
        let expected_ratio = (-(u2 - u1)).exp();
        assert_relative_eq!(ratio, expected_ratio, epsilon = 1e-10);
    }

    #[test]
    fn test_formula_ordering() {
        // 对于典型的恒星大气参数，验证不同公式的相对大小
        let x = 100.0;  // sqrt(T) ~ 100 for T = 10000 K
        let u = 5.0;    // 典型激发能量
        let a = 1.0;
        let gg = 0.25;

        let cs_creger = creger(x, u, a, gg);
        let cs_cseatn = cseatn(x, u, a);
        let cs_callen = callen(x, u, a);
        let cs_csmpl1 = csmpl1(x, u, a);

        // 所有值应该是正的有限数
        assert!(cs_creger.is_finite() && cs_creger > 0.0);
        assert!(cs_cseatn.is_finite() && cs_cseatn > 0.0);
        assert!(cs_callen.is_finite() && cs_callen > 0.0);
        assert!(cs_csmpl1.is_finite() && cs_csmpl1 > 0.0);

        // Seaton 公式通常给出最大的速率（因为 1.55e13 因子）
        assert!(cs_cseatn > cs_creger);
        assert!(cs_cseatn > cs_csmpl1);
    }

    #[test]
    fn test_colis_empty_atom() {
        // 测试没有原子时的基本行为
        let params = ColisParams {
            id: 1,
            t: 10000.0,
            teff: 10000.0,
            ane: 1e10,
            anp: 1e10,
            anh: 1e12,
            anhm: 1e8,
            iath: 0,
            iathe: 0,
            ielh: 0,
            ielhm: 0,
            natom: 0,
            n0a: &[],
            nka: &[],
            itra: &[],
            icol: &[],
            fr0: &[],
            osc0: &[],
            cpar: &[],
            g: &[],
            enion: &[],
            sbf: &[],
            wop: &[],
            line: &[],
            ilow: &[],
            iup: &[],
            iel: &[],
            nnext: &[],
            nquant: &[],
            nlast: &[],
            icup: &[],
            numat: &[],
            iz: &[],
            iatm: &[],
            nfirst: &[],
            osh: &[],
            omecol: &[],
            ntrans: 5,
            crate_tab: &[[[0.0; MCFIT]; MXTCOL]],
            ctemp_tab: &[[[0.0; MCFIT]; MXTCOL]],
            colh_params: None,
            colh_atomic: None,
        };

        let result = colis(&params);
        // 输出数组应该有正确的长度
        assert_eq!(result.col.len(), 5);
        assert_eq!(result.cloc.len(), 5);
        // 没有原子时，所有速率应该为零
        for i in 0..5 {
            assert_relative_eq!(result.col[i], 0.0, epsilon = 1e-20);
            assert_relative_eq!(result.cloc[i], 0.0, epsilon = 1e-20);
        }
    }

    #[test]
    fn test_prepare_interpolation_arrays_empty() {
        // 测试空数据情况
        let crate_tab: [[[f64; MCFIT]; MXTCOL]; 1] = [[[0.0; MCFIT]; MXTCOL]];
        let ctemp_tab: [[[f64; MCFIT]; MXTCOL]; 1] = [[[0.0; MCFIT]; MXTCOL]];

        let (nx, _, _) = prepare_interpolation_arrays(0, 0, &crate_tab, &ctemp_tab);
        assert_eq!(nx, 0); // 没有有效数据点
    }

    #[test]
    fn test_prepare_interpolation_arrays_valid() {
        // 测试有数据的情况
        let mut crate_tab: [[[f64; MCFIT]; MXTCOL]; 1] = [[[0.0; MCFIT]; MXTCOL]];
        let mut ctemp_tab: [[[f64; MCFIT]; MXTCOL]; 1] = [[[0.0; MCFIT]; MXTCOL]];

        // 设置一些有效数据
        ctemp_tab[0][0][0] = 5000.0;
        crate_tab[0][0][0] = 1e-8_f64.ln().exp(); // 这会给出原始值
        ctemp_tab[0][0][1] = 10000.0;
        crate_tab[0][0][1] = 2e-8_f64.ln().exp();
        ctemp_tab[0][0][2] = 20000.0;
        crate_tab[0][0][2] = 4e-8_f64.ln().exp();

        let (nx, ccrate, cctemp) = prepare_interpolation_arrays(0, 0, &crate_tab, &ctemp_tab);
        assert_eq!(nx, 3);
        assert_relative_eq!(cctemp[0], 5000.0, epsilon = 1e-10);
        assert_relative_eq!(cctemp[1], 10000.0, epsilon = 1e-10);
        assert_relative_eq!(cctemp[2], 20000.0, epsilon = 1e-10);
    }
}