jwspecfit

jwspecfit — JWST NIRSpec emission-line fitting.

Functions

fit_lines(spectrum, z, *[, grating, R, ...])

Fit emission lines in a spectrum.

class jwspecfit.BroadFitResult(best_fit, selected_model, bic_narrow, bic_broad1, bic_broad2, bic_both, all_fits, bic_bootstrap=<factory>, oiii_selected='off', bic_oiii_off=nan, bic_oiii_broad1=nan, bic_oiii_broad2=nan, bic_oiii_both=nan, bic_oiii_bootstrap=<factory>, hei_selected='off', bic_hei_off=nan, bic_hei_broad1=nan, bic_hei_broad2=nan, bic_hei_both=nan, bic_hei_bootstrap=<factory>)[source]

Bases: object

Result of broad component model selection.

Parameters:

best_fit (FitResult)
selected_model (str)
bic_narrow (float)
bic_broad1 (float)
bic_broad2 (float)
bic_both (float)
all_fits (dict[str, FitResult])
bic_bootstrap (dict[str, ndarray])
oiii_selected (str)
bic_oiii_off (float)
bic_oiii_broad1 (float)
bic_oiii_broad2 (float)
bic_oiii_both (float)
bic_oiii_bootstrap (dict[str, ndarray])
hei_selected (str)
bic_hei_off (float)
bic_hei_broad1 (float)
bic_hei_broad2 (float)
bic_hei_both (float)
bic_hei_bootstrap (dict[str, ndarray])

best_fit

The selected best-fit result.

Type:: FitResult

selected_model

Model name: "narrow", "broad1", "broad2", or "both".

Type:: str

bic_narrow

Median BIC for narrow-only model (from bootstrap when used).

Type:: float

bic_broad1

Median BIC for narrow + BROAD1 model (NaN if not attempted).

Type:: float

bic_broad2

Median BIC for narrow + BROAD2 model (NaN if not attempted).

Type:: float

bic_both

Median BIC for narrow + BROAD1 + BROAD2 model (NaN if not attempted).

Type:: float

all_fits

All fitted models for inspection.

Type:: dict[str, FitResult]

bic_bootstrap

Full BIC distributions from bootstrap iterations. Empty if bootstrap BIC was not run.

Type:: dict[str, np.ndarray]

all_fits: dict[str, FitResult]

best_fit: FitResult

bic_bootstrap: dict[str, ndarray]

bic_both: float

bic_broad1: float

bic_broad2: float

bic_hei_bootstrap: dict[str, ndarray]

bic_hei_both: float = nan

bic_hei_broad1: float = nan

bic_hei_broad2: float = nan

bic_hei_off: float = nan

bic_narrow: float

bic_oiii_bootstrap: dict[str, ndarray]

bic_oiii_both: float = nan

bic_oiii_broad1: float = nan

bic_oiii_broad2: float = nan

bic_oiii_off: float = nan

property chi2: float: Reduced chi-squared.

property constraints: Applied constraints.

property continuum: ndarray: Best-fit continuum.

property hei_broad_selected: bool

True if any HeI broad component was selected.

Type:: Convenience

hei_selected: str = 'off'

property line_names: list[str]: Ordered line names.

property lines: dict: Per-line results (delegates to best_fit.lines).

property model_flux: ndarray: Best-fit emission-line model.

property oiii_broad_selected: bool

True if any OIII broad component was selected.

Type:: Convenience

oiii_selected: str = 'off'

property params: ndarray: Best-fit parameter vector.

property residuals: ndarray: Fit residuals.

selected_model: str

property spectrum: Input spectrum.

property success: bool: Whether the optimiser converged.

class jwspecfit.DLAResult(log_NHI, log_NHI_err, beta_UV, beta_UV_err, log_F0, log_F0_err, Sigma_HI, samples, wave_fit, flux_fit, flux_err_fit, model_best, z=0.0, Av=0.0, log_evidence=0.0, _lya_subtracted=False, lya_params=None, x_HI=0.0, x_HI_err=(0.0, 0.0), fit_x_HI=False, igm_z_min=5.3, log_NHI_upper95=0.0, is_upper_limit=False, lam_pivot_A=None)[source]

Bases: object

Result of a DLA column density fit.

Parameters:

log_NHI (float)
log_NHI_err (tuple[float, float])
beta_UV (float)
beta_UV_err (tuple[float, float])
log_F0 (float)
log_F0_err (tuple[float, float])
Sigma_HI (float)
samples (dict[str, ndarray])
wave_fit (ndarray)
flux_fit (ndarray)
flux_err_fit (ndarray)
model_best (ndarray)
z (float)
Av (float)
log_evidence (float)
_lya_subtracted (bool)
lya_params (ndarray | None)
x_HI (float)
x_HI_err (tuple[float, float])
fit_x_HI (bool)
igm_z_min (float)
log_NHI_upper95 (float)
is_upper_limit (bool)
lam_pivot_A (float | None)

log_NHI

Median posterior log10(N_HI / cm^-2).

Type:: float

log_NHI_err

(lower, upper) 68% CI half-widths.

Type:: tuple of float

beta_UV

Median posterior UV spectral slope.

Type:: float

beta_UV_err

(lower, upper) 68% CI half-widths.

Type:: tuple of float

log_F0

Median posterior log10(F0) continuum normalisation.

Type:: float

log_F0_err

(lower, upper) 68% CI half-widths.

Type:: tuple of float

Sigma_HI

H I gas surface density in M_sun pc^-2 (from log_NHI).

Type:: float

samples

Full posterior samples: {“log_NHI”: arr, “beta_UV”: arr, “log_F0”: arr}.

Type:: dict

wave_fit

Wavelengths used in fit (after masking), observed frame.

Type:: np.ndarray

flux_fit

Dust-corrected fluxes used in fit.

Type:: np.ndarray

flux_err_fit

Dust-corrected flux errors used in fit.

Type:: np.ndarray

model_best

Best-fit model evaluated on wave_fit.

Type:: np.ndarray

z

Redshift used in fit.

Type:: float

Av

Dust correction applied.

Type:: float

log_evidence

Log-evidence (logZ) from dynesty.

Type:: float

Av: float = 0.0

Sigma_HI: float

beta_UV: float

beta_UV_err: tuple[float, float]

corner(**kwargs)[source]

Plot a corner plot of the posterior samples.

Requires the corner package.

Parameters:: **kwargs (Any) – Passed to corner.corner().
Returns:: The corner plot figure.
Return type:: matplotlib Figure

fit_x_HI: bool = False

flux_err_fit: ndarray

flux_fit: ndarray

igm_z_min: float = 5.3

is_upper_limit: bool = False

lam_pivot_A: float | None = None

log_F0: float

log_F0_err: tuple[float, float]

log_NHI: float

log_NHI_err: tuple[float, float]

log_NHI_upper95: float = 0.0

log_evidence: float = 0.0

lya_params: ndarray | None = None

model_best: ndarray

plot(ax=None, show_residuals=True, flux_unit='fnu', **kwargs)[source]

Plot the DLA fit over the data.

Parameters:

ax (matplotlib Axes, optional) – Axes to plot on. If None, creates a new figure. If show_residuals is True, this is ignored and a new figure with two panels is created.
show_residuals (bool) – If True, show a residual panel below the main plot.
flux_unit (str) – "fnu" for F_nu (default, same units as input) or "flam" for F_lambda (converted via F_lam = F_nu * c / lam^2).
**kwargs (Any) – Passed to the data plot (e.g. color, alpha).

Returns:

The figure object.

Return type:

matplotlib Figure

samples: dict[str, ndarray]

summary()[source]

Return a formatted summary string.

Return type:: str

wave_fit: ndarray

x_HI: float = 0.0

x_HI_err: tuple[float, float] = (0.0, 0.0)

z: float = 0.0

class jwspecfit.FitResult(lines, params, model_flux, continuum, residuals, chi2, spectrum, line_names=<factory>, constraints=None, success=True, lya_params=None)[source]

Bases: object

Container for a complete line-fit result.

Parameters:

lines (dict[str, LineResult])
params (ndarray)
model_flux (ndarray)
continuum (ndarray)
residuals (ndarray)
chi2 (float)
spectrum (Spectrum)
line_names (list[str])
constraints (ConstraintSet | None)
success (bool)
lya_params (ndarray | None)

lines

Per-line results, keyed by line name.

Type:: dict of LineResult

params

Full best-fit parameter vector.

Type:: np.ndarray

model_flux

Best-fit emission-line model (continuum-subtracted flux density).

Type:: np.ndarray

continuum

Best-fit continuum (µJy).

Type:: np.ndarray

residuals

Fit residuals (data − continuum − model), in µJy.

Type:: np.ndarray

chi2

Reduced χ² of the fit.

Type:: float

spectrum

Input spectrum.

Type:: Spectrum

line_names

Ordered line names.

Type:: list of str

constraints

Applied constraints.

Type:: ConstraintSet

success

Whether the optimiser converged.

Type:: bool

chi2: float

constraints: ConstraintSet | None = None

continuum: ndarray

flux_upper_limit(line_name, n_sigma=3.0)[source]

Compute a noise-based flux upper limit for a line.

Uses the local RMS of the residuals near the line position, multiplied by n_sigma and the line width.

Parameters:

line_name (str) – Line name.
n_sigma (float) – Number of sigma for the upper limit (default 3).

Returns:

Integrated flux upper limit in f_lam units, or None if there are insufficient pixels near the line.

Return type:

float or None

flux_upper_limits(line_names=None, n_sigma=3.0, snr_threshold=3.0)[source]

Compute noise-based upper limits for low-SNR lines.

Parameters:

line_names (list of str, optional) – Lines to check. If None, checks all fitted lines.
n_sigma (float) – Number of sigma for the upper limit (default 3).
snr_threshold (float) – Only compute upper limits for lines with SNR below this (default 3.0).

Returns:

dict of {str – {line_name: flux_upper_limit} for each line below the SNR threshold.

Return type:

float}

line_names: list[str]

lines: dict[str, LineResult]

lya_params: ndarray | None = None

model_flux: ndarray

params: ndarray

residuals: ndarray

spectrum: Spectrum

success: bool = True

class jwspecfit.LineResult(name, rest_wave_A, amplitude, centroid_A, sigma_A, flux, flux_err, ew_A, snr_int_err, snr_int_cont, snr_peak_err, snr_peak_cont)[source]

Bases: object

Fit result for a single emission line.

Parameters:

name (str)
rest_wave_A (float)
amplitude (float)
centroid_A (float)
sigma_A (float)
flux (float)
flux_err (float)
ew_A (float)
snr_int_err (float)
snr_int_cont (float)
snr_peak_err (float)
snr_peak_cont (float)

name

Line name.

Type:: str

rest_wave_A

Rest-frame wavelength (Å).

Type:: float

amplitude

Best-fit amplitude (flux × Å, in flux-density units × Å).

Type:: float

centroid_A

Best-fit observed centroid (Å).

Type:: float

sigma_A

Best-fit observed Gaussian σ (Å).

Type:: float

flux

Integrated line flux (amplitude, since profile is area-normalised).

Type:: float

flux_err

Bootstrap uncertainty on flux.

Type:: float

ew_A

Rest-frame equivalent width (Å).

Type:: float

snr_int_err

Integrated SNR: flux / bootstrap flux error.

Type:: float

snr_int_cont

Integrated SNR: flux / (local continuum RMS × √N_pix).

Type:: float

snr_peak_err

Peak SNR: model peak / uncertainty array at peak pixel.

Type:: float

snr_peak_cont

Peak SNR: model peak / local continuum RMS.

Type:: float

amplitude: float

centroid_A: float

ew_A: float

flux: float

flux_err: float

name: str

rest_wave_A: float

sigma_A: float

property snr: float: Backward-compatible alias for snr_int_err.

snr_int_cont: float

snr_int_err: float

snr_peak_cont: float

snr_peak_err: float

upper_limit(n_sigma=3.0)[source]

Return n-sigma flux upper limit.

Return type:: float
Parameters:: n_sigma (float)

class jwspecfit.Peak(z, prob, dchi2, ci68, ci95, n_lines_used, lines_used)[source]

Bases: object

A single local minimum of the chi^2(z) curve, refined and ranked.

Parameters:

z (float)
prob (float)
dchi2 (float)
ci68 (tuple[float, float])
ci95 (tuple[float, float])
n_lines_used (int)
lines_used (list[str])

ci68: tuple[float, float]

ci95: tuple[float, float]

dchi2: float

lines_used: list[str]

n_lines_used: int

prob: float

z: float

jwspecfit.R_from_pixels(lam_um)[source]

Estimate resolving power from the pixel spacing.

Assumes the spectrum is Nyquist-sampled, so that the FWHM of the line-spread function spans ~2 pixels: R ≈ λ / (2 Δλ).

Returns a callable R(lam_um) that can be passed directly to resolve_R(), sigma_inst_A(), or fit_lines().

This is a rough estimate and should only be used as a fallback when neither a grating name nor an explicit R is available.

Parameters:: lam_um (array_like) – Wavelength in microns (must be sorted).
Returns:: Function R(lam_um) -> np.ndarray that interpolates the estimated resolving power.
Return type:: callable

jwspecfit.R_prism(lam_um, post_launch=True)[source]

Resolving power R(λ) for NIRSpec PRISM/CLEAR (Jakobsen+22 Fig. 6).

Tabulated values from Jakobsen et al. 2022 (A&A 661, A80, Fig. 6 — pre-launch instrument-model curve), interpolated linearly in λ. When post_launch=True (default), the curve is multiplied by 1.3, the correction factor adopted by Pollock+26 / de Graaff+25 to reflect the post-launch on-orbit performance.

Parameters:

lam_um (array_like) – Wavelength in microns.
post_launch (bool) – If True (default), apply the 1.3× post-launch correction.

Returns:

Resolving power R at each wavelength.

Return type:

np.ndarray

class jwspecfit.RedshiftResult(z_best, z_ci68, z_ci95, peaks, is_decisive, z_grid_coarse, chi2_coarse, P_z_coarse, z_grid_fine, chi2_fine, P_z_fine, lines_used, grating, spec)[source]

Bases: object

Output of fit_redshift().

Parameters:

z_best (float)
z_ci68 (tuple[float, float])
z_ci95 (tuple[float, float])
peaks (list[Peak])
is_decisive (bool)
z_grid_coarse (ndarray)
chi2_coarse (ndarray)
P_z_coarse (ndarray)
z_grid_fine (ndarray)
chi2_fine (ndarray)
P_z_fine (ndarray)
lines_used (list[str])
grating (str | None)
spec (Spectrum)

P_z_coarse: ndarray

P_z_fine: ndarray

chi2_coarse: ndarray

chi2_fine: ndarray

grating: str | None

is_decisive: bool

lines_used: list[str]

peaks: list[Peak]

plot()[source]: Two-panel diagnostic: Delta chi^2(z) and the spectrum at z_best.

spec: Spectrum

z_best: float

z_ci68: tuple[float, float]

z_ci95: tuple[float, float]

z_grid_coarse: ndarray

z_grid_fine: ndarray

class jwspecfit.Spectrum(wave_um, flux_ujy, err_ujy, grating=None, z=None, R=None, meta=<factory>, sci_2d=None)[source]

Bases: object

Container for a 1-D spectrum.

Parameters:

wave_um (ndarray)
flux_ujy (ndarray)
err_ujy (ndarray)
grating (str | None)
z (float | None)
R (float | None)
meta (dict[str, Any])
sci_2d (ndarray | None)

wave_um

Observed wavelength in microns.

Type:: np.ndarray

flux_ujy

Flux density in micro-Jansky.

Type:: np.ndarray

err_ujy

1-sigma uncertainty in micro-Jansky.

Type:: np.ndarray

grating

Grating name (e.g. "PRISM", "G395M"). None for stacked spectra.

Type:: str or None

z

Source redshift (set by user, not from header).

Type:: float or None

R

Spectral resolving power. Overrides grating when set (useful for stacks).

Type:: float or callable or None

meta

Arbitrary metadata from the FITS header or user.

Type:: dict

sci_2d

Optional 2-D rectified spectrum image with shape (n_spatial, n_pix) where n_pix == len(wave_um). Populated by read_fits() when an SCI ImageHDU is present and its wavelength axis matches the 1-D extraction; otherwise None. Read by plot_spectrum_interactive() to render a 2-D + 1-D stacked view; ignored by every other code path.

Type:: np.ndarray or None

R: float | None = None

copy()[source]

Return a shallow copy with copied arrays.

Return type:: Spectrum

property dlam_A: ndarray: Pixel widths in Angstroms.

property err_flam: ndarray: Error in erg/s/cm²/Å.

err_ujy: ndarray

property flux_flam: ndarray: Flux density in erg/s/cm²/Å.

flux_ujy: ndarray

grating: str | None = None

mask_valid()[source]

Boolean mask: True where flux and error are finite and err > 0.

Return type:: ndarray

meta: dict[str, Any]

property n_pix: int

sci_2d: ndarray | None = None

property wave_A: ndarray: Wavelength in Angstroms.

property wave_edges_A: ndarray: Pixel-edge wavelengths in Angstroms (length n_pix + 1).

wave_um: ndarray

z: float | None = None

jwspecfit.compute_D_Lya(wave_A, flux, flux_err, z, *, cont_range_A=(1250.0, 2600.0), int_range_A=(1180.0, 1350.0), mask_lines=True, mask_width_A=20.0, n_mc=1000, seed=42)[source]

Heintz+25 D_Lyα equivalent-width statistic (Eq. 1).

Fits a power-law continuum F_cont(λ) = F_1550 (λ/1550)^β over cont_range_A (lines masked) then integrates

D_Lyα = ∫ (1 − F_λ / F_cont(λ)) dλ / (1 + z)

over int_range_A (rest-frame). D_Lyα < 0 indicates a Lyα emitter; ≳50 Å indicates a strong damped absorber.

Parameters:

wave_A (np.ndarray) – Spectrum in observed-frame Å and arbitrary linear flux units.
flux (np.ndarray) – Spectrum in observed-frame Å and arbitrary linear flux units.
flux_err (np.ndarray) – Spectrum in observed-frame Å and arbitrary linear flux units.
z (float) – Source redshift.
cont_range_A (tuple) – Rest-frame range used to fit the power-law continuum (default 1250–2600 Å, matching Heintz+25).
int_range_A (tuple) – Rest-frame range over which the EW integral is taken (default 1180–1350 Å, matching Heintz+25 Eq. 1).
mask_lines (bool) – Whether to mask known UV emission lines from the continuum fit (default True). The integral is always taken over the full int_range_A — emission/absorption inside it is the signal.
mask_width_A (float) – Half-width of line masks for the continuum fit (rest-frame Å).
n_mc (int) – Monte-Carlo iterations for D_Lyα uncertainty propagation.
seed (int) – RNG seed.

Returns:

{"D_Lya": float, "D_Lya_err": float, "beta_UV": float, "beta_UV_err": float, "F_1550": float, "F_1550_err": float}

Return type:

dict

jwspecfit.export_lines_txt(result, path, z=None)[source]

Export per-line measurements to a text file.

Columns: name, rest_wave_A, centroid_A, flux, flux_err, ew_A, sigma_v_kms, snr_integrated, snr_peak.

Parameters:

result (FitResult) – Fit result.
path (str or Path) – Output text file path.
z (float, optional) – Redshift (for velocity calculation). If None, uses result.spectrum.z.

Return type:

None

jwspecfit.fit_NHI(wave_A, flux, flux_err, z=0.0, *, fit_x_HI=False, igm_z_min=5.3, b_kms=30.0, continuum=None, Av=0.0, dust_law='cardelli', Rv=3.1, R=None, mask_lines=True, mask_width_A=20.0, mask_lya_emission_width_A=8.0, fit_range_A=(1216.0, 3000.0), mask_regions_A=None, fit_lya=False, lya_params=None, emission_lines=None, lam_pivot_A=None, prior_log_NHI=(18.0, 24.0), prior_x_HI=(0.0, 1.0), prior_beta_UV=(-4.0, 0.0), prior_log_F0=None, prior_log_A_lya=(-22.0, -15.0), prior_sigma_lya=(1.0, 30.0), prior_alpha_lya=(0.0, 5.0), n_live=500, dlogz=0.5, seed=42)[source]

Joint Bayesian fit of the DLA + IGM damping wing (Pollock+26-style).

Implements the Pollock et al. (2026) joint nested-sampling fit over (log_F0, β_UV, log_NHI[, x_HI[, Lyα]]) with an exact Voigt-Hjerting profile and the Miralda-Escudé (1998) IGM damping wing in the Totani+06 closed form.

Parameters:

wave_A (np.ndarray) – Wavelength (observed-frame Å) and flux density arrays.
flux (np.ndarray) – Wavelength (observed-frame Å) and flux density arrays.
flux_err (np.ndarray) – Wavelength (observed-frame Å) and flux density arrays.
z (float) – Source redshift (use 0 for rest-frame stacks).
fit_x_HI (bool) – If True, sample over IGM neutral fraction x_HI in addition to the local DLA. Recommended at z > 6.
igm_z_min (float) – Lower-redshift bound of the IGM integration (default 5.3, Bosman+22 reionisation end).
b_kms (float) – Doppler parameter (km/s) for the Voigt profile. The damping wing is insensitive to b for log N_HI > 20.3.
continuum (np.ndarray, optional) – Pre-smoothed continuum estimate. When provided the model is fit against this smoothed array instead of the raw flux — emission lines are then implicitly removed.
Av (float, str, float) – Optional foreground dust correction applied before fitting.
dust_law (float, str, float) – Optional foreground dust correction applied before fitting.
Rv (float, str, float) – Optional foreground dust correction applied before fitting.
R (float or None) – Spectral resolving power for LSF convolution.
mask_lines (bool) – If True (default), mask known UV emission lines. Lyα, the NV λλ1238.8, 1242.8 doublet, and any abs_* low-ionisation absorption features are excluded from this mask: Lyα and NV sit inside the damping wing (1216–1263 Å rest) where the N_HI constraint lives, and any DLA strong enough to fit absorbs NV at τ ≫ 1 (no emission left to contaminate); abs_* features are foreground-gas signatures aligned with the DLA.
mask_width_A (float) – Half-width of the general emission-line mask (rest-frame Å). Default 20 Å is appropriate for NIRSpec PRISM resolution.
mask_lya_emission_width_A (float) – Half-width of the narrow Lyα-emission core mask in rest-frame Å (default 8 Å). Set this small enough not to eat the damping wing. Ignored when fit_lya=True.
fit_range_A (tuple) – Rest-frame fit window (default 1216–3000 Å, matching Pollock+26 — Lyα to just below the Balmer break).
mask_regions_A (list of (float, float), optional) – Extra rest-frame wavelength ranges to ignore.
fit_lya (bool) – Jointly fit a skewed Lyα emission profile.
lya_params (array of length 4, optional) – Pre-determined Lyα profile (A_peak, μ, σ, α) added to the intrinsic spectrum before attenuation.
prior_log_NHI (tuple) – Uniform-prior bounds for the named parameter (Pollock+26 Sect. 3.2 defaults: log_NHI ∈ [18,24], x_HI ∈ [0,1], β_UV ∈ [-4,0]).
prior_x_HI (tuple) – Uniform-prior bounds for the named parameter (Pollock+26 Sect. 3.2 defaults: log_NHI ∈ [18,24], x_HI ∈ [0,1], β_UV ∈ [-4,0]).
prior_beta_UV (tuple) – Uniform-prior bounds for the named parameter (Pollock+26 Sect. 3.2 defaults: log_NHI ∈ [18,24], x_HI ∈ [0,1], β_UV ∈ [-4,0]).
prior_log_F0 (tuple, optional) – Uniform prior bounds on log10(F0) in the same flux units as the input data. Defaults to ±3 dex around the median continuum flux derived from the input spectrum.
prior_log_A_lya (tuple) – Priors used when fit_lya=True.
prior_sigma_lya (tuple) – Priors used when fit_lya=True.
prior_alpha_lya (tuple) – Priors used when fit_lya=True.
n_live (int) – Number of dynesty live points (default 500).
dlogz (float) – dynesty convergence threshold on log Z (default 0.5).
seed (int) – RNG seed.
emission_lines (list[tuple[float, float, float]] | None)
lam_pivot_A (float | None)

Returns:

Posterior samples for every free parameter, median ± 16/84 percentiles, log-evidence, and a 95th-percentile upper limit on log N_HI when the posterior is unconstrained.

Return type:

DLAResult

Notes

Pollock+26 Eq. 4: F = (F0·(λ/λ_pivot)^β + Lyα) · exp(−τ_DLA) · exp(−τ_IGM).
Lyα emission is added before attenuation (Pollock+26 convention).
Reports an upper limit when the posterior touches the prior lower bound (Pollock+26 do this for 21/48 sources).

jwspecfit.fit_lines(spectrum, z, *, grating=None, R=None, lines=None, wave_range_A=None, wave_windows_A=None, deg=2, n_boot=1000, clip_sigma=2.5, n_jobs=-1, save_path=None, fit_balmer_broad=False, fit_oiii_broad=False, fit_hei_broad=False, n_boot_bic=100, snr_threshold=5.0, oiii_snr_threshold=5.0, hei_snr_threshold=5.0, bic_delta=6.0, sigma_factor=1.0, centroid_vmax=500.0, centroid_max_sigma=1.0, moving_average=False, tie_uv_doublets=True, tie_uv_centroids=True, tie_uv_widths=True, sigma_overrides=None, centroid_overrides=None, lya_break=False, niv_doublet_ratio=None, ciii_doublet_ratio=None)[source]

Fit emission lines in a spectrum.

Narrow-only fit by default. Set any flag to opt in to a BIC-based broad component selection (via fit_with_broad()):

fit_balmer_broad=True — Balmer broad component selection.
fit_oiii_broad=True — independent [OIII] outflow selection.
fit_hei_broad=True — independent He I broad selection (kinematics shared across HeI lines within each tier).

The tests are independent — enable any combination.

Parameters:

spectrum (Spectrum) – Input spectrum.
z (float) – Source redshift.
grating (str, optional) – Grating name.
R (float or callable, optional) – Resolving power.
lines (list of str, optional) – Lines to fit.
wave_range_A (tuple, optional) – Observed wavelength range (Angstrom). Mutually exclusive with wave_windows_A.
wave_windows_A (list of (float, float), optional) – Fit the spectrum in multiple independent wavelength windows, each with its own continuum subtraction. Useful for stacked spectra where the continuum shape varies across the full wavelength range (e.g. UV-normalised stacks). Each tuple is (lo, hi) in observed-frame Angstroms. Results from all windows are merged into a single FitResult. Mutually exclusive with wave_range_A.
deg (int) – Continuum polynomial degree (default 2).
n_boot (int) – Bootstrap iterations for flux uncertainties (default 1000).
clip_sigma (float) – Continuum sigma-clipping threshold (default 2.5).
n_jobs (int) – Parallel jobs for bootstrap (default -1).
save_path (str or Path, optional) – Path to save the result.
fit_balmer_broad (bool) – If True, run BIC selection for Balmer broad components. Default False (narrow-only).
fit_oiii_broad (bool) – If True, run independent BIC test for [OIII] outflow broad components. Default False (narrow-only).
fit_hei_broad (bool) – If True, run independent BIC test for He I broad components (shared kinematics across HeI lines). Default False (narrow-only).
n_boot_bic (int) – Bootstrap iterations for BIC model selection (default 100). Only used when at least one broad flag is True.
snr_threshold (float) – Minimum Hα SNR to attempt Balmer broad fitting (default 5.0).
oiii_snr_threshold (float) – Minimum [OIII] 5007/4959 SNR to attempt OIII broad fitting (default 5.0).
bic_delta (float) – ΔBIC threshold for model selection (default 6.0).
sigma_factor (float) – Multiplicative factor on the upper line-width bound. Use values > 1 for stacked spectra (default 1.0).
centroid_vmax (float) – Maximum centroid offset in km/s (default 500). Sets the velocity-space ceiling on how far any line’s centroid can wander from systemic. Increase for stacked spectra with larger velocity offsets between lines.
centroid_max_sigma (float) – Resolution-aware cap on narrow-line centroids: each narrow line’s centroid can drift at most centroid_max_sigma × σ_inst from systemic, where σ_inst is the instrumental Gaussian σ evaluated at the line. The effective bound is the tighter of this and the velocity cap. Broad components (_BROAD/_BROAD2) bypass this so real outflow blueshifts aren’t clipped. Default 1.0.
moving_average (bool or int) – If False (default), use polynomial continuum. If True, use a median filter with a default window of 75 pixels. If an int, use that as the median-filter window size.
tie_uv_doublets (bool) – Tie UV doublet kinematics and fix resonance-line flux ratios. Recommended for stacked spectra where doublets are poorly resolved (default True).
hei_snr_threshold (float)
tie_uv_centroids (bool)
tie_uv_widths (bool)
sigma_overrides (dict[str, tuple[float, float]] | None)
centroid_overrides (dict[str, tuple[float, float]] | None)
lya_break (bool)
niv_doublet_ratio (float | None)
ciii_doublet_ratio (float | None)

Return type:

FitResult | BroadFitResult

Returns:

BroadFitResult – When at least one broad flag is True and wave_windows_A is None. Delegates all FitResult attributes (lines, params, model_flux, etc.) via properties, so it can be used as a drop-in replacement.
FitResult – When both broad flags are False, or when wave_windows_A is specified (merged result).

jwspecfit.fit_redshift(spec, *, lines=None, z_min=0.0, z_max=20.0, dz_coarse=None, prior=None, n_peaks=5, refine_dchi2=9.0, continuum_order=3, nonneg=True, min_lines_in_range=2, verbose=True)[source]

Fit the redshift of spec by joint strong-line fitting on a z grid.

Parameters:

spec (Spectrum) – Input spectrum. Must have valid wave_um, flux_ujy, err_ujy.
lines (list of str, optional) – Line names (keys of REST_LINES_A) to include. Defaults to DEFAULT_LINES.
z_min (float) – Bounds of the redshift search (default 0 to 20).
z_max (float) – Bounds of the redshift search (default 0 to 20).
dz_coarse (float, optional) – Coarse-grid step. Defaults to a value tied to the spectral resolving power (see _default_dz_coarse()).
prior (None | (lo, hi) | list[(lo, hi)] | callable) – Optional non-Gaussian prior over z. None means flat.
n_peaks (int) – Minimum number of peaks to refine and return.
refine_dchi2 (float) – All local minima within this Δχ² of the global minimum are refined, in addition to the top n_peaks.
continuum_order (int) – Order of the Legendre-polynomial continuum that is added to the per-z linear fit. Default 3.
nonneg (bool) – If True, line amplitudes are constrained to be non-negative (emission only). Strongly recommended.
min_lines_in_range (int) – Trial redshifts where fewer than this many default lines fall inside the spectrum get the continuum-only chi^2 baseline.
verbose (bool) – If True, log a few timing milestones at INFO level.

Return type:

RedshiftResult

jwspecfit.fit_with_broad(spectrum, z, *, grating=None, R=None, lines=None, deg=2, fit_balmer_broad=False, fit_oiii_broad=False, fit_hei_broad=False, n_boot=1000, n_boot_bic=100, n_jobs=-1, snr_threshold=5.0, oiii_snr_threshold=5.0, hei_snr_threshold=5.0, bic_delta=6.0, sigma_factor=1.0, centroid_vmax=500.0, centroid_max_sigma=1.0, moving_average=False, tie_uv_doublets=True, tie_uv_centroids=True, tie_uv_widths=True, sigma_overrides=None, centroid_overrides=None, niv_doublet_ratio=None, ciii_doublet_ratio=None, _print_R=True)[source]

Fit emission lines with optional BIC-selected broad components.

Two independent BIC tests can be enabled:

Balmer broad (fit_balmer_broad): tests narrow vs. narrow + intermediate broad (FWHM 500–2000 km/s) vs. narrow + very broad (FWHM 2000–5000 km/s) vs. narrow + both. Gated by snr_threshold on Hα.
[OIII] broad (fit_oiii_broad): tests baseline vs. baseline + broad [OIII] 5007/4959 (outflow signature). Gated by oiii_snr_threshold on [OIII] 5007 / 4959.

The two tests are independent — both can be selected, only one, or neither.

Parameters:

spectrum (Spectrum) – Input spectrum.
z (float) – Source redshift.
grating (str, optional) – Grating name.
R (float or callable, optional) – Resolving power.
lines (list of str, optional) – Narrow line list. If None, auto-detected.
deg (int) – Continuum polynomial degree.
fit_balmer_broad (bool) – If True, run BIC model selection for a broad Balmer component (intermediate / very-broad / both). Default False — opt in only when you want to test for a BLR or outflow signature on Hα/Hβ.
fit_oiii_broad (bool) – If True, run an independent BIC test for a broad component on [OIII] 5007/4959. Decoupled from the Balmer test. Default False — opt in only when looking for outflow signatures.
fit_hei_broad (bool) – If True, run an independent BIC test for a broad component on all observable He I lines (5877/6680/4472/…). Within each tier all HeI broads share kinematics (anchored on first present). Default False.
n_boot (int) – Number of bootstrap iterations for flux uncertainties (default 1000).
n_boot_bic (int) – Number of bootstrap iterations for BIC model selection (default 100). Set to 0 for single-point BIC (legacy behaviour).
n_jobs (int) – Number of parallel jobs for bootstrap. -1 uses all cores, 1 runs sequentially. Default -1.
snr_threshold (float) – Minimum Hα SNR to attempt the Balmer broad fit (default 5.0).
oiii_snr_threshold (float) – Minimum [OIII] 5007 / 4959 SNR to attempt the OIII broad fit (default 5.0).
bic_delta (float) – ΔBIC threshold for accepting a more complex model (default 6.0).
hei_snr_threshold (float)
sigma_factor (float)
centroid_vmax (float)
centroid_max_sigma (float)
moving_average (bool | int)
tie_uv_doublets (bool)
tie_uv_centroids (bool)
tie_uv_widths (bool)
sigma_overrides (dict[str, tuple[float, float]] | None)
centroid_overrides (dict[str, tuple[float, float]] | None)
niv_doublet_ratio (float | None)
ciii_doublet_ratio (float | None)
_print_R (bool)

Return type:

BroadFitResult

jwspecfit.get_line_list(grating='prism')[source]

Return default line names for a given grating.

For prism (R ~ 30–300), merged doublet entries are used for pairs that cannot be resolved. For medium and high-resolution gratings (R ≥ 1000), individual doublet components are returned.

The user can always override by passing an explicit lines= argument to fit_lines().

Parameters:: grating (str) – One of "prism", "medium" / "g140m" / "g235m" / "g395m", "high" / "g140h" / "g235h" / "g395h", or "grating" (generic resolved mode).
Returns:: Line names present in REST_LINES_A.
Return type:: list of str

jwspecfit.igm_transmission(wave_obs_A, z_source)[source]

Mean IGM transmission using Inoue et al. (2014).

Computes the average transmission through the intergalactic medium for photons observed at wave_obs_A from a source at z_source, accounting for Lyman-series line absorption (LAF) and damped Lyman-alpha systems (DLA).

Parameters:

wave_obs_A (np.ndarray) – Observed wavelength in Angstroms.
z_source (float) – Source redshift.

Returns:

Transmission T(λ) ∈ [0, 1] at each wavelength.

Return type:

np.ndarray

jwspecfit.load_result(path)[source]

Load a FitResult from a .npz file saved by save_result().

Parameters:: path (str or Path) – Path to the .npz file.
Return type:: FitResult

jwspecfit.lya_model(lam_left_A, lam_right_A, z, amplitude, mu_A, sigma_A, skew)[source]

Lyman-alpha emission model: skewed Gaussian × IGM transmission.

Parameters:

lam_left_A (np.ndarray) – Bin edges (Å).
lam_right_A (np.ndarray) – Bin edges (Å).
z (float) – Source redshift.
amplitude (float) – Integrated Lyα flux (before IGM absorption).
mu_A (float) – Observed Lyα centroid (Å).
sigma_A (float) – Observed Lyα width (Å).
skew (float) – Skewness parameter (positive = red-asymmetric, typical for Lyα).

Returns:

Model flux density in each bin (Å⁻¹ × flux units).

Return type:

np.ndarray

jwspecfit.observable_lines(line_names, z, wave_min_um, wave_max_um, *, margin_sigma=3.0, sigma_um=0.005)[source]

Filter lines to those observable in the wavelength range.

Parameters:

line_names (list of str) – Candidate line names (keys of REST_LINES_A).
z (float) – Source redshift.
wave_min_um (float) – Observed wavelength range in microns.
wave_max_um (float) – Observed wavelength range in microns.
margin_sigma (float) – Number of sigma margin from the edges.
sigma_um (float) – Approximate line width in microns (for margin calculation).

Returns:

Lines whose observed wavelength falls within the range.

Return type:

list of str

jwspecfit.plot_2d_1d(path, z=None, *, sci_ext='SCI', hdu=None, wave_col=None, flux_col=None, err_col=None, flux_scale=1000.0, flux_label='$F_\\\\nu$ [nJy]', xlabel='$\\\\lambda_{\\\\rm obs}\\\\,[\\\\mu{\\\\rm m}]$', cmap='Blues', vmin_pct=5.0, vmax_pct=99.5, y_crop=(0.25, 0.75), xlim=None, ylim=None, line_colour='steelblue', line_width=0.5, err_colour='grey', err_alpha=0.35, figsize=(8, 4), dpi=300, title=None, lines=None, add_lines=None, height_ratios=(0.35, 1.0))[source]

Plot a JWST/NIRSpec 2D + 1D spectrum from a .spec.fits file.

Reads the 2D image from sci_ext and the 1D extraction from SPEC1D via jwspecfit.read_fits(), then renders a stacked figure (2D pcolormesh above, 1D flux below) with optional emission- line markers at the supplied redshift.

Parameters:

path (str | Path) – Path to the .spec.fits file.
z (float | None) – Redshift used to place observed-frame line markers. If None, no markers are drawn.
sci_ext (str) – Name of the 2D image extension (default "SCI").
hdu (str | int | None) – Forwarded to read_fits() for 1D extraction (override the default SPEC1D / column auto-detection).
wave_col (str | None) – Forwarded to read_fits() for 1D extraction (override the default SPEC1D / column auto-detection).
flux_col (str | None) – Forwarded to read_fits() for 1D extraction (override the default SPEC1D / column auto-detection).
err_col (str | None) – Forwarded to read_fits() for 1D extraction (override the default SPEC1D / column auto-detection).
flux_scale (float) – Multiplier and y-axis label applied to the 1D flux (default converts µJy → nJy).
flux_label (str) – Multiplier and y-axis label applied to the 1D flux (default converts µJy → nJy).
xlabel (str) – Shared x-axis label (default observed wavelength in µm).
cmap (str) – Colormap for the 2D panel.
vmin_pct (float) – Percentile clip for the 2D colour scale.
vmax_pct (float) – Percentile clip for the 2D colour scale.
y_crop (tuple[float, float]) – (y_lo_frac, y_hi_frac) fractional crop of the 2D vertical extent (default keeps the middle 50 % rows).
xlim (tuple[float, float] | None) – Optional axis limits. xlim applies to both panels.
ylim (tuple[float, float] | None) – Optional axis limits. xlim applies to both panels.
line_colour (str) – Colour and line width of the 1D flux trace.
line_width (float) – Colour and line width of the 1D flux trace.
err_colour (str) – Fill colour and alpha for the ±1σ error band.
err_alpha (float) – Fill colour and alpha for the ±1σ error band.
figsize (tuple[float, float]) – Matplotlib figure size, resolution (default 300), and 2D-vs-1D height ratio.
dpi (float) – Matplotlib figure size, resolution (default 300), and 2D-vs-1D height ratio.
height_ratios (tuple[float, float]) – Matplotlib figure size, resolution (default 300), and 2D-vs-1D height ratio.
title (str | None) – Optional title placed above the 2D panel.
lines (list[str] | bool | None) – Default-marker control: None uses the package defaults; a list of REST_LINES_A keys restricts to those names; False disables defaults entirely.
add_lines (list[str] | dict[str, float] | None) – Extra markers: either REST_LINES_A keys, or a {label: rest_wavelength_A} dict.

Returns:

Matplotlib figure and the two axes.

Return type:

fig, (ax_2d, ax_1d)

jwspecfit.plot_fit(result, *, fig=None, wave_unit='A', flux_unit='fnu', show_residuals=True, show_components=True, label_lines=True, y_pad=1.3, exclude_wave_A=None, rest_frame=False, save_path=None)[source]

Plot a spectral fit with data, model, continuum, and residuals.

The y-axis upper limit is set to the peak of the tallest emission line (above continuum) times y_pad, so the plot is scaled to the lines rather than noise spikes.

Parameters:

result (FitResult) – Output of fit_lines().
fig (Figure, optional) – Matplotlib figure to draw on. If None, creates a new one.
wave_unit (str) – "A" for Angstroms (default) or "um" for microns.
flux_unit (str) – "fnu" for µJy (default) or "flam" for erg/s/cm²/Å.
show_residuals (bool) – Show residual panel below the main plot (default True).
show_components (bool) – Show individual Gaussian components as filled curves (default True). Broad components are drawn with hatching for clarity.
label_lines (bool) – Annotate line identifications (default True).
y_pad (float) – Multiplicative padding above the tallest line peak (default 1.3).
exclude_wave_A (list of (float, float), optional) – Wavelength ranges in Angstroms to hide from the plot. Each tuple is (lo, hi). Useful for masking noisy detector regions.
rest_frame (bool) – If True, plot wavelengths in the rest frame by dividing by (1 + z) using the redshift stored in the spectrum. Default False (observed frame).
save_path (str, optional) – If given, save the figure to this file path (e.g. "fit.pdf").

Returns:

The matplotlib figure.

Return type:

Figure

jwspecfit.plot_fit_interactive(result, *, wave_unit='A', flux_unit='fnu', show_components=True, show_residuals=True, y_pad=1.3, exclude_wave_A=None, rest_frame=False, z=None, lines=False, add_lines=None, line_color='darkred', show_2d='auto', cmap_2d='Blues', vmin_pct=5.0, vmax_pct=99.5, y_crop=(0.25, 0.75))[source]

Interactive plotly plot of a spectral fit with zoom and hover.

Optionally stacks a 2-D rectified-spectrum panel above the main fit panel (when result.spectrum.sci_2d is populated by read_fits()) and a residual panel below. By default no curated emission-line markers are drawn — every fitted line is already labelled at its peak above the model. Pass lines=None to add the package-default markers, or an explicit list of jwspecfit.lines.REST_LINES_A keys to mark only those.

Parameters:

result (FitResult) – Output of fit_lines().
wave_unit (str) – "A" for Angstroms (default) or "um" for microns.
flux_unit (str) – "fnu" for µJy (default) or "flam" for erg/s/cm²/Å.
show_components (bool) – Show individual line components (default True).
show_residuals (bool) – Show residual panel below the main plot (default True).
y_pad (float) – Multiplicative padding above tallest line (default 1.3).
exclude_wave_A (list of (float, float), optional) – Wavelength ranges in Angstroms to hide from the plot.
rest_frame (bool) – If True, plot wavelengths in the rest frame by dividing by (1 + z). Default False (observed frame). Emission-line markers are placed at rest wavelengths when True.
z (float, optional) – Redshift override used for rest-frame conversion and marker placement. When None (default), result.spectrum.z is used. Raises ValueError if rest_frame=True and neither source provides a redshift.
lines (sequence of str, bool, or None) – Curated emission-line markers (vertical dashed lines + labels). False (default) draws none — the fit-component peak annotations already identify every fitted line. None opts in to the package-default marker set; an explicit list of jwspecfit.lines.REST_LINES_A keys marks only those names.
add_lines (dict[str, float] or sequence of str, optional) – Extra markers to overlay on top of lines. Same semantics as in plot_spectrum_interactive().
line_color (str) – Colour for the emission-line markers and their labels (default "darkred").
show_2d (bool or {"auto"}) – Whether to stack the 2-D rectified spectrum image above the fit panel. "auto" (default) shows it iff result.spectrum.sci_2d is populated. True forces it (silently no-ops when no 2-D is available); False always hides it.
cmap_2d (str) – Plotly colorscale for the 2-D panel (default "Blues").
vmin_pct (float) – Percentile clip for the 2-D colour scale (default 5 / 99.5).
vmax_pct (float) – Percentile clip for the 2-D colour scale (default 5 / 99.5).
y_crop ((float, float)) – Fractional crop (lo, hi) of the spatial axis on the 2-D panel (default (0.25, 0.75)).

Return type:

plotly.graph_objects.Figure

jwspecfit.plot_spectrum_interactive(source, *, z=None, wave_unit='A', flux_unit='fnu', rest_frame=False, exclude_wave_A=None, title=None, labels=None, lines=None, add_lines=None, line_color='darkred', show_zero=True, show_2d='auto', cmap_2d='Blues', vmin_pct=5.0, vmax_pct=99.5, y_crop=(0.25, 0.75), **read_kwargs)[source]

Open and interactively plot one or more 1-D spectra.

Accepts a Spectrum object, a path to a .fits / .npz file, or a list / tuple of such items. When given a path, the file is read via read_fits() (or read_npz() for .npz) and any extra read_kwargs are forwarded to the reader (e.g. hdu=, wave_col= for FITS overrides).

Parameters:

source (Spectrum, str, Path, or sequence of these) – A single spectrum / file path, or a list / tuple of them to overplot.
z (float, optional) – Source redshift. Used only when a source item is a path; for an existing Spectrum, its own z is preserved. Forwarded to every reader call.
wave_unit (str) – "A" for Angstroms (default) or "um" for microns.
flux_unit (str) – "fnu" for µJy (default) or "flam" for erg/s/cm²/Å.
rest_frame (bool) – If True and a spectrum has a redshift, divide wavelengths by (1 + z). Default False (observed frame). Applied per spectrum.
exclude_wave_A (list of (float, float), optional) – Wavelength ranges in Angstroms to hide from the plot.
title (str, optional) – Figure title. When a single spectrum is supplied, defaults to filename + redshift + grating if available; when multiple spectra are supplied, defaults to no title.
labels (str or sequence of str, optional) – Legend label(s). When None (default), a single spectrum uses "Data" and multiple spectra use each spectrum’s filename (or "Spectrum {i}" as a fallback).
lines (sequence of str, bool, or None) – Emission lines to mark as vertical dashed lines at the supplied redshift. None (default) draws a curated list of common UV/optical lines. Pass an explicit list of keys from jwspecfit.lines.REST_LINES_A to override, or False to disable. The effective redshift is taken from z if given, else from a single spectrum’s own spec.z. In rest-frame mode the markers sit at the rest wavelengths.
add_lines (dict[str, float] or sequence of str, optional) –
Extra lines to overlay on top of lines. Two accepted forms:
- dict — {label: rest_wavelength_A}. Free-form labels with explicit rest-frame wavelengths in Angstroms. Use this for lines not in jwspecfit.lines.REST_LINES_A, e.g. add_lines={"Mg II 2796": 2796.352}.
- list of str — names from REST_LINES_A (e.g. add_lines=["H8", "HEPSILON", "FeII_2382"]). The rest wavelength is looked up automatically. Call jwspecfit.show_lines() to see what’s available.
Each entry is redshifted by (1 + z) and staggered alongside the default markers.
line_color (str) – Colour for the emission-line markers and their labels (default "darkred").
show_zero (bool) – Draw a light-grey dashed horizontal line at y = 0 to make continuum detection easier to gauge by eye (default True).
show_2d (bool or {"auto"}) – Whether to render the 2-D rectified spectrum image above the 1-D panel. "auto" (default) shows it iff a single spectrum is supplied and its sci_2d is populated (e.g. when read from an msaexp .spec.fits file). True forces it (silently no-ops when multiple spectra are supplied or no 2-D is available); False always shows the 1-D only. The two panels share the wavelength axis; emission-line markers span both.
cmap_2d (str) – Plotly colorscale for the 2-D panel (default "Blues" to mirror plot_2d_1d()).
vmin_pct (float) – Percentile clip for the 2-D colour scale (default 5 / 99.5).
vmax_pct (float) – Percentile clip for the 2-D colour scale (default 5 / 99.5).
y_crop ((float, float)) – Fractional crop (lo, hi) of the spatial axis on the 2-D panel (default (0.25, 0.75) keeps the middle 50 % of rows).
**read_kwargs – Forwarded to the file reader when a source item is a path.

Return type:

plotly.graph_objects.Figure

jwspecfit.read_dict(data, z=None, grating=None, R=None)[source]

Create a Spectrum from a dict with keys wave/lam, flux, err.

Wavelength assumed in microns.

Parameters:

data (dict) – Must contain "wave" or "lam" (µm), "flux" (µJy), "err" (µJy).
z (float, optional) – Source redshift.
grating (str, optional) – Grating name.
R (float, optional) – Resolving power.

Return type:

Spectrum

jwspecfit.read_fits(path, z=None, *, hdu=None, wave_col=None, flux_col=None, err_col=None)[source]

Read a 1-D spectrum from a FITS file.

By default tries the JWST/NIRSpec SPEC1D BinTable convention with columns wave (µm), flux (µJy), err (µJy). Falls through to auto-detection across all extensions when SPEC1D is absent or when hdu is given:

BinTable HDUs — the first table containing wavelength- and flux-like columns wins. Recognised column-name aliases include wave/wavelength/lam/lambda/loglam and flux/fnu/flam and err/error/sigma/noise (or ivar). Units are read from TUNITn keywords; common aliases (µm/Å/nm, µJy/mJy/Jy, erg/s/cm²/Å) are converted automatically.
Image HDUs — a 1-D image is read with the WCS keywords CRVAL1, CDELT1/CD1_1, CRPIX1, CTYPE1, CUNIT1; flux unit is BUNIT. Errors are not available from an image and are set to zero.

Parameters:

path (str or Path) – Path to the FITS file.
z (float, optional) – Source redshift to attach to the returned Spectrum.
hdu (str or int, optional) – Force a specific HDU to read (name or index). When None (default), tries SPEC1D first, then auto-detects.
wave_col (str, optional) – Force specific column names instead of auto-detecting. Only used for BinTable HDUs.
flux_col (str, optional) – Force specific column names instead of auto-detecting. Only used for BinTable HDUs.
err_col (str, optional) – Force specific column names instead of auto-detecting. Only used for BinTable HDUs.

Return type:

Spectrum

jwspecfit.read_npz(path, z=None, R=None)[source]

Read a stacked spectrum from a NumPy .npz file.

Expected keys: wave_angstrom, flux, err. Optionally n_stacked.

Parameters:

path (str or Path) – Path to the .npz file.
z (float, optional) – Source redshift.
R (float, optional) – Effective spectral resolving power of the stack.

Return type:

Spectrum

jwspecfit.resolve_R(lam_um, grating=None, R=None)[source]

Return R(λ) array from either a grating name or user-supplied R.

Parameters:

lam_um (array_like) – Wavelength in microns.
grating (str, optional) – Grating name ("PRISM", "G395M", etc.).
R (float or callable, optional) – Numeric resolving power (constant) or callable R(lam_um). Takes precedence over grating.

Returns:

Resolving power at each wavelength.

Return type:

np.ndarray

Raises:

ValueError – If neither grating nor R is specified.

jwspecfit.save_result(result, path)[source]

Save a FitResult to a .npz file for later replotting.

Parameters:

result (FitResult) – Fit result to save.
path (str or Path) – Output file path (should end in .npz).

Return type:

None

jwspecfit.show_lines(*, rest_min_A=None, rest_max_A=None, search=None)[source]

Print the available emission/absorption lines in REST_LINES_A.

The printed names are the keys you can pass to the lines= or add_lines= arguments of plot_spectrum_interactive() (and to the fitter), grouped by wavelength region.

Parameters:

rest_min_A (float, optional) – Restrict to lines with rest wavelengths in this range (Angstroms).
rest_max_A (float, optional) – Restrict to lines with rest wavelengths in this range (Angstroms).
search (str, optional) – Case-insensitive substring filter on the line name.

Return type:

None

Examples

>>> import jwspecfit
>>> jwspecfit.show_lines(rest_min_A=4000, rest_max_A=5100)
>>> jwspecfit.show_lines(search="Fe")

jwspecfit.sigma_inst_A(lam_um, grating=None, R=None)[source]

Instrumental Gaussian σ in Angstroms.

σ_inst = λ / (R × 2.3548)

Parameters:

lam_um (array_like) – Wavelength in microns.
grating (str, optional) – Grating name.
R (float or callable, optional) – Resolving power (overrides grating).

Returns:

Instrumental σ in Angstroms at each wavelength.

Return type:

np.ndarray

jwspecfit.tau_DLA(wave_A, log_NHI, z=0.0, b_kms=30.0)[source]

Compute DLA optical depth at each wavelength.

Implements Eq. 1 of Pollock et al. (2026):: tau_DLA = N_HI * sigma_0 * H(a, u)

where sigma_0 = sqrt(pi) * e^2 * f_alpha / (m_e * c * Delta_nu_D), a = Gamma_alpha / (4 pi Delta_nu_D), and u is the dimensionless frequency offset from Lya.

Parameters:

wave_A (np.ndarray) – Observed-frame wavelength in Angstrom.
log_NHI (float) – log10(N_HI / cm^-2).
z (float) – Source redshift (0 for rest-frame spectra).
b_kms (float) – Doppler parameter in km/s (default 30).

Returns:

Optical depth tau(lambda).

Return type:

np.ndarray

jwspecfit.tau_IGM_DW(wave_A, x_HI, z_source, z_min=5.3, cosmo=None)[source]

IGM damping-wing optical depth from a uniformly neutral IGM.

Implements the closed-form solution of Miralda-Escudé (1998), as written in Totani et al. (2006) Eq. 3 and used in Pollock et al. (2026) Eq. 3. Computes the redward Lyα damping-wing absorption by neutral hydrogen of mean fraction x_HI distributed uniformly between redshift z_min and z_source.

Parameters:

wave_A (np.ndarray) – Observed-frame wavelength array in Angstrom.
x_HI (float) – Volume-averaged neutral hydrogen fraction of the IGM.
z_source (float) – Source redshift (upper bound of the IGM integration).
z_min (float) – Lower bound of the IGM integration; default 5.3 (Bosman+22).
cosmo (astropy.cosmology, optional) – Cosmology used to compute the line-centre Gunn-Peterson optical depth. Defaults to Planck18.

Returns:

Optical depth τ_IGM(λ). Zero for pixels redward of the source’s Lyα (the formula returns zero by construction outside the integration range), and a finite damping-wing contribution between λ_α(1+z_source) and longer wavelengths.

Return type:

np.ndarray

Notes

For pixels blueward of λ_α(1+z_source) the photon would have been absorbed inside the GP forest already; the IGM damping-wing formula is not the right description there. Callers should enforce a separate model = 0 cut blueward of source Lyα.

Modules

`broad`	Broad Balmer component fitting with BIC model selection.
`constraints`	Parameter constraints: tied kinematics and fixed flux ratios.
`continuum`	Continuum fitting via iterative σ-clipped polynomial or median filter.
`data`
`dla`	DLA column density fitter — measure N_HI from Lya damping wings.
`fitter`	Core emission-line fitting engine.
`io`	Spectrum I/O: FITS and NPZ readers, Spectrum container.
`lines`	Emission-line database and line-list helpers.
`lyman_alpha`	Lyman-alpha fitting: skewed Gaussian profile × IGM absorption.
`models`	Gaussian emission-line profiles with bin-averaged evaluation.
`plotting`	Publication-quality visualisation of spectral fits.
`redshift`	Strong-line redshift fitting via brute-force grid evaluation.
`resolution`	Spectral resolution models for JWST NIRSpec gratings.