//! 对流温度计算模块。
//!
//! 重构自 TLUSTY `contmp.f`
//!
//! # 功能
//!
//! LTEGR 的辅助过程，确定对流不稳定层中的温度：
//! - 求解能量平衡方程 F(rad) + F(conv) = F(mech)
//! - 这会产生关于对数温度梯度的三次方程
//! - 迭代计算温度分布、电子密度、平均不透明度

use crate::tlusty::state::constants::{BOLK, HALF, SIG4P, UN, MDEPTH};
// f2r_depends: CONOUT, CONVEC, CUBIC, ELDENS, MEANOP, MEANOPT, OPACF0, STEQEQ, WNSTOR

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

/// 温度相对误差容限
const ERRT: f64 = 1e-3;

// ============================================================================
// 配置结构体
// ============================================================================

/// CONTMP 配置参数。
#[derive(Debug, Clone)]
pub struct ContmpConfig {
    /// 深度点数 (ND)
    pub nd: usize,
    /// 频率点数 (NFREQ)
    pub nfreq: usize,
    /// 有效温度 (TEFF)
    pub teff: f64,
    /// 表面重力 (GRAV)
    pub grav: f64,
    /// 引力参数 (QGRAV) - 盘模式使用
    pub qgrav: f64,

    // 控制标志
    /// 盘模式标志 (IDISK)
    pub idisk: i32,
    /// 不透明度表标志 (IOPTAB)
    pub ioptab: i32,
    /// 对数梯度标志 (ILGDER)
    pub ilgder: i32,
    /// 辐射压标志 (IFPRAD)
    pub ifprad: i32,
    /// 对流迭代次数 (NCONIT)
    pub nconit: i32,
    /// 对流起始深度 (IDCONZ)
    pub idconz: usize,

    // 对流参数
    /// 混合长度参数 (HMIX0)
    pub hmix0: f64,
    /// 对流常数 A (ACONML)
    pub aconml: f64,
    /// 对流常数 B (BCONML)
    pub bconml: f64,
    /// 对流常数 C (CCONML)
    pub cconml: f64,

    // 打印控制
    /// 打印标志 (IPRING)
    pub ipring: i32,

    // 通道模式
    /// 通道模式标志 (ICHANM)
    pub ichanm: i32,
}

impl Default for ContmpConfig {
    fn default() -> Self {
        Self {
            nd: 50,
            nfreq: 1000,
            teff: 10000.0,
            grav: 1e4,
            qgrav: 0.0,
            idisk: 0,
            ioptab: 0,
            ilgder: 0,
            ifprad: 1,
            nconit: 3,
            idconz: 1,
            hmix0: 1.0,
            aconml: 1.0,
            bconml: 1.0,
            cconml: 1.0,
            ipring: 0,
            ichanm: 0,
        }
    }
}

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

/// CONTMP 输入参数（可变引用）。
pub struct ContmpParams<'a> {
    /// 配置
    pub config: ContmpConfig,

    // 深度相关数组 (nd)
    /// 温度 (TEMP) - 会被修改
    pub temp: &'a mut [f64],
    /// 电子密度 (ELEC) - 会被修改
    pub elec: &'a mut [f64],
    /// 总粒子密度 (DENS) - 会被修改
    pub dens: &'a mut [f64],
    /// 分子质量 (WMM)
    pub wmm: &'a [f64],
    /// 深度 (柱质量密度, DM) - 可能被修改
    pub dm: &'a mut [f64],
    /// 深度变量 (ZD) - 盘模式使用
    pub zd: &'a [f64],
    /// 角度 (THETA) - 盘模式使用
    pub theta: &'a [f64],
    /// 总压力 (PTOTAL) - 可能被修改
    pub ptotal: &'a mut [f64],
    /// 气压 (PGS) - 会被修改
    pub pgs: &'a mut [f64],
    /// 辐射压 (PRADT) - 会被修改
    pub pradt: &'a mut [f64],
    /// 湍流速度 (VTURB)
    pub vturb: &'a [f64],
    /// Rosseland 光学深度 (TAUROS)
    pub tauros: &'a [f64],
    /// Rosseland 不透明度/密度 (ABROSD) - 会被修改
    pub abrosd: &'a mut [f64],
    /// Planck 不透明度/密度 (ABPLAD) - 会被修改
    pub abplad: &'a mut [f64],
    /// 氢分子压力 (PHMOL) - 会被修改
    pub phmol: &'a mut [f64],
}

/// CONTMP 输出结果。
#[derive(Debug, Clone)]
pub struct ContmpOutput {
    /// 最大温度变化
    pub chantm: f64,
    /// 迭代次数
    pub iconit: i32,
}

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

/// 简化的对流通量计算（避免复杂依赖）
fn compute_convec_simplified(
    t: f64,
    ptot: f64,
    pg: f64,
    prad: f64,
    delta: f64,
    abros: f64,
    hmix0: f64,
    grav: f64,
) -> (f64, f64, f64, f64, f64) {
    // 简化的热力学导数
    let rho = ptot / (BOLK * t * (1.0 + ptot / pg));
    let grdadb = 0.4;  // 简化的绝热梯度
    let heatcp = 1.5 * BOLK;  // 简化的定压比热
    let dlrdlt = 1.0;  // 简化的密度导数

    let ddel = delta - grdadb;
    if ddel < 0.0 || hmix0 < 0.0 {
        return (0.0, 0.0, grdadb, rho, heatcp);
    }

    // 压力标高
    let hscale = ptot / rho / grav;
    let hmix = if hmix0 == 0.0 { 1.0 } else { hmix0 };

    // 对流参数
    let vco = hmix * (ptot / rho * dlrdlt).abs().sqrt();
    let flco = rho * heatcp * t * hmix / 12.566370614359172;

    // 光学厚度
    let taue = hmix * abros * rho * hscale;
    let fac = taue / (UN + HALF * taue * taue);

    // 简化的辐射耗散
    let b = 5.67e-5 * t.powi(3) / (rho * heatcp * vco) * fac * HALF;

    let d = b * b / 2.0;
    let dlt = if d / 2.0 + ddel >= 0.0 {
        d + ddel - b * (d / 2.0 + ddel).sqrt()
    } else {
        d + ddel
    };

    let dlt = if dlt < 0.0 { 0.0 } else { dlt };

    let vconv = vco * dlt.sqrt();
    let flxcnv = flco * vconv * dlt;

    (flxcnv, vconv, grdadb, rho, heatcp)
}

/// 简化的三次方程求解
fn solve_cubic(a: f64, b: f64, del: f64, grdadb: f64) -> f64 {
    const THIRD: f64 = 1.0 / 3.0;

    if a.abs() < 1e-30 {
        return grdadb + del;
    }

    let aa = THIRD / a;
    let bb = b / a;
    let cc = -del / a;

    let p = bb * THIRD - aa * aa;
    let q = aa.powi(3) - (bb * aa - cc) / 2.0;
    let d = q * q + p * p * p;

    let sol = if d > 0.0 {
        let d_sqrt = d.sqrt();
        if (d_sqrt - q.abs()) < 1e-14 * d_sqrt {
            (2.0 * d_sqrt).powf(THIRD) - aa
        } else {
            let d1 = (d_sqrt - q).abs();
            let d2 = (d_sqrt + q).abs();
            d1 / (d_sqrt - q) * d1.powf(THIRD) - d2 / (d_sqrt + q) * d2.powf(THIRD) - aa
        }
    } else {
        let cosf = -q / (p * p * p).abs().sqrt();
        let tanf = (UN - cosf * cosf).sqrt() / cosf;
        let fi = tanf.atan() * THIRD;
        2.0 * p.abs().sqrt() * fi.cos() - aa
    };

    // Newton-Raphson 修正
    let delda = sol * (b + sol);
    if delda > del || delda < 0.0 {
        let mut x0 = sol;
        for _ in 0..50 {
            let delx = (del - x0 * (b + x0 + a * x0 * x0))
                / (3.0 * a * x0 * x0 + 2.0 * x0 + b);
            x0 += delx;
            if (delx / x0).abs() < 1e-6 {
                break;
            }
        }
        grdadb + b * x0 + x0 * x0
    } else {
        grdadb + b * sol + sol * sol
    }
}

/// 简化的电子密度计算
fn compute_eldens_simplified(t: f64, pg: f64, wmm: f64) -> (f64, f64, f64) {
    // Saha 方程简化
    let chi_h = 13.6;  // 氢电离能 (eV)
    let k_t = 8.617e-5 * t;  // kT in eV

    let saha_factor = 2.415e15 * t.powf(1.5) * (-chi_h / k_t).exp();
    let an = pg / (BOLK * t);

    // 求解二次方程
    let b = saha_factor;
    let c = -an;
    let discriminant = b * b + 4.0 * c.abs();
    let ane = if discriminant > 0.0 {
        0.5 * (-b + discriminant.sqrt())
    } else {
        0.0
    };
    let ane = ane.max(0.0).min(an);

    let dens = wmm * (an - ane);
    let ahmol = 0.0;  // 简化的氢分子丰度

    (ane, dens, ahmol)
}

/// 简化的平均不透明度计算
fn compute_mean_opacity_simplified(t: f64, rho: f64) -> (f64, f64) {
    // Kramers 不透明度近似
    let abros = 4.34e24 * (1.0 + 0.7) * (rho + 1e-30).sqrt() * t.powf(-3.5);
    let abpla = abros * 1.2;  // Planck 平均略高于 Rosseland
    (abros, abpla)
}

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

/// 计算对流不稳定层中的温度 (CONTMP)。
///
/// 通过求解能量平衡方程 F(rad) + F(conv) = F(mech) 来确定温度，
/// 这会产生关于对数温度梯度的三次方程。
///
/// # 参数
///
/// * `params` - 输入参数（包含可变引用）
///
/// # 返回值
///
/// 返回 `ContmpOutput`，包含最大温度变化和迭代次数。
pub fn contmp(params: &mut ContmpParams) -> ContmpOutput {
    let nd = params.config.nd;
    let teff = params.config.teff;

    // 创建临时数组
    let mut deltr = vec![0.0; MDEPTH];
    let mut tempr = vec![0.0; MDEPTH];
    let mut icon0 = vec![0i32; MDEPTH];

    // 初始化常量
    let t4 = teff.powi(4);
    let flxto0 = SIG4P * t4;
    let mut dprad = 1.891204931e-15 * t4;
    if params.config.ifprad == 0 {
        dprad = 0.0;
    }
    let prad0 = dprad / 1.732;

    // 存储纯辐射平衡模型的温度和梯度
    for id in 0..nd {
        tempr[id] = params.temp[id];
        if id == 0 {
            deltr[id] = 0.0;
        } else {
            let ptotal_id = params.ptotal[id];
            let ptotal_idm1 = params.ptotal[id - 1];
            let temp_id = params.temp[id];
            let temp_idm1 = params.temp[id - 1];

            if params.config.ilgder == 0 {
                // 线性梯度
                let denom = ptotal_id - ptotal_idm1;
                if denom.abs() > 1e-30 {
                    deltr[id] = (temp_id - temp_idm1) / denom
                        * (ptotal_id + ptotal_idm1)
                        / (temp_id + temp_idm1);
                }
            } else {
                // 对数梯度
                let denom = (ptotal_id / ptotal_idm1).ln();
                if denom.abs() > 1e-30 {
                    deltr[id] = (temp_id / temp_idm1).ln() / denom;
                }
            }
        }
    }

    let mut iconit = 0;
    let mut chantm = 0.0;
    let mut deltc = 0.0;

    // 全局迭代循环
    loop {
        iconit += 1;
        let mut iconbe = 0;
        chantm = 0.0;

        let mut pgm = 0.0;
        let mut pradm = 0.0;
        let mut deltt0 = 0.0;

        for id in 0..nd {
            let t = params.temp[id];
            let ptot = params.ptotal[id];
            let pgas = params.pgs[id];
            let pturb = HALF * params.dens[id] * params.vturb[id] * params.vturb[id];
            let mut prad = ptot - pgas - pturb;
            params.pradt[id] = prad;

            let mut flxtot = flxto0;
            let mut gravd = 0.0;
            if params.config.idisk == 1 {
                flxtot = flxto0 * (UN - params.theta[id]);
                gravd = params.zd[id] * params.config.qgrav;
            }

            icon0[id] = 0;

            if id == 0 {
                deltt0 = params.temp[id] - t;
                pgm = pgas;
                pradm = prad;
                continue;
            }

            let mut j = 0;
            let mut tm_val = t;
            if iconit == 1 {
                tm_val = t - tempr[id - 1] + params.temp[id - 1];
            }
            let tm = params.temp[id - 1];
            let mut t_val = tm_val;
            if t_val < 0.0 {
                t_val = tm;
            }
            let ptotm = params.ptotal[id - 1];

            let pt0 = if params.config.ilgder == 0 {
                HALF * (ptot + ptotm)
            } else {
                (ptot * ptotm).sqrt()
            };

            let delr = deltr[id];

            // 内层迭代循环
            loop {
                j += 1;
                let told = t_val;

                let (t0, pg0, pr0, ab0) = if params.config.ilgder == 0 {
                    (
                        HALF * (t_val + tm),
                        HALF * (pgas + pgm),
                        HALF * (prad + pradm),
                        HALF * (params.abrosd[id] + params.abrosd[id - 1]),
                    )
                } else {
                    (
                        (t_val * tm).sqrt(),
                        (pgas * pgm).sqrt(),
                        (prad * pradm).sqrt(),
                        (params.abrosd[id] * params.abrosd[id - 1]).sqrt(),
                    )
                };

                let mut pgas_new = pgas;

                if id >= nd - 2 && iconbe == 0 {
                    // 跳过对流计算
                } else {
                    // 调用简化的对流计算
                    let grav_use = if params.config.idisk == 1 {
                        if gravd == 0.0 {
                            params.config.grav  // 使用默认引力
                        } else {
                            gravd
                        }
                    } else {
                        params.config.grav
                    };

                    let (flxcnv, vconv, grdadb, rho, _heatcp) = compute_convec_simplified(
                        t0, pt0, pg0, pr0, delr, ab0,
                        params.config.hmix0, grav_use,
                    );

                    if flxcnv == 0.0 || id < params.config.idconz {
                        // 没有对流
                    } else {
                        icon0[id] = 1;
                        iconbe = 1;

                        // 构建三次方程参数
                        let a = if flxtot > 0.0 {
                            flxcnv * vconv / flxtot * delr
                        } else {
                            0.0
                        };
                        let b = 5.67e-5 * t0.powi(3) / (rho + 1e-30) / vconv * HALF;

                        let delta0 = solve_cubic(a, b, grdadb - delr, grdadb);

                        let reff = if delr.abs() > 1e-30 {
                            delta0 / delr
                        } else {
                            UN
                        };
                        prad = pradm + (params.tauros[id] - params.tauros[id - 1]) * dprad * reff;
                        params.pradt[id] = prad;
                        pgas_new = ptot - prad - pturb;

                        if params.config.ilgder == 0 {
                            let reff_clamped = reff.clamp(0.0, UN);
                            let fac = delta0 * (ptot - ptotm) / (ptot + ptotm);
                            t_val = tm * (UN + fac) / (UN - fac);
                            if t_val < tm {
                                t_val = tm;
                            }
                        } else {
                            t_val = tm * (ptot / ptotm).powf(delta0);
                            if t_val < tm {
                                t_val = tm * 1.0001;
                            }
                        }

                        if ((UN - t_val / told).abs() > ERRT) && (j < 10) {
                            continue;
                        }
                    }
                }

                // 存储最终量
                if id > 0 && icon0[id] == 0 && icon0[id - 1] == 1 {
                    deltc = deltt0;
                }
                deltt0 = params.temp[id] - t_val;

                let mut chant0 = 0.0;
                if params.temp[id] != 0.0 {
                    chant0 = ((t_val - params.temp[id]) / params.temp[id]).abs();
                }
                if chant0 > chantm {
                    chantm = chant0;
                }

                params.temp[id] = t_val;
                if iconit > 1 && icon0[id] == 0 && iconbe == 1 {
                    params.temp[id] = t_val - deltc;
                }

                pgm = pgas_new;
                pradm = prad;
                params.pgs[id] = pgas_new;
                break;
            }
        }

        // 更新电子密度、密度和平均不透明度
        // FORMAT 600: diagnostic output for CONTMP iteration
        if params.config.ipring == 2 {
            eprintln!("\n CONVECTIVE FLUX: AT CONTMP, ITER={:2}", iconit);
        }

        for id in 0..nd {
            let t = params.temp[id];
            let p = params.ptotal[id];
            let mut itint = 0;

            loop {
                itint += 1;

                if params.config.ioptab >= -1 {
                    // 使用简化的电子密度计算
                    let (ane, dens, ahmol) =
                        compute_eldens_simplified(t, params.pgs[id], params.wmm[id]);

                    params.elec[id] = ane;
                    params.dens[id] = dens;
                    params.phmol[id] = ahmol;

                    // 使用简化的平均不透明度计算
                    let (abros, abpla) = compute_mean_opacity_simplified(t, dens);
                    params.abrosd[id] = abros;
                    params.abplad[id] = abpla;
                } else {
                    // 简化的状态方程
                    let rho = p / (BOLK * t * (1.0 + 0.7));
                    params.dens[id] = rho;

                    let (abros, abpla) = compute_mean_opacity_simplified(t, rho);
                    params.abrosd[id] = abros;
                    params.abplad[id] = abpla;
                }

                // 更新柱质量
                let ptold = params.ptotal[id];

                if params.config.idisk == 0 && params.config.ichanm > 0 {
                    if id == 0 {
                        params.dm[id] = params.tauros[id] / params.abrosd[id].max(1e-30);
                        params.ptotal[id] = params.dm[id] * params.config.grav + prad0;
                    } else {
                        params.dm[id] = params.dm[id - 1]
                            + (params.tauros[id] - params.tauros[id - 1])
                                / HALF
                                / (params.abrosd[id - 1] + params.abrosd[id]).max(1e-30);
                        params.ptotal[id] = params.dm[id] * params.config.grav + prad0;
                    }

                    let pturb = HALF * params.dens[id] * params.vturb[id] * params.vturb[id];
                    params.pgs[id] = params.ptotal[id] - params.pradt[id] - pturb;
                }

                if ptold.abs() > 1e-30 && (params.ptotal[id] - ptold) / ptold >= 1e-3 {
                    if itint > 5 {
                        // FORMAT 601: slow convergence warning
                        eprintln!("\n SLOW CONVERGENCE OF INTERNAL ITERATIONS IN CONTMP: ID, PTOT(OLD), PTOT(NEW) =\n{:3}{:10.2e}{:10.2e}", id + 1, ptold, params.ptotal[id]);
                        break;
                    } else {
                        continue;
                    }
                }
                break;
            }
        }

        if iconit >= params.config.nconit {
            break;
        }
    }

    ContmpOutput { chantm, iconit }
}

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

    fn create_test_arrays(nd: usize) -> (
        Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>,
        Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>,
        Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>,
    ) {
        let temp: Vec<f64> = (0..nd).map(|i| 10000.0 - i as f64 * 100.0).collect();
        let elec: Vec<f64> = vec![1e12; nd];
        let dens: Vec<f64> = vec![1e15; nd];
        let wmm: Vec<f64> = vec![1.0; nd];
        let dm: Vec<f64> = (0..nd).map(|i| 0.01 * (i + 1) as f64).collect();
        let zd: Vec<f64> = vec![0.0; nd];
        let theta: Vec<f64> = vec![0.0; nd];
        let ptotal: Vec<f64> = (0..nd).map(|i| 1e4 + i as f64 * 1e3).collect();
        let pgs: Vec<f64> = ptotal.clone();
        let vturb: Vec<f64> = vec![0.0; nd];
        let tauros: Vec<f64> = (0..nd).map(|i| 0.1 * (i + 1) as f64).collect();
        let abrosd: Vec<f64> = vec![0.1; nd];
        let abplad: Vec<f64> = vec![0.1; nd];
        let phmol: Vec<f64> = vec![0.0; nd];
        let pradt: Vec<f64> = vec![0.0; nd];

        (temp, elec, dens, wmm, dm, zd, theta, ptotal, pgs, vturb,
         tauros, abrosd, abplad, phmol, pradt)
    }

    #[test]
    fn test_contmp_basic() {
        let nd = 10;
        let (
            mut temp, mut elec, mut dens, wmm, mut dm,
            zd, theta, mut ptotal, mut pgs, vturb,
            tauros, mut abrosd, mut abplad, mut phmol, mut pradt,
        ) = create_test_arrays(nd);

        let config = ContmpConfig {
            nd,
            nfreq: 100,
            teff: 10000.0,
            grav: 1e4,
            nconit: 1,
            ..Default::default()
        };

        let mut params = ContmpParams {
            config,
            temp: &mut temp,
            elec: &mut elec,
            dens: &mut dens,
            wmm: &wmm,
            dm: &mut dm,
            zd: &zd,
            theta: &theta,
            ptotal: &mut ptotal,
            pgs: &mut pgs,
            pradt: &mut pradt,
            vturb: &vturb,
            tauros: &tauros,
            abrosd: &mut abrosd,
            abplad: &mut abplad,
            phmol: &mut phmol,
        };

        let output = contmp(&mut params);

        assert!(output.iconit >= 1);
        assert!(output.chantm.is_finite());
    }

    #[test]
    fn test_contmp_no_convection() {
        let nd = 5;
        let (
            mut temp, mut elec, mut dens, wmm, mut dm,
            zd, theta, mut ptotal, mut pgs, vturb,
            tauros, mut abrosd, mut abplad, mut phmol, mut pradt,
        ) = create_test_arrays(nd);

        let config = ContmpConfig {
            nd,
            hmix0: -1.0,  // 禁用对流
            nconit: 1,
            ..Default::default()
        };

        let mut params = ContmpParams {
            config,
            temp: &mut temp,
            elec: &mut elec,
            dens: &mut dens,
            wmm: &wmm,
            dm: &mut dm,
            zd: &zd,
            theta: &theta,
            ptotal: &mut ptotal,
            pgs: &mut pgs,
            pradt: &mut pradt,
            vturb: &vturb,
            tauros: &tauros,
            abrosd: &mut abrosd,
            abplad: &mut abplad,
            phmol: &mut phmol,
        };

        let output = contmp(&mut params);

        assert!(output.iconit >= 1);
    }

    #[test]
    fn test_contmp_log_gradient() {
        let nd = 5;
        let (
            mut temp, mut elec, mut dens, wmm, mut dm,
            zd, theta, mut ptotal, mut pgs, vturb,
            tauros, mut abrosd, mut abplad, mut phmol, mut pradt,
        ) = create_test_arrays(nd);

        let config = ContmpConfig {
            nd,
            ilgder: 1,
            nconit: 1,
            ..Default::default()
        };

        let mut params = ContmpParams {
            config,
            temp: &mut temp,
            elec: &mut elec,
            dens: &mut dens,
            wmm: &wmm,
            dm: &mut dm,
            zd: &zd,
            theta: &theta,
            ptotal: &mut ptotal,
            pgs: &mut pgs,
            pradt: &mut pradt,
            vturb: &vturb,
            tauros: &tauros,
            abrosd: &mut abrosd,
            abplad: &mut abplad,
            phmol: &mut phmol,
        };

        let output = contmp(&mut params);

        assert!(output.iconit >= 1);
        assert!(output.chantm.is_finite());
    }

    #[test]
    fn test_compute_convec_simplified() {
        let (flxcnv, vconv, grdadb, rho, heatcp) = compute_convec_simplified(
            10000.0,  // t
            1e5,      // ptot
            1e5,      // pg
            0.0,      // prad
            0.5,      // delta
            0.1,      // abros
            1.0,      // hmix0
            1e4,      // grav
        );

        assert!(flxcnv.is_finite());
        assert!(vconv.is_finite());
        assert!(grdadb > 0.0);
        assert!(rho > 0.0);
        assert!(heatcp > 0.0);
    }

    #[test]
    fn test_solve_cubic() {
        let delta = solve_cubic(1.0, 0.5, 0.1, 0.4);
        assert!(delta.is_finite());
        assert!(delta > 0.4);  // 应该大于绝热梯度
    }

    #[test]
    fn test_compute_eldens_simplified() {
        let (ane, dens, ahmol) = compute_eldens_simplified(
            10000.0,  // t
            1e5,      // pg
            1.0,      // wmm
        );

        assert!(ane >= 0.0);
        assert!(dens >= 0.0);
    }

    #[test]
    fn test_compute_mean_opacity_simplified() {
        let (abros, abpla) = compute_mean_opacity_simplified(
            10000.0,  // t
            1e-7,     // rho
        );

        assert!(abros > 0.0);
        assert!(abpla > 0.0);
    }
}