//! Opacf0State — OPACF0 连续谱不透明度状态。
//!
//! 预计算 BF 截面参数和 FF 离子数据，用于在 RESOLV 频率循环中
//! 计算不透明度时使用模型种群 (POPUL) 而非重新计算 Saha 平衡。
//!
//! # 与 Fortran OPACF0 的对应
//!
//! Fortran OPACF0(ID, NFRQ) 对每个深度 ID 计算所有频率的不透明度。
//! 使用预计算的 CROSS(IBFT, IJ) 截面和 POPUL(II, ID) 模型种群。
//!
//! 本模块将静态原子数据（BF 跃迁列表、FF 离子数据）预计算到 Opacf0State，
//! 然后在 RESOLV 中使用 model.levpop.popul 计算不透明度。

use crate::tlusty::math::hydrogen::hephot;
use crate::tlusty::math::atomic::gfree1;
use crate::tlusty::math::special::gaunt;
use crate::tlusty::math::hydrogen::sgmer1;
use crate::tlusty::state::constants::{HK, SIGE, H, UN};
use crate::tlusty::state::model::{GffPar, ModelState};
use crate::tlusty::state::atomic::AtomicData;

// ============================================================================
// 常量 (from Fortran OPACF0)
// ============================================================================

/// SIGK 常数: SIH0 = 2.815e29 (氢原子 bf 截面基准, cm² × Hz³)
const SIH0: f64 = 2.815e29;

/// FF 基准常数: SGFF0 = 3.694e8
const SGFF0: f64 = 3.694e8;

/// C14 = 2.99793e14 (用于 Gaunt 因子计算)
const C14: f64 = 2.99793e14;

// ============================================================================
// BF 跃迁数据
// ============================================================================

/// 束缚-自由跃迁信息。
#[derive(Debug, Clone)]
pub struct BfTransition {
    /// 下能级索引 (0-based)
    pub ilow: usize,
    /// 上能级 (连续区) 索引 (0-based)
    pub iup: usize,
    /// 离子索引 (0-based)
    pub ion_idx: usize,
    /// 有效主量子数
    pub nquant: f64,
    /// Z² (原子序数的平方)
    pub iz2: f64,
    /// 电离频率 (Hz), ν_n = ENION(n)/H
    pub nu_edge: f64,
    /// He I 量子数 (s, l, n) — 使用 HEPHOT OP 截面而非氢原子 SIGK
    pub he1_quant: Option<(i32, i32, i32)>,
}

// ============================================================================
// FF 离子数据
// ============================================================================

/// 自由-自由离子信息。
#[derive(Debug, Clone)]
pub struct FfIon {
    /// 连续区能级索引 NNEXT (0-based)
    pub nnext: usize,
    /// Z² (CHARG2)
    pub charg2: f64,
    /// FF 模式: 1=hydrogenic Gaunt=1, 2=exact Gaunt (GFREE1)
    pub itra: i32,
    /// FF 边缘频率 (Hz), 用于 stimulated emission switch
    pub ff_edge: f64,
}

// ============================================================================
// Opacf0State
// ============================================================================

/// OPACF0 预计算状态。
///
/// 存储从原子数据预计算的 BF 跃迁和 FF 离子信息，
/// 在 RESOLV 中与 model.levpop.popul 一起使用计算不透明度。
#[derive(Debug, Clone)]
pub struct Opacf0State {
    /// BF 跃迁列表 (对应 Fortran NTRANC 个 ITRBF 跃迁)
    pub bf_transitions: Vec<BfTransition>,
    /// FF 离子列表 (对应 Fortran NION 个离子)
    pub ff_ions: Vec<FfIon>,
    /// 是否使用 HEPHOT OP 截面 (替代氢原子 SIGK) 计算 He I bf
    /// 注意: 启用后需要完整 OPACF0+SOLVES 管线才能收敛
    pub use_hephot: bool,
}

impl Opacf0State {
    /// 创建 HHe 模型的 Opacf0State。
    ///
    /// HHe 模型结构 (Rust 0-based):
    /// - Ion 0 (H I):   levels 0-8 (n=1..9), continuum=9
    /// - Ion 1 (He I):  levels 10-23 (14 levels), continuum=24
    /// - Ion 2 (He II): levels 24-37 (n=1..14), continuum=38
    /// - NLEVEL=39, NION=3
    pub fn new_hhe() -> Self {
        let h_planck: f64 = HK * 1.3806e-16; // H = HK × k = 6.6256e-27 erg·s
        let eh: f64 = 2.17853041e-11; // H 电离能 (erg)

        let mut bf_transitions = Vec::new();

        // === H I bf 跃迁 (n=1..9 → continuum level 9) ===
        // Hydrogenic: ν_n = EH/(n² × H)
        for n in 1..=9_usize {
            let nn = n as f64;
            let nu_n = eh / (nn * nn * h_planck); // ionization frequency of level n
            bf_transitions.push(BfTransition {
                ilow: n - 1,
                iup: 9,
                ion_idx: 0,
                nquant: nn,
                iz2: 1.0, // Z=1, Z²=1
                nu_edge: nu_n,
                he1_quant: None, // H I: use hydrogenic SIGK
            });
        }

        // === He I bf 跃迁 (14 levels → continuum level 24) ===
        // 使用实际频率数据 (from he1.dat) + HEPHOT OP 截面
        // 量子数 (s, l, n) from he1.dat level definitions
        let he1_freq: [f64; 14] = [
            5.94503520e15, 1.15267210e15, 9.60145430e14, 8.75933720e14,
            8.14536220e14, 4.51727350e14, 4.02921120e14, 3.81935640e14,
            3.65836790e14, 3.65746870e14, 3.62599020e14, 2.40043860e14,
            2.20797190e14, 2.12492940e14,
        ];
        let he1_nq: [f64; 14] = [
            1.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 4.0, 4.0, 4.0,
        ];
        // He I quantum numbers (s, l, n) from he1.dat / sbfhe1.f
        let he1_sln: [(i32, i32, i32); 14] = [
            (1, 0, 1), // 1 sing S   (G=1)
            (3, 0, 2), // 2 trip S   (G=3)
            (1, 0, 2), // 2 sing S   (G=1)
            (3, 1, 2), // 2 trip P   (G=9)
            (1, 1, 2), // 2 sing P   (G=3)
            (3, 0, 3), // 3 trip S   (G=3)
            (1, 0, 3), // 3 sing S   (G=1)
            (3, 1, 3), // 3 trip P   (G=9)
            (3, 2, 3), // 3 trip D   (G=15)
            (1, 2, 3), // 3 sing D   (G=5)
            (1, 1, 3), // 3 sing P   (G=3)
            (3, 0, 4), // 4 trip S   (G=3)
            (1, 0, 4), // 4 sing S   (G=1)
            (3, 1, 4), // 4 trip P   (G=9)
        ];
        for ilev in 0..14 {
            bf_transitions.push(BfTransition {
                ilow: 10 + ilev,
                iup: 24,
                ion_idx: 1,
                nquant: he1_nq[ilev],
                iz2: 4.0,
                nu_edge: he1_freq[ilev],
                he1_quant: Some(he1_sln[ilev]),
            });
        }

        // === He II bf 跃迁 (n=1..14 → continuum level 38) ===
        // Hydrogenic Z=2: ν_n = 4×EH/(n² × H)
        let he2_freq: [f64; 14] = [
            1.31575980e16, 3.28939940e15, 1.46195530e15, 8.22349860e14,
            5.26303910e14, 3.65488820e14, 2.68522400e14, 2.05587460e14,
            1.62439480e14, 1.31575980e14, 1.08740480e14, 9.13722060e13,
            7.78556080e13, 6.71306000e13,
        ];
        for n in 1..=14_usize {
            let nn = n as f64;
            bf_transitions.push(BfTransition {
                ilow: 24 + n - 1,
                iup: 38,
                ion_idx: 2,
                nquant: nn,
                iz2: 4.0, // Z=2, Z²=4
                nu_edge: he2_freq[n - 1],
                he1_quant: None, // He II: hydrogenic, use SIGK
            });
        }

        // === FF 离子 ===
        let ff_ions = vec![
            FfIon {
                nnext: 9,   // H II (continuum)
                charg2: 1.0, // Z²=1 for H⁺
                itra: 2,     // exact Gaunt (GFREE1)
                ff_edge: 0.0, // H: FF(ION)=0 → no edge
            },
            FfIon {
                nnext: 24,  // He II (continuum)
                charg2: 1.0, // Z²=1 for He⁺
                itra: 2,
                ff_edge: 0.0,
            },
            FfIon {
                nnext: 38,  // He III (continuum)
                charg2: 4.0, // Z²=4 for He²⁺
                itra: 2,
                ff_edge: 0.0,
            },
        ];

        Self {
            bf_transitions,
            ff_ions,
            use_hephot: true, // He I: use Opacity Project cross-sections (HEPHOT)
        }
    }

    /// 计算给定 (depth, frequency) 的连续谱不透明度。
    ///
    /// 对应 Fortran OPACF0 的频率循环体 (行 168-346)。
    /// 使用模型种群 POPUL(level, depth) 而非重新计算 Saha 平衡。
    ///
    /// # 参数
    ///
    /// * `fr` - 频率 (Hz)
    /// * `t` - 温度 (K)
    /// * `ne` - 电子密度 (cm⁻³)
    /// * `popul` - 种群数组 popul[level][depth] (cm⁻³)
    /// * `id` - 深度索引 (0-based)
    /// * `gffpar` - Gaunt 因子预计算参数 (由 gfree0 填充)
    ///
    /// # 返回
    ///
    /// Returns (true_abs, scat, emis_pre_stim, rayleigh) where:
    /// - true_abs: true absorption after stimulated emission correction (cm⁻¹)
    /// - scat: total scattering = elscat + rayleigh (cm⁻¹)
    /// - emis_pre_stim: accumulated emission before stimulated emission correction (cm⁻¹)
    ///   Needed for OPACFL finalization: EMIS1L = emis_pre_stim * XKF * Bν
    /// - rayleigh: Rayleigh scattering only (cm⁻¹), needed for SCAT1 in Lucy
    pub fn compute_opacity(
        &self,
        fr: f64,
        t: f64,
        ne: f64,
        popul: &[Vec<f64>],
        id: usize,
        gffpar: &GffPar,
    ) -> (f64, f64, f64, f64) {
        let hkt = HK / t;
        let sqrt_t = t.sqrt();
        let _tk1 = 1.0 / (H * t); // 1/(kT) in erg⁻¹

        // 电子散射 (对应 Fortran: ELSCAT = ELEC(ID) × SIGE)
        let elscat = SIGE * ne;

        // 预计算频率相关量
        let fr_inv = 1.0 / fr;
        let fr3_inv = fr_inv * fr_inv * fr_inv;

        // 受激发射因子 (XKF = exp(-hν/kT))
        let xkf = if hkt * fr < 150.0 {
            (-hkt * fr).exp()
        } else {
            0.0
        };

        let mut abso = elscat;
        let mut emis = 0.0;

        // ================================================================
        // 1. 束缚-自由贡献
        // ================================================================
        // 对应 Fortran: DO 30 IBFT=1,NTRANC → ABSO += CROSS × ABTRA; EMIS += CROSS × EMTRA
        for bt in &self.bf_transitions {
            if fr < bt.nu_edge {
                continue;
            }

            let pop_low = popul[bt.ilow][id];
            if pop_low <= 0.0 {
                continue;
            }

            // 截面计算:
            // He I: HEPHOT (OP cross-sections) if enabled
            // H I / He II: hydrogenic SIGK with exact Gaunt factor (IBF=1)
            //   SIGK = SIH0 * Z^4 / (nu^3 * n^5) * GAUNT(n, nu/Z^2)
            let sigma = if self.use_hephot {
                if let Some((s, l, n)) = bt.he1_quant {
                    hephot(s, l, n, fr)
                } else {
                    let sigma_base = SIH0 * bt.iz2 * bt.iz2 * fr3_inv / bt.nquant.powi(5);
                    let g_bf = gaunt(bt.nquant as usize, fr / bt.iz2);
                    sigma_base * g_bf
                }
            } else {
                let sigma_base = SIH0 * bt.iz2 * bt.iz2 * fr3_inv / bt.nquant.powi(5);
                let g_bf = gaunt(bt.nquant as usize, fr / bt.iz2);
                sigma_base * g_bf
            };
            if sigma <= 0.0 {
                continue;
            }

            // ABTRA = POPUL(II, ID) — 下能级（束缚态）种群
            abso += sigma * pop_low;

            // EMTRA — LTE 中 Kirchhoff 定律: EMTRA = ABTRA = pop_low × σ
            // 对应 Fortran OPACF0 行 70-71:
            //   ABTRA(ITR,ID) = POPUL(II,ID)
            //   EMTRA(ITR,ID) = POPUL(JJ,ID)*ANE*SBF(II)*WOP(II,ID)*CORR
            // LTE 中 (Saha-Boltzmann 平衡, WOP=1, CORR=1): EMTRA = POPUL(II) = pop_low
            emis += sigma * pop_low;
        }

        // ================================================================
        // 2. 自由-自由贡献
        // ================================================================
        // 对应 Fortran: DO 40 ION=1,NION
        let sgff = SGFF0 / sqrt_t * ne;

        for ion in &self.ff_ions {
            let pop_cont = popul[ion.nnext][id];
            if pop_cont <= 0.0 {
                continue;
            }

            // SF1 = SFF3(ION,ID) × FR3INV
            let sf1 = sgff * ion.charg2 * pop_cont * fr3_inv;

            // SF2 = SFF2(ION,ID) = exp(FF(ION) × HKT1(ID))
            // FF(ION) = 0 for our model → SF2 = 1
            let sf2 = 1.0_f64; // FF=0 → exp(0) = 1

            let mut absoff = sf1 * sf2;

            // Exact Gaunt factor (ITRA=2)
            // Fortran OPACF0 line ~210: SFF3(ION,ID) = GFREE1(ID,X) * SFF2(ION,ID)
            // where X = C14 * CHARG2(ION) / FR
            if ion.itra == 2 {
                let x = C14 * ion.charg2 / fr;
                let gf1 = gfree1(id, x, gffpar);
                // SFF3 = GF1 * SFF2 → absoff = SF1 * GF1 (replacing the default G=1)
                absoff = sf1 * gf1;
            }

            abso += absoff;
            emis += absoff; // ff: Kirchhoff (emission = absorption for thermal ff)
        }

        // ================================================================
        // 受激发射修正
        // ================================================================
        // 对应 Fortran: ABSO(IJ) = ABSO(IJ) - EMIS(IJ) × XKF(ID)
        abso -= emis * xkf;

        // ================================================================
        // 3. OPADD0 — 额外不透明度来源 (Rayleigh 散射)
        // ================================================================
        // 对应 Fortran OPADD0: HI Rayleigh, HeI Rayleigh
        let mut rayleigh = 0.0;

        // HI Rayleigh scattering (IRSCT=1 in Fortran)
        // σ_Ray(HI) = (CR0 + (CR1 + CR2/X)/X) / X², X = (c/λ)²
        // where λ = c/ν, so X = (c/ν)² = c²/ν²
        // clamped at FRRAY = 2.463e15 Hz
        {
            const CLS: f64 = 2.997925e18;
            const CR0: f64 = 5.799e-13;
            const CR1: f64 = 1.422e-6;
            const CR2: f64 = 2.784;
            const FRRAY: f64 = 2.463e15;
            let frm = fr.min(FRRAY);
            let x = (CLS / frm) * (CLS / frm); // X = (c/λ)²
            let cs_ray_hi = (CR0 + (CR1 + CR2 / x) / x) / (x * x);
            // POPUL(0, id) = H I ground state population
            rayleigh += popul[0][id] * cs_ray_hi;
        }

        // HeI Rayleigh scattering (IRSCHE=1 in Fortran)
        // σ_Ray(HeI) = 5.484e-14 / X² * (1 + (2.44e5 + 5.94e10/(X-2.90e5))/X)²
        // clamped at FRAYHe = 5.150e15 Hz
        {
            const CLS: f64 = 2.997925e18;
            const FRAYHE: f64 = 5.150e15;
            let frm = fr.min(FRAYHE);
            let x = (CLS / frm) * (CLS / frm);
            let cs_ray_he1 =
                5.484e-14 / (x * x) * (1.0 + (2.44e5 + 5.94e10 / (x - 2.90e5)) / x).powi(2);
            // He I ground state is level 10 (index 10)
            rayleigh += popul[10][id] * cs_ray_he1;
        }

        // 分离散射和真吸收
        // abso = elscat + true_absorption (after stimulated emission)
        // true_abs = abso - elscat = emis × (1 - XKF)
        let true_abs = (abso - elscat).max(0.0);
        let scat = elscat + rayleigh;

        (true_abs, scat, emis, rayleigh)
    }

    /// 使用真实模型数据（BFCS 和 SGMER）计算连续谱不透明度，完全匹配 Fortran OPACF0
    pub fn compute_opacity_from_model(
        fr: f64,
        ij: usize,
        t: f64,
        ne: f64,
        popul: &[Vec<f64>],
        id: usize,
        gffpar: &GffPar,
        model: &ModelState,
        atomic: &AtomicData,
    ) -> (f64, f64, f64, f64) {
        let hkt = HK / t;
        let sqrt_t = t.sqrt();
        let tk1 = 1.0 / (H * t);

        let elscat = SIGE * ne;

        let fr_inv = 1.0 / fr;
        let fr3_inv = fr_inv * fr_inv * fr_inv;

        let xkf = if hkt * fr < 150.0 {
            (-hkt * fr).exp()
        } else {
            0.0
        };

        let mut abso = elscat;
        let mut emis = 0.0;

        // ================================================================
        // 1. 束缚-自由贡献
        // ================================================================
        let ntranc = model.obfpar.itrbf.iter().filter(|&&x| x > 0).count();
        for ibft in 0..ntranc {
            let itr = model.obfpar.itrbf[ibft] as usize - 1;
            if atomic.trapar.indexp[itr] == 0 { continue; }
            
            let ii = atomic.trapar.ilow[itr] as usize - 1;
            let jj = atomic.trapar.iup[itr] as usize - 1;
            let it = atomic.trapar.itra[ii][jj] as usize;
            if it == 0 { continue; }

            let pop_low = popul[ii][id];
            
            // 交叉截面 CROSS
            let ij0 = model.frqall.ijbf[ij] as usize;
            let a1 = model.phoexp.aijbf[ij];
            let sig0 = model.phoexp.bfcs[ibft][ij0] as f64;
            let sig1 = model.phoexp.bfcs[ibft][ij0 + 1] as f64;
            let mut sigma = a1 * sig0 + (UN - a1) * sig1;

            // SGMER 修正 (合并能级)
            let imer = model.mrgpar.imrg[ii];
            if imer > 0 {
                sigma = sgmer1(fr_inv, fr3_inv, imer, id + 1, &model.mrgpar);
            }

            if sigma <= 0.0 { continue; }

            // CORR 计算
            let ie = atomic.levpar.iel[ii] as usize - 1;
            let nke = atomic.ionpar.nnext[ie] as usize - 1;
            let corr = if nke != jj {
                let g_ratio = atomic.levpar.g[nke] / atomic.levpar.g[jj];
                let delta_e = atomic.levpar.enion[nke] - atomic.levpar.enion[jj];
                g_ratio * (delta_e * tk1).exp()
            } else {
                UN
            };

            let pop_up = popul[jj][id];
            let wop_ii = model.wmcomp.wop[ii][id];
            let emis_val = pop_up * ne * model.levpop.sbf[ii] * wop_ii * corr;

            abso += sigma * pop_low;
            emis += sigma * emis_val;
        }

        // ================================================================
        // 2. 自由-自由贡献
        // ================================================================
        let sgff = SGFF0 / sqrt_t * ne;

        let nion = atomic.ionpar.iz.iter().filter(|&&x| x > 0).count();
        for ion in 0..nion {
            let nnext = atomic.ionpar.nnext[ion] as usize - 1;
            let pop_cont = popul[nnext][id];
            if pop_cont <= 0.0 { continue; }

            let charg2 = atomic.ionpar.charg2[ion];
            let sf1 = sgff * charg2 * pop_cont * fr3_inv;
            
            // sf2 = exp(ff * hkt1)
            let ff_val = atomic.ionpar.ff[ion];
            let mut sf2 = (ff_val * hkt).exp();
            if fr < ff_val {
                sf2 = UN / xkf;
            }

            let mut absoff = sf1 * sf2;

            // Gaunt factor
            let itra = atomic.trapar.itra[nnext][nnext];
            if itra == 2 {
                let x = C14 * charg2 / fr;
                let gf1 = gfree1(id, x, gffpar);
                sf2 = sf2 - UN + gf1;
                absoff = sf1 * sf2;
            }

            abso += absoff;
            emis += absoff;
        }

        abso -= emis * xkf;

        // ================================================================
        // 3. OPADD0 — 额外不透明度来源 (Rayleigh 散射)
        // ================================================================
        let mut rayleigh = 0.0;

        // HI Rayleigh scattering
        if nion > 0 {
            const CLS: f64 = 2.997925e18;
            const CR0: f64 = 5.799e-13;
            const CR1: f64 = 1.422e-6;
            const CR2: f64 = 2.784;
            const FRRAY: f64 = 2.463e15;
            let frm = fr.min(FRRAY);
            let x = (CLS / frm) * (CLS / frm);
            let cs_ray_hi = (CR0 + (CR1 + CR2 / x) / x) / (x * x);
            
            // 假设 H I 基态为能级 0
            if atomic.levpar.iel[0] == 1 {
                rayleigh += popul[0][id] * cs_ray_hi;
            }
        }

        // HeI Rayleigh scattering
        if nion > 1 {
            const CLS: f64 = 2.997925e18;
            const FRAYHE: f64 = 5.150e15;
            let frm = fr.min(FRAYHE);
            let x = (CLS / frm) * (CLS / frm);
            let cs_ray_he1 =
                5.484e-14 / (x * x) * (1.0 + (2.44e5 + 5.94e10 / (x - 2.90e5)) / x).powi(2);
            
            // 寻找 He I 基态
            let nlevel = atomic.levpar.iel.iter().filter(|&&x| x > 0).count();
            for i in 0..nlevel {
                if atomic.levpar.iel[i] == 2 && atomic.ionpar.iz[1] == 1 { // He I
                    rayleigh += popul[i][id] * cs_ray_he1;
                    break;
                }
            }
        }

        let true_abs = (abso - elscat).max(0.0);
        let scat = elscat + rayleigh;

        (true_abs, scat, emis, rayleigh)
    }


    /// 计算给定 (frequency, temperature, populations) 的连续谱不透明度。
    ///
    /// 单深度版本，直接接收 `popul[level]` 扁平数组（无需 2D 数组 + depth 索引）。
    /// 用于有限差分温度导数计算（需要扰动种群但不修改全局模型）。
    pub fn compute_opacity_1d(
        &self,
        fr: f64,
        t: f64,
        ne: f64,
        popul: &[f64],
        gffpar: &GffPar,
    ) -> (f64, f64, f64, f64) {
        let hkt = HK / t;
        let sqrt_t = t.sqrt();

        let elscat = SIGE * ne;

        let fr_inv = 1.0 / fr;
        let fr3_inv = fr_inv * fr_inv * fr_inv;

        let xkf = if hkt * fr < 150.0 {
            (-hkt * fr).exp()
        } else {
            0.0
        };

        let mut abso = elscat;
        let mut emis = 0.0;

        // BF contributions
        for bt in &self.bf_transitions {
            if fr < bt.nu_edge { continue; }
            let pop_low = popul[bt.ilow];
            if pop_low <= 0.0 { continue; }

            let sigma = if self.use_hephot {
                if let Some((s, l, n)) = bt.he1_quant {
                    hephot(s, l, n, fr)
                } else {
                    let sigma_base = SIH0 * bt.iz2 * bt.iz2 * fr3_inv / bt.nquant.powi(5);
                    let g_bf = gaunt(bt.nquant as usize, fr / bt.iz2);
                    sigma_base * g_bf
                }
            } else {
                let sigma_base = SIH0 * bt.iz2 * bt.iz2 * fr3_inv / bt.nquant.powi(5);
                let g_bf = gaunt(bt.nquant as usize, fr / bt.iz2);
                sigma_base * g_bf
            };
            if sigma <= 0.0 { continue; }

            abso += sigma * pop_low;
            emis += sigma * pop_low;
        }

        // FF contributions
        let sgff = SGFF0 / sqrt_t * ne;
        for ion in &self.ff_ions {
            let pop_cont = popul[ion.nnext];
            if pop_cont <= 0.0 { continue; }

            let sf1 = sgff * ion.charg2 * pop_cont * fr3_inv;
            let mut absoff = sf1; // SF2=1 (FF=0)
            if ion.itra == 2 {
                let x = C14 * ion.charg2 / fr;
                let gf1 = gfree1(0, x, gffpar);
                absoff = sf1 * gf1;
            }
            abso += absoff;
            emis += absoff;
        }

        abso -= emis * xkf;
        let true_abs = (abso - elscat).max(0.0);

        // Rayleigh scattering (1D version — same formula as compute_opacity)
        let mut rayleigh = 0.0;
        {
            const CLS: f64 = 2.997925e18;
            const CR0: f64 = 5.799e-13;
            const CR1: f64 = 1.422e-6;
            const CR2: f64 = 2.784;
            const FRRAY: f64 = 2.463e15;
            let frm = fr.min(FRRAY);
            let x = (CLS / frm) * (CLS / frm);
            let cs_ray_hi = (CR0 + (CR1 + CR2 / x) / x) / (x * x);
            rayleigh += popul[0] * cs_ray_hi;
        }
        {
            const CLS: f64 = 2.997925e18;
            const FRAYHE: f64 = 5.150e15;
            let frm = fr.min(FRAYHE);
            let x = (CLS / frm) * (CLS / frm);
            let cs_ray_he1 =
                5.484e-14 / (x * x) * (1.0 + (2.44e5 + 5.94e10 / (x - 2.90e5)) / x).powi(2);
            rayleigh += popul[10] * cs_ray_he1;
        }

        (true_abs, elscat + rayleigh, emis, rayleigh)
    }
}

impl Default for Opacf0State {
    fn default() -> Self {
        Self::new_hhe()
    }
}

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

#[cfg(test)]
mod tests {
    use super::*;
    use crate::tlusty::state::constants::MDEPTH;
    use crate::tlusty::math::atomic::gfree0;

    #[test]
    fn test_opacf0_state_hhe() {
        let state = Opacf0State::new_hhe();
        // H I: 9 bf, He I: 14 bf, He II: 14 bf → total 37
        assert_eq!(state.bf_transitions.len(), 37);
        assert_eq!(state.ff_ions.len(), 3);

        // Check first H I bf transition
        assert_eq!(state.bf_transitions[0].ilow, 0); // H I n=1
        assert_eq!(state.bf_transitions[0].iup, 9);  // H II
        assert_eq!(state.bf_transitions[0].iz2, 1.0);

        // Check first He I bf transition
        assert_eq!(state.bf_transitions[9].ilow, 10); // He I level 0
        assert_eq!(state.bf_transitions[9].iup, 24);  // He II

        // Check first He II bf transition
        assert_eq!(state.bf_transitions[23].ilow, 24); // He II n=1
        assert_eq!(state.bf_transitions[23].iup, 38);  // He III
        assert_eq!(state.bf_transitions[23].iz2, 4.0);
    }

    #[test]
    fn test_compute_opacity_basic() {
        let state = Opacf0State::new_hhe();
        let nlevel = 39;
        let nd = 1;

        // Create simple populations
        let mut popul = vec![vec![0.0_f64; nd]; nlevel];
        popul[0][0] = 1e12;  // H I n=1
        popul[9][0] = 1e14;  // H II
        popul[10][0] = 1e10; // He I n=1
        popul[24][0] = 1e10; // He II n=1
        popul[38][0] = 1e12; // He III

        let fr = 3.28805e15; // H Lyman edge
        let t = 30000.0;
        let ne = 1e14;

        // Precompute Gaunt factor coefficients
        let mut gffpar = GffPar::new();
        let temp = [30000.0_f64; MDEPTH];
        gfree0(0, &temp, &mut gffpar);

        let (ab, sct, _emis, ray) = state.compute_opacity(fr, t, ne, &popul, 0, &gffpar);

        assert!(ab > 0.0, "absorption should be positive");
        assert!(sct > 0.0, "scattering should be positive");
        assert!(ab > sct, "at H Lyman edge, absorption should dominate scattering");
        assert!(ray > 0.0, "Rayleigh scattering should be positive");
    }

    #[test]
    fn test_compute_opacity_consistency() {
        // Test that opacity is consistent across nearby frequencies
        let state = Opacf0State::new_hhe();
        let nlevel = 39;
        let nd = 1;
        let mut popul = vec![vec![0.0_f64; nd]; nlevel];
        popul[0][0] = 1e12;
        popul[9][0] = 1e14;

        let t = 30000.0;
        let ne = 1e14;

        let mut gffpar = GffPar::new();
        let temp = [30000.0_f64; MDEPTH];
        gfree0(0, &temp, &mut gffpar);

        // Just above H Lyman edge
        let fr1 = 3.3e15;
        // Just below H Lyman edge
        let fr2 = 3.2e15;

        let (ab1, _, _, _) = state.compute_opacity(fr1, t, ne, &popul, 0, &gffpar);
        let (ab2, _, _, _) = state.compute_opacity(fr2, t, ne, &popul, 0, &gffpar);

        // Above edge: bf contributes → higher opacity
        // Below edge: no H I bf → lower opacity
        assert!(ab1 > ab2, "opacity above Lyman edge should be higher");
    }
}