SpectraRust/src/io/settrm.rs

//! 热力学表读取和初始化。
//!
//! 重构自 TLUSTY `SETTRM` 子程序。
//!
//! # 功能
//!
//! - 读取熵表 (stab.dat) 和压力表 (ptab.dat)
//! - 计算边缘热力学量

use super::{FortranReader, IoError, Result};
use crate::math::{prsent, PrsentParams, ThermTables};
use std::path::Path;

// 气体常数 (erg/mol/K)
const RCON: f64 = 8.31434e7;

/// 读取热力学表并填充 ThermTables 结构体。
///
/// # 参数
/// * `tables` - ThermTables 结构体可变引用
/// * `stab_path` - 熵表文件路径 (stab.dat)
/// * `ptab_path` - 压力表文件路径 (ptab.dat)
///
/// # 返回值
/// 成功返回 Ok(())，失败返回错误。
pub fn settrm<P: AsRef<Path>>(
    tables: &mut ThermTables,
    stab_path: P,
    ptab_path: P,
) -> Result<()> {
    // 读取熵表
    read_stab(tables, stab_path)?;

    // 读取压力表
    read_ptab(tables, ptab_path)?;

    // 计算边缘数组
    compute_edge_arrays(tables)?;

    Ok(())
}

/// 读取熵表 (stab.dat)
fn read_stab<P: AsRef<Path>>(tables: &mut ThermTables, path: P) -> Result<()> {
    let mut reader = FortranReader::from_file(path)?;

    // 读取头部参数
    let _yhea: f64 = reader.read_value()?;
    let index: usize = reader.read_value()?;
    let r1: f64 = reader.read_value()?;
    let r2: f64 = reader.read_value()?;
    let t1: f64 = reader.read_value()?;
    let t2: f64 = reader.read_value()?;
    let t12: f64 = reader.read_value()?;
    let t22: f64 = reader.read_value()?;

    tables.index = index;
    tables.r1 = r1;
    tables.r2 = r2;
    tables.t1 = t1;
    tables.t2 = t2;
    tables.t12 = t12;
    tables.t22 = t22;

    // 读取熵表数据 SL(INDEX, 100)
    // Fortran 格式: 10F8.5，每行 10 个值
    // 注意：Fortran 列优先，但这里直接读取为行优先
    for jr in 1..=index {
        for jqs in 1..=10 {
            let jl = 1 + (jqs - 1) * 10;
            let ju = jl + 9;
            for jq in jl..=ju {
                let value: f64 = reader.read_value()?;
                tables.sl[jr - 1][jq - 1] = value;
            }
        }
    }

    Ok(())
}

/// 读取压力表 (ptab.dat)
fn read_ptab<P: AsRef<Path>>(tables: &mut ThermTables, path: P) -> Result<()> {
    let mut reader = FortranReader::from_file(path)?;

    // 读取头部参数（与熵表相同格式）
    let _yhea: f64 = reader.read_value()?;
    let _index: usize = reader.read_value()?;
    let _r1: f64 = reader.read_value()?;
    let _r2: f64 = reader.read_value()?;
    let _t1: f64 = reader.read_value()?;
    let _t2: f64 = reader.read_value()?;
    let _t12: f64 = reader.read_value()?;
    let _t22: f64 = reader.read_value()?;

    // 读取压力表数据 PL(INDEX, 100)
    for jr in 1..=tables.index {
        for jqp in 1..=10 {
            let jl = 1 + (jqp - 1) * 10;
            let ju = jl + 9;
            for jq in jl..=ju {
                let value: f64 = reader.read_value()?;
                tables.pl[jr - 1][jq - 1] = value;
            }
        }
    }

    Ok(())
}

/// 计算边缘热力学数组
fn compute_edge_arrays(tables: &mut ThermTables) -> Result<()> {
    // 计算参考密度和温度范围
    let r = 1.5 * 10.0_f64.powf(tables.r1);
    let tmin = 1.5 * 10.0_f64.powf(tables.t1);
    let tmax = 0.9 * 10.0_f64.powf(tables.t2);

    tables.redge = r;

    for i in 1..=100 {
        // 计算温度
        let t_log = tables.t1 + (tables.t2 - tables.t1) * (i - 1) as f64 / 99.0;
        let mut t = 10.0_f64.powf(t_log);
        t = tmax.min(tmin.max(t));

        tables.tedge[i - 1] = t;

        let rho = r;

        // 计算数值导数
        let params1 = PrsentParams {
            r: rho * 1.1,
            t,
            tables: tables as &ThermTables,
        };
        let out1 = prsent(&params1);
        let p1 = out1.fp;
        let s1 = out1.fs;

        let params2 = PrsentParams {
            r: rho,
            t: t * 1.1,
            tables: tables as &ThermTables,
        };
        let out2 = prsent(&params2);
        let p2 = out2.fp;
        let s2 = out2.fs;

        let params0 = PrsentParams {
            r: rho,
            t,
            tables: tables as &ThermTables,
        };
        let out0 = prsent(&params0);
        let p0 = out0.fp;
        let s0 = out0.fs;

        // 转换单位
        let s1 = RCON * s1;
        let s2 = RCON * s2;
        let s0 = RCON * s0;

        // 计算导数
        let dpdr = (p1 - p0) / (0.1 * rho);
        let dpdt = (p2 - p0) / (0.1 * t);
        let dsdt = (s2 - s0) / (0.1 * t);
        let dsdr = (s1 - s0) / (0.1 * rho);

        // 计算热力学量
        let den = dpdr * dsdt - dpdt * dsdr;
        let p = p0;
        let s = s0 / RCON;
        let _cv = t * dsdt;
        let cp = t * den / dpdr;
        let dq = dsdt * p / (den * rho);
        let gamma = 1.0 / dq;

        // 存储结果
        tables.pedge[i - 1] = p;
        tables.sedge[i - 1] = s;
        tables.gammaedge[i - 1] = gamma;
    }

    Ok(())
}

/// 仅初始化表为默认值（不读取文件）。
///
/// 当数据文件不可用时使用。
pub fn settrm_default(tables: &mut ThermTables) {
    *tables = ThermTables::default();
}

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

    #[test]
    fn test_settrm_default() {
        let mut tables = ThermTables::default();
        settrm_default(&mut tables);

        assert_eq!(tables.index, 330);
        assert_eq!(tables.r1, -10.0);
        assert_eq!(tables.r2, 0.0);
    }

    #[test]
    fn test_settrm_file_not_found() {
        let mut tables = ThermTables::default();
        let result = settrm(&mut tables, "/nonexistent/stab.dat", "/nonexistent/ptab.dat");
        assert!(result.is_err());
    }

    #[test]
    fn test_compute_edge_arrays() {
        let mut tables = ThermTables::default();

        // 填充一些合理的测试数据
        for i in 0..330 {
            for j in 0..100 {
                tables.sl[i][j] = 10.0 + 0.01 * i as f64 + 0.1 * j as f64;
                tables.pl[i][j] = 5.0 + 0.01 * i as f64 + 0.05 * j as f64;
            }
        }

        // 计算边缘数组
        let result = compute_edge_arrays(&mut tables);
        assert!(result.is_ok());

        // 检查边缘值是否有限（测试数据可能不真实，所以只检查有限性）
        for i in 0..100 {
            assert!(tables.tedge[i] > 0.0, "tedge[{}] = {} should be positive", i, tables.tedge[i]);
            assert!(tables.pedge[i].is_finite(), "pedge[{}] = {} should be finite", i, tables.pedge[i]);
            assert!(tables.sedge[i].is_finite(), "sedge[{}] = {} should be finite", i, tables.sedge[i]);
            assert!(tables.gammaedge[i].is_finite(), "gammaedge[{}] = {} should be finite", i, tables.gammaedge[i]);
        }
    }
}