#include <PhotoIonizationSolver.hpp>

Classes
struct	CellResult

Public Member Functions
	PhotoIonizationSolver ()=default
void	computePhotoionizationRates (const Array &Jv, double gamma[]) const
int	effectiveStages (int e) const
void	initialize (const Array &lambdav, const Array &dlambdav, bool coolingTable=false)
void	initializeIonFractions (double ionFracs[]) const
int	numBins () const
double	opacityAbs (double lambda, const double ionFracs[totalStages], double nH, double yHe, const double metalAbundances[8]) const
void	setCosmicRayHeating (bool enabled, double scaleFactor=1.0)
void	setMaxIonizationEnergy (double maxEV)
CellResult	solve (const Array &Jv, double nH, double yHe, const double metalAbundances[8], double Tprior, const double *ionFracsPrior, double TdampFactor=0.) const
CellResult	solveIonizationAtFixedT (const Array &Jv, double nH, double yHe, const double metalAbundances[8], double T, const double *ionFracsPrior, bool useEffectiveStages=false) const
int	totalEffectiveStages () const

Static Public Attributes
static constexpr int	numElements
static constexpr int	numStages []
static constexpr int	stageOffset []
static constexpr int	totalStages

Private Member Functions
double	computeCoolingRate (double T, const double ionFracs[], double nH, double yHe, const double metalAbundances[8], double ne, const int effNS[]) const
double	computeHeatingRate (const Array &Jv, const double ionFracs[], double nH, double yHe, const double metalAbundances[8], double ne, const int effNS[]) const
double	interpolateMetalCooling (int ionFracsIdx, double T, double ne) const
void	loadCoolingTable ()
double	solveIonizationBalance (double T, double ne, const double gamma[], double nH, double yHe, const double metalAbundances[8], double ionFracs[], const int effNS[]) const

Private Attributes
std::vector< double >	_coolData
std::vector< double >	_coolLogNe
std::vector< double >	_coolLogT
int	_coolNNe
int	_coolNT
bool	_cosmicRayHeating
double	_cosmicRayScale
std::vector< double >	_dlambda
std::vector< double >	_energy
std::vector< double >	_energyErg
bool	_hasCoolingTable
std::vector< double >	_lambda
int	_nStagesEff [numElements]
int	_numBins
std::vector< int >	_stabIdxForIon

Static Private Attributes
static constexpr int	_metalOffset
static constexpr int	_numMetalIons

Detailed Description

The PhotoIonizationSolver class implements a first-principles photoionization solver for computing the ionization state, equilibrium temperature, opacity, and line emissivities of ionized gas given a radiation field J_lambda (mean intensity per unit wavelength). It handles H, He, C, N, O, Ne, Mg, Si, S, and Fe through a unified pipeline.

Note: in SKIRT, the array conventionally called "Jv" stores J_lambda [W m^-3 sr^-1], not J_nu. The solver converts from wavelength to frequency/energy internally.

The solver implements:

Ionization balance via the Katz, Weinberg & Hernquist (1996) recursion
Temperature from explicit heating = cooling equilibrium (bisection)
Heating from full J_lambda integration (no 5-bin compression)
Cooling from analytic H/He channels plus pre-tabulated metal line cooling
Opacity computed analytically from ion fractions and Verner+ cross-sections

The solver is designed as a standalone utility class (not a SimulationItem) so it can be embedded as an inline ion-balance helper by any material mix that needs first-principles photoionisation. The interface accepts a prior state (T_prior, ion fractions) as the initial estimate, which accelerates convergence. All internal computations use CGS units (cm^-3 for densities, erg/s/cm^3 for rates, eV for energies).

Constructor & Destructor Documentation

◆ PhotoIonizationSolver()

PhotoIonizationSolver::PhotoIonizationSolver ( )

default

Default constructor.

Member Function Documentation

◆ computeCoolingRate()

double PhotoIonizationSolver::computeCoolingRate	(	double	T,
		const double	ionFracs[],
		double	nH,
		double	yHe,
		const double	metalAbundances[8],
		double	ne,
		const int	effNS[] ) const

private

Computes the total cooling rate [erg s^-1 cm^-3] from analytic H/He cooling channels (recombination, free-free, collisional ionization, collisional excitation) plus metal line cooling. The effNS[] array specifies effective stages per element.

◆ computeHeatingRate()

double PhotoIonizationSolver::computeHeatingRate	(	const Array &	Jv,
		const double	ionFracs[],
		double	nH,
		double	yHe,
		const double	metalAbundances[8],
		double	ne,
		const int	effNS[] ) const

private

Computes the photoheating rate [erg s^-1 cm^-3] from the radiation field and ion fractions. Integrates 4*pi*J_nu*sigma(nu)*(h*nu - IP)/(h*nu) over the wavelength grid, weighted by ion number densities n_ion. The effNS[] array specifies effective stages per element.

◆ computePhotoionizationRates()

void PhotoIonizationSolver::computePhotoionizationRates	(	const Array &	Jv,
		double	gamma[] ) const

Computes photoionization rates [s^-1] for all ions by integrating 4*pi*J_nu*sigma(nu)/(h*nu) over the wavelength grid. Results stored in gamma[].

◆ effectiveStages()

int PhotoIonizationSolver::effectiveStages ( int e ) const

inline

Returns the effective number of stages for element e.

◆ initialize()

void PhotoIonizationSolver::initialize	(	const Array &	lambdav,
		const Array &	dlambdav,
		bool	coolingTable = false )

Initializes the solver for a given radiation field wavelength grid.

Parameters

lambdav	wavelength grid bin centers [m], must be sorted ascending
dlambdav	wavelength grid bin widths [m]
coolingTable	true to load cooling table, false to skip

◆ initializeIonFractions()

void PhotoIonizationSolver::initializeIonFractions ( double ionFracs[] ) const

Initializes ion fractions to a default state (fully ionized for H/He, singly ionized for metals).

◆ interpolateMetalCooling()

double PhotoIonizationSolver::interpolateMetalCooling	(	int	ionFracsIdx,
		double	T,
		double	ne ) const

private

Returns the metal line cooling rate [erg/s per ion] for the given ionFracs index by bilinear interpolation in log T and log ne from the pre-loaded table. Returns 0 if no table is loaded or the ion has no cooling data.

◆ loadCoolingTable()

void PhotoIonizationSolver::loadCoolingTable ( )

private

Loads the CHIANTI metal line cooling table from file. Parses ion names to build the mapping from ionFracs index to table quantity index.

◆ numBins()

int PhotoIonizationSolver::numBins ( ) const

inline

Returns the number of wavelength bins in the solver's grid.

◆ opacityAbs()

double PhotoIonizationSolver::opacityAbs	(	double	lambda,
		const double	ionFracs[totalStages],
		double	nH,
		double	yHe,
		const double	metalAbundances[8] ) const

Returns the absorption opacity [cm^-1] at wavelength lambda [m] given the ion fractions and hydrogen number density.

◆ setCosmicRayHeating()

void PhotoIonizationSolver::setCosmicRayHeating	(	bool	enabled,
		double	scaleFactor = 1.0 )

Enables or disables the cosmic-ray background heating term (Indriolo+ 2007; Wolfire+ 1995). When enabled, adds a heating rate of scaleFactor * zeta_CR * 35 eV * nHI per cell. The default CR ionization rate is zeta_CR = 2e-16 s^-1 (Indriolo+ 2007).

◆ setMaxIonizationEnergy()

void PhotoIonizationSolver::setMaxIonizationEnergy ( double maxEV )

Sets the maximum ionization energy threshold [eV]. Ion stages with ionization potential above this threshold are excluded from the PIO solve (capped). CIE mode always uses full stages. Negative value means no cap (all stages active). Must be called after initialize().

◆ solve()

CellResult PhotoIonizationSolver::solve	(	const Array &	Jv,
		double	nH,
		double	yHe,
		const double	metalAbundances[8],
		double	Tprior,
		const double *	ionFracsPrior,
		double	TdampFactor = 0. ) const

Solves the ionization balance and thermal equilibrium for a single cell.

Parameters

Jv	mean intensity J_lambda at each wavelength bin [W m^-3 sr^-1]
nH	hydrogen number density [cm^-3]
yHe	helium abundance by number relative to H
metalAbundances	number abundances relative to H for {C, N, O, Ne, Mg, Si, S, Fe}
Tprior	prior temperature estimate [K]
ionFracsPrior	prior ion fractions (can be nullptr for default initialization)
TdampFactor	maximum factor by which T can change per iteration (0 or <=1 = no damping)

◆ solveIonizationAtFixedT()

CellResult PhotoIonizationSolver::solveIonizationAtFixedT	(	const Array &	Jv,
		double	nH,
		double	yHe,
		const double	metalAbundances[8],
		double	T,
		const double *	ionFracsPrior,
		bool	useEffectiveStages = false ) const

Solves the ionization balance at a fixed temperature. The temperature is kept fixed; only the ionization balance is solved (no T bisection). Photoionization rates from Jv are included, so this is a photo+collisional hybrid. Heating and cooling rates are computed for diagnostics but do not affect T. Used for CIE cells (shock-heated gas retaining snapshot T) and for re-solving ion fractions after T sign-flip damping.

Parameters

Jv	mean intensity J_lambda at each wavelength bin [W m^-3 sr^-1]
nH	hydrogen number density [cm^-3]
yHe	helium abundance by number relative to H
metalAbundances	number abundances relative to H for {C, N, O, Ne, Mg, Si, S, Fe}
T	fixed temperature [K]
ionFracsPrior	prior ion fractions (can be nullptr for default initialization)
useEffectiveStages	if true, use maxIonizationEnergy-capped stages (for sign-flip re-solve consistency); if false, use full stages (for CIE cells)

◆ solveIonizationBalance()

double PhotoIonizationSolver::solveIonizationBalance	(	double	T,
		double	ne,
		const double	gamma[],
		double	nH,
		double	yHe,
		const double	metalAbundances[8],
		double	ionFracs[],
		const int	effNS[] ) const

private

Solves the ionization balance for all species at temperature T and electron density ne, given photoionization rates gamma[]. Uses the Katz+ (1996) recursion for metals and iterative solution for H/He. Updates ionFracs[] and returns the new electron density. The effNS[] array specifies effective stages per element.

◆ totalEffectiveStages()

int PhotoIonizationSolver::totalEffectiveStages ( ) const

Returns the total number of effective (active) ion stages across all elements.

Member Data Documentation

◆ numElements

int PhotoIonizationSolver::numElements

staticconstexpr

Number of elements tracked by the solver.

◆ numStages

int PhotoIonizationSolver::numStages[]

staticconstexpr

Number of ionization stages per element (includes fully ionized stage).

◆ stageOffset

int PhotoIonizationSolver::stageOffset[]

staticconstexpr

Index of first ionization stage for each element in the ion fraction array.

◆ totalStages

int PhotoIonizationSolver::totalStages

staticconstexpr

Total number of ion fractions stored (sum of numStages).

The documentation for this class was generated from the following file:

PhotoIonizationSolver.hpp

Classes

Public Member Functions

Static Public Attributes

Private Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ PhotoIonizationSolver()

Member Function Documentation

◆ computeCoolingRate()

◆ computeHeatingRate()

◆ computePhotoionizationRates()

◆ effectiveStages()

◆ initialize()

◆ initializeIonFractions()

◆ interpolateMetalCooling()

◆ loadCoolingTable()

◆ numBins()

◆ opacityAbs()

◆ setCosmicRayHeating()

◆ setMaxIonizationEnergy()

◆ solve()

◆ solveIonizationAtFixedT()

◆ solveIonizationBalance()

◆ totalEffectiveStages()

Member Data Documentation

◆ numElements

◆ numStages

◆ stageOffset

◆ totalStages