SpectraRust/src/math/matcon.rs

//! 对流对线性化矩阵的贡献。
//!
//! 重构自 TLUSTY `matcon.f`。
//!
//! 功能：
//! - 计算对流对矩阵 A 和 B 的贡献
//! - 修改能量平衡行 (NRE) 和新的 DELTA 行 (NDEL)
//! - 参考 Grenfell, Astr.Ap. 20, 293 (1972)

use crate::math::convec::{convec, ConvecConfig, ConvecOutput, ConvecParams};
use crate::state::constants::{BOLK, HALF, UN};

/// MATCON 配置参数
#[derive(Debug, Clone)]
pub struct MatconConfig {
    /// 混合长度参数 (HMIX0)
    pub hmix0: f64,
    /// 压力模式标志 (IPRESS)
    pub ipress: i32,
    /// 对数导数标志 (ILGDER)
    pub ilgder: i32,
    /// 对流模式标志 (ICONV)
    pub iconv: i32,
    /// 盘模式标志 (IDISK)
    pub idisk: i32,
    /// DELTA 方程标志 (INDL)
    pub indl: i32,
    /// 氦方程标志 (INHE)
    pub inhe: i32,
    /// 能量方程行号偏移 (INRE)
    pub inre: i32,
    /// 压力方程行号偏移 (INPC)
    pub inpc: i32,
    /// 频率数 (NFREQE)
    pub nfreqe: usize,
}

impl Default for MatconConfig {
    fn default() -> Self {
        Self {
            hmix0: 1.0,
            ipress: 0,
            ilgder: 0,
            iconv: 0,
            idisk: 0,
            indl: 0,
            inhe: 0,
            inre: 0,
            inpc: 0,
            nfreqe: 0,
        }
    }
}

/// MATCON 输入参数
pub struct MatconParams<'a> {
    /// 深度点索引 (1-indexed)
    pub id: usize,
    /// 总深度点数
    pub nd: usize,
    /// 温度数组 [K]
    pub temp: &'a [f64],
    /// 总压力数组
    pub ptotal: &'a [f64],
    /// 气体压力数组
    pub pgs: &'a [f64],
    /// 密度数组 [g/cm³]
    pub dens: &'a [f64],
    /// 电子密度数组
    pub elec: &'a [f64],
    /// 平均分子量数组
    pub wmm: &'a [f64],
    /// 湍流速度数组
    pub vturb: &'a [f64],
    /// Rosseland 不透明度数组
    pub abrosd: &'a [f64],
    /// 柱密度数组
    pub dm: &'a [f64],
    /// DELTA 参数数组 (dlnT/dlnP)
    pub delta: &'a mut [f64],
    /// 对流通量数组
    pub flxc: &'a mut [f64],
    /// 微分方程权重数组
    pub redif: &'a [f64],
    /// 积分方程权重数组
    pub reint: &'a [f64],
    /// 几何因子数组 (盘模型用)
    pub zd: &'a [f64],
    /// 重力加速度缩放因子
    pub qgrav: f64,
    /// 配置参数
    pub config: MatconConfig,
    /// CONVEC 配置
    pub convec_config: ConvecConfig,
}

/// MATCON 矩阵元素
pub struct MatconMatrices<'a> {
    /// 矩阵 A (三对角，a[i] = 上一行对角线左边)
    pub a: &'a mut [f64],
    /// 矩阵 B (对角线)
    pub b: &'a mut [f64],
    /// 矩阵 C (下一行对角线右边)
    pub c: &'a mut [f64],
    /// 右端向量
    pub vecl: &'a mut [f64],
}

/// MATCON 输出
#[derive(Debug, Clone)]
pub struct MatconOutput {
    /// 是否执行了计算
    pub computed: bool,
    /// 更新后的 DELTA 值
    pub delta: f64,
    /// 对流通量
    pub flxc: f64,
}

/// 计算对流对矩阵的贡献。
///
/// # 参数
/// * `params` - 输入参数
/// * `matrices` - 矩阵元素（会被修改）
///
/// # 返回值
/// 输出结构体，包含更新后的 DELTA 和对流通量
pub fn matcon(params: &mut MatconParams, matrices: &mut MatconMatrices) -> MatconOutput {
    // 检查是否启用对流
    if params.config.hmix0 <= 0.0 {
        return MatconOutput {
            computed: false,
            delta: 0.0,
            flxc: 0.0,
        };
    }

    let id = params.id;
    let cfg = &params.config;

    // 计算行索引 (0-indexed in Rust)
    let nhe = cfg.nfreqe + cfg.inhe as usize;
    let nre = cfg.nfreqe + cfg.inre as usize;
    let npc = cfg.nfreqe + cfg.inpc as usize;
    let ndel = cfg.nfreqe + cfg.indl as usize;

    // 计算电子相对密度
    let anerel = params.elec[0] / (params.dens[0] / params.wmm[0] + params.elec[0]);

    // 上边界条件 (ID = 1)
    if id == 1 {
        params.delta[0] = 0.0;
        params.flxc[0] = 0.0;
        if cfg.indl > 0 {
            // B(NDEL, NDEL) = 1
            let idx = ndel * ndel;
            if idx < matrices.b.len() {
                matrices.b[idx] = UN;
            }
        }
        return MatconOutput {
            computed: true,
            delta: 0.0,
            flxc: 0.0,
        };
    }

    // 正常深度点 1 < ID < ND
    let t = params.temp[id - 1];
    let p = params.ptotal[id - 1];
    let pg = params.pgs[id - 1];
    let prad = p - pg - HALF * params.dens[id - 1] * params.vturb[id - 1].powi(2);

    let tm = params.temp[id - 2];
    let pm = params.ptotal[id - 2];
    let pgm = params.pgs[id - 2];
    let pradm = pm - pgm - HALF * params.dens[id - 2] * params.vturb[id - 2].powi(2);

    let (t0, p0, pg0, pr0, ab0, dlt, ddt0, ddtm, dlp) = if cfg.ilgder == 0 {
        // 算术平均
        let t0 = HALF * (t + tm);
        let p0 = HALF * (p + pm);
        let pg0 = HALF * (pg + pgm);
        let pr0 = HALF * (prad + pradm);
        let ab0 = HALF * (params.abrosd[id - 1] + params.abrosd[id - 2]);
        let dlt = (t - tm) / (p - pm) * p0 / t0;
        let tt = t * t - tm * tm;
        let ddt0 = dlt / HALF * tm / tt;
        let ddtm = -ddt0 * t / tm;
        (t0, p0, pg0, pr0, ab0, dlt, ddt0, ddtm, 0.0)
    } else {
        // 几何平均
        let t0 = (t * tm).sqrt();
        let p0 = (p * pm).sqrt();
        let pg0 = (pg * pgm).sqrt();
        let pr0 = (prad * pradm).sqrt();
        let ab0 = (params.abrosd[id - 1] * params.abrosd[id - 2]).sqrt();
        let dlp = UN / (p / pm).ln();
        let dlt = (t / tm).ln() * dlp;
        let ddt0 = dlp / t;
        let ddtm = -dlp / tm;
        (t0, p0, pg0, pr0, ab0, dlt, ddt0, ddtm, dlp)
    };

    params.delta[id - 1] = dlt;

    // DELTA 方程的矩阵元素
    if cfg.indl > 0 {
        // B(NDEL, NDEL) = -1
        let idx = ndel * ndel;
        if idx < matrices.b.len() {
            matrices.b[idx] = -UN;
        }
        // VECL(NDEL) = DELTA(ID) - DLT
        if ndel < matrices.vecl.len() {
            matrices.vecl[ndel] = params.delta[id - 1] - dlt;
        }

        // 压力导数项
        let (ddp0, ddpm) = if cfg.ipress > 0 {
            if cfg.ilgder == 0 {
                let pp0 = p * p - pm * pm;
                let ddp0 = -dlt / HALF * pm / pp0;
                let ddpm = -ddp0 * p / pm;
                (ddp0, ddpm)
            } else {
                let pp0 = (t / tm).ln() * dlp * dlp;
                let ddp0 = -pp0 / p;
                let ddpm = pp0 / pm;
                (ddp0, ddpm)
            }
        } else {
            (0.0, 0.0)
        };

        // A 矩阵 DELTA 行
        if cfg.inhe > 0 {
            let idx_a = ndel * nhe;
            if idx_a < matrices.a.len() {
                matrices.a[idx_a] = BOLK * tm * ddpm;
            }
        }
        let idx_a = ndel * nre;
        if idx_a < matrices.a.len() {
            matrices.a[idx_a] = pgm / tm * ddpm + ddtm;
        }

        // B 矩阵 DELTA 行
        if cfg.inhe > 0 {
            let idx_b = ndel * nhe;
            if idx_b < matrices.b.len() {
                matrices.b[idx_b] = BOLK * t * ddp0;
            }
        }
        let idx_b = ndel * nre;
        if idx_b < matrices.b.len() {
            matrices.b[idx_b] = pg / t * ddp0 + ddt0;
        }
    }

    // 计算对流通量及其导数
    let gravd = if cfg.idisk == 1 {
        params.zd[id - 1] * params.qgrav
    } else {
        0.0
    };

    let convec_params = ConvecParams {
        id,
        t: t0,
        ptot: p0,
        pg: pg0,
        prad: pr0,
        abros: ab0,
        delta: dlt,
        taurs: 0.0, // 需要从模型获取
        config: params.convec_config.clone(),
        trmder_config: None,
        therm_tables: None,
    };

    let convec_out = convec(&convec_params);
    let flxcnv = convec_out.flxcnv;
    let vcon = convec_out.vconv;
    params.flxc[id - 1] = flxcnv;

    // 计算对流通量导数（数值微分）
    let delmde = 0.0; // 需要从 CONVEC 输出获取
    let dhcdd = if delmde > 0.0 {
        1.5 / delmde * flxcnv
    } else {
        0.0
    };

    // T 导数
    let t1 = 1.001 * t0;
    let convec_params_t = ConvecParams {
        id,
        t: t1,
        ptot: p0,
        pg: pg0,
        prad: pr0,
        abros: ab0,
        delta: dlt,
        taurs: 0.0,
        config: params.convec_config.clone(),
        trmder_config: None,
        therm_tables: None,
    };
    let convec_out_t = convec(&convec_params_t);
    let flxc1 = convec_out_t.flxcnv;

    let dhcdt0 = (flxc1 - flxcnv) * 1e3 * HALF;
    let mut dhcdt = dhcdt0 / t;
    let mut dhcdtm = dhcdt0 / tm;

    // P 导数
    let mut dhcdp = 0.0;
    if cfg.ipress > 0 {
        let pg1 = 1.001 * pg0;
        let convec_params_p = ConvecParams {
            id,
            t: t0,
            ptot: p0,
            pg: pg1,
            prad: pr0,
            abros: ab0,
            delta: dlt,
            taurs: 0.0,
            config: params.convec_config.clone(),
            trmder_config: None,
            therm_tables: None,
        };
        let convec_out_p = convec(&convec_params_p);
        let flxc2 = convec_out_p.flxcnv;

        dhcdp = (flxc2 - flxcnv) * 1e3 / pg0 * HALF;

        if cfg.ipress > 1 {
            let p1 = 1.001 * p0;
            let convec_params_pt = ConvecParams {
                id,
                t: t0,
                ptot: p1,
                pg: pg0,
                prad: pr0,
                abros: ab0,
                delta: dlt,
                taurs: 0.0,
                config: params.convec_config.clone(),
                trmder_config: None,
                therm_tables: None,
            };
            let convec_out_pt = convec(&convec_params_pt);
            let flxc3 = convec_out_pt.flxcnv;

            let dhcdpt = (flxc3 - flxcnv) * 1e3 / p0 * HALF;
            dhcdp += dhcdpt;
            dhcdt += dhcdpt * 4.0 * pr0 / t0;
        }
    }

    if cfg.indl == 0 {
        dhcdt += dhcdd * ddt0;
        dhcdtm += dhcdd * ddtm;
    }

    // 微分方程形式的矩阵贡献
    let redif_val = params.redif[id - 1];
    if redif_val > 0.0 {
        if cfg.iconv > 0 {
            if cfg.inhe > 0 {
                let idx_a = nre * nhe;
                if idx_a < matrices.a.len() {
                    matrices.a[idx_a] -= dhcdp * BOLK * tm * redif_val;
                }
                let idx_b = nre * nhe;
                if idx_b < matrices.b.len() {
                    matrices.b[idx_b] += dhcdp * BOLK * t * redif_val;
                }
            }
            let idx_a = nre * nre;
            if idx_a < matrices.a.len() {
                matrices.a[idx_a] -= (dhcdp * pgm / tm + dhcdtm) * redif_val;
            }
            let idx_b = nre * nre;
            if idx_b < matrices.b.len() {
                matrices.b[idx_b] += (dhcdp * pg / t + dhcdt) * redif_val;
            }
            if cfg.indl > 0 {
                let idx_b = nre * ndel;
                if idx_b < matrices.b.len() {
                    matrices.b[idx_b] += dhcdd * redif_val;
                }
            }
        }
        if nre < matrices.vecl.len() {
            matrices.vecl[nre] -= flxcnv * redif_val;
        }
    }

    // 积分方程形式 - 简化实现，完整实现需要处理 ID+1 点
    let reint_val = params.reint[id - 1];
    if reint_val > 0.0 && cfg.iconv <= 2 && id < params.nd {
        // 完整实现需要计算与 ID+1 点相关的项
        // 这里简化处理
    }

    MatconOutput {
        computed: true,
        delta: dlt,
        flxc: flxcnv,
    }
}

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

    fn create_test_params<'a>(
        id: usize,
        temp: &'a [f64],
        ptotal: &'a [f64],
        pgs: &'a [f64],
        dens: &'a [f64],
        elec: &'a [f64],
        wmm: &'a [f64],
        vturb: &'a [f64],
        abrosd: &'a [f64],
        dm: &'a [f64],
        delta: &'a mut [f64],
        flxc: &'a mut [f64],
        redif: &'a [f64],
        reint: &'a [f64],
        zd: &'a [f64],
    ) -> MatconParams<'a> {
        let config = MatconConfig {
            hmix0: 1.0,
            ipress: 0,
            ilgder: 0,
            iconv: 1,
            idisk: 0,
            indl: 0,
            inhe: 0,
            inre: 1,
            inpc: 0,
            nfreqe: 10,
        };

        MatconParams {
            id,
            nd: temp.len(),
            temp,
            ptotal,
            pgs,
            dens,
            elec,
            wmm,
            vturb,
            abrosd,
            dm,
            delta,
            flxc,
            redif,
            reint,
            zd,
            qgrav: 1e4,
            config,
            convec_config: ConvecConfig::default(),
        }
    }

    #[test]
    fn test_matcon_disabled() {
        // 对流禁用时 (hmix0 <= 0)
        let temp = vec![10000.0, 9500.0, 9000.0];
        let ptotal = vec![1e5, 2e5, 3e5];
        let pgs = vec![0.9e5, 1.9e5, 2.9e5];
        let dens = vec![1e-7, 2e-7, 3e-7];
        let elec = vec![1e-8, 2e-8, 3e-8];
        let wmm = vec![1.4e-24; 3];
        let vturb = vec![2e5; 3];
        let abrosd = vec![0.4; 3];
        let dm = vec![1e2, 1e2, 1e2];
        let mut delta = vec![0.0; 3];
        let mut flxc = vec![0.0; 3];
        let redif = vec![1.0; 3];
        let reint = vec![0.0; 3];
        let zd = vec![1.0; 3];

        let mut params = create_test_params(
            2, &temp, &ptotal, &pgs, &dens, &elec, &wmm, &vturb,
            &abrosd, &dm, &mut delta, &mut flxc, &redif, &reint, &zd,
        );
        params.config.hmix0 = -1.0; // 禁用对流

        let mut a = vec![0.0; 100];
        let mut b = vec![0.0; 100];
        let mut c = vec![0.0; 100];
        let mut vecl = vec![0.0; 100];

        let mut matrices = MatconMatrices {
            a: &mut a,
            b: &mut b,
            c: &mut c,
            vecl: &mut vecl,
        };

        let result = matcon(&mut params, &mut matrices);

        assert!(!result.computed);
    }

    #[test]
    fn test_matcon_upper_boundary() {
        // 上边界条件 (ID = 1)
        let temp = vec![10000.0, 9500.0, 9000.0];
        let ptotal = vec![1e5, 2e5, 3e5];
        let pgs = vec![0.9e5, 1.9e5, 2.9e5];
        let dens = vec![1e-7, 2e-7, 3e-7];
        let elec = vec![1e-8, 2e-8, 3e-8];
        let wmm = vec![1.4e-24; 3];
        let vturb = vec![2e5; 3];
        let abrosd = vec![0.4; 3];
        let dm = vec![1e2, 1e2, 1e2];
        let mut delta = vec![0.5; 3];
        let mut flxc = vec![0.0; 3];
        let redif = vec![1.0; 3];
        let reint = vec![0.0; 3];
        let zd = vec![1.0; 3];

        let mut params = create_test_params(
            1, &temp, &ptotal, &pgs, &dens, &elec, &wmm, &vturb,
            &abrosd, &dm, &mut delta, &mut flxc, &redif, &reint, &zd,
        );

        let mut a = vec![0.0; 100];
        let mut b = vec![0.0; 100];
        let mut c = vec![0.0; 100];
        let mut vecl = vec![0.0; 100];

        let mut matrices = MatconMatrices {
            a: &mut a,
            b: &mut b,
            c: &mut c,
            vecl: &mut vecl,
        };

        let result = matcon(&mut params, &mut matrices);

        assert!(result.computed);
        assert_eq!(result.delta, 0.0);
        assert_eq!(result.flxc, 0.0);
        assert_eq!(params.delta[0], 0.0);
    }
}