SpectraRust/src/synspec/math/opac.rs
fmq 0dfe6facd6 feat: 完善 SYNSPEC 不透明度调用链 + 新增旋转卷积和状态结构体
- 新增 rotin.rs: ROTINS 旋转卷积函数(含 kernel、interpolate_at 等辅助函数)
- 新增 synspec/state/: COMMON 块翻译(constants, model, params, wind)
- 重构 opac.rs: 连接 LINOP、MOLOP、HYDLIN、HE2LIN、PHTION、PHTX 调用
  - 新增 OpacLinopData、OpacHydlinData、OpacHe2linData 等可选数据结构
  - 支持 IHYL=0(插值模式)和 IHYL>0(详细模式)的氢线处理
  - 支持分子线不透明度(MOLOP)和光致电离(PHTION/PHTX)
- 更新 resolv.rs: 适配新的 OpacParams 签名
- 更新 he2lin.rs: 修复 minor import
- 更新 mod.rs: 导出 rotin 模块

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
2026-06-07 13:22:15 +08:00

637 lines
19 KiB
Rust

//! Opacity and emissivity calculation for SYNSPEC.
//!
//! Translated from SYNSPEC `OPAC` subroutine (synspec54.f:4541).
//!
//! Calculates absorption, emission, and scattering coefficients
//! at a given depth for multiple frequencies.
//!
//! This routine calls downstream functions for line opacity (LINOP),
//! molecular opacity (MOLOP), hydrogen lines (HYDLIN), He II lines (HE2LIN),
//! and photoionization (PHTION, PHTX).
use crate::synspec::math::{
gfree, he2lin, hydlin, linop, lymlin, molop, phtion, phtx,
He2linParams, HydlinParams, LinopParams, LymlinParams,
MolLineData, MolopConfig, MolopFreqParams, MolopModelState, MolopOutput,
PhtionParams, PhtionOutput, PhtxParams, PhtxOutput,
};
use crate::synspec::math::voigtk::MVOI;
use crate::tlusty::math::hydrogen::sffhmi;
/// Physical constants
const UN: f64 = 1.0;
const TEN15: f64 = 1.0e-15;
const CSB: f64 = 2.0706e-16;
const CFF: f64 = 3.694e8;
/// Line opacity data for LINOP call.
///
/// Corresponds to Fortran COMMON/LINDAT/ data needed by OPAC → LINOP.
pub struct OpacLinopData<'a> {
/// Temperature at depth ID
pub temp: f64,
/// Total frequency count (NFREQS)
pub nfreqs: usize,
/// Frequency spacing factor (DFRCON)
pub dfrcon: f64,
/// Planck function at depth ID
pub plan: f64,
/// Stimulated emission correction
pub stim: f64,
/// Number of lines
pub nlin: usize,
/// Line data array
pub lines: &'a [crate::synspec::math::LineData],
/// RRR(ID,ION,IAT) — Saha/Boltzmann factors
pub rrr: &'a [f64],
/// RRR dimensions: (natom, nion)
pub rrr_dims: (usize, usize),
/// Doppler width DOPA1(IAT, ID)
pub dopa1: &'a [f64],
/// DOPA1 number of atoms
pub dopa1_nat: usize,
/// Level statistical weights G(level)
pub g: &'a [f64],
/// Level populations POPUL(level, ID)
pub popul: &'a [f64],
/// POPUL number of levels
pub popul_nlev: usize,
/// NLTE populations PNLT(IAT, ION, ID)
pub pnlt: &'a [f64],
/// PNLT dimensions: (natom, nion)
pub pnlt_dims: (usize, usize),
/// Ionization energies ENEV(IAT, ION)
pub enev: &'a [f64],
/// ENEV number of atoms
pub enev_nat: usize,
/// Level energies ENION(level)
pub enion: &'a [f64],
/// Laser deletion flag
pub lasdel: bool,
/// He II special line count
pub nsp: usize,
/// He II special line indices
pub isp0: &'a [usize],
/// Voigt table H0
pub h0tab: &'a [f64; MVOI],
/// Voigt table H1
pub h1tab: &'a [f64; MVOI],
/// Voigt table H2
pub h2tab: &'a [f64; MVOI],
}
/// Hydrogen line data for HYDLIN call.
pub struct OpacHydlinData {
/// Lower level for H lines
pub ilowh: i32,
/// Upper level limits
pub m10: usize,
pub m20: usize,
/// Turbulent velocity (cm/s)
pub vturb: f64,
/// wnHint factors
pub wn_hint: Vec<Vec<f64>>,
}
/// He II line data for HE2LIN call.
pub struct OpacHe2linData {
/// He II common data
pub common: crate::synspec::math::He2Common<'static>,
/// He II lowest series index
pub ilwhe2: usize,
/// Max principal quantum number
pub mhe10: usize,
/// Upper limit for He II lines
pub mhe20: usize,
}
/// Molecular opacity data for MOLOP call.
pub struct OpacMolopData<'a> {
/// MOLOP configuration
pub config: &'a MolopConfig,
/// Molecular line data
pub line_data: &'a MolLineData<'a>,
/// Voigt tables (H0, H1, H2)
pub voigt_tables: (&'a [f64; MVOI], &'a [f64; MVOI], &'a [f64; MVOI]),
}
/// Photoionization data for PHTION call.
pub struct OpacPhtionData<'a> {
/// Photoionization cross-section data
pub photcs: &'a crate::synspec::math::PhotcsData,
/// Atomic number density RRR
pub rrr: &'a [f64],
/// Level populations POPUL
pub popul: &'a [f64],
}
/// Extended photoionization data for PHTX call.
pub struct OpacPhtxData<'a> {
/// Level photoionization data
pub level_photo: &'a crate::synspec::math::LevelPhotoData,
/// H/He photoionization data
pub hhe_photo: &'a crate::synspec::math::HhePhotoData,
/// Level populations POPUL
pub popul: &'a [f64],
/// Atomic number density RRR
pub rrr: &'a [f64],
/// Total frequency count
pub nfreq_total: usize,
/// Continuous frequency count
pub nfreqc: usize,
/// IASV flag
pub iasv: i32,
/// NQHT flag
pub nqht: i32,
/// PHTX state for interpolation caching
pub state: &'a mut crate::synspec::math::PhtxState,
}
/// Parameters for opacity calculation.
///
/// Includes optional data for downstream functions (LINOP, MOLOP, HYDLIN,
/// HE2LIN, PHTION, PHTX). When `None`, the corresponding call is skipped.
pub struct OpacParams<'a> {
/// Depth index
pub id: usize,
/// Temperature at depth ID (K)
pub t: f64,
/// Electron density at depth ID
pub ane: f64,
/// Number of frequencies
pub nfreq: usize,
/// Frequency array (Hz)
pub freq: Vec<f64>,
/// Wavelength array (Å)
pub wlam: Vec<f64>,
/// Frequency interpolation weights
pub frx1: Vec<f64>,
pub frx2: Vec<f64>,
/// Mode flag
pub imode: i32,
/// Standard depth index
pub idstd: usize,
/// Continuum opacity switch
pub icontl: i32,
/// Lyman line treatment switch
pub iophli: i32,
/// Hydrogen line mode (0=interpolated, >0=detailed)
pub ihyll: i32,
/// He II line switch
pub ihe2l: i32,
/// Molecular line switch
pub ifmol: i32,
/// Number of molecular lists
pub nmlist: usize,
/// H atom index
pub iath: i32,
/// H ground level population
pub pop_h: f64,
/// H continuum level population
pub pop_h_cont: f64,
/// Electron scattering coefficient
pub sce: f64,
/// Planck function at depth ID
pub plan: f64,
/// H/kT
pub hkt: f64,
/// h/k
pub hk: f64,
/// Boltzmann constant * temperature
pub bn: f64,
/// wnHint factors for Lyman lines
pub wn_hint: [f64; 5],
/// RELOP factor for AVAB calculation
pub relop: f64,
// --- Downstream function data (optional) ---
/// Line opacity data for LINOP
pub linop_data: Option<OpacLinopData<'a>>,
/// Hydrogen line data for HYDLIN
pub hydlin_data: Option<OpacHydlinData>,
/// He II line data for HE2LIN
pub he2lin_data: Option<OpacHe2linData>,
/// Molecular opacity data for MOLOP
pub molop_data: Option<OpacMolopData<'a>>,
/// Photoionization data for PHTION
pub phtion_data: Option<OpacPhtionData<'a>>,
/// Extended photoionization data for PHTX
pub phtx_data: Option<OpacPhtxData<'a>>,
}
/// Result of opacity calculation
pub struct OpacResult {
/// Absorption coefficient array
pub abso: Vec<f64>,
/// Emission coefficient array
pub emis: Vec<f64>,
/// Scattering coefficient array
pub scat: Vec<f64>,
/// Average absorption for line selection
pub avab: f64,
}
/// Calculate opacity and emissivity.
///
/// This is the main opacity calculation routine for SYNSPEC.
/// It computes absorption, emission, and scattering coefficients
/// including continuum, line, and molecular contributions.
///
/// # Arguments
/// * `params` - Input parameters (including optional downstream data)
///
/// # Returns
/// Opacity coefficients for all frequencies
pub fn opac(params: &mut OpacParams) -> OpacResult {
let nfreq = params.nfreq;
let mut abso = vec![0.0; nfreq];
let mut emis = vec![0.0; nfreq];
let mut scat = vec![0.0; nfreq];
// Skip if not needed
if params.imode == -1 && params.id != params.idstd {
return OpacResult {
abso,
emis,
scat,
avab: 0.0,
};
}
let t = params.t;
let ane = params.ane;
let t1 = UN / t;
let hkt = params.hkt;
let _tk = hkt / params.hk;
let srt = UN / t.sqrt();
let sgff = CFF * srt;
let con = CSB * t1 * srt;
let conts = 1.0e-36 / con;
// Continuum opacity (first and last frequency)
let ij0 = if nfreq == 1 {
1
} else if params.imode == 2 {
nfreq
} else {
2
};
let mut ably1 = 0.0;
let mut emly1 = 0.0;
let mut scly1 = 0.0;
for ij in 0..ij0 {
let fr = params.freq[ij];
let fr15 = fr * TEN15;
let bnu = params.bn * fr15 * fr15 * fr15;
let hkf = hkt * fr;
// Bound-free and free-free continuum
let mut abf = 0.0;
let mut ebf = 0.0;
let mut aff = 0.0;
// Free-free Gaunt factor contribution
let gfree_val = gfree(t, fr);
let sf1 = sgff / (fr * fr * fr);
let hktm = hkt * fr.min(fr);
let sf2 = hktm.exp() + gfree_val - UN;
let sff = sf1 * sf2;
aff += sff;
// H- opacity
let sffhmi_val = sffhmi(params.pop_h, fr, t);
aff += sffhmi_val;
// Additional opacities (OPADD placeholder)
let abad = 0.0;
let emad = 0.0;
let scad = 0.0;
// Lyman line wings
let (ably_ij, emly_ij, scly_ij) = if params.iophli != 0 {
let lymlin_params = LymlinParams {
id: params.id,
freq: fr,
pop_h: params.pop_h,
pop_excited: [params.pop_h_cont; 4],
t,
ane,
iath: params.iath,
iophli: params.iophli,
wn_hint: params.wn_hint,
};
let result = lymlin(&lymlin_params);
(result.ably, result.emly, result.scly)
} else {
(0.0, 0.0, 0.0)
};
// Total opacity and emissivity
let x = (-hkf).exp();
let x1 = UN - x;
let bne = bnu * x * ane;
abso[ij] = abf + ane * (x1 * aff - x * ebf) + abad;
emis[ij] = bne * (aff + ebf) + emad + emly_ij;
scat[ij] = scad + scly_ij + params.sce;
if ij == 0 {
ably1 = ably_ij;
emly1 = emly_ij;
scly1 = scly_ij;
}
}
let avab = (abso[0] + abso[1] + scat[0] + scat[1]) * 0.5 * params.relop;
if nfreq <= 2 || params.imode == -1 {
return OpacResult {
abso,
emis,
scat,
avab,
};
}
// Determine M (start index for line accumulation)
let m: usize = if params.icontl == 1 { 1 } else { 3 };
if params.imode != 2 {
// Interpolated continuum for all frequencies
for ij in 2..nfreq {
abso[ij] = params.frx1[ij] * abso[1] + params.frx2[ij] * abso[0];
emis[ij] = params.frx1[ij] * emis[1] + params.frx2[ij] * emis[0];
scat[ij] = params.frx1[ij] * scat[1] + params.frx2[ij] * scat[0];
}
// Hydrogen lines for IHYL=0 (interpolated mode)
if params.ihyll == 0 {
if let Some(ref hyd_data) = params.hydlin_data {
let hyd_params = HydlinParams {
id: params.id,
i0: 0,
i1: 1,
nfreq,
wlam: params.wlam.clone(),
freq: params.freq.clone(),
t,
ane,
iath: params.iath,
ilowh: hyd_data.ilowh,
m10: hyd_data.m10,
m20: hyd_data.m20,
pop_h: params.pop_h,
pop_h_cont: params.pop_h_cont,
vturb: hyd_data.vturb,
wn_hint: hyd_data.wn_hint.clone(),
};
let result = hydlin(&hyd_params);
for ij in m..nfreq {
abso[ij] += params.frx1[ij] * result.absoh[1]
+ params.frx2[ij] * result.absoh[0];
emis[ij] += params.frx1[ij] * result.emish[1]
+ params.frx2[ij] * result.emish[0];
}
}
}
// Line opacity (LINOP)
if let Some(ref linop_d) = params.linop_data {
let linop_params = LinopParams {
id: params.id,
temp: linop_d.temp,
nfreq,
freq: &params.freq,
nfreqs: linop_d.nfreqs,
dfrcon: linop_d.dfrcon,
avab,
plan: linop_d.plan,
stim: linop_d.stim,
nlin: linop_d.nlin,
lines: linop_d.lines,
rrr: linop_d.rrr,
rrr_dims: linop_d.rrr_dims,
dopa1: linop_d.dopa1,
dopa1_nat: linop_d.dopa1_nat,
g: linop_d.g,
popul: linop_d.popul,
popul_nlev: linop_d.popul_nlev,
pnlt: linop_d.pnlt,
pnlt_dims: linop_d.pnlt_dims,
enev: linop_d.enev,
enev_nat: linop_d.enev_nat,
enion: linop_d.enion,
lasdel: linop_d.lasdel,
nsp: linop_d.nsp,
isp0: linop_d.isp0,
h0tab: linop_d.h0tab,
h1tab: linop_d.h1tab,
h2tab: linop_d.h2tab,
phe1_data: None,
phe2_common: None,
};
let result = linop(&linop_params);
for ij in 3..nfreq {
abso[ij] += result.ablin[ij];
emis[ij] += result.emlin[ij];
}
}
// Molecular line opacity (MOLOP)
if params.ifmol > 0 {
if let Some(ref mol_d) = params.molop_data {
for ilist in 0..params.nmlist {
let mol_params = MolopFreqParams {
nfreq,
nfreqs: nfreq,
freq: &params.freq,
};
let result = molop(
mol_d.config,
&MolopModelState {
temp: t,
elec: ane,
stim: 1.0,
plan: params.plan,
mol: crate::synspec::math::MolModelState {
rrmol: &[],
dopmol: &[],
},
},
mol_d.line_data,
&mol_params,
avab,
mol_d.voigt_tables,
);
for ij in 3..nfreq {
abso[ij] += result.ablin[ij];
emis[ij] += result.emlin[ij];
}
let _ = ilist; // suppress unused warning
}
}
}
}
// Detailed hydrogen line opacity (IHYL>0 or IMODE=2)
if params.ihyll > 0 || params.imode == 2 {
if let Some(ref hyd_data) = params.hydlin_data {
let hyd_params = HydlinParams {
id: params.id,
i0: m.saturating_sub(1),
i1: nfreq.saturating_sub(1),
nfreq,
wlam: params.wlam.clone(),
freq: params.freq.clone(),
t,
ane,
iath: params.iath,
ilowh: hyd_data.ilowh,
m10: hyd_data.m10,
m20: hyd_data.m20,
pop_h: params.pop_h,
pop_h_cont: params.pop_h_cont,
vturb: hyd_data.vturb,
wn_hint: hyd_data.wn_hint.clone(),
};
let result = hydlin(&hyd_params);
for ij in m..nfreq {
abso[ij] += result.absoh[ij];
emis[ij] += result.emish[ij];
}
}
}
// Detailed He II line opacity (IHE2L>0)
if params.ihe2l > 0 {
if let Some(ref he2_d) = params.he2lin_data {
let he2_params = He2linParams {
common: he2_d.common.clone(),
i0: m.saturating_sub(1),
i1: nfreq.saturating_sub(1),
ilwhe2: he2_d.ilwhe2,
mhe10: he2_d.mhe10,
mhe20: he2_d.mhe20,
};
let result = he2lin(&he2_params);
for ij in m..nfreq {
abso[ij] += result.absoh[ij];
emis[ij] += result.emish[ij];
}
}
}
// Photoionization opacity (PHTION)
if let Some(ref phtion_d) = params.phtion_data {
let mut phtion_out = PhtionOutput {
abso: &mut abso,
emis: &mut emis,
};
let phtion_params = PhtionParams {
temp: t,
fre: &params.freq,
nfre: nfreq,
rrr: phtion_d.rrr,
popul: phtion_d.popul,
photcs: phtion_d.photcs,
};
phtion(&phtion_params, &mut phtion_out);
}
// Photoionization opacity extended (PHTX)
if let Some(ref mut phtx_d) = params.phtx_data {
let mut phtx_out = PhtxOutput {
abso: &mut abso,
emis: &mut emis,
};
let phtx_params = PhtxParams {
temp: t,
depth_id: params.id,
fre: &params.freq,
nfre: nfreq,
nfreq,
nfreqc: phtx_d.nfreqc,
icon: 0,
iasv: phtx_d.iasv,
nqht: phtx_d.nqht,
popul: phtx_d.popul,
rrr: phtx_d.rrr,
level_photo: phtx_d.level_photo,
hhe_photo: phtx_d.hhe_photo,
};
phtx(&phtx_params, &mut phtx_out, phtx_d.state);
}
// Add scattering to absorption for positive imode
if params.imode >= 0 {
for ij in 0..nfreq {
abso[ij] += scat[ij];
}
}
// Correct for Lyman line wings
if params.icontl != 1 {
abso[0] -= ably1;
emis[0] -= emly1;
scat[0] -= scly1;
}
OpacResult {
abso,
emis,
scat,
avab,
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_opac_basic() {
let mut params = OpacParams {
id: 0,
t: 6000.0,
ane: 1.0e13,
nfreq: 5,
freq: vec![1.0e14, 2.0e14, 3.0e14, 4.0e14, 5.0e14],
wlam: vec![30000.0, 15000.0, 10000.0, 7500.0, 6000.0],
frx1: vec![0.0, 0.0, 0.5, 0.5, 0.5],
frx2: vec![0.0, 0.0, 0.5, 0.5, 0.5],
imode: 0,
idstd: 0,
icontl: 0,
iophli: 0,
ihyll: 0,
ihe2l: 0,
ifmol: 0,
nmlist: 0,
iath: 1,
pop_h: 1.0e16,
pop_h_cont: 1.0e10,
sce: 1.0e-3,
plan: 1.0,
hkt: 4.79928144e-11 / 6000.0,
hk: 4.79928144e-11,
bn: 1.0,
wn_hint: [1.0, 1.0, 1.0, 1.0, 1.0],
relop: 1.0,
linop_data: None,
hydlin_data: None,
he2lin_data: None,
molop_data: None,
phtion_data: None,
phtx_data: None,
};
let result = opac(&mut params);
assert!(result.abso.iter().all(|&x| x.is_finite()));
assert!(result.emis.iter().all(|&x| x.is_finite()));
assert!(result.scat.iter().all(|&x| x.is_finite()));
assert!(result.avab.is_finite());
}
}