jwspecabund.forward
Bayesian forward model for emission-line abundance determination.
Implements the forward modelling approach of Cullen et al. (2025): free parameters (ionic abundances, T_e, n_e) are sampled via MCMC or nested sampling. For each parameter vector, CEL emissivities from PyNEB and the Aller (1984) Hbeta recombination formula predict line flux ratios that are compared to the observed ratios via a Gaussian likelihood.
References
Aller L. H., 1984, Physics of Thermal Gaseous Nebulae (Reidel)
Cullen et al. (2025) — forward modelling for EXCELS galaxies
Functions
|
Run a Bayesian forward model on emission-line fluxes. |
Hbeta volume emissivity coefficient using the Aller (1984) Case B formula. |
- jwspecabund.forward.forward_model(line_fluxes, line_errors, *, sampler='emcee', n_walkers=32, n_steps=5000, n_burn=1000, n_live=500, seed=42, progress=True)[source]
Run a Bayesian forward model on emission-line fluxes.
Free parameters are log(T_e), log(n_e), and one log(ionic abundance) per detected ion. The model predicts line/Hbeta flux ratios using PyNEB CEL emissivities and the Aller (1984) Hbeta formula, and the Gaussian likelihood compares predictions to observations.
- Parameters:
line_fluxes (dict) –
{line_name: flux}. Must include"HBETA".line_errors (dict) –
{line_name: flux_err}.sampler (str) –
"emcee"(default) or"dynesty".n_walkers (int) – Number of walkers for emcee (default 32).
n_steps (int) – Total MCMC steps for emcee (default 5000).
n_burn (int) – Burn-in steps to discard for emcee (default 1000).
n_live (int) – Number of live points for dynesty (default 500).
seed (int) – Random seed (default 42).
progress (bool) – Show progress bar (default
True).
- Returns:
Result dictionary with keys:
param_names,param_labels— parameter metadatasamples— raw posterior samples, shape(n, n_dim)Te,Te_err,Te_posteriorne,ne_err,ne_posteriorOH,OH_err,OH_posteriorNO,NO_err,NO_posterior(if nitrogen detected)ionic— dict of median log-ionic abundancesionic_posteriors— dict of log-ionic posterior arraysratio_names— observable labels
- Return type:
Examples
>>> import jwspecabund >>> fluxes = { ... "HBETA": 1.0, "OIII_5007": 5.0, ... "OIII_4363": 0.05, "NeIII_3869": 0.3, ... } >>> errors = {k: 0.05 * v for k, v in fluxes.items()} >>> res = jwspecabund.forward_model(fluxes, errors, n_steps=2000) >>> print(f"12+log(O/H) = {res['OH']:.2f}")
- jwspecabund.forward.hbeta_emissivity_aller84(Te)[source]
Hbeta volume emissivity coefficient using the Aller (1984) Case B formula.
\[\varepsilon_{H\beta} = \alpha_{H\beta}^{\rm eff}(T)\;h\nu_{H\beta}\]where \(\alpha_{H\beta}^{\rm eff} = 3.03\times10^{-14}\,(T/10^4)^{-0.874}\) cm3 s-1.
Unlike PyNEB’s recombination tables, this formula is valid at all temperatures — no upper-bound limitation at 30,000 K.