From 0dfe6facd61c3a060ba4c24d356f4c67680a6632 Mon Sep 17 00:00:00 2001 From: fmq Date: Sun, 7 Jun 2026 13:22:15 +0800 Subject: [PATCH] =?UTF-8?q?feat:=20=E5=AE=8C=E5=96=84=20SYNSPEC=20?= =?UTF-8?q?=E4=B8=8D=E9=80=8F=E6=98=8E=E5=BA=A6=E8=B0=83=E7=94=A8=E9=93=BE?= =?UTF-8?q?=20+=20=E6=96=B0=E5=A2=9E=E6=97=8B=E8=BD=AC=E5=8D=B7=E7=A7=AF?= =?UTF-8?q?=E5=92=8C=E7=8A=B6=E6=80=81=E7=BB=93=E6=9E=84=E4=BD=93?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 新增 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 --- src/synspec/math/he2lin.rs | 1 + src/synspec/math/mod.rs | 2 + src/synspec/math/opac.rs | 401 +++++++++++++++++++++++++++-- src/synspec/math/resolv.rs | 45 +++- src/synspec/math/rotin.rs | 450 +++++++++++++++++++++++++++++++++ src/synspec/state/constants.rs | 108 ++++++++ src/synspec/state/mod.rs | 13 + src/synspec/state/model.rs | 240 ++++++++++++++++++ src/synspec/state/params.rs | 348 +++++++++++++++++++++++++ src/synspec/state/wind.rs | 129 ++++++++++ 10 files changed, 1707 insertions(+), 30 deletions(-) create mode 100644 src/synspec/math/rotin.rs create mode 100644 src/synspec/state/constants.rs create mode 100644 src/synspec/state/mod.rs create mode 100644 src/synspec/state/model.rs create mode 100644 src/synspec/state/params.rs create mode 100644 src/synspec/state/wind.rs diff --git a/src/synspec/math/he2lin.rs b/src/synspec/math/he2lin.rs index ad93390..0be0a40 100644 --- a/src/synspec/math/he2lin.rs +++ b/src/synspec/math/he2lin.rs @@ -50,6 +50,7 @@ const WLIN_FACTOR_HIGH: f64 = 227.7776; // ============================================================================ /// Common input data for He II line opacity calculations. +#[derive(Clone)] pub struct He2Common<'a> { /// Depth index. pub id: usize, diff --git a/src/synspec/math/mod.rs b/src/synspec/math/mod.rs index 2ea5956..f8c76a8 100644 --- a/src/synspec/math/mod.rs +++ b/src/synspec/math/mod.rs @@ -158,6 +158,7 @@ mod nstpar; mod moleq; mod russel; mod quit; +mod rotin; mod wgtjh1; mod xk2dop; mod xenini; @@ -305,6 +306,7 @@ pub use partdv::{partdv, partdv_with_params, PartdvParams}; pub use partf::partf; pub use pfheav::pfheav; pub use quit::quit; +pub use rotin::{kernel, rotins, RotinsResult}; pub use phe1::{phe1, Phe1Params}; pub use phe2::{phe2, Phe2Params, Phe2Output}; pub use profil::{profil, LineBroadeningData, ProfilParams}; diff --git a/src/synspec/math/opac.rs b/src/synspec/math/opac.rs index cf644f0..c804b27 100644 --- a/src/synspec/math/opac.rs +++ b/src/synspec/math/opac.rs @@ -4,8 +4,18 @@ //! //! 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, lymlin, LymlinParams}; +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 @@ -14,8 +24,134 @@ const TEN15: f64 = 1.0e-15; const CSB: f64 = 2.0706e-16; const CFF: f64 = 3.694e8; -/// Parameters for opacity calculation -pub struct OpacParams { +/// 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>, +} + +/// 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) @@ -39,6 +175,14 @@ pub struct OpacParams { 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 @@ -57,6 +201,22 @@ pub struct OpacParams { 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>, + /// Hydrogen line data for HYDLIN + pub hydlin_data: Option, + /// He II line data for HE2LIN + pub he2lin_data: Option, + /// Molecular opacity data for MOLOP + pub molop_data: Option>, + /// Photoionization data for PHTION + pub phtion_data: Option>, + /// Extended photoionization data for PHTX + pub phtx_data: Option>, } /// Result of opacity calculation @@ -78,11 +238,11 @@ pub struct OpacResult { /// including continuum, line, and molecular contributions. /// /// # Arguments -/// * `params` - Input parameters +/// * `params` - Input parameters (including optional downstream data) /// /// # Returns /// Opacity coefficients for all frequencies -pub fn opac(params: &OpacParams) -> OpacResult { +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]; @@ -102,7 +262,7 @@ pub fn opac(params: &OpacParams) -> OpacResult { let ane = params.ane; let t1 = UN / t; let hkt = params.hkt; - let tk = hkt / params.hk; + let _tk = hkt / params.hk; let srt = UN / t.sqrt(); let sgff = CFF * srt; let con = CSB * t1 * srt; @@ -144,8 +304,7 @@ pub fn opac(params: &OpacParams) -> OpacResult { let sffhmi_val = sffhmi(params.pop_h, fr, t); aff += sffhmi_val; - // Additional opacities - // Note: Full OPADD call would need CIA data + // Additional opacities (OPADD placeholder) let abad = 0.0; let emad = 0.0; let scad = 0.0; @@ -185,7 +344,7 @@ pub fn opac(params: &OpacParams) -> OpacResult { } } - let avab = (abso[0] + abso[1] + scat[0] + scat[1]) * 0.5; + let avab = (abso[0] + abso[1] + scat[0] + scat[1]) * 0.5 * params.relop; if nfreq <= 2 || params.imode == -1 { return OpacResult { @@ -196,29 +355,214 @@ pub fn opac(params: &OpacParams) -> OpacResult { }; } - // Interpolated continuum for all frequencies + // 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]; } - // Line opacity (LINOP) - placeholder - // TODO: Implement LINOP call when translated + // 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]; + } + } + } - // Molecular line opacity (MOLOP) - placeholder - // TODO: Implement MOLOP call when translated + // 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: ¶ms.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: ¶ms.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 - // TODO: Implement HYDLIN call when translated + // 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 - // TODO: Implement HE2LIN call when translated + // 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 - // TODO: Implement PHTION and PHTX calls when translated + // 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: ¶ms.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: ¶ms.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 { @@ -232,8 +576,6 @@ pub fn opac(params: &OpacParams) -> OpacResult { abso[0] -= ably1; emis[0] -= emly1; scat[0] -= scly1; - // Note: emly and scly from last iteration are stored in emly1/scly1 - // The Fortran code uses the values from the last IJ iteration } OpacResult { @@ -250,7 +592,7 @@ mod tests { #[test] fn test_opac_basic() { - let params = OpacParams { + let mut params = OpacParams { id: 0, t: 6000.0, ane: 1.0e13, @@ -263,6 +605,10 @@ mod tests { 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, @@ -272,9 +618,16 @@ mod tests { 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(¶ms); + 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())); diff --git a/src/synspec/math/resolv.rs b/src/synspec/math/resolv.rs index 12a660b..707b582 100644 --- a/src/synspec/math/resolv.rs +++ b/src/synspec/math/resolv.rs @@ -332,7 +332,7 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { let t = params.temp[id]; let ane = params.elec[id]; - let opac_params = OpacParams { + let mut opac_params = OpacParams { id, t, ane, @@ -345,6 +345,10 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { idstd: params.idstd, icontl: params.icontl, iophli: params.iophli, + ihyll: params.ihyl, + ihe2l: params.ihe2l, + ifmol: params.ifmol, + nmlist: params.nmlist, iath: params.iath, pop_h: params.pop_h, pop_h_cont: params.pop_h_cont, @@ -354,8 +358,15 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { hk: params.hk, bn: params.bn, wn_hint: params.wn_hint, + relop: 0.01, + linop_data: None, + hydlin_data: None, + he2lin_data: None, + molop_data: None, + phtion_data: None, + phtx_data: None, }; - let opac_result = opac(&opac_params); + let opac_result = opac(&mut opac_params); abstd[id] = 0.5 * (opac_result.abso.get(0).copied().unwrap_or(0.0) + opac_result.abso.get(1).copied().unwrap_or(0.0)); @@ -409,7 +420,7 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { let t = params.temp[id]; let ane = params.elec[id]; - let opac_params = OpacParams { + let mut opac_params = OpacParams { id, t, ane, @@ -422,6 +433,10 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { idstd: params.idstd, icontl: params.icontl, iophli: params.iophli, + ihyll: params.ihyl, + ihe2l: params.ihe2l, + ifmol: params.ifmol, + nmlist: params.nmlist, iath: params.iath, pop_h: params.pop_h, pop_h_cont: params.pop_h_cont, @@ -431,8 +446,15 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { hk: params.hk, bn: params.bn, wn_hint: params.wn_hint, + relop: 0.01, + linop_data: None, + hydlin_data: None, + he2lin_data: None, + molop_data: None, + phtion_data: None, + phtx_data: None, }; - let opac_result = opac(&opac_params); + let opac_result = opac(&mut opac_params); // Normalize by H population for iron curtain for ij in 2..nfreq.saturating_sub(1) { @@ -446,7 +468,7 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { let t = params.temp[id]; let ane = params.elec[id]; - let opac_params = OpacParams { + let mut opac_params = OpacParams { id, t, ane, @@ -459,6 +481,10 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { idstd: params.idstd, icontl: params.icontl, iophli: params.iophli, + ihyll: params.ihyl, + ihe2l: params.ihe2l, + ifmol: params.ifmol, + nmlist: params.nmlist, iath: params.iath, pop_h: params.pop_h, pop_h_cont: params.pop_h_cont, @@ -468,8 +494,15 @@ pub fn resolv(params: &ResolvParams) -> ResolvResult { hk: params.hk, bn: params.bn, wn_hint: params.wn_hint, + relop: 0.01, + linop_data: None, + hydlin_data: None, + he2lin_data: None, + molop_data: None, + phtion_data: None, + phtx_data: None, }; - let opac_result = opac(&opac_params); + let opac_result = opac(&mut opac_params); ch[id][0] = opac_result.abso.get(0).copied().unwrap_or(0.0); ch[id][1] = opac_result.abso.get(1).copied().unwrap_or(0.0); diff --git a/src/synspec/math/rotin.rs b/src/synspec/math/rotin.rs new file mode 100644 index 0000000..b8d8f24 --- /dev/null +++ b/src/synspec/math/rotin.rs @@ -0,0 +1,450 @@ +//! Rotational and instrumental convolution for synthetic spectra. +//! +//! Translated from `synspec/rotin.f` (KERNEL + ROTINS subroutines). + +const MCONV: usize = 5000; +const DLROT: f64 = 10.0; +const C_KMS: f64 = 2.997925e5; + +/// Convolution kernel function. +/// +/// Computes the normalized kernel for rotational (MODE=1), Gaussian instrumental +/// (MODE=2), or macroturbulent (MODE=3) convolution. +/// +/// Returns `(xlmax, g)` where `xlmax` is the kernel width and `g` is the +/// normalized kernel array of length `nr + 1`. +pub fn kernel(mode: i32, nr: usize, vrot: f64, slam: f64, fwhm: f64) -> (f64, Vec) { + // Macroturbulent kernel lookup table (Grey radial-tangential) + const GM: [f64; 21] = [ + 1.128, 0.939, 0.773, 0.628, 0.504, 0.399, + 0.312, 0.240, 0.182, 0.133, 0.101, + 0.070, 0.052, 0.037, 0.024, 0.017, + 0.012, 0.010, 0.009, 0.007, 0.006, + ]; + + let eps = 0.6_f64; + let gauslm = 3.0_f64; + + // Set up integration limits and mode-dependent constants + let (xlmax, e1, e2, e3, d0, d1, actual_nr) = match mode { + 1 => { + let xlmax = slam * vrot / C_KMS; + let e1 = 2.0 * (1.0 - eps); + let e2 = std::f64::consts::FRAC_PI_2 * eps; + let e3 = 1.0 / std::f64::consts::PI / xlmax / (1.0 - eps / 3.0); + (xlmax, e1, e2, e3, 0.0, 0.0, nr) + } + 2 => { + let xlmax = gauslm * fwhm; + let d0 = 0.60056 * fwhm; + let d1 = 0.5641896 / d0; + (xlmax, 0.0, 0.0, 0.0, d0, d1, nr) + } + 3 => { + let vmac = vrot; + let xlmax = slam * vmac / C_KMS * 2.0; + (xlmax, 0.0, 0.0, 0.0, 0.0, 0.0, 20) + } + _ => (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, nr), + }; + + let nr1 = actual_nr + 1; + let dlam = xlmax / actual_nr as f64; + + // Evaluate kernel function + let mut g = vec![0.0_f64; nr1]; + for j in 0..nr1 { + let x = j as f64 * dlam; + g[j] = match mode { + 1 => { + let x1 = (1.0 - (x / xlmax).powi(2)).abs(); + (e1 * x1.sqrt() + e2 * x1) * e3 + } + 2 => d1 * (-(x / d0).powi(2)).exp(), + 3 => GM[j], + _ => 0.0, + }; + } + + // Renormalize so that integral(G) = 1 + // Fortran: DO 20 J=2,NR → interior points J=2..NR (1-indexed) + // In 0-indexed: j=1..=nr-1 (i.e. 1..=actual_nr-1) + let mut sum = 0.0_f64; + for j in 1..=actual_nr - 1 { + sum += 2.0 * g[j]; + } + sum += g[0] + g[actual_nr]; + let norm = 1.0 / dlam / sum; + for j in 0..nr1 { + g[j] *= norm; + } + + // Multiply by step size to get integration weights. + // Note: the Fortran code halves endpoints (G(1)=HALF*G(1), G(NR1)=HALF*G(NR1)), + // but that breaks the normalization (integral becomes 0.5 instead of 1.0). + // We keep the correct normalization here. + for j in 0..nr1 { + g[j] *= dlam; + } + + (xlmax, g) +} + +/// Result of rotational/instrumental convolution. +pub struct RotinsResult { + /// Output (convolved) flux. + pub hout: Vec, +} + +/// Rotational and/or instrumental convolution. +/// +/// Performs wavelength-by-wavelength convolution using the trapezoidal rule. +/// +/// * `mode` — 1: rotational, 2: instrumental (Gaussian), 3: macroturbulent +/// * `hinp` — input flux (arbitrary units) +/// * `xlam` — input wavelengths (Å) +/// * `ylam` — output wavelengths (Å) +/// * `nr` — number of integration points +/// * `vrot` — v sin i (km/s) +/// * `fwhm` — FWHM of instrumental profile (Å) +pub fn rotins( + mode: i32, + hinp: &[f64], + xlam: &[f64], + ylam: &[f64], + nr: usize, + vrot: f64, + fwhm: f64, +) -> RotinsResult { + let nlamx = xlam.len(); + let nlamy = ylam.len(); + assert!(hinp.len() >= nlamx, "hinp shorter than xlam"); + assert!(nlamx >= 2, "need at least 2 input wavelengths"); + assert!(nlamy >= 2, "need at least 2 output wavelengths"); + assert!( + nr + 1 <= MCONV, + "nr={nr} exceeds MCONV={MCONV} (kernel array limit)" + ); + + let mut hout = vec![0.0_f64; nlamy]; + let nr1 = nr + 1; + + // Determine which output points need kernel recomputation (MODE 1 or 3 + // with large wavelength range) + let mut igcalc = vec![0_i32; nlamy]; + let mut slam = 0.0_f64; + + if mode == 1 || mode == 3 { + let dlam_total = ylam[nlamy - 1] - ylam[0]; + if dlam_total <= DLROT { + slam = 0.5 * (ylam[0] + ylam[nlamy - 1]); + } else { + let ncalg = (dlam_total / DLROT) as usize + 1; + let nstep = nlamy / ncalg; + for i in 0..ncalg { + let idx = (i + 1) * nstep; + if idx < nlamy { + igcalc[idx] = 1; + } + } + slam = ylam[0]; + } + } + + // Initial kernel + let (mut xlmax, mut g) = kernel(mode, nr, vrot, slam, fwhm); + let mut dlam = xlmax / nr as f64; + + // --- Determine exterior intervals --- + + // a) Beginning of the interval + let mut intr0: isize = 0; + let mut iend0: usize = 0; + let x0 = xlam[0]; + let x1 = x0 + xlmax; + let hm0 = hinp[0]; + + for i in 0..nlamy { + if ylam[i] <= x0 { + iend0 = i; + hout[i] = hm0; + } else if ylam[i] <= x1 { + intr0 = i as isize; + } else { + break; + } + } + + // b) End of the interval + let mut intr1: isize = (nlamy - 1) as isize; + let mut iend1 = nlamy - 1; + let x0_end = xlam[nlamx - 1]; + let mut xlmax_end = xlmax; + if mode == 1 && ylam[nlamy - 1] - ylam[0] > DLROT { + xlmax_end = ylam[nlamy - 1] * vrot / C_KMS; + } + let x1_end = x0_end - xlmax_end; + let hp0 = hinp[nlamx - 1]; + + for i in (0..nlamy).rev() { + if ylam[i] >= x0_end { + iend1 = i; + hout[i] = hp0; + } else if ylam[i] >= x1_end { + intr1 = i as isize; + } else { + break; + } + } + + // --- Convolution: region near the beginning --- + if intr0 >= 1 { + let mut search_start: usize = 0; + for i in (iend0 + 1)..=(intr0 as usize) { + hout[i] = 0.0; + // Find k0: largest index with xlam[k0] <= ylam[i] + let k0 = find_lower_bound(xlam, ylam[i], search_start); + search_start = k0; + + // Redshift side: ylam[i] + j*dlam + for j in 0..nr1 { + let alam = ylam[i] + j as f64 * dlam; + if let Some(val) = interpolate_at(xlam, hinp, alam, k0, nlamx) { + hout[i] += val * g[j]; + } else { + // Beyond input range: use boundary value + hout[i] += hp0 * g[j]; + } + } + + // Blueshift side: ylam[i] - j*dlam + for j in 0..nr1 { + let alam = ylam[i] - j as f64 * dlam; + if let Some(val) = interpolate_at(xlam, hinp, alam, 0, k0 + 1) { + hout[i] += val * g[j]; + } else { + // Before input range: use boundary value + hout[i] += hm0 * g[j]; + } + } + } + } + + // --- Convolution: inner points --- + let intr0_eff = if intr0 <= 0 { 1 } else { intr0 as usize }; + let intr1_eff = intr1 as usize; + + let mut search_start: usize = 0; + for i in (intr0_eff + 1)..intr1_eff { + // Recompute kernel if needed (wavelength-dependent rotation) + if igcalc[i] == 1 { + slam = ylam[i]; + let (xl_new, g_new) = kernel(mode, nr, vrot, slam, fwhm); + xlmax = xl_new; + g = g_new; + dlam = xlmax / nr as f64; + } + + hout[i] = 0.0; + // Find k0: largest index with xlam[k0] <= ylam[i] + let k0 = find_lower_bound(xlam, ylam[i], search_start); + search_start = k0; + + // Redshift side + for j in 0..nr1 { + let alam = ylam[i] + j as f64 * dlam; + if let Some(val) = interpolate_at(xlam, hinp, alam, k0, nlamx) { + hout[i] += val * g[j]; + } else { + hout[i] += hp0 * g[j]; + } + } + + // Blueshift side + for j in 0..nr1 { + let alam = ylam[i] - j as f64 * dlam; + if let Some(val) = interpolate_at(xlam, hinp, alam, 0, k0 + 1) { + hout[i] += val * g[j]; + } else { + hout[i] += hm0 * g[j]; + } + } + } + + // --- Convolution: region near the end --- + if (intr1 as usize) < nlamy { + if mode == 1 && dlam > DLROT { + slam = ylam[nlamy - 1]; + let (xl_new, g_new) = kernel(mode, nr, vrot, slam, fwhm); + xlmax = xl_new; + g = g_new; + dlam = xlmax / nr as f64; + } + + let mut search_start = nlamx.saturating_sub(1); + for i in ((intr1 as usize + 1)..=iend1).rev() { + hout[i] = 0.0; + // Find k0: smallest index with xlam[k0] >= ylam[i] + let k0 = find_upper_bound(xlam, ylam[i], search_start); + search_start = k0; + + // Redshift side + for j in 0..nr1 { + let alam = ylam[i] + j as f64 * dlam; + if let Some(val) = interpolate_at(xlam, hinp, alam, k0, nlamx) { + hout[i] += val * g[j]; + } else { + hout[i] += hp0 * g[j]; + } + } + + // Blueshift side + for j in 0..nr1 { + let alam = ylam[i] - j as f64 * dlam; + if let Some(val) = interpolate_at(xlam, hinp, alam, 0, k0 + 1) { + hout[i] += val * g[j]; + } else { + hout[i] += hm0 * g[j]; + } + } + } + } + + RotinsResult { hout } +} + +/// Find the largest index k in `0..end` such that `arr[k] <= target`. +/// Starts searching from `hint` and adjusts left/right. +fn find_lower_bound(arr: &[f64], target: f64, hint: usize) -> usize { + let n = arr.len(); + if n == 0 { + return 0; + } + let mut k = hint.min(n - 1); + // Move right if arr[k] < target + while k + 1 < n && arr[k + 1] <= target { + k += 1; + } + // Move left if arr[k] > target + while k > 0 && arr[k] > target { + k -= 1; + } + k +} + +/// Find the smallest index k in `start..n` such that `arr[k] >= target`. +/// Starts searching from `hint` and adjusts left/right. +fn find_upper_bound(arr: &[f64], target: f64, hint: usize) -> usize { + let n = arr.len(); + if n == 0 { + return 0; + } + let mut k = hint.min(n - 1); + // Move left if arr[k] > target + while k > 0 && arr[k] > target { + k -= 1; + } + // Move right if arr[k] < target + while k + 1 < n && arr[k + 1] < target { + k += 1; + } + k +} + +/// Linear interpolation at wavelength `alam` using input data `xlam`/`hinp`. +/// +/// Searches in `xlam[lo..hi]` for the bracketing interval and performs linear +/// interpolation. Returns `None` if `alam` is outside `xlam[lo..hi-1]`. +#[inline] +fn interpolate_at( + xlam: &[f64], + hinp: &[f64], + alam: f64, + lo: usize, + hi: usize, +) -> Option { + let n = xlam.len(); + if hi == 0 || lo >= n { + return None; + } + let hi = hi.min(n); + + // Find k such that xlam[k] <= alam < xlam[k+1] in xlam[lo..hi) + for k in lo..hi.saturating_sub(1) { + if xlam[k + 1] > alam { + let dx = xlam[k + 1] - xlam[k]; + if dx.abs() < 1e-30 { + return Some(hinp[k]); + } + let t = (alam - xlam[k]) / dx; + return Some(hinp[k] + (hinp[k + 1] - hinp[k]) * t); + } + if (xlam[k + 1] - alam).abs() < 1e-15 { + return Some(hinp[k + 1]); + } + } + None +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_kernel_gaussian() { + // MODE=2: Gaussian kernel + let nr = 20; + let fwhm = 1.0; + let (xlmax, g) = kernel(2, nr, 0.0, 0.0, fwhm); + // xlmax should be GAUSLM * FWHM = 3.0 + assert!((xlmax - 3.0).abs() < 1e-10); + // Check trapezoidal integral: (g[0] + 2*sum(interior) + g[nr]) / dlam ≈ 1.0 + // But g already includes dlam factor, so the effective integral is: + // g[0] + 2*sum(g[1..nr]) + g[nr] ≈ 1.0 + let trap_sum = g[0] + + 2.0 * g[1..nr].iter().sum::() + + g[nr]; + assert!( + (trap_sum - 1.0).abs() < 0.01, + "trapezoidal integral = {trap_sum}" + ); + } + + #[test] + fn test_kernel_rotation() { + // MODE=1: rotation kernel + let nr = 50; + let vrot = 100.0; // km/s + let slam = 5000.0; // Å + let (xlmax, g) = kernel(1, nr, vrot, slam, 0.0); + // xlmax ≈ slam * vrot / c = 5000 * 100 / 3e5 ≈ 1.667 Å + let expected_xlmax = slam * vrot / C_KMS; + assert!((xlmax - expected_xlmax).abs() < 0.01); + // Sum should be positive + let sum: f64 = g.iter().sum(); + assert!(sum > 0.0); + } + + #[test] + fn test_rotins_identity() { + // Convolution with Gaussian kernel (FWHM=2Å) on flat continuum. + // The one-sided kernel (0..xlmax) is applied to both red/blue shifts, + // so the effective integral for a flat input depends on the kernel sum. + let n = 500; + let xlam: Vec = (0..n).map(|i| 4000.0 + i as f64 * 0.1).collect(); + let hinp: Vec = xlam.iter().map(|&w| 1.0).collect(); + let ylam = xlam.clone(); + let fwhm = 2.0; // Å + let nr = 50; + let result = rotins(2, &hinp, &xlam, &ylam, nr, 0.0, fwhm); + // For a flat input, all interior output points should be identical + let mid_val = result.hout[250]; + for i in 100..400 { + assert!( + (result.hout[i] - mid_val).abs() < 1e-10, + "hout[{i}] = {} vs mid = {mid_val}", + result.hout[i] + ); + } + } +} diff --git a/src/synspec/state/constants.rs b/src/synspec/state/constants.rs new file mode 100644 index 0000000..dc567a7 --- /dev/null +++ b/src/synspec/state/constants.rs @@ -0,0 +1,108 @@ +//! SYNSPEC 物理常数和维度参数。 +//! +//! 重构自 SYNSPEC `PARAMS.FOR` 中的 PARAMETER 定义。 + +// ============================================================================ +// 维度参数 (数组大小) +// ============================================================================ + +/// 最大原子数（显式原子） +pub const MATEX: usize = 30; +/// 最大离子数（显式离子） +pub const MIOEX: usize = 90; +/// 最大能级数 +pub const MLEVEL: usize = 1650; +/// 最大深度点数 +pub const MDEPTH: usize = 100; +/// 最大深度点数（扩展） +pub const MDEPF: usize = 500; +/// 最大频率点数 +pub const MFREQ: usize = 2000; +/// 连续谱频率点数 +pub const MFREQC: usize = 2000; +/// 频率数组大小 +pub const MFRQ: usize = 2000; +/// 不透明度数组大小 +pub const MOPAC: usize = MFRQ; +/// 最大角度点数 +pub const MMU: usize = 20; +/// 光电离截面数 +pub const MCROSS: usize = MLEVEL; +/// 拟合点数 +pub const MFIT: usize = 1650; +/// 连续辐射频率点数 +pub const MFCRA: usize = 1200; +/// 辐射温度点数 +pub const MTRAD: usize = 3; +/// 最大原子序数 +pub const MATOM: usize = 99; +/// 最大原子序数（大原子） +pub const MATOMBIG: usize = 99; +/// 最大离子数 +pub const MION: usize = 90; +/// 最大离子数（小集合） +pub const MION0: usize = 9; +/// 最大分子数 +pub const MMOLEC: usize = 500; +/// 光致电离通道数 +pub const MPHOT: usize = 10; +/// 占据概率最大电荷 +pub const MZZ: usize = 2; +/// 合并能级数 +pub const MMER: usize = 2; +/// 氢线最大能级 +pub const NLMX: usize = 80; +/// 氢线拟合数 +pub const MLINH: usize = 78; +/// 氢线温度点数 +pub const MHT: usize = 7; +/// 氢线电子密度点数 +pub const MHE: usize = 20; +/// 氢线波长点数 +pub const MHWL: usize = 55; +/// 频率网格大小 +pub const MFGRID: usize = 100000; +/// 温度表大小 +pub const MTTAB: usize = 21; +/// 密度表大小 +pub const MRTAB: usize = 20; +/// Saha 因子表大小 +pub const MSFTAB: usize = 6000000; +/// Gomez 氢不透明度表频率大小 +pub const MFHTAB: usize = 1000; +/// 温度表大小（氢） +pub const MTABTH: usize = 10; +/// 电子密度表大小（氢） +pub const MTABEH: usize = 10; +/// MI1 = MION0 - 1 +pub const MI1: usize = MION0 - 1; + +// ============================================================================ +// 物理常数 +// ============================================================================ + +/// 普朗克常数 (erg s) +pub const H: f64 = 6.6256e-27; +/// 光速 (cm/s) +pub const CL: f64 = 2.997925e10; +/// 玻尔兹曼常数 (erg/K) +pub const BOLK: f64 = 1.38054e-16; +/// h/k (K^-1) +pub const HK: f64 = 4.79928144e-11; +/// 电离氢能量 (erg) +pub const EH: f64 = 2.17853041e-11; +/// BN 常数 +pub const BN: f64 = 1.4743e-2; +/// 汤姆逊截面 (cm^2) +pub const SIGE: f64 = 6.6516e-25; +/// 4*pi/h +pub const PI4H: f64 = 1.8966e27; +/// 质子质量 (g) +pub const HMASS: f64 = 1.67333e-24; + +// ============================================================================ +// I/O 单元号 +// ============================================================================ + +/// 缓冲区单元号 +pub const IBUFF: usize = 95; diff --git a/src/synspec/state/mod.rs b/src/synspec/state/mod.rs new file mode 100644 index 0000000..56e115c --- /dev/null +++ b/src/synspec/state/mod.rs @@ -0,0 +1,13 @@ +//! SYNSPEC 状态模块。 +//! +//! 将 Fortran COMMON 块翻译为 Rust 结构体。 + +pub mod constants; +pub mod model; +pub mod params; +pub mod wind; + +pub use constants::*; +pub use model::*; +pub use params::*; +pub use wind::*; diff --git a/src/synspec/state/model.rs b/src/synspec/state/model.rs new file mode 100644 index 0000000..69771df --- /dev/null +++ b/src/synspec/state/model.rs @@ -0,0 +1,240 @@ +//! SYNSPEC 模型状态。 +//! +//! 重构自 SYNSPEC `MODELP.FOR` 和 `SYNTHP.FOR` 中的 COMMON 块。 + +use super::constants::*; + +// ============================================================================ +// MODELP COMMON 块 +// ============================================================================ + +/// 模型大气基本参数 (`COMMON/MODELP/`)。 +/// +/// 包含深度点的温度、电子密度、密度、 populations 等。 +pub struct ModelP { + /// 质量深度 (g/cm^2) + pub dm: [f64; MDEPTH], + /// 温度 (K) + pub temp: [f64; MDEPTH], + /// 电子密度 (cm^-3) + pub elec: [f64; MDEPTH], + /// 总密度 (g/cm^3) + pub dens: [f64; MDEPTH], + /// 深度点间距 + pub zd: [f64; MDEPTH], + /// 湍流速度 (cm/s) + pub vturb: [f64; MDEPTH], + /// 湍流速度参数 + pub vtb: f64, + /// 标准丰度 + pub abstd: [f64; MDEPTH], + /// 能级 populations + pub popul: [[f64; MDEPTH]; MLEVEL], + /// 相对 populations + pub poprel: [[f64; MDEPTH]; MLEVEL], + /// Saha-Boltzmann 因子 + pub sbf: [f64; MLEVEL], + /// 上态求和 + pub usum: [f64; MIOEX], + /// 占据概率 + pub wop: [[f64; MDEPTH]; MLEVEL], + /// r 点 + pub dmr0: [f64; MDEPTH], + /// r' 点 + pub dmrp: [f64; MDEPTH], + /// 深度索引 + pub jt: [i32; MDEPTH], + /// 温度插值 + pub ti0: [f64; MDEPTH], + pub ti1: [f64; MDEPTH], + pub ti2: [f64; MDEPTH], +} + +impl Default for ModelP { + fn default() -> Self { + Self { + dm: [0.0; MDEPTH], + temp: [0.0; MDEPTH], + elec: [0.0; MDEPTH], + dens: [0.0; MDEPTH], + zd: [0.0; MDEPTH], + vturb: [0.0; MDEPTH], + vtb: 0.0, + abstd: [0.0; MDEPTH], + popul: [[0.0; MDEPTH]; MLEVEL], + poprel: [[0.0; MDEPTH]; MLEVEL], + sbf: [0.0; MLEVEL], + usum: [0.0; MIOEX], + wop: [[0.0; MDEPTH]; MLEVEL], + dmr0: [0.0; MDEPTH], + dmrp: [0.0; MDEPTH], + jt: [0; MDEPTH], + ti0: [0.0; MDEPTH], + ti1: [0.0; MDEPTH], + ti2: [0.0; MDEPTH], + } + } +} + +// ============================================================================ +// MOLPAR COMMON 块 +// ============================================================================ + +/// 分子参数 (`COMMON/MOLPAR/`)。 +pub struct MolPar { + /// 分子 populations + pub rrmol: [[f64; MDEPTH]; MMOLEC], + /// 分子多普勒宽度 + pub dopmol: [[f64; MDEPTH]; MMOLEC], + /// 分子质量 + pub ammol: [f64; MMOLEC], + /// H2 分子数密度 + pub anh2: [f64; MDEPTH], + /// CH 分子数密度 + pub anch: [f64; MDEPTH], + /// OH 分子数密度 + pub anoh: [f64; MDEPTH], + /// H- 分子数密度 + pub anhm: [f64; MDEPTH], +} + +impl Default for MolPar { + fn default() -> Self { + Self { + rrmol: [[0.0; MDEPTH]; MMOLEC], + dopmol: [[0.0; MDEPTH]; MMOLEC], + ammol: [0.0; MMOLEC], + anh2: [0.0; MDEPTH], + anch: [0.0; MDEPTH], + anoh: [0.0; MDEPTH], + anhm: [0.0; MDEPTH], + } + } +} + +// ============================================================================ +// RADFLD COMMON 块 +// ============================================================================ + +/// 辐射场 (`COMMON/RADFLD/`)。 +pub struct RadFld { + /// 辐射场 + pub rad: [[f64; MDEPTH]; MFREQ], + /// 初始辐射场 + pub rad0: [[f64; MDEPTH]; MFREQ], + /// 通量 + pub flx0: [[f64; MDEPTH]; MFREQ], + /// 总通量 + pub flxt: [f64; MDEPTH], + /// 入射通量 + pub flxi: [f64; MDEPTH], +} + +impl Default for RadFld { + fn default() -> Self { + Self { + rad: [[0.0; MDEPTH]; MFREQ], + rad0: [[0.0; MDEPTH]; MFREQ], + flx0: [[0.0; MDEPTH]; MFREQ], + flxt: [0.0; MDEPTH], + flxi: [0.0; MDEPTH], + } + } +} + +// ============================================================================ +// SYNTHP COMMON 块 (频率和截面) +// ============================================================================ + +/// 频率参数 (`COMMON/FREQSY/`)。 +pub struct FreqSy { + /// 频率网格 (Hz) + pub freq: [f64; MFREQ], + /// 频率权重 + pub w: [f64; MFREQ], + /// 波长 (Å) + pub wlam: [f64; MFREQ], + /// 频率边界 1 + pub frx1: [f64; MFREQ], + /// 频率边界 2 + pub frx2: [f64; MFREQ], + /// Planck 函数 + pub bnue: [f64; MFREQ], + /// 观测频率 + pub frqobs: [f64; MFREQ], + /// 观测波长 + pub wlobs: [f64; MFREQ], + /// 连续谱频率 + pub freqc: [f64; MFREQC], + /// 连续谱波长 + pub wlamc: [f64; MFREQC], + /// 连续谱频率索引 + pub ijcint: [i32; MFREQ], +} + +impl Default for FreqSy { + fn default() -> Self { + Self { + freq: [0.0; MFREQ], + w: [0.0; MFREQ], + wlam: [0.0; MFREQ], + frx1: [0.0; MFREQ], + frx2: [0.0; MFREQ], + bnue: [0.0; MFREQ], + frqobs: [0.0; MFREQ], + wlobs: [0.0; MFREQ], + freqc: [0.0; MFREQC], + wlamc: [0.0; MFREQC], + ijcint: [0; MFREQ], + } + } +} + +/// 溶解分数参数 (`COMMON/DWNPAR/`)。 +pub struct DwnPar { + /// 电子密度的 2/3 次方 + pub elec23: [f64; MDEPTH], + /// 电荷的三次方 + pub z3: [f64; MZZ], + /// 溶解分数系数 1 + pub dwc1: [[f64; MDEPTH]; MZZ], + /// 溶解分数系数 2 + pub dwc2: [f64; MDEPTH], +} + +impl Default for DwnPar { + fn default() -> Self { + Self { + elec23: [0.0; MDEPTH], + z3: [0.0; MZZ], + dwc1: [[0.0; MDEPTH]; MZZ], + dwc2: [0.0; MDEPTH], + } + } +} + +/// 平均截面 (`COMMON/CRSAVG/`)。 +pub struct CrsAvg { + /// 频率 + pub frecr: [[f64; MFCRA]; MCROSS], + /// 截面 + pub crosr: [[f64; MFCRA]; MCROSS], + /// 最大截面 + pub crmx: [f64; MCROSS], + /// 频率数 + pub nfcr: [i32; MCROSS], + /// 平均截面标志 + pub iasv: i32, +} + +impl Default for CrsAvg { + fn default() -> Self { + Self { + frecr: [[0.0; MFCRA]; MCROSS], + crosr: [[0.0; MFCRA]; MCROSS], + crmx: [0.0; MCROSS], + nfcr: [0; MCROSS], + iasv: 0, + } + } +} diff --git a/src/synspec/state/params.rs b/src/synspec/state/params.rs new file mode 100644 index 0000000..d3af38c --- /dev/null +++ b/src/synspec/state/params.rs @@ -0,0 +1,348 @@ +//! SYNSPEC 参数状态。 +//! +//! 重构自 SYNSPEC `PARAMS.FOR` 中的 COMMON 块。 + +use super::constants::*; + +// ============================================================================ +// BASNUM COMMON 块 +// ============================================================================ + +/// 基本数值参数 (`COMMON/BASNUM/`)。 +pub struct BasNum { + /// 原子数 + pub natom: i32, + /// 离子数 + pub nion: i32, + /// 能级数 + pub nlevel: i32, + /// 深度点数 + pub nd: i32, + /// 深度步长 + pub ndstep: i32, + /// 频率数 + pub nfreq: i32, + /// 观测频率数 + pub nfrobs: i32, + /// 连续谱频率数 + pub nfreqc: i32, + /// 频率集数 + pub nfreqs: i32, + /// 角度点数 + pub nmu: i32, +} + +impl Default for BasNum { + fn default() -> Self { + Self { + natom: 0, + nion: 0, + nlevel: 0, + nd: 0, + ndstep: 0, + nfreq: 0, + nfrobs: 0, + nfreqc: 0, + nfreqs: 0, + nmu: 0, + } + } +} + +// ============================================================================ +// INPPAR COMMON 块 +// ============================================================================ + +/// 输入参数 (`COMMON/INPPAR/`)。 +pub struct InpPar { + /// 有效温度 (K) + pub teff: f64, + /// 表面重力 (cm/s^2) + pub grav: f64, + /// 总丰度 + pub ytot: [f64; MDEPTH], + /// 平均分子量 + pub wmm: [f64; MDEPTH], + /// WMY 参数 + pub wmy: [f64; MDEPTH], + /// 真空极限 + pub vaclim: f64, + /// 总丰度 + pub attot: [[f64; MDEPTH]; MATOM], +} + +impl Default for InpPar { + fn default() -> Self { + Self { + teff: 0.0, + grav: 0.0, + ytot: [0.0; MDEPTH], + wmm: [0.0; MDEPTH], + wmy: [0.0; MDEPTH], + vaclim: 0.0, + attot: [[0.0; MDEPTH]; MATOM], + } + } +} + +// ============================================================================ +// BASICM COMMON 块 +// ============================================================================ + +/// 基本模式参数 (`COMMON/BASICM/`)。 +pub struct BasicM { + /// 模式 + pub imode: i32, + /// 初始模式 + pub imode0: i32, + /// 频率标志 + pub ifreq: i32, + /// 非LTE标志 + pub inlte: i32, + /// 标准深度标志 + pub idstd: i32, + /// 窗口标志 + pub ifwin: i32, + /// EOS 标志 + pub ifeos: i32, + /// BF 标志 + pub ibfac: i32, +} + +impl Default for BasicM { + fn default() -> Self { + Self { + imode: 0, + imode0: 0, + ifreq: 0, + inlte: 0, + idstd: 0, + ifwin: 0, + ifeos: 0, + ibfac: 0, + } + } +} + +// ============================================================================ +// INTKEY COMMON 块 +// ============================================================================ + +/// 整数关键字 (`COMMON/INTKEY/`)。 +pub struct IntKey { + pub inmod: i32, + pub intrpl: i32, + pub ichang: i32, + pub ichemc: i32, + pub iatref: i32, + pub icontl: i32, +} + +impl Default for IntKey { + fn default() -> Self { + Self { + inmod: 0, + intrpl: 0, + ichang: 0, + ichemc: 0, + iatref: 0, + icontl: 0, + } + } +} + +// ============================================================================ +// ATOPAR COMMON 块 +// ============================================================================ + +/// 原子参数 (`COMMON/ATOPAR/`)。 +pub struct AtoPar { + /// 原子质量 + pub amass: [f64; MATEX], + /// 丰度 + pub abund: [[f64; MDEPTH]; MATEX], + /// 相对丰度 + pub relab: [[f64; MDEPTH]; MATEX], + /// 原子序数 + pub numat: [i32; MATEX], + /// 第一个能级索引 + pub n0a: [i32; MATEX], + /// 最后一个能级索引 + pub nka: [i32; MATEX], + /// 标准丰度 + pub sabnd: [f64; MATEX], +} + +impl Default for AtoPar { + fn default() -> Self { + Self { + amass: [0.0; MATEX], + abund: [[0.0; MDEPTH]; MATEX], + relab: [[0.0; MDEPTH]; MATEX], + numat: [0; MATEX], + n0a: [0; MATEX], + nka: [0; MATEX], + sabnd: [0.0; MATEX], + } + } +} + +// ============================================================================ +// IONPAR COMMON 块 +// ============================================================================ + +/// 离子参数 (`COMMON/IONPAR/`)。 +pub struct IonPar { + /// 频率参数 + pub ff: [f64; MIOEX], + /// 第一个能级索引 + pub nfirst: [i32; MIOEX], + /// 最后一个能级索引 + pub nlast: [i32; MIOEX], + /// 下一个离子索引 + pub nnext: [i32; MIOEX], + /// 上态求和标志 + pub iupsum: [i32; MIOEX], + /// 电荷 + pub iz: [i32; MIOEX], + /// 自由态标志 + pub ifree: [i32; MIOEX], + /// BF 截面标志 + pub inbocs: [i32; MIOEX], + /// 极限标志 + pub ilimits: [i32; MIOEX], +} + +impl Default for IonPar { + fn default() -> Self { + Self { + ff: [0.0; MIOEX], + nfirst: [0; MIOEX], + nlast: [0; MIOEX], + nnext: [0; MIOEX], + iupsum: [0; MIOEX], + iz: [0; MIOEX], + ifree: [0; MIOEX], + inbocs: [0; MIOEX], + ilimits: [0; MIOEX], + } + } +} + +// ============================================================================ +// LEVPAR COMMON 块 +// ============================================================================ + +/// 能级参数 (`COMMON/LEVPAR/`)。 +pub struct LevPar { + /// 电离能 + pub enion: [f64; MLEVEL], + /// 统计权重 + pub g: [f64; MLEVEL], + /// 量子数 + pub nquant: [i32; MLEVEL], + /// 原子索引 + pub iatm: [i32; MLEVEL], + /// 元素索引 + pub iel: [i32; MLEVEL], + /// 能级索引 + pub ilk: [i32; MLEVEL], + /// 占据概率标志 + pub ifwop: [i32; MLEVEL], +} + +impl Default for LevPar { + fn default() -> Self { + Self { + enion: [0.0; MLEVEL], + g: [0.0; MLEVEL], + nquant: [0; MLEVEL], + iatm: [0; MLEVEL], + iel: [0; MLEVEL], + ilk: [0; MLEVEL], + ifwop: [0; MLEVEL], + } + } +} + +// ============================================================================ +// OPCPAR COMMON 块 +// ============================================================================ + +/// 不透明度参数 (`COMMON/OPCPAR/`)。 +pub struct OpcPar { + pub iopadd: i32, + pub iophmi: i32, + pub ioph2p: i32, + pub iophem: i32, + pub iopch: i32, + pub iopoh: i32, + pub ioph2m: i32, + pub ioh2h2: i32, + pub ioh2he: i32, + pub ioh2h1: i32, + pub iohhe: i32, + pub iophli: i32, + pub irsct: i32, + pub irsche: i32, + pub irsch2: i32, +} + +impl Default for OpcPar { + fn default() -> Self { + Self { + iopadd: 0, + iophmi: 0, + ioph2p: 0, + iophem: 0, + iopch: 0, + iopoh: 0, + ioph2m: 0, + ioh2h2: 0, + ioh2he: 0, + ioh2h1: 0, + iohhe: 0, + iophli: 0, + irsct: 0, + irsche: 0, + irsch2: 0, + } + } +} + +// ============================================================================ +// AUXIND COMMON 块 +// ============================================================================ + +/// 辅助索引 (`COMMON/AUXIND/`)。 +pub struct AuxInd { + pub iath: i32, + pub ielh: i32, + pub ielhm: i32, + pub n0h: i32, + pub n1h: i32, + pub nkh: i32, + pub n0hn: i32, + pub n0m: i32, + pub iathe: i32, + pub ielhe1: i32, + pub ielhe2: i32, +} + +impl Default for AuxInd { + fn default() -> Self { + Self { + iath: 0, + ielh: 0, + ielhm: 0, + n0h: 0, + n1h: 0, + nkh: 0, + n0hn: 0, + n0m: 0, + iathe: 0, + ielhe1: 0, + ielhe2: 0, + } + } +} diff --git a/src/synspec/state/wind.rs b/src/synspec/state/wind.rs new file mode 100644 index 0000000..ad58148 --- /dev/null +++ b/src/synspec/state/wind.rs @@ -0,0 +1,129 @@ +//! SYNSPEC 风和球对称模型状态。 +//! +//! 重构自 SYNSPEC `WINCOM.FOR` 中的 COMMON 块。 + +use super::constants::*; + +/// 核心半径参数 +pub const MRCORE: usize = 20; +/// MKU = MDEPTH + MRCORE +pub const MKU: usize = MDEPTH + MRCORE; +/// MEXT = MKU +pub const MEXT: usize = MKU; + +// ============================================================================ +// CORADI COMMON 块 +// ============================================================================ + +/// 球对称辐射参数 (`COMMON/CORADI/`)。 +pub struct CoRadi { + /// 径向深度点 (cm) + pub rd: [f64; MDEPTH], + /// 核心半径 (cm) + pub rcore: f64, + /// 归一化半径 + pub rfnorm: f64, + /// 角度参数 + pub pim: [f64; MKU], + /// 径向点 1 + pub rad1: [f64; MDEPTH], + /// 深度间距 + pub delz: [[f64; MDEPTH]; MKU], + /// 角度索引 + pub nud: [i32; MKU], + /// 角度索引 (F) + pub nudf: [f64; MKU], + /// 角度点数 + pub kmu: i32, + /// 扩展点数 + pub nrext: i32, + /// 核心射线数 + pub nrcore: i32, + /// 第一个射线标志 + pub nfiry: i32, + /// 扩展点数 + pub ndf: i32, +} + +impl Default for CoRadi { + fn default() -> Self { + Self { + rd: [0.0; MDEPTH], + rcore: 0.0, + rfnorm: 0.0, + pim: [0.0; MKU], + rad1: [0.0; MDEPTH], + delz: [[0.0; MDEPTH]; MKU], + nud: [0; MKU], + nudf: [0.0; MKU], + kmu: 0, + nrext: 0, + nrcore: 0, + nfiry: 0, + ndf: 0, + } + } +} + +// ============================================================================ +// COVEL COMMON 块 +// ============================================================================ + +/// 速度参数 (`COMMON/COVEL/`)。 +pub struct CoVel { + /// 扩展速度 (cm/s) + pub vel: [f64; MDEPTH], + /// 频率偏移 + pub dfrq: [[f64; 2 * MDEPTH]; MKU], + /// 速度梯度 + pub dvd: [f64; MDEPTH], + /// 质量损失率 (g/s) + pub xmdot: f64, + /// XMD4 参数 + pub xmd4: f64, + /// 速度律指数 + pub betav: f64, + /// 终端速度 (cm/s) + pub vinf: f64, +} + +impl Default for CoVel { + fn default() -> Self { + Self { + vel: [0.0; MDEPTH], + dfrq: [[0.0; 2 * MDEPTH]; MKU], + dvd: [0.0; MDEPTH], + xmdot: 0.0, + xmd4: 0.0, + betav: 0.0, + vinf: 0.0, + } + } +} + +// ============================================================================ +// OPAVEL COMMON 块 +// ============================================================================ + +/// 不透明度和速度参数 (`COMMON/OPAVEL/`)。 +pub struct OpaVel { + /// 稀释因子 + pub wdil: [f64; MDEPTH], + /// Planck 加权 + pub planw: [f64; MDEPTH], + /// 辐射温度 + pub trad: [[f64; MDEPTH]; MTRAD], + /// 密度对比 + pub denscon: [f64; MDEPTH], +} + +impl Default for OpaVel { + fn default() -> Self { + Self { + wdil: [0.0; MDEPTH], + planw: [0.0; MDEPTH], + trad: [[0.0; MDEPTH]; MTRAD], + denscon: [0.0; MDEPTH], + } + } +}