API Reference

The PST public API is organized by functionality. The sections below follow the current package layout under pst and replace the older documentation layout that grouped most classes under pst.models.

Simple Stellar Population Libraries

class pst.SSP.SSPBase

Bases: object

Base class that represents a model of Simple Stellar Populations.

Description

This class is meant for representing a model of Simple Stellar Populations that is mainly composed by a discreted grid of ages and metallicities and their associated spectral energy distributions.

ages

Ages of the SSPs.

Type:

astropy.units.Quantity

metallicities

Metallicities of the SSPs.

Type:

astropy.units.Quantity

L_lambda

Spectral energy distribution of each SSP. Each dimension correspond to (metallicity, ages, wavelength).

Type:

astropy.units.Quantity

wavelength

Wavelength array associated to the SED of the SSPs.

Type:

astroy.units.Quantity

default_path = '/home/docs/checkouts/readthedocs.org/user_builds/population-synthesis-toolkit/checkouts/latest/src/pst/data/ssp'

property name

Name of the SSP model.

property isochrone

Isochrone model associated to the SSP model.

property stellar_library

Stellar atmospheres model associated to the SSP model.

property imf

Return the IMF associated to the SSP model.

property ages

Ages of the SSPs.

property metallicities

Metallicities of the SSPs.

property n_ages

Number of ages in the SSP model.

property n_metallicities

Number of metallicities in the SSP model.

property n_ssps

Number of SSPs in the model.

property n_wavelengths

Number of wavelength points in the SSP model.

property L_lambda

Spectral energy distribution of each SSP. Each dimension correspond to (metallicity, ages, wavelength). Units must be luminosity per wavelength unit per initial mass unit (e.g. Lsun / Angstrom / Msun).

property wavelength

Wavelength array associated to the SED of the SSPs.

property remnant_mass_frac

Return the remnant mass of the SSP model.

property returned_mass_frac

Return the returned mass of the SSP model.

property current_mass

Return the current mass of the SSP model.

property log_ionising_HI_photons

Number of ionising photons of the SSP model in dex(s^-1).

property log_ionising_HeI_photons

Number of ionising photons of the SSP model in dex(s^-1).

property log_ionising_HeII_photons

Number of ionising photons of the SSP model in dex(s^-1).

property supernova_rate

Return the supernova rate of the SSP model.

get_weights(ages, metallicities, masses=None)

2D interpolation of a list of ages and metallicities.

Parameters:

• ages (np.array or astropy.units.Quantity) – SSP ages to interpolate.

• metallicites (np.array or astropy.units.Quantity) – Metallicity associated to each age.

• masses (np.array, astropy.units.Quantity or None, optional) – Stellar mass corresponding to each SSP.

get_ssp_l_lambda(age, metallicity)

Compute the SED associated to an SSP of a given age and metallicity.

Parameters:

• age (float or u.Quantity) – SSP age.

• metallicity (float or u.Quantity) – SSP metallicity

Returns:

sed (u.Quantity) – Spectra energy distribution associated to the SSP.

get_ssp_logedges()

Get the edges of the SSP metallicities and ages.

regrid(age_bins, metallicity_bins, verbose=True)

Interpolate the SSP model to a new grid of input ages and metallicities.

Parameters:

age_bins (np.array or astropy.units.Quantity)

cut_wavelength(wl_min=None, wl_max=None, verbose=True)

Cut model wavelength edges.

Parameters:

• wl_min (float or astropy.units.Quantity, optional) – Minimum wavelength value.

• wl_max (float or astropy.units.Quantity, optional) – Maximum wavelength value.

interpolate_sed(new_wl, verbose=True, log=False, **interp_kwargs)

Flux-conserving interpolation.

Parameters:

- new_wl (bin centers of the new interpolated points.)

get_mass_lum_ratio(wl_range)

Compute the mass-to-light ratio within a wavelength range.

get_specific_mass_lum_ratio(wl_range)

Compute the mass-to-light ratio per wavelength unit within a wavelength range.

get_ionising_photon_rate(species='HI')

Compute the ionising photon rate of the SSP model for a given species.

Parameters:

species (str, optional) – Species for which to compute the ionising photon rate. Must be one of ‘HI’, ‘HeI’ or ‘HeII’. Default is ‘HI’.

Returns:

photon_rate (astropy.units.Quantity) – Ionising photon rate of the SSP model for the given species in units of s^-1 / Msun.

compute_photometry(filter_list, z_obs=0.0, verbose=True)

Compute the SSP synthetic photometry of a set of filters.

Parameters:

filter_list (list of pst.observable.Filter) – A list of photometric filters.

Returns:

photometry (np.ndarray) – Array storing the photometry. The dimensions correspond to filter, metallicity and age.

copy()

Return a copy of the SSP model.

class pst.SSP.PopStar(IMF, nebular=False, path=None, verbose=True)

Bases: SSPBase

PopStar SSP models (Mollá et al. 2009).

This class represents the PopStar evolutionary synthesis models for simple stellar populations as presented in Mollá et al. (2009). It provides access to spectral energy distributions (SEDs) across a range of metallicities and ages, optionally including nebular emission.

Parameters:

• IMF (str) – Initial Mass Function name. Must be one of the supported IMFs listed below:

  • ‘sal’ : Salpeter (1955)

  • ‘fer’ : Ferrini, Penco & Palla (1990)

  • ‘kro’ : Kroupa (2002)

  • ‘cha’ : Chabrier (2003)

  This keyword is used to select the corresponding PopStar model files.

• nebular (bool, optional) – If True, include nebular emission in the SEDs. Default is False, meaning only the stellar continuum is considered.

• path (str or None, optional) – Filesystem path to the PopStar model data. If None (default), the package default path plus ‘PopStar’ subdirectory is used.

• verbose (bool, optional) – If True (default), print informational messages during initialization.

Example

Initialize with the desired IMF and nebular option:

>>> from pst.SSP import PopStar
>>> help(PopStar) # Retrieve the class documentation
>>> popstar = PopStar(IMF="cha", nebular=True) # Initialise the model
>>> print(popstar.wavelength)
<Quantity [91. , 94. , 96. , ..., 1200000., 1400000., 1600000.] Angstrom>
>>> oldest_sed = popstar.L_lambda[-1, -1] # Get the SED of the oldest SSP with highest metallicity

Notes

• The PopStar models cover metallicities from 0.0001 to 0.05 and ages spanning from 10^5 to ~10^10 years.

• Nebular emission is included only if nebular=True is set.

• Model files are expected to be in compressed FITS format with specific naming conventions as per Mollá et al. 2009.

For more details see the SSPBase base class documentation.

References

Mollá, M., García-Vargas, M. L., & Bressan, A. (2009). PopStar I: evolutionary synthesis model description. MNRAS, 398(1), 451-470. Fractal data repository

class pst.SSP.PyPopStar(IMF, nebular=False, path=None, verbose=True)

Bases: SSPBase

PyPopStar SSP models (Millán-Irigoyen et al. 2021).

This class represents the PyPopStar simple stellar population (SSP) models as described in Millán-Irigoyen et al. (2021). It provides spectral energy distributions (SEDs) for a range of metallicities and ages, optionally including nebular emission.

Parameters:

• IMF (str) – Initial Mass Function identifier. Must be one of the supported IMFs:

  • ‘KRO’ : Kroupa (2002)

  • ‘SAL’ : Salpeter (1955)

  • ‘CHA’ : Chabrier (2003)

  This keyword selects the corresponding model subdirectory and file set.

• nebular (bool, optional) – If True, nebular emission is included in the SEDs. Default is False, meaning only the stellar continuum is included.

• path (str or None, optional) – Filesystem path to the PyPopStar model data. If None (default), the package default path plus ‘PyPopStar’ and IMF subdirectories are used.

• verbose (bool, optional) – If True (default), prints informational messages during initialization.

Example

Initialize the PyPopStar model with a chosen IMF and nebular emission option:

>>> from pst.SSP import PyPopStar
>>> psp = PyPopStar(IMF='KRO', nebular=False)
>>> print(psp.wavelength)
<Quantity [ ... ] Angstrom>
>>> sed_example = psp.L_lambda[0, 0]  # SED for lowest metallicity, youngest age

Notes

• PyPopStar models cover metallicities from 0.004 to 0.05 and ages spanning approximately 10^5 to 10^10 years.

• Nebular emission is included only if nebular=True.

• Model files must follow the naming convention: ‘SSP-{IMF}_Z{metallicity}_logt{log_age}.fits’.

For additional details, see the SSPBase base class documentation.

References

Millán-Irigoyen, I.; Mollá, M.; Cerviño, M.; et al. (2021). HR-PYPOPSTAR: high-wavelength-resolution stellar populations evolutionary synthesis model MNRAS, 506 (4), 4781-4799.

Fractal data repository

class pst.SSP.BC03_2003(isochrone='Padova1994', model='BaSeL', imf='Chabrier', path=None, verbose=True)

Bases: SSPBase

BC03 SSP models (Bruzual & Charlot 2003, original version).

This class implements the 2003 version of the Bruzual & Charlot (BC03) evolutionary synthesis models. It supports both the STELIB and BaSeL stellar spectral libraries and various IMFs.

Parameters:

• isochrone (str, optional) – Isochrone set to use:’Padova1994’ (default) or ‘Padova2000’.

• model (str, optional) – Stellar spectral library. Must be one of:

  • ‘BaSeL’ : Low-resolution spectra

  • ‘STELIB’ : High-resolution spectra

  Case-insensitive. Default is ‘BaSeL’.

• imf (str, optional) – Initial Mass Function. Must be one of:

  • ‘chabrier’ or ‘cha’ : Chabrier (2003)

  • ‘salpeter’ or ‘sal’ : Salpeter (1955)

  • ‘kroupa’ or ‘kro’ : Kroupa (2002)

  Case-insensitive. Default is ‘chabrier’.

• path (str or None, optional) – Filesystem path to the BC03 model data. If None (default), the package default path is used.

• verbose (bool, optional) – If True (default), print informational messages during initialization.

Example

Create a BC03 SSP object using the Chabrier IMF and BaSeL library:

>>> from pst.SSP import BC03_2003
>>> ssp = BC03_2003(model='BaSeL', imf='Chabrier')
>>> print(ssp.wavelength)
<Quantity [ ... ] Angstrom>
>>> sed = ssp.L_lambda[0, -1]  # SED for lowest metallicity and oldest age

Notes

• This implementation loads the original 2003 release of BC03.

• Supported metallicities (Padova1994 isochrones):

  • ‘m22’: Z = 0.0001

  • ‘m32’: Z = 0.0004

  • ‘m42’: Z = 0.004

  • ‘m52’: Z = 0.008

  • ‘m62’: Z = 0.02

  • ‘m72’: Z = 0.05

• SEDs are read from FITS files with filenames of the form: bc2003_{res}_{metallicity}_{imf}_ssp.fits

• Units for L_lambda are Lsun / Angstrom / Msun.

• Ages are loaded from a standard time grid provided in TIME_SCALE.DAT.

See also

SSPBase

Base class for all SSP models.

References

Bruzual, G. & Charlot, S. (2003). Stellar population synthesis at the resolution of 2003. MNRAS, 344(4), 1000–1028.

Gustavo Bruzual page

class pst.SSP.BC03_2013(isochrone='Padova1994', model='BaSeL', imf='Chabrier', path=None, verbose=True)

Bases: SSPBase

BC03 SSP models (Bruzual & Charlot 2003, 2013 updated version).

This class implements the updated 2013 version of the Bruzual & Charlot (BC03) simple stellar population (SSP) models, using improved coverage of high metallicities and updated isochrone input. It supports both the STELIB and BaSeL stellar spectral libraries and several initial mass functions (IMFs).

Parameters:

• isochrone (str, optional) – Isochrone set to use. Currently only ‘Padova1994’ is supported (default).

• model (str, optional) – Stellar spectral library. Must be one of:

  • ‘BaSeL’ : Low-resolution spectra

  • ‘STELIB’ : High-resolution spectra

  Case-insensitive. Default is ‘BaSeL’.

• imf (str, optional) – Initial Mass Function. Must be one of:

  • ‘chabrier’ or ‘cha’ : Chabrier (2003)

  • ‘salpeter’ or ‘sal’ : Salpeter (1955)

  • ‘kroupa’ or ‘kro’ : Kroupa (2002)

  Case-insensitive. Default is ‘chabrier’.

• path (str or None, optional) – Filesystem path to the BC03 model data. If None (default), the package default path is used.

• verbose (bool, optional) – If True (default), print informational messages during initialization.

Example

Create a BC03 SSP object using the Chabrier IMF and BaSeL library:

>>> from pst.SSP import BC03_2013
>>> ssp = BC03_2013(model='BaSeL', imf='Chabrier')
>>> print(ssp.wavelength)
<Quantity [ ... ] Angstrom>
>>> sed = ssp.L_lambda[0, -1]  # SED for lowest metallicity and oldest age

Notes

• This version extends the original 2003 release by including an additional high-metallicity track (Z = 0.1).

• Supported metallicities (Padova1994 isochrones):

  • ‘m22’: Z = 0.0001

  • ‘m32’: Z = 0.0004

  • ‘m42’: Z = 0.004

  • ‘m52’: Z = 0.008

  • ‘m62’: Z = 0.02

  • ‘m72’: Z = 0.05

  • ‘m82’: Z = 0.1

• File naming convention: bc2003_{res}_{model}_{metallicity}_{imf}_ssp.fits

• Units for L_lambda are Lsun / Angstrom / Msun.

• Ages are taken from the same TIME_SCALE.DAT file as BC03_2003.

See also

SSPBase

Common functionality for SSP models.

References

Bruzual, G. & Charlot, S. (2003). Stellar population synthesis at the resolution of 2003. MNRAS, 344(4), 1000–1028.

Gustavo Bruzual page

class pst.SSP.BC03_2016(model='BaSeL', imf='Kroupa', path=None, load_properties=True, verbose=True)

Bases: SSPBase

BC03 SSP models (Bruzual & Charlot 2003, 2016 updated version).

This class implements the 2016 version of the Bruzual & Charlot (BC03) simple stellar population (SSP) models. It includes higher metallicity coverage and new stellar libraries, including MILES, STELIB, and BaSeL. Stellar tracks and isochrones are based on the PARSEC code (Bressan et al. 2012). This version uses updated file naming conventions and directory structure.

Parameters:

• model (str, optional) – Stellar spectral library. Must be one of:

  • ‘BaSeL’ : BaSeL 3.1 + ATLAS (low resolution)

  • ‘STELIB’ : STELIB + ATLAS (high resolution)

  • ‘MILES’ : MILES + STELIB (high resolution)

  Case-insensitive. Default is ‘BaSeL’.

• imf (str, optional) – Initial Mass Function. Must be one of:

  • ‘chabrier’ or ‘cha’ : Chabrier IMF

  • ‘salpeter’ or ‘sal’ : Salpeter IMF

  • ‘kroupa’ or ‘kro’ : Kroupa IMF

  Case-insensitive. Default is ‘Kroupa’.

• path (str or None, optional) – Filesystem path to the BC03 model data. If None (default), the package default path is used.

• load_properties (bool, optional) – If True (default), load the model properties during initialization. The properties include remnant mass fraction, returned mass fraction, supernova rate, and ionizing photon rates (Q_HI, Q_HeI, Q_HeII) as a function of age and metallicity (see SSPBase).

• verbose (bool, optional) – If True (default), print informational messages during initialization.

Example

Create a BC03 SSP object using the Chabrier IMF and BaSeL library:

>>> from pst.SSP import BC03_2016
>>> ssp = BC03_2016(model='BaSeL', imf='Chabrier')
>>> print(ssp.wavelength)
<Quantity [ ... ] Angstrom>
>>> sed = ssp.L_lambda[0, -1]  # SED for lowest metallicity and oldest age

Notes

• This version supports extended metallicity coverage up to Z = 0.1.

• The directory and file naming conventions follow:

  bc2003_{res}_{model_key}_{metallicity}_{imf_key}_ssp.fits

• Supported metallicities (Padova1994 isochrones):

  • ‘m22’: Z = 0.0001

  • ‘m32’: Z = 0.0004

  • ‘m42’: Z = 0.004

  • ‘m52’: Z = 0.008

  • ‘m62’: Z = 0.02

  • ‘m72’: Z = 0.05

  • ‘m82’: Z = 0.1

• MILES input is selected via model=’MILES’ (resolved as xmiless internally).

See also

SSPBase

Base class for all SSP models.

References

Bruzual, G., & Charlot, S. (2003). Stellar population synthesis at the resolution of 2003. MNRAS, 344(4), 1000-1028.

Gustavo Bruzual page

Gutkin, J., Charlot, S., & Bruzual, G. (2016). Modelling the nebular emission from primeval to present-day star-forming galaxies MNRAS, 462(2), 1757-1774.

isochrone = 'PARSEC2012'

class pst.SSP.BaseGM(path=None, verbose=True)

Bases: SSPBase

SSP model class for the GM (Granada-MILES) base used in the CALIFA survey.

This class loads a combined grid of simple stellar population (SSP) models:

• The young stellar populations (ages < 63 Myr) are taken from GRANADA models developed by González Delgado et al. (2005).

• The older populations (ages ≥ 63 Myr) are based on the MILES spectral library, specifically using the models from Vazdekis et al. (2010).

Parameters:

• path (str or None, optional) – Path to the gsd01_156.fits file containing the GRANADA SSP model cube. If None, a default path is used. The companion fits_like_properties.dat file containing the age and metallicity grid must reside in the same directory.

• verbose (bool, optional) – If True, print initialization messages. Default is True.

Notes

• The Granada models are provided as 156 SEDs, reshaped to a (39 × 4) grid of ages and metallicities.

• All models assume a Salpeter (1955) initial mass function (IMF).

• Spectral range: 3000 Å to 7000 Å, with 0.3 Å resolution.

Example

>>> from pst.SSP import BaseGM
>>> gm = BaseGM(verbose=False)
>>> print(gm.wavelength[0], gm.wavelength[-1])
3000.0 Angstrom 7000.0 Angstrom
>>> print(gm.ages.min(), gm.ages.max())
1e6 yr 1.3e10 yr

References

González Delgado, R. M., et al. (2005). “Evolutionary stellar population synthesis at high spectral resolution: optical wavelengths .” MNRAS, 357, 945. https://ui.adsabs.harvard.edu/abs/2005MNRAS.357..945G

Sánchez, S. F., et al. (2012). “CALIFA, the Calar Alto Legacy Integral Field Area survey: I. Survey presentation.” A&A, 538, A8. https://ui.adsabs.harvard.edu/abs/2012A%26A…538A…8S

class pst.SSP.XSL(IMF, ISO, path=None, verbose=True)

Bases: SSPBase

X-Shooter SSP empirical models.

These models provide spectral energy distributions for simple stellar populations observed with the X-Shooter instrument. They cover a wide range of ages and metallicities, using two sets of isochrones and two IMFs.

Parameters:

• IMF (str) – Initial Mass Function to use. Options:

  • ‘Kroupa’ : Kroupa (2001) IMF

  • ‘Salpeter’ : Salpeter (1955) IMF

• ISO (str) – Set of isochrones to use. Options:

  • ‘P00’ : Padova 2000 isochrones

  • ‘PC’ : PARSEC-COLIBRI isochrones

• path (str or None, optional) – Path to directory containing the SSP FITS files. If None, defaults to default_path/XSL/IMF.

• verbose (bool, optional) – If True, prints initialization messages. Default is True.

Notes

• Model files must be available in the specified path.

References

Verro, K., Trager, S.C., Peletier, R.F., et al. 2022, A&A, 661, A50.

Project home page: http://xsl.u-strasbg.fr/index.html

See also

SSPBase

Base class for SSP models.

Example

>>> from pst.SSP import XSL
>>> xsl = XSL(IMF='Kroupa', ISO='P00', verbose=False)
>>> print(xsl.wavelength[0], xsl.wavelength[-1])
3000.0 Angstrom 25000.0 Angstrom
>>> print(xsl.ages.min(), xsl.ages.max())
5e7 yr 1.6e10 yr
>>> print(xsl.metallicities)
[0.0004 0.001 0.004 0.008 0.019 0.03]
>>> print(xsl.L_lambda.shape)
(6, 26, 58642)

class pst.SSP.EMILES(iso, imf, path=None, verbose=True)

Bases: SSPBase

E-MILES simple stellar population (SSP) models.

These models provide high-quality empirical SSP spectra covering the spectral range 1680−50000 Å at moderately high resolution. The E-MILES library combines multiple stellar libraries and theoretical isochrones to deliver spectra of single-age, single-metallicity populations.

The UV spectral range (1680–3540 Å) is computed using the NGSL space-based stellar library, offering a significant improvement over earlier space-based models. The optical spectra use the MILES empirical library, and redder wavelengths are covered by Indo-US, CaT, and IRTF empirical libraries, all computed with consistent methods.

The SSP spectra span metallicities from −1.79 < [M/H] < +0.26 and ages greater than 30 Myr, across several IMF types with varying slopes.

Spectral resolution:

• UV range (λ < 3060 Å): FWHM ≈ 3 Å

• UV range (3060 Å < λ < 3540 Å): FWHM ≈ 5 Å

• Optical range (3540 Å to 8950 Å): FWHM ≈ 2.5 Å

• Infrared (longer wavelengths): σ = 60 km/s

Parameters:

• iso (str) – Isochrone set to use. Options:

  • ‘BASTI’ : BASTI isochrones (Pietrinferni et al. 2004)

  • ‘PADOVA00’ : Padova 2000 isochrones (Girardi et al. 2000)

• imf (str) – Initial mass function to use. Options:

  • ‘KROUPA_UNIVERSAL’ : Kroupa universal IMF (Kroupa 2001)

  • ‘UNIMODAL’ : Unimodal IMF (single power-law, Vazdekis et al. 1996)

  • ‘BIMODAL’ : Bimodal IMF (Vazdekis et al. 1996)

  • ‘CHABRIER’ : Chabrier IMF (Chabrier 2003)

• path (str or None, optional) – Path to the directory containing the E-MILES FITS files. Defaults to default_path/EMILES.

• verbose (bool, optional) – If True, prints initialization messages. Default is True.

Notes

• The spectral resolution and wavelength coverage follow Vazdekis et al. (2016) and Röck et al. (2016).

• The models cover a wavelength range of approximately 1680–50000 Å with spectral resolution ~2.5 Å (FWHM).

• The metallicities are given as logarithmic values relative to solar, with solar metallicity fixed at Z=0.019.

References

Vazdekis, A., Koleva, M., Ricciardelli, E., et al. 2016, MNRAS, 463, 3409 Röck, B., Vazdekis, A., Ricciardelli, E., et al. 2016, A&A, 589, A73, 8

See also

SSPBase

Base class for simple stellar population models.

Example

>>> from pst.SSP import EMILES
>>> emiles = EMILES(iso='BASTI', imf='KROUPA_UNIVERSAL', verbose=False)
>>> print(emiles.wavelength[0], emiles.wavelength[-1])
1680.0 Angstrom 50000.0 Angstrom
>>> print(emiles.ages.min(), emiles.ages.max())
30.0 Myr 14.0 Gyr
>>> print(emiles.metallicities)
[0.0001 0.00016 ... 0.03]  # actual values depend on the isochrone
>>> print(emiles.L_lambda.shape)
(12, 52, 53689)

_get_models_prefix(iso, imf)

Get the file prefix used to load the model files.

_load_models(model_name)

Load the SSP models from individual FITS files.

Parameterized Model Infrastructure

class pst.model.Parameter(value: int | float | number | Quantity, unit: Unit | None = None, vrange: Tuple[int | float | number | Quantity, int | float | number | Quantity] | None = None, fixed: bool = False, doc: str = '')

Bases: object

A model parameter that behaves like an astropy Quantity.

property vrange

Allowed interval as (vmin, vmax), or None.

property q: Quantity

Return the underlying Quantity.

property value_raw

Raw numeric value of the underlying Quantity (unit-stripped).

property unit_raw: Unit

Unit of the underlying Quantity.

to(unit: Unit, equivalencies=None) → Quantity

Return a converted Quantity (does not modify the Parameter).

to_value(unit: Unit | None = None, equivalencies=None)

Return numeric value in the requested unit.

convert_to(unit: Unit, equivalencies=None) → Parameter

Convert the parameter in place and return self.

set(new_value: int | float | number | Quantity, *, validate: bool = True) → None

Update the parameter value.

Parameters:

• new_value (number or astropy.units.Quantity) – New parameter value.

• validate (bool, optional) – If True, enforce vrange constraints when defined.

as_quantity() → Quantity

Return the current parameter value as an astropy Quantity.

property size

Number of scalar elements in the underlying quantity.

__array_ufunc__(ufunc, method, *inputs, **kwargs)

Intercept numpy ufuncs and delegate to Quantity.

class pst.model.ModelBase

Bases: object

Base class for models with named parameters.

Subclasses typically declare Parameter attributes directly, for example:

class MyModel(ModelBase):

name: str = “my_model” a_v: Parameter = Parameter(0.2, vrange=(0.0, 5.0)) r_v: Parameter = Parameter(3.1, vrange=(2.0, 6.0), fixed=True)

Notes

This base class provides: - discovery of Parameter attributes - getting and setting values by name - filtering fixed or free parameters - conversion to plain dictionaries for IO

name: str = 'model'

models_recursive(*, prefix: str = '', max_depth: int | None = None) → Dict[str, ModelBase]

Return nested models including self.

Returns:

models (dict) – Mapping from dotted path to model instance.

parameters_recursive(*, prefix: str = '', max_depth: int | None = None, include_fixed: bool = True) → Dict[str, Parameter]

Return parameters from this model and nested models.

Returns:

params (dict) – Mapping from dotted parameter path to Parameter.

parameters() → Dict[str, Parameter]

Return direct Parameter attributes defined on this model.

Returns:

params (dict) – Mapping from attribute name to Parameter.

parameter_names(*, include_fixed: bool = True) → List[str]

Return parameter names.

Parameters:

include_fixed (bool, optional) – If False, only return free parameters.

Returns:

names (list of str)

get(name: str) → Parameter

Get a Parameter object by name.

Raises:

KeyError – If the parameter does not exist.

get_values(*, include_fixed: bool = True, as_quantity: bool = False) → Dict[str, Any]

Get parameter values as a dict.

Parameters:

• include_fixed (bool, optional) – If False, only returns free parameters.

• as_quantity (bool, optional) – If True, returns values as Quantities when possible.

Returns:

values (dict) – Mapping from parameter name to value.

set_values(values: Mapping[str, Any], *, validate: bool = True, strict: bool = True) → None

Set parameter values from a mapping.

Parameters:

• values (mapping) – Mapping from parameter name to new value.

• validate (bool, optional) – If True, checks fixed and vrange constraints. Default is True.

• strict (bool, optional) – If True, unknown keys raise KeyError. If False, unknown keys are ignored.

freeze(names: Sequence[str] | None = None) → None

Freeze parameters.

Parameters:

names (sequence of str or None) – If None, freeze all parameters. Otherwise freeze selected ones.

unfreeze(names: Sequence[str] | None = None) → None

Unfreeze parameters.

Parameters:

names (sequence of str or None) – If None, unfreeze all parameters. Otherwise unfreeze selected ones.

search(text: str, *, in_docs: bool = True, in_names: bool = True) → List[str]

Search parameters by substring.

Parameters:

• text (str) – Substring to search for, case-insensitive.

• in_docs (bool, optional) – If True, search in Parameter.doc strings.

• in_names (bool, optional) – If True, search in parameter names.

Returns:

matches (list of str) – Matching parameter names.

to_dict(*, include_fixed: bool = True) → Dict[str, Any]

Serialize the model configuration to a plain dict.

Notes

This returns a simple structure suitable for JSON/YAML. Quantities are represented as dicts with value and unit.

Chemical Evolution Models (CEM)

pst.cem._check_time_dec(func)

Decorator that normalizes a time argument to an astropy.units.Quantity.

The wrapped method may be called with:

• a bare number / numpy array (assumed to be in Gyr),

• an astropy.units.Quantity convertible to Gyr,

• a pst.model.Parameter whose .q is convertible to Gyr.

The wrapped function is always called with a Quantity in Gyr.

Notes

This decorator only converts the first positional argument after self (conventionally named time or times). Any additional positional and keyword arguments are passed through unchanged.

class pst.cem.MassPropMetallicityMixin

Bases: object

Mixin implementing a mass-dependent metallicity evolution.

This mixin provides a default implementation of ism_metallicity() where the ISM metallicity scales with the cumulative stellar mass formed:

\[Z(t) = Z_{\rm today}\,\left(\frac{M_\star(t)}{M_\star({\rm today})}\right)^\alpha\]

The intention is to allow any star-formation-history (SFH) model to be combined with a simple analytic chemical enrichment prescription.

Required host attributes

The host class must define:

• mass_todayQuantity or pst.model.Parameter

  Total stellar mass formed at the observing time (same units as stellar_mass_formed()).

• ism_metallicity_todayQuantity or Parameter

  Present-day ISM metallicity (dimensionless mass fraction).

• alpha_powerlawQuantity/Parameter or float

  Power-law exponent controlling the enrichment rate.

• stellar_mass_formed(times)()

  Cumulative stellar mass formed as a function of cosmic time.

Notes

• The returned metallicity is dimensionless (mass fraction).

• For physical results, host models should ensure stellar_mass_formed() is non-decreasing and non-negative.

property ism_metallicity_today: Parameter

ISM metals mass fraction at present.

property alpha_powerlaw: Parameter

Stellar mass power-law exponent.

ism_metallicity(times: array | Quantity) → Quantity

Evaluate ISM metallicity at one or more cosmic times.

Parameters:

times (array-like or astropy.units.Quantity) – Cosmic times since the Big Bang. If a bare array is provided, it is interpreted as Gyr.

Returns:

z_t (astropy.units.Quantity) – Dimensionless ISM metallicity (mass fraction) evaluated at each input time. Shape matches times.

pst.cem.sfh_quenching_decorator(stellar_mass_formed)

Decorator that enforces a hard quenching event in a cumulative SFH.

The decorated stellar_mass_formed() behaves normally for time < quenching_time. For later times it returns the constant value \(M_\star(\mathrm{quenching\_time})\), i.e. no further stellar mass is formed.

The host instance is expected to provide an attribute quenching_time, which may be a pst.model.Parameter, an astropy.units.Quantity, or a bare number (assumed Gyr).

Notes

• This decorator assumes the wrapped function returns the cumulative stellar mass formed.

• The quenching is instantaneous and permanent (a step in the cumulative SFH).

pst.cem.weights_cache_decorator(func)

Decorator that caches SSP weights computed by a CEM instance.

The decorated method is expected to compute SSP weights for a given SSP model and observing time. The cache is stored in the instance attribute _ssp_weights_cache as a dictionary mapping (ssp_name, t_obs) tuples to weight arrays.

Notes

• This decorator assumes the wrapped method has signature func(self, ssp, t_obs, *args, **kwargs) where ssp is an instance of pst.SSP.SSPBase and t_obs is an astropy.units.Quantity.

  • The cache keys are based only on the SSP name and the observing time value in Gyr.

    Any other argument that changes the interpolation result, such as oversample_factor, is ignored by the cache key.

• The cache values are the computed weights arrays. The decorator does not attempt to verify the consistency of the cached weights with the input

  arguments beyond the SSP name and observing time.

  • Different SSP objects sharing the same ssp.name will collide in the

    cache.

class pst.cem.ChemicalEvolutionModel(today: Parameter | Quantity | float = None, cache_interp_ssp_mass: bool = False)

Bases: ModelBase, ABC

Base class for chemical evolution models (CEMs).

A Chemical Evolution Model defines two core functions of cosmic time:

• stellar_mass_formed(t)(): cumulative stellar mass formed up to time t

• ism_metallicity(t)(): ISM metallicity (mass fraction) at time t

The class also provides utilities to synthesize spectra and photometry from a pst.SSP.SSPBase grid by integrating the SFH over SSP age bins.

Time convention

All “time” arguments are cosmic time since the Big Bang, not lookback time. The optional attribute today represents the cosmic age of the Universe at the observing epoch and is used by helper methods that require a normalization point (e.g. formation-time percentiles).

Caching

The method interpolate_ssp_masses() can be decorated with the weights_cache_decorator() to enable caching of computed SSP weights for given SSP models and observing times. This can speed up repeated SED synthesis at the cost of increased memory usage. The cache is stored in the instance attribute _ssp_weights_cache as a dictionary mapping (ssp_name, t_obs) tuples to weight arrays.

The cache keys are based only on the SSP name and the observing time value in Gyr, so cached weights can be reused even when other inputs that affect the interpolation change. In particular, changing oversample_factor while caching is enabled can silently return weights computed with a different oversampling setup.

Parameters:

• today (float, astropy.units.Quantity, or pst.model.Parameter, optional) – Cosmic time of observation. If a bare number is provided it is assumed to be in Gyr.

• cache_interp_ssp_mass (bool, optional) – If True, cache the results of interpolate_ssp_masses() by SSP name and observing time. This can accelerate repeated evaluations with identical inputs, but it is only safe when all other inputs affecting the interpolation, such as oversample_factor, remain unchanged.

Notes

Implementations should strive to ensure:

• stellar_mass_formed() is non-decreasing with time,

• stellar_mass_formed(t < 0) = 0(),

• behavior for t > today is well-defined (commonly clipped to \(M_\star(\mathrm{today})\)).

name: str = 'base_cem'

property today

abstractmethod stellar_mass_formed(time: Quantity) → Quantity

Cumulative stellar mass formed as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

mstar (astropy.units.Quantity) – Cumulative stellar mass formed up to each input time. Units should be mass (typically Msun). Shape matches time.

abstractmethod ism_metallicity(time: Quantity) → Quantity

ISM metallicity (mass fraction) as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

z_t (astropy.units.Quantity) – Dimensionless metallicity (mass fraction). Shape matches time.

interpolate_ssp_masses(ssp: SSPBase, t_obs: Quantity, oversample_factor=10) → Quantity

Compute SSP mass weights by integrating the SFH over SSP age bins.

This routine discretizes the interval [0, t_obs] into SSP age bins and computes the stellar mass formed in each bin via differences of the cumulative SFH, then maps each bin to the SSP grid using pst.SSP.SSPBase.get_weights(). When cache_interp_ssp_mass is enabled, the cached result depends only on ssp.name and t_obs.

Parameters:

• ssp (pst.SSP.SSPBase) – SSP model providing an age grid ssp.ages (shape (A,)), metallicity grid, and a mapping function get_weights.

• t_obs (astropy.units.Quantity) – Cosmic time of observation. Only SSP ages younger than t_obs contribute.

• oversample_factor (int, optional) – Factor by which SSP age bins are oversampled for smoother integration. Default is 10.

Returns:

weights (astropy.units.Quantity) – Mass weights on the SSP grid. The expected shape is (Z, A) (metallicity, age). Units are mass (typically Msun).

Notes

The method uses midpoints of SSP age bins and evaluates stellar_mass_formed() and ism_metallicity() at corresponding cosmic times (t_obs - age).

If caching is enabled, callers should keep oversample_factor fixed for repeated evaluations at the same ssp and t_obs. Changing oversample_factor does change the interpolation, but it is not part of the cache key.

surviving_stellar_mass(ssp: SSPBase, t_obs: Quantity)

Compute the surviving stellar mass at an observing time.

This method integrates the SFH over SSP age bins as in interpolate_ssp_masses() but applies SSP mass loss to compute the surviving mass rather than the formed mass.

Parameters:

• ssp (pst.SSP.SSPBase) – SSP model providing an age grid ssp.ages (shape (A,)), metallicity grid, and a mapping function get_weights.

• t_obs (astropy.units.Quantity) – Cosmic time of observation. Only SSP ages younger than t_obs contribute.

Returns:

m_surviving (astropy.units.Quantity) – Total surviving stellar mass at t_obs. Units are mass (typically Msun).

supernova_rate(ssp: SSPBase, t_obs: Quantity, cc=True, ia=True)

Compute the supernova rate at an observing time.

This method integrates the SFH over SSP age bins as in interpolate_ssp_masses() but applies the SSP supernova rate to compute the total SN rate rather than the formed mass.

Parameters:

• ssp (pst.SSP.SSPBase) – SSP model providing an age grid ssp.ages (shape (A,)), metallicity grid, and a mapping function get_weights.

• t_obs (astropy.units.Quantity) – Cosmic time of observation. Only SSP ages younger than t_obs contribute.

• cc (bool, optional) – If True (default), include core-collapse supernovae in the rate.

• ia (bool, optional) – If True (default), include Type Ia supernovae in the rate.

Returns:

sn_rate (astropy.units.Quantity) – Total supernova rate at t_obs. Units are typically SN/yr.

ionising_photon_rate_hi(ssp: SSPBase, t_obs: Quantity, species='HI')

Compute the ionising photon production rate at an observing time.

This method integrates the SFH over SSP age bins as in interpolate_ssp_masses() but applies the SSP ionising photon rate to compute the total production rate rather than the formed mass.

Parameters:

• ssp (pst.SSP.SSPBase) – SSP model providing an age grid ssp.ages (shape (A,)), metallicity grid, and a mapping function get_weights.

• t_obs (astropy.units.Quantity) – Cosmic time of observation. Only SSP ages younger than t_obs contribute.

• species ({‘HI’, ‘HeI’, ‘HeII’}, optional) – Ion species for which to compute the photon rate. Default is ‘HI’.

Returns:

q_ion (astropy.units.Quantity) – Total ionising photon production rate at t_obs in dex(photons/s).

time_at_stellar_mass_frac(frac, today=None, time_res=<Quantity 30. Myr>)

Time at which a fraction of the final stellar mass has formed.

This function samples the cumulative SFH on a regular cosmic-time grid and linearly interpolates the time(s) at which:

\[\frac{M_\star(t)}{M_\star(\mathrm{today})} = f\]

Parameters:

• frac (float or array-like) – Target mass fraction(s) in the range [0, 1].

• today (float, astropy.units.Quantity, or pst.model.Parameter, optional) – Cosmic observing time. If not provided, uses self.today.

• time_res (float or astropy.units.Quantity, optional) – Sampling resolution of the time grid. If a bare number is provided it is assumed to be in Gyr. Default is 30 Myr.

Returns:

t_frac (astropy.units.Quantity) – Cosmic time(s) at which the specified fraction(s) are reached. Shape matches frac.

Raises:

ValueError – If neither today nor self.today is set.

mean_stellar_age(ssp: SSPBase, t_obs: Quantity, log: bool = False, surviving_mass: bool = False)

Compute the mean stellar age at an observing time.

Parameters:

• ssp (pst.SSP.SSPBase) – SSP model providing an age grid ssp.ages (shape (A,)), metallicity grid, and a mapping function get_weights.

• t_obs (astropy.units.Quantity) – Cosmic time of observation. Only SSP ages younger than t_obs contribute.

• log (bool, optional) – If True, compute the logarithmic mean age (base 10) rather than the linear mean. Default is False.

• surviving_mass (bool, optional) – If True, weight SSP ages by the surviving stellar mass rather than the formed mass. Default is False.

Returns:

mean_age (astropy.units.Quantity) – Mean stellar age at t_obs. Units are time (typically Gyr).

compute_SED(ssp: SSPBase, t_obs: Quantity, *, allow_negative: bool = False, age_bin_edges=None) → Quantity

Synthesize the spectral energy distribution (SED) at an observing time.

The SED is obtained by summing SSP spectra weighted by the stellar mass formed in each SSP age/metallicity bin.

Parameters:

• ssp (pst.SSP.SSPBase) – SSP model providing ssp.L_lambda with shape (Z, A, W) where W is wavelength.

• t_obs (astropy.units.Quantity) – Cosmic time of observation.

• allow_negative (bool, optional) – If False (default), any negative SSP weights are clipped to zero.

• age_bin_edges (array-like or astropy.units.Quantity, optional) – If provided, the SED is returned in coarse age bins. Edges are in the same units as ssp.ages and define Nbin = len(edges) - 1 bins.

Returns:

sed (astropy.units.Quantity) – If age_bin_edges is None: shape (W,). If age_bin_edges is provided: shape (Nbin, W). Units are (mass unit) * (ssp.L_lambda unit).

Notes

age_bin_edges bins SSP ages (not cosmic times). Internally the mapping is performed using numpy.digitize() on ssp.ages.

compute_photometry(ssp, t_obs, *, photometry: Quantity = None, allow_negative: bool = False, age_bin_edges=None) → Quantity

Synthesize broadband photometry at an observing time.

Parameters:

• ssp (pst.SSP.SSPBase) – SSP model defining the SSP grid used to compute weights.

• t_obs (astropy.units.Quantity) – Cosmic time of observation.

• photometry (astropy.units.Quantity or ndarray, optional) – Photometry grid on the SSP axes with shape (B, Z, A), where B is the number of bands. If an ndarray is provided, it is assumed to be in Jy / Msun. If None, uses ssp.photometry.

• allow_negative (bool, optional) – If False (default), negative SSP weights are clipped to zero.

• age_bin_edges (array-like or astropy.units.Quantity, optional) – If provided, photometry is returned per coarse age bin defined over SSP ages. Output will have Nbin = len(edges) - 1 bins.

Returns:

model_photometry (astropy.units.Quantity) – If age_bin_edges is None: shape (B,). If age_bin_edges is provided: shape (Nbin, B). Units are photometry.unit * weights.unit.

Notes

This method assumes the photometry grid is expressed per unit formed mass (e.g. Jy / Msun). The output is therefore an absolute flux/luminosity in the same system as the input grid.

class pst.cem.SingleBurstCEM(*, mass_burst: Parameter | Quantity | float = <Quantity 1. solMass>, time_burst: Parameter | Quantity | float = <Quantity 0. Gyr>, burst_metallicity: Parameter | Quantity | float = <Quantity 0.02>, today: Parameter | Quantity | float | None = None)

Bases: ChemicalEvolutionModel

Instantaneous single-burst star formation history.

The model forms a total mass mass_burst at cosmic time time_burst and forms no additional stars thereafter. The cumulative formed mass is:

• 0 for t < time_burst

• mass_burst for t >= time_burst

Parameters:

• mass_burst (float, astropy.units.Quantity, or pst.model.Parameter, optional) – Total stellar mass formed in the burst. Bare numbers are assumed in Msun.

• time_burst (float, astropy.units.Quantity, or pst.model.Parameter, optional) – Cosmic time of the burst. Bare numbers are assumed in Gyr.

• burst_metallicity (float, astropy.units.Quantity, or pst.model.Parameter, optional) – Constant ISM metallicity assigned to all stars (dimensionless mass fraction).

Notes

This is a cumulative SFH model; the instantaneous SFR is a delta function.

name: str = 'single_burst_cem'

stellar_mass_formed(time: Unit('Gyr')) → Quantity

Total stellar mass formed at a given time.

ism_metallicity(time: Unit('Gyr')) → Quantity

ISM metals mass fraction at a given time.

class pst.cem.ExponentialCEM(*, stellar_mass_inf: Parameter | Quantity | float = <Quantity 1. solMass>, tau: Parameter | Quantity | float = <Quantity 1. Gyr>, metallicity: Parameter | Quantity | float = <Quantity 0.02>, today: Parameter | Quantity | float | None = None)

Bases: ChemicalEvolutionModel

Exponentially-declining cumulative SFH.

The cumulative formed mass follows:

\[M_\star(t) = M_{\infty}\,\left(1 - e^{-t/\tau}\right)\]

where M_infty is the asymptotic formed mass and tau is the e-folding time.

Parameters:

• stellar_mass_inf (float, Quantity, or Parameter, optional) – Asymptotic formed stellar mass M_infty. Bare numbers assumed in Msun.

• tau (float, Quantity, or Parameter, optional) – E-folding time in Gyr.

• metallicity (float, Quantity, or Parameter, optional) – Constant ISM metallicity (dimensionless mass fraction).

Notes

The metallicity returned by ism_metallicity() is constant in time.

name: str = 'exponential_cem'

stellar_mass_formed(time: Unit('Gyr')) → Quantity

Cumulative stellar mass formed as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

mstar (astropy.units.Quantity) – Cumulative stellar mass formed up to each input time. Units should be mass (typically Msun). Shape matches time.

ism_metallicity(time: Unit('Gyr')) → Quantity

ISM metallicity (mass fraction) as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

z_t (astropy.units.Quantity) – Dimensionless metallicity (mass fraction). Shape matches time.

class pst.cem.ExponentialQuenchedCEM(*, quenching_time: Parameter | Quantity | float | None = None, **kwargs)

Bases: ExponentialCEM

Exponentially declining CEM model including a quenching event.

Examples

Create an instance of the ExponentialQuenchedCEM model and compute the stellar mass formed:

>>> from astropy import units as u
>>> from pst.models import ExponentialQuenchedCEM
>>> model = ExponentialQuenchedCEM(stellar_mass_inf=1e10 * u.Msun, tau=2 * u.Gyr,
...                                 metallicity=0.02, quenching_time=1.5 * u.Gyr)
>>> time = [0.5, 1.0, 2.0] * u.Gyr
>>> model.stellar_mass_formed(time)
<Quantity [2.21199217e+09, 3.93469340e+09, 5.27633447e+09] solMass>

See also

ExponentialCEM

name: str = 'exponential_quenched_cem'

stellar_mass_formed(time: Unit('Gyr'))

Cumulative stellar mass formed as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

mstar (astropy.units.Quantity) – Cumulative stellar mass formed up to each input time. Units should be mass (typically Msun). Shape matches time.

class pst.cem.ExponentialDelayedCEM(*, tau: Parameter | Quantity | float = <Quantity 2. Gyr>, mass_today: Parameter | Quantity | float = <Quantity 1. solMass>, ism_metallicity_today: Parameter | Quantity | float = 0.02, today: Parameter | Quantity | float | None = None)

Bases: ChemicalEvolutionModel

Delayed-exponential cumulative SFH normalized at today.

The cumulative formed mass is:

\[M_\star(t) \propto 1 - e^{-t/\tau}\,\frac{t+\tau}{\tau}\]

This class rescales the proportional form so that \(M_\star(\\mathrm{today}) = \\mathrm{mass\\_today}\).

Parameters:

• tau (float, Quantity, or Parameter, optional) – Timescale parameter (Gyr).

• mass_today (float, Quantity, or Parameter, optional) – Total formed stellar mass at today (Msun).

• ism_metallicity_today (float, Quantity, or Parameter, optional) – Constant metallicity (dimensionless mass fraction).

• today (float, Quantity, or Parameter) – Cosmic observing time required for normalization.

Raises:

ValueError – If today is not set.

name: str = 'exponential_delayed_cem'

stellar_mass_formed(time)

Cumulative stellar mass formed as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

mstar (astropy.units.Quantity) – Cumulative stellar mass formed up to each input time. Units should be mass (typically Msun). Shape matches time.

ism_metallicity(time: Unit('Gyr'))

ISM metallicity (mass fraction) as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

z_t (astropy.units.Quantity) – Dimensionless metallicity (mass fraction). Shape matches time.

class pst.cem.ExponentialDelayedZPowerLawCEM(*, ism_metallicity_today: Parameter | Quantity | float = <Quantity 0.02>, alpha_powerlaw: Parameter | Quantity | float = <Quantity 1.>, **kwargs)

Bases: MassPropMetallicityMixin, ExponentialDelayedCEM

A ExponentialDelayedCEM with a Mass-dependent Metallicity chemical model.

See also

MassPropMetallicityMixin

name: str = 'exponential_delayed_zpowlaw_cem'

class pst.cem.ExponentialDelayedQuenchedCEM(*, quenching_time: Parameter | Quantity | float | None = None, **kwargs)

Bases: ExponentialDelayedZPowerLawCEM

A ExponentialDelayedZPowerLawCEM with a quenching event.

name: str = 'exponential_delayed_quenched_cem'

stellar_mass_formed(times: Quantity)

Cumulative stellar mass formed as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

mstar (astropy.units.Quantity) – Cumulative stellar mass formed up to each input time. Units should be mass (typically Msun). Shape matches time.

class pst.cem.GaussianBurstCEM(*, mass_burst: Parameter | Quantity | float, time_burst: Parameter | Quantity | float, sigma_burst: Parameter | Quantity | float, burst_metallicity: Parameter | Quantity | float, today: Parameter | Quantity | float | None = None)

Bases: ChemicalEvolutionModel

Gaussian burst star formation model.

This model represents a galaxy’s star formation history as a single Gaussian burst of star formation. The burst occurs at a specific time and is characterized by a Gaussian distribution with a given standard deviation.

mass_burst

Total stellar mass formed in the burst.

Type:

float or astropy.Quantity

time_burst

Time of the starburst in cosmic time.

Type:

float or astropy.Quantity

sigma_burst

Span time of the burst in terms of the standard deviation.

Type:

float or astropy.Quantity

burst_metallicity

Metallicity of the burst.

Type:

float or astropy.Quantity

Examples

Create an instance of the GaussianBurstCEM model and compute the stellar mass formed:

>>> from astropy import units as u
>>> from pst.models import GaussianBurstCEM
>>> model = GaussianBurstCEM(mass_burst=1e10 * u.Msun, time_burst=2 * u.Gyr,
...                          sigma_burst=0.5 * u.Gyr, burst_metallicity=0.02)
>>> time = [1.0, 2.0, 3.0] * u.Gyr
>>> model.stellar_mass_formed(time)
<Quantity [2.27501319e+08, 5.00000000e+09, 9.77249868e+09] solMass>

name: str = 'gaussian_burst_cem'

stellar_mass_formed(time)

Cumulative stellar mass formed as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

mstar (astropy.units.Quantity) – Cumulative stellar mass formed up to each input time. Units should be mass (typically Msun). Shape matches time.

ism_metallicity(time: Unit('Gyr'))

ISM metallicity (mass fraction) as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

z_t (astropy.units.Quantity) – Dimensionless metallicity (mass fraction). Shape matches time.

class pst.cem.LogNormalCEM(*, t0: Parameter | Quantity | float = <Quantity 1. Gyr>, scale: Parameter | Quantity | float = <Quantity 1.>, mass_today: Parameter | Quantity | float = <Quantity 1. solMass>, ism_metallicity_today: Parameter | Quantity | float = <Quantity 0.02>, today: Parameter | Quantity | float | None = None)

Bases: ChemicalEvolutionModel

Log-normal star formation history model.

This CEM models a galaxy’s star formation rate as a log-normal function of time (e.g. Gladders et al 2013), where the SFR rises initially and then decays:

\[M_\star(t) = \frac{M_{\rm today}}{2} \cdot \left(1 - {\rm erf}\left(\frac{{\rm ln}(t) - {\rm ln}(t_0)}{\sigma \sqrt{2}}\right) \right)\]

lnt0

Mean of the lognormal SFH.

Type:

float

tau

Standard deviation of the lognormal SFH.

Type:

astropy.Quantity

mass_today

Total stellar mass formed at present.

Type:

float or astropy.Quantity

metallicity

Metallicity of the gas (constant).

Type:

float

name: str = 'lognormal_cem'

stellar_mass_formed(times: Quantity)

Cumulative stellar mass formed as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

mstar (astropy.units.Quantity) – Cumulative stellar mass formed up to each input time. Units should be mass (typically Msun). Shape matches time.

ism_metallicity(time: Unit('Gyr'))

ISM metallicity (mass fraction) as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

z_t (astropy.units.Quantity) – Dimensionless metallicity (mass fraction). Shape matches time.

class pst.cem.LogNormalZPowerLawCEM(*, ism_metallicity_today: Parameter | Quantity | float = <Quantity 0.02>, alpha_powerlaw: Parameter | Quantity | float = <Quantity 1.>, **kwargs)

Bases: MassPropMetallicityMixin, LogNormalCEM

A LogNormalCEM with a Mass-dependent Metallicity chemical model.

See also

MassPropMetallicityMixin

name: str = 'lognormal_zpowlaw_cem'

class pst.cem.LogNormalQuenchedCEM(*, quenching_time: Parameter | Quantity | float | None = None, **kwargs)

Bases: LogNormalZPowerLawCEM

A LogNormalCEM with a quenching event.

name: str = 'lognormal_quenched_cem'

stellar_mass_formed(times: Quantity)

Cumulative stellar mass formed as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

mstar (astropy.units.Quantity) – Cumulative stellar mass formed up to each input time. Units should be mass (typically Msun). Shape matches time.

class pst.cem.BetaCEM(*, mass_today: Parameter | Quantity | float = <Quantity 1. solMass>, alpha: Parameter | Quantity | float | None = None, beta: Parameter | Quantity | float | None = None, t_mean: Parameter | Quantity | float | None = None, kappa: Parameter | Quantity | float | None = None, t_start: Parameter | Quantity | float = <Quantity 0. Gyr>, t_end: Parameter | Quantity | float | None = None, ism_metallicity_today: Parameter | Quantity | float = <Quantity 0.02>, today: Parameter | Quantity | float | None = None)

Bases: ChemicalEvolutionModel

Beta-CMF chemical evolution model.

This model defines the cumulative formed stellar mass using the regularised incomplete Beta function:

\[f_\star(t) = I_x(\alpha, \beta), \qquad x = \frac{t - t_{\rm start}}{t_{\rm end} - t_{\rm start}} \in [0,1]\]

and

\[M_{\rm formed}(t) = M_{\rm today}\, f_\star(t)\]

where \(I_x\) is the regularised incomplete Beta function.

Parameterisation

You can specify the Beta shape parameters directly (alpha, beta), or provide the mean formation time and concentration (t_mean, kappa):

• kappa = alpha + beta

• m = t_mean / t_end (fractional mean in (0,1))

• alpha = m * kappa, beta = (1-m) * kappa

Parameters:

• mass_today (float, astropy.units.Quantity, or pst.model.Parameter) – Total formed stellar mass by t_end. Bare numbers are assumed in Msun.

• alpha, beta (float, Quantity, or Parameter, optional) – Beta-CMF shape parameters (>0). Provide both if using the direct form. If provided as Quantity/Parameter they must be dimensionless.

• t_mean (float, Quantity, or Parameter, optional) – Mean formation time. If a Quantity/Parameter, interpreted as cosmic time (Gyr). If a bare float, interpreted as fractional mean m in (0,1). Must be provided together with kappa.

• kappa (float, Quantity, or Parameter, optional) – Concentration parameter (>0). Dimensionless. Used with t_mean.

• t_start (float, Quantity, or Parameter, optional) – Cosmic time at which star formation begins. Default is 0 Gyr.

• t_end (float, Quantity, or Parameter, optional) – Cosmic time at which the model stops. If None, uses today. Must satisfy t_end > t_start.

• ism_metallicity_today (float, Quantity, or Parameter, optional) – Constant ISM metallicity (dimensionless mass fraction). Default 0.02.

Notes

• All time arguments are cosmic time since the Big Bang (not lookback).

• ism_metallicity() returns a constant metallicity (no enrichment model).

name: str = 'beta_cem'

_resolve_shape_params() → None

Set internal dimensionless Parameters _alpha/_beta from either: - direct alpha/beta, or - mean+concentration (t_mean,kappa)

property alpha_val: float

Beta-CMF early-time shape parameter (alpha > 0).

property beta_val: float

Beta-CMF late-time shape parameter (beta > 0).

property kappa_val: float

Concentration parameter kappa = alpha + beta.

property t_mean_q: Quantity

Mean formation time E[t] implied by alpha/beta (in Gyr).

_alpha_beta_from_mean_concentration(t_mean: Parameter | Quantity | float, kappa: Parameter | Quantity | float) → tuple[float, float]

Map mean+concentration to alpha/beta.

Parameters:

• t_mean (float or Quantity or Parameter) – If float: fractional mean m in (0,1). If Quantity/Parameter: cosmic mean time in Gyr, converted to m=t_mean/t_end.

• kappa (float or Quantity or Parameter) – Concentration (>0), dimensionless.

Returns:

alpha, beta (float)

_rescaled_x(time: Quantity) → ndarray

Rescale cosmic time to x in [0,1] over [t_start, t_end], clipping outside.

stellar_mass_formed(time: Unit('Gyr')) → Quantity

Cumulative formed stellar mass M_formed(t).

ism_metallicity(time: Unit('Gyr')) → Quantity

Constant ISM metallicity (mass fraction).

sfr(time: Unit('Gyr')) → Quantity

Instantaneous SFR(t) (analytic Beta PDF mapped onto [t_start,t_end]).

Notes

• Outside (t_start, t_end) the SFR is set to 0 by construction.

• Integrating this over time yields mass_today.

t_peak() → Quantity | None

Cosmic time of peak SFR (exists only if alpha>1 and beta>1).

class pst.cem.TabularCEM(*, times: Parameter | u.Quantity | np.array, masses: Parameter | u.Quantity | np.array, metallicities: Parameter | u.Quantity | np.array, **kwargs)

Bases: ChemicalEvolutionModel

Chemical evolution model based on a grid of times and metallicities.

Description

This model represents the chemical evolution of a galaxy by means of a discrete grid of ages and metallicities

table_t

Tabulated cosmic time.

Type:

astropy.Quantity

table_M

Total stellar mass at each cosmic time step.

Type:

astropy.Quantity

table_metallicity

ISM metallicity at each cosmic time step.

Type:

astropy.Quantity

See also

pst.models.ChemicalEvolutionModel

name: str = 'tabular_cem'

property times

property masses

property metallicities

property table_t

property table_mass

property table_metallicity

stellar_mass_formed(times: Quantity)

Evaluate the integral of the SFR over a given set of times.

Description

This method evaluates the integral:

\[\int_{0}^{t} {\rm SFR}(t') dt'\]

at each time input time \(t\).

Parameters:

times (astropy.units.Quantity) – Array of cosmic times at which the integral will be evaluated.

returns:

integral (astropy.units.Quantity) – The cumulative stellar mass formed at each input time.

ism_metallicity(times: Quantity)

Evaluate the integral of the SFR over a given set of times.

Description

Return the metallicity Z(t) of the interstellar medium (gas and dust) at a certain set of cosmic times (i.e. since the Big Bang).

Parameters:

times (astropy.units.Quantity) – Cosmic times at which the metallicity will be evaluated.

returns:

z_t (astropy.units.Quantity) – Vector with the ISM metallicity at each input time.

class pst.cem.TabularCEM_ZPowerLaw(*, mass_today: Parameter, ism_metallicity_today: Parameter, alpha_powerlaw: Parameter, **kwargs)

Bases: MassPropMetallicityMixin, TabularCEM

Chemical evolution model based on a grid of times that assumes an analytic chemical enrichment history.

See also

TabularCEM, MassPropMetallicityMixin

name: str = 'tabular_zpowlaw_cem'

class pst.cem.CC25TabularCEM(*, tau_ssfr, ssfr, **kwargs)

Bases: TabularCEM_ZPowerLaw

Chemical evolution model based on the SFH parameterization from Corcho-Caballero et al. (2025).

The SFH is defined in terms of the average specific star formation rate (ssfr) several timescales (tau_ssfr), expressed in lookback time with respect to the age of the Universe at the time of observation, related by the expression:

\[M_\star(t) = M_{\rm today} \cdot \left(1 - \tau_{\rm ssfr} \cdot ssfr \right)\]

See also

TabularCEM_ZPowerLaw, MassPropMetallicityMixin

property ssfr

property tau_ssfr

class pst.cem.TabularMassFracCEM(*, mass_frac, **kwargs)

Bases: TabularCEM_ZPowerLaw

Chemical evolution model based on mass percentiles.

The SFH is defined in terms of the time at which a given fraction of the total stellar mass seen today was formed and normalized by the total mass at present.

Parameters:

• mass_frac (u.Quantity) – Cumulative mass fraction at a given time.

• times (u.Quantity) – Times associated to `mass_frac.

See also

TabularCEM_ZPowerLaw, MassPropMetallicityMixin

property mass_frac

property times

class pst.cem.ParticleListCEM(*, time_form: Parameter | Quantity | ndarray, metallicities: Parameter | Quantity | ndarray, masses: Parameter | Quantity | ndarray, today: Parameter | Quantity | float | None = None)

Bases: ChemicalEvolutionModel

Chemical Evolution Model using individual Simple Stellar Population (SSP) data.

This model represents the chemical evolution of a galaxy by reconstructing a composite stellar population (CSP) using individual stellar population (SSP) particles, each of which is defined by its formation time, metallicity, and mass.

Parameters:

• time_form (numpy.array or astropy.units.Quantity) – Array representing the formation times of each SSP particle. If the input is a numpy.array, it is assumed to be in Gyr. Otherwise, an astropy.units.Quantity with appropriate units is required.

• metallicities (numpy.array or astropy.units.Quantity) – Array representing the metallicities of each SSP particle. If the input is a numpy.array, it is assumed to be dimensionless (i.e., no units).

• masses (numpy.array or astropy.units.Quantity) – Array representing the masses of each SSP particle. If the input is a numpy.array, it is assumed to be in solar masses.

name: str = 'particle_list_cem'

property time_form

property metallicities

property masses

interpolate_ssp_masses(ssp: SSPBase, t_obs: Quantity)

Interpolate the SSP particles onto an SSP base model at the observed time.

This method interpolates the stellar masses of the SSP particles based on their ages and metallicities at the observed cosmic time t_obs using an SSP model grid.

Parameters:

• ssp (pst.SSP.SSPBase) – SSP model providing the ages and metallicities for interpolation.

• t_obs (astropy.units.Quantity) – The age of the Universe at the time of the observation.

Returns:

ssp_weights (astropy.units.Quantity) – A 2D array representing the stellar mass associated with each SSP particle in the base SSP grid. Units are in solar masses.

stellar_mass_formed(time)

Cumulative stellar mass formed as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

mstar (astropy.units.Quantity) – Cumulative stellar mass formed up to each input time. Units should be mass (typically Msun). Shape matches time.

ism_metallicity(time)

ISM metallicity (mass fraction) as a function of cosmic time.

Parameters:

time (astropy.units.Quantity) – Cosmic time(s) since the Big Bang.

Returns:

z_t (astropy.units.Quantity) – Dimensionless metallicity (mass fraction). Shape matches time.

Spectral Energy Distribution (SED) Components

class pst.sed.SedComponent

Bases: ModelBase, ABC

Base class for spectral energy distribution components.

Subclasses implement emission_spectrum() and return a spectrum sampled on the requested wavelength grid.

abstractmethod emission_spectrum(wavelength: Quantity, **params) → Quantity

Returns emission spectrum as Quantity.

The units are determined by your chosen normalization parameterization (e.g., L_IR scaling or amplitude in L_lambda). Keep this consistent with PST.

static integrate_sed(wavelength, sed, wl_min: Quantity = None, wl_max: Quantity = None)

Integrate a spectrum over wavelength with optional bounds.

Parameters:

• wavelength (astropy.units.Quantity) – Wavelength grid associated with sed.

• sed (astropy.units.Quantity) – Spectral density values sampled on wavelength.

• wl_min (astropy.units.Quantity, optional) – Lower integration limit.

• wl_max (astropy.units.Quantity, optional) – Upper integration limit.

Returns:

integral (astropy.units.Quantity) – Integrated quantity with units sed.unit * wavelength.unit.

classmethod q_ionising_photons(wavelength, sed)

Compute ionising photon production integrated below 912 Angstrom.

Parameters:

• wavelength (astropy.units.Quantity) – Wavelength grid.

• sed (astropy.units.Quantity) – Spectral energy density.

Returns:

q_ion (astropy.units.Quantity) – Integrated ionising photon rate-like quantity.

class pst.sed.TabularSedComponent

Bases: SedComponent

Base interface for SED components backed by tabulated data.

abstract property default_unit

Return the native output spectral unit for this component.

abstractmethod load_table()

Load the underlying tabulated model data.

abstractmethod emission_spectrum(wavelength: Quantity, **params) → Quantity

Returns emission spectrum as Quantity.

The units are determined by your chosen normalization parameterization (e.g., L_IR scaling or amplitude in L_lambda). Keep this consistent with PST.

class pst.sed.StellarComponent(ssp: SSPBase, sfh: ChemicalEvolutionModel, name: str = 'stellar_emission')

Bases: SedComponent

Stellar emission component built from an SSP model and an SFH/CEM.

ssp: SSPBase

sfh: ChemicalEvolutionModel

name: str = 'stellar_emission'

default_unit = Unit("solLum / Angstrom")

emission_spectrum(wavelength: Quantity, **params)

Compute stellar emission on the requested wavelength grid.

Parameters:

• wavelength (astropy.units.Quantity) – Target wavelength grid.

• **params – Parameters forwarded to sfh.compute_SED.

Returns:

sed (astropy.units.Quantity) – Stellar SED in default_unit.

Galaxy Assembly and Composite SEDs

class pst.galaxy.GalaxySED(*, stellar_model: SedComponent = None, dust_attenuation_model: AttenuationModel = None, dust_model: SedComponent = None, target_wavelength: Quantity = None, redshift: Parameter | float = 0.0, cosmology=None, filters: FilterList = None)

Bases: SedComponent

Composite galaxy SED model (stellar emission + attenuation + dust emission).

This component combines: - a stellar emission component (typically driven by a CEM/SSP model), - an optional attenuation model applied to the stellar emission, and - an optional dust emission component (optionally calorimetric / energy-balance).

Parameters are handled via pst.model.Parameter objects and may be updated at runtime using dotted parameter paths (see build_param_index() and update_parameters()), enabling integration with samplers/fitters such as BESTA.

Notes

• The internal spectral grid is target_wavelength in the rest frame.

• If to_obs_frame=True, spectra are converted to observed-frame fluxes using luminosity distance and redshift and resampled onto lambda_obs and target_wavelength is treated as the observed frame.

name: str = 'galaxy_sed'

rest_unit = Unit("solLum / Angstrom")

obs_unit = Unit("erg / (Angstrom s cm2)")

property redshift: Parameter

Cosmological redshift parameter used for observation-frame mapping.

build_param_index(*, include_fixed: bool = True, prefix: str = '') → Dict[str, Parameter]

Build and cache a mapping from parameter path to Parameter object.

Parameters:

• include_fixed (bool, optional) – If True include fixed parameters in the index.

• prefix (str, optional) – Optional path prefix.

Returns:

index (dict) – Mapping from dotted parameter path to Parameter.

update_parameters(values: Mapping[str, Any], *, validate: bool = True, strict: bool = True) → None

Update parameters from a dict.

Parameters:

• values (mapping) – Mapping from dotted parameter path to value. Values may be floats or Quantities.

• validate (bool, optional) – If True, checks fixed and vrange constraints.

• strict (bool, optional) – If True, unknown keys raise KeyError. If False, unknown keys are ignored.

Notes

This uses a cached parameter index when available. If the cache has not been built, it is built on the fly.

emission_components(stars_params=None, dust_att_params=None, dust_em_params=None)

Compute individual SED components on the rest-frame wavelength grid.

Returns a dict with (when available): - stellar_sed_unatt : intrinsic stellar emission (no attenuation) - stellar_sed : attenuated stellar emission (if attenuation model is set) - dust_sed : dust emission component (if present)

If self.dust is calorimetric (energy-balance), the dust component is computed consistently from the absorbed stellar energy.

emission_spectrum(stars_params=None, dust_att_params=None, dust_em_params=None, to_obs_frame=False)

Compute the composite SED (stellar + dust).

Parameters:

• stars_params (dict, optional) – Parameters forwarded to the stellar component.

• dust_att_params (dict, optional) – Parameters forwarded to the attenuation model.

• dust_em_params (dict, optional) – Parameters forwarded to the dust emission component.

• to_obs_frame (bool) – If True, returns observed-frame flux density and shifts/resamples the spectrum to lambda_obs = (1+z) lambda_rest.

Returns:

sed (astropy.units.Quantity) – Composite spectrum in rest_unit (rest frame) or obs_unit (observed frame when to_obs_frame=True).

emission_photometry(stars_params=None, dust_att_params=None, dust_em_params=None, to_obs_frame=False)

Compute synthetic photometry from the composite spectrum.

Parameters:

• stars_params (dict, optional) – Parameters forwarded to the stellar component.

• dust_att_params (dict, optional) – Parameters forwarded to the attenuation model.

• dust_em_params (dict, optional) – Parameters forwarded to the dust emission component.

• to_obs_frame (bool, optional) – If True, compute photometry from observed-frame fluxes.

Returns:

photometry (astropy.units.Quantity) – Band-integrated synthetic flux densities for configured filters.

Notes

Requires filters to be provided at initialization (or set later). The filter container must support evaluation of fluxes from a spectrum sampled on target_wavelength.

Dust Attenuation and Emission

Dust Extinction/Emission Models for Stellar Populations

This module implements dust extinction and emission models for stellar population synthesis, including a base class for general dust models and specific implementations such as a dust screen and the Charlot & Fall (2000) extinction model.

Usage

This module is intended for applying dust extinction or emission to stellar spectra, either to synthetic stellar populations or other types of spectra.

pst.dust._exp_bounded(exponent: ndarray) → ndarray

Exponentiate with clipping to float64 finite range.

pst.dust._pow_positive_bounded(base: ndarray, exponent: float) → ndarray

Compute base**exponent for positive bases using a log-space bounded form.

pst.dust.modified_blackbody(lam: Quantity, T: Quantity, beta: float, lam_0: Quantity | None = None, lam_ref: Quantity = <Quantity 100. um>, per_freq: bool = False) → Quantity

Compute a modified blackbody spectrum.

This function returns a spectral shape given by a blackbody multiplied by an emissivity term. Two emissivity parameterizations are supported.

If lam_0 is provided, the emissivity uses a finite optical depth form

\[S_{nu} \propto (1 - \exp(-\tau_{nu})) B_{nu}(T)\]

with

\[\tau_{nu} = \left(\frac{\nu}{\nu_0}\right)^{\beta} = \left(\frac{\lambda_0}{\lambda}\right)^{\beta}\]

If lam_0 is None, the emissivity uses an optically thin approximation

\[S_{nu} \propto \left(\frac{\nu}{\nu_{ref}}\right)^{\beta} B_{nu}(T)\]

Parameters:

• lam (astropy.units.Quantity) – Wavelength grid. Internally converted to microns.

• T (astropy.units.Quantity) – Dust temperature. Internally converted to Kelvin.

• beta (float) – Emissivity index.

• lam_0 (astropy.units.Quantity, optional) – Turnover wavelength where tau equals 1. If provided, the finite optical depth emissivity term is used.

• lam_ref (astropy.units.Quantity, optional) – Reference wavelength used to normalize the optically thin emissivity term. Only used when lam_0 is None. Default is 100 micron.

• per_freq (bool, optional) – If True, return a per-frequency spectrum proportional to B_nu. If False, return a per-wavelength spectrum using a spectral density equivalency. Default is False.

Returns:

mbb (astropy.units.Quantity) – Modified blackbody spectrum evaluated on lam. The normalization is not set by this function. If per_freq is True the result is proportional to B_nu(T). If per_freq is False the result is returned as a per-wavelength density.

Notes

This function is typically used to build a template shape that is later normalized to an integrated luminosity.

The per-wavelength conversion uses astropy spectral density equivalencies.

class pst.dust.AttenuationCurve

Bases: ModelBase

Base class for attenuation curve shapes.

This class defines the wavelength dependence of an attenuation law via a dimensionless curve k_lambda. The normalization (for example an A_V value) is applied separately by scaling k_lambda.

Subclasses must implement k_lambda.

Notes

This interface intentionally separates shape from normalization. The helper methods provide common conversions:

• a_lambda: convert k_lambda to attenuation in magnitudes given a_v

• tau_lambda: convert magnitudes to optical depth using A = 1.086 tau

• attenuation_factor: compute 10**(-0.4 A_lambda)

name

Curve identifier.

Type:

str

name: str = 'attenuation_curve'

abstractmethod k_lambda(wavelength: Quantity, **params) → Quantity

Return the dimensionless attenuation curve shape.

Parameters:

• wavelength (astropy.units.Quantity) – Wavelength grid where the curve is evaluated.

• **params – Optional call time parameters. Most implementations store curve parameters as Parameter attributes, so this is often unused.

Returns:

k (astropy.units.Quantity) – Dimensionless curve with the same shape as wavelength.

a_lambda(wavelength: Quantity, *, a_v: Quantity, **params) → Quantity

Convert curve shape to attenuation in magnitudes.

The attenuation is computed as

\[A_{\lambda} = A_{v} k_{\lambda}\]

Parameters:

• wavelength (astropy.units.Quantity) – Wavelength grid where the attenuation is evaluated.

• a_v (astropy.units.Quantity) – Attenuation normalization in magnitudes. The interpretation depends on the curve convention. Commonly, a_v corresponds to attenuation at the pivot wavelength used to normalize k_lambda.

• **params – Parameters forwarded to k_lambda.

Returns:

A_lambda (astropy.units.Quantity) – Attenuation in magnitudes with the same shape as wavelength.

tau_lambda(wavelength: Quantity, *, a_v: float = None, **params) → Quantity

Convert attenuation in magnitudes to optical depth.

Uses the relation

\[A_{\lambda} = 1.086 \tau_{\lambda}\]

Parameters:

• wavelength (astropy.units.Quantity) – Wavelength grid where the optical depth is evaluated.

• a_v (astropy.units.Quantity) – Attenuation normalization in magnitudes.

• **params – Parameters forwarded to a_lambda.

Returns:

tau (astropy.units.Quantity) – Optical depth as a dimensionless quantity with the same shape as wavelength.

attenuation_factor(wavelength: Quantity, *, a_v: Quantity, **params) → Quantity

Compute the multiplicative attenuation factor.

The factor is defined as

\[f_{\lambda} = 10^{-0.4 A_{\lambda}}\]

Parameters:

• wavelength (astropy.units.Quantity) – Wavelength grid where the attenuation is evaluated.

• a_v (astropy.units.Quantity) – Attenuation normalization in magnitudes.

• **params – Parameters forwarded to a_lambda.

Returns:

factor (astropy.units.Quantity) – Dimensionless attenuation factor with the same shape as wavelength.

class pst.dust.PowerLawAttenuationCurve(*, name: str = 'powerlaw_attenuation_curve', alpha: Parameter = <factory>, pivot: Parameter = <factory>)

Bases: AttenuationCurve

Power law attenuation curve shape.

This curve defines the dimensionless shape

\[k_{\lambda} = \left(\frac{\lambda}{\lambda_{pivot}}\right)^{\alpha}\]

where lambda_pivot is the pivot wavelength and alpha is the power law exponent.

Parameters:

• name (str, optional) – Curve identifier.

• alpha (pst.model.Parameter) – Power law exponent. Commonly negative for dust attenuation.

• pivot (pst.model.Parameter) – Pivot wavelength used to normalize the curve so that k_lambda at pivot is 1.

Notes

This class returns k_lambda. The normalization A_v is applied by AttenuationCurve.a_lambda or by an AttenuationModel.

See also

AttenuationCurve

Base interface for attenuation curve shapes.

name: str = 'powerlaw_attenuation_curve'

alpha: Parameter

pivot: Parameter

__post_init__()

Validate curve parameters.

Raises:

ValueError – If alpha is positive.

k_lambda(wavelength: Quantity, **params) → Quantity

Evaluate the curve shape.

Parameters:

• wavelength (astropy.units.Quantity) – Wavelength grid where the curve is evaluated.

• **params – Optional call time parameters. Reserved for future extensions.

Returns:

k (astropy.units.Quantity) – Dimensionless curve with the same shape as wavelength.

class pst.dust.ExtinctionLibCurve(*, name: str = 'extinction_lib', law: str = 'ccm89', r_v: Parameter = <factory>)

Bases: AttenuationCurve

Attenuation curve shape based on the extinction package.

This class wraps an extinction law implemented by the extinction package and returns a dimensionless curve k_lambda computed from the law evaluated at A_v equals 1 magnitude.

Parameters:

• name (str, optional) – Curve identifier.

• law (str, optional) – Name of the extinction law function in the extinction package.

• r_v (pst.model.Parameter) – Total to selective extinction ratio.

Raises:

ValueError – If the requested law is not present in the extinction package.

Notes

The extinction package functions typically accept wavelength, a_v, and r_v and return A_lambda in magnitudes. This class calls the selected law with a_v set to 1 and returns the resulting A_lambda as the dimensionless curve shape.

If you want k_lambda normalized such that k_lambda at a pivot wavelength is 1, normalize externally using a_v or implement a pivot normalization in this class.

See also

AttenuationCurve

Base interface for attenuation curve shapes.

name: str = 'extinction_lib'

law: str = 'ccm89'

r_v: Parameter

__post_init__()

Resolve the extinction law from the extinction package.

Raises:

ValueError – If the law is not found.

k_lambda(wavelength: Quantity, **params) → Quantity

Evaluate the curve shape using the selected extinction law.

Parameters:

• wavelength (astropy.units.Quantity) – Wavelength grid where the curve is evaluated.

• **params – Optional call time parameters. Reserved for future extensions.

Returns:

k (astropy.units.Quantity) – Dimensionless curve with the same shape as wavelength.

class pst.dust.AttenuationModel

Bases: ModelBase, ABC

Base class for geometry dependent attenuation models.

An attenuation model returns a multiplicative attenuation factor that can be applied to spectra. Models represent different geometries such as a foreground screen or multi component prescriptions.

Subclasses must implement attenuation_factor.

Notes

The apply method multiplies the input spectra by the attenuation factor, broadcasting the factor along the wavelength axis.

name

Model identifier.

Type:

str

name: str = 'attenuation_model'

abstractmethod attenuation_factor(wavelength: Quantity, **params) → Quantity

Compute the multiplicative attenuation factor.

Parameters:

• wavelength (astropy.units.Quantity) – Wavelength grid where the attenuation factor is evaluated.

• **params – Optional call time parameters. Most implementations store parameters as Parameter attributes, so this is often unused.

Returns:

factor (astropy.units.Quantity) – Dimensionless attenuation factor. The shape is model dependent. For a single screen this is (n_wave,). For a multi component model this may be (n_comp, n_wave).

apply(wavelength: Quantity, spectra, axis: int = -1, **params)

Apply attenuation to spectra.

Parameters:

• wavelength (astropy.units.Quantity) – Wavelength grid corresponding to the spectra.

• spectra (array-like or astropy.units.Quantity) – Input spectra to be attenuated.

• axis (int, optional) – Axis in spectra corresponding to wavelength. Default is -1.

• **params – Parameters forwarded to attenuation_factor.

Returns:

attenuated (same type as spectra) – Spectra multiplied by the attenuation factor, broadcast to the selected axis.

Notes

If attenuation_factor returns multiple components, the returned factor cannot be directly applied to a single spectrum. In that case, the caller should select or combine components before calling apply.

class pst.dust.DustScreenAttenuation(name: str = 'dust_screen', curve: AttenuationCurve = <factory>, a_v: Parameter = <factory>)

Bases: AttenuationModel

Foreground dust screen attenuation model.

This model applies a single attenuation curve with a single normalization parameter A_v.

Parameters:

• name (str, optional) – Model identifier.

• curve (AttenuationCurve, optional) – Attenuation curve shape used by the model.

• a_v (pst.model.Parameter) – Attenuation normalization in magnitudes.

Notes

The attenuation factor is computed as

\[f_{\lambda} = 10^{-0.4 A_{v} k_{\lambda}}\]

where k_lambda is provided by the curve.

See also

AttenuationCurve

Curve interface used to compute k_lambda.

name: str = 'dust_screen'

curve: AttenuationCurve

a_v: Parameter

__post_init__()

Normalize the curve input.

Notes

If curve is provided as a string, it is interpreted as an extinction package law name and converted into an ExtinctionLibCurve.

attenuation_factor(wavelength: Quantity, **params) → Quantity

Compute the attenuation factor for a foreground screen.

Parameters:

• wavelength (astropy.units.Quantity) – Wavelength grid where the attenuation is evaluated.

• a_v (float, optional) – Unused. The model uses the internal Parameter a_v. This argument is kept for backward compatibility with older call sites.

• **params – Additional parameters forwarded to the curve. Reserved for future extensions.

Returns:

factor (astropy.units.Quantity) – Dimensionless attenuation factor with the same shape as wavelength.

class pst.dust.CharlotFall00Attenuation(*, name: str = 'CF00', curve_young: AttenuationCurve = <factory>, curve_old: AttenuationCurve = <factory>, a_v_young: Parameter = <factory>, a_v_old: Parameter = <factory>, young_age: Parameter = <factory>)

Bases: AttenuationModel

Two component attenuation model inspired by Charlot and Fall 2000.

The model returns two attenuation factors corresponding to young and old stellar populations.

Parameters:

• name (str, optional) – Model identifier.

• curve_young (AttenuationCurve) – Curve used for young component.

• curve_old (AttenuationCurve) – Curve used for old component.

• a_v_young (pst.model.Parameter) – Normalization for young component in magnitudes.

• a_v_old (pst.model.Parameter) – Normalization for old component in magnitudes.

• young_age (pst.model.Parameter) – Age threshold used by callers that split stellar emission into components.

Returns:

factor (astropy.units.Quantity) – Dimensionless array with shape (2, n_wave). The first row corresponds to young component and the second row corresponds to old component.

Notes

This method returns factors only. The logic that splits stellar emission into young and old components is implemented elsewhere, for example in the stellar emission component or in a calorimetric dust model.

See also

CalorimetricDustComponent

Uses the returned factors to compute absorbed luminosity by component.

name: str = 'CF00'

curve_young: AttenuationCurve

curve_old: AttenuationCurve

a_v_young: Parameter

a_v_old: Parameter

young_age: Parameter

attenuation_factor(wavelength: Quantity, **params) → Quantity

Compute attenuation factors for young and old components.

Parameters:

• wavelength (astropy.units.Quantity) – Wavelength grid where the attenuation is evaluated.

• a_v_young (float, optional) – V band attenuation in magnitudes for the young component. Default is 1.0.

• a_v_old (float, optional) – V band attenuation in magnitudes for the old component. Default is 0.3.

• **params – Additional parameters forwarded to each attenuation curve, such as r_v.

Returns:

factor (astropy.units.Quantity) – Dimensionless attenuation factors with shape (2, n_wave). The first row corresponds to the young component and the second row corresponds to the old component.

class pst.dust.Casey2012DustComponent(*, name: str = 'Casey2012', t_dust: Parameter = <factory>, beta: Parameter = <factory>, alpha: Parameter = <factory>, lam0: Parameter = <factory>, min_wavelength: Quantity = <Quantity 1. um>, ir_range: Quantity, ~astropy.units.quantity.Quantity]=(<Quantity 8. um>, <Quantity 1000. um>))

Bases: SedComponent

Dust emission component inspired by Casey 2012.

This component builds an unnormalized template shape as the sum of a modified blackbody and a mid infrared power law. The shape is then normalized so that the integrated luminosity over a configurable infrared wavelength range matches lum_ir.

Parameters:

• name (str, optional) – Component identifier.

• t_dust (pst.model.Parameter) – Dust temperature.

• beta (pst.model.Parameter) – Emissivity index.

• alpha (pst.model.Parameter) – Mid infrared power law exponent.

• lam0 (pst.model.Parameter) – Optical depth scale wavelength.

• min_wavelength (astropy.units.Quantity, optional) – Emission below this wavelength is set to zero.

• ir_range (tuple of astropy.units.Quantity, optional) – Integration bounds used to normalize the template.

• default_unit (astropy.units.Unit, optional) – Output spectral density unit.

Notes

The output is a spectral density in default_unit. The overall normalization is set by lum_ir.

See also

modified_blackbody

Function used to compute the modified blackbody component.

name: str = 'Casey2012'

t_dust: Parameter

beta: Parameter

alpha: Parameter

lam0: Parameter

min_wavelength: Quantity = <Quantity 1. um>

ir_range: Tuple[Quantity, Quantity] = (<Quantity 8. um>, <Quantity 1000. um>)

default_unit = Unit("solLum / Angstrom")

emission_spectrum(wavelength: Quantity, *, lum_ir: Quantity, lam_pivot: Quantity | None = None, **kwargs) → Quantity

Compute the dust emission spectrum.

Parameters:

• wavelength (astropy.units.Quantity) – Wavelength grid where the emission is evaluated.

• lum_ir (astropy.units.Quantity) – Target integrated infrared luminosity used for normalization. Expected units are luminosity, for example Lsun.

• **kwargs – Additional unused parameters for compatibility.

Returns:

l_lambda (astropy.units.Quantity) – Dust emission spectrum in default_unit, sampled on wavelength.

Raises:

ValueError – If the template has zero integrated luminosity in the IR range and cannot be normalized.

_lambda_pivot(t_dust: Quantity, alpha: float) → Quantity

Compute the pivot wavelength used to connect template components.

Parameters:

• t_dust (astropy.units.Quantity) – Dust temperature.

• alpha (float) – Mid infrared power law slope.

Returns:

lam_pivot (astropy.units.Quantity) – Pivot wavelength in microns.

class pst.dust.CalorimetricDustComponent(attenuation: AttenuationModel, dust_sed_component: SedComponent, name: str = 'calorimetric_dust_sed')

Bases: SedComponent

Calorimetric dust emission component.

This component couples an attenuation model and a dust emission component using energy balance. It computes the stellar luminosity absorbed by dust and uses it as the infrared luminosity normalization for the dust emission component.

Parameters:

• attenuation (AttenuationModel) – Attenuation model used to compute attenuation factors.

• dust_sed_component (SedComponent) – Dust emission component that accepts lum_ir as normalization.

• default_unit (astropy.units.Unit, optional) – Output spectral density unit.

• name (str, optional) – Component identifier.

Returns:

• Lsrc (astropy.units.Quantity) – Intrinsic stellar spectrum.

• Latt (astropy.units.Quantity) – Attenuated stellar spectrum.

• Ldust (astropy.units.Quantity) – Dust emission spectrum normalized to absorbed luminosity.

Notes

If the attenuation model returns multiple factors (for example two component CF00), the stellar source must return a matching binned output.

See also

CharlotFall00Attenuation

Two component attenuation model.

Casey2012DustComponent

Dust emission model normalized by lum_ir.

attenuation: AttenuationModel

dust_sed_component: SedComponent

default_unit = Unit("solLum / Angstrom")

name: str = 'calorimetric_dust_sed'

emission_spectrum(wavelength: Quantity, *, source: StellarComponent, source_params: dict = None, **params) → Tuple[Quantity, Quantity, Quantity]

Compute intrinsic, attenuated, and dust emission spectra.

Parameters:

• wavelength (astropy.units.Quantity) – Wavelength grid where spectra are evaluated.

• source (StellarComponent) – Stellar emission source.

• source_params (dict, optional) – Parameters forwarded to the stellar source emission_spectrum.

• **params – Parameters forwarded to the attenuation model and the dust emission component. For example a_v or a_v_young and a_v_old.

Returns:

• Lsrc (astropy.units.Quantity) – Intrinsic stellar spectrum in default_unit.

• Latt (astropy.units.Quantity) – Attenuated stellar spectrum in default_unit.

• Ldust (astropy.units.Quantity) – Dust emission spectrum in default_unit normalized to absorbed energy.

Raises:

ValueError – If CharlotFall00Attenuation is used and the stellar source does not return a 2D binned spectrum with shape (n_bins, n_wave).

class pst.dust.DustModelBase

Bases: ABC

Abstract base class for dust extinction and emission models.

Description

This class provides the framework for implementing dust models. Subclasses should define methods to compute the extinction and emission due to dust.

extinction_law

The name of the extinction law to be used. This is retrieved from the extinction library.

Type:

str

abstractmethod get_extinction(*args, **kwargs)

Compute the dust extinction for a given set of parameters.

This method must be implemented by subclasses.

abstractmethod get_emission(*args, **kwargs)

Compute the dust emission for a given set of parameters.

This method must be implemented by subclasses.

apply_extinction(wavelength, spectra, axis=-1, **kwargs)

Apply the dust extinction model to a given spectra.

Parameters:

• wavelength (np.ndarray or astropy.Quantity) – Wavelength array. Can be either a numpy array of floats or an astropy.Quantity with associated units (e.g., Angstroms).

• spectra (np.ndarray or astropy.Quantity) – Array of spectra to which the extinction will be applied.

• axis (int, optional) – The axis of the spectra array corresponding to the wavelength dimension. Default is -1.

• **kwargs – Additional keyword arguments passed to the get_extinction method.

Returns:

reddened_spectra (np.ndarray) – The spectra array with dust extinction applied.

apply_emission(wavelength, spectra, axis=-1, **kwargs)

Add the predicted dust emission to a given spectra.

Parameters:

• wavelength (np.ndarray or astropy.Quantity) – Wavelength array. Can be either a numpy array of floats or an astropy.Quantity with associated units.

• spectra (np.ndarray or astropy.Quantity) – Array of spectra to which the dust emission will be added.

• axis (int, optional) – The axis of the spectra array corresponding to the wavelength dimension. Default is -1.

• **kwargs – Additional keyword arguments passed to the get_emission method.

Returns:

spectra_with_emission (np.ndarray) – The spectra array with dust emission added.

redden_ssp_model(ssp_model, **kwargs)

Apply extinction to a simple stellar population (SSP) model.

Parameters:

• ssp_model (pst.SSPBase object) – A simple stellar population (SSP) model instance.

• **kwargs – Additional keyword arguments passed to the apply_extinction method.

Returns:

reddened_ssp_model (pst.SSPBase object) – The SSP model with dust extinction applied.

class pst.dust.DustScreen(extinction_law_name, r_extinction=3.1)

Bases: DustModelBase

Dust screen extinction model.

Implements a simple dust screen model where dust extinction is applied to spectra using a specified extinction law and R_V parameter.

extinction_law_name

The name of the extinction law from the extinction library (e.g., ‘ccm89’, ‘odonnell94’).

Type:

str

r_extinction

The R_V value for the extinction law. Default is 3.1.

Type:

float

get_extinction(wavelength, a_v=1.0)

Compute the dust extinction.

Parameters:

• wavelength (np.ndarray or astropy.Quantity) – Wavelength array in Angstroms.

• a_v (float, optional) – The V-band extinction (in magnitudes). Default is 1.0.

Returns:

extinction_curve (np.ndarray) – Dimensionless extinction factor to be applied to the spectra.

get_emission(wavelength)

Compute the dust emission.

For this model, no dust emission is included, so this method returns zeros.

Parameters:

wavelength (np.ndarray or astropy.Quantity) – Wavelength array.

Returns:

emission (np.ndarray) – An array of zeros with the same shape as wavelength.

class pst.dust.CF03DustScreen(extinction_law_name, young_ssp_age, r_extinction=3.1)

Bases: DustScreen

Charlot & Fall (2000) dust screen model for young and old stellar populations.

This model applies different extinction curves to young and old populations based on their ages.

Parameters:

• extinction_law_name (str) – The name of the extinction law from the extinction library.

• young_ssp_age (astropy.Quantity) – The age threshold for defining young populations (in years).

• r_extinction (float, optional) – The R_V value for the extinction law. Default is 3.1.

get_extinction(wavelength, age, a_v_young=1.0, a_v_old=0.3)

Compute the dust extinction for young and old stellar populations.

Parameters:

• wavelength (np.ndarray or astropy.Quantity) – Wavelength array.

• age (np.ndarray or astropy.Quantity) – Array of stellar population ages.

• a_v_young (float, optional) – V-band extinction for young populations. Default is 1.0.

• a_v_old (float, optional) – V-band extinction for old populations. Default is 0.3.

Returns:

extinction_curve (np.ndarray) – 2D array of extinction factors with shape (age.size, wavelength.size).

Observables, Filters, and Spectral Indices

This module contains some tools for computing observable quantities (e.g. photometry, equivalent widths) from spectra.

pst.observables.list_of_available_filters()

List the currently available filters in the default directory.

pst.observables.load_photometric_filters(filters, to_filter_list=False)

Convenience function for constructing a list of photometric filters.

Parameters:

filters (list of str) – List of filters to load. The list might contain the absolute path to a filter response file, or just the filter name following the SVO convention.

Returns:

filters_out (list of pst.observables.Filter) – List of filters.

pst.observables.download_svo_filter(name: str, dest_dir: str, verbose=True, retry=3)

Download a filter from the Spanish Virtual Observatory (SVO) Filter Profile Service.

Parameters:

• name (str) – SVO-compliant filename. The naming convention for SVO filters is TELESC_INSTRUMENT.BAND (e.g. WISE_WISE.W1, Subaru_HSC.g)

• dest_dir (str) – Path to the directory where to store the data.

Returns:

file_path (str) – Path to the downloaded filter file.

class pst.observables.Filter(wavelength=None, response=None, filter_wavelength=None, filter_response=None, name=None)

Bases: object

A photometric filter.

filter_resp

Original photometric passband response curve.

Type:

np.ndarray

filter_wavelength

Original wavelength associated to filter_resp. If type is numpy.ndarray, the value is converted in to an ``astropy.units.Quantity` in Angstrom.

Type:

numpy.ndarray or astropy.units.Quantity

response

Filter passband response curve after interpolation.

Type:

np.ndarray

wavelength

Wavelength vector associated to response. If type is numpy.ndarray, the value is converted in to an ``astropy.units.Quantity` in Angstrom.

Type:

numpy.ndarray or astropy.units.Quantity

default_dir

Default directory containing filter files.

Type:

str

Example

>>> from pst.observables import Filter
>>> ps_r_filter = Filter("PANSTARRS_PS1.r")
>>> wl = np.linspace(5000, 7000, 1000) * u.angstrom
>>> ps_r_filter.interpolate(wl)
>>> ps_r_filter.plot(add_props=True, show=True)

default_dir = '/home/docs/checkouts/readthedocs.org/user_builds/population-synthesis-toolkit/checkouts/latest/src/pst/data/filters'

property wavelength

Interpolated wavelength grid used by response.

property filter_wavelength

Native wavelength grid associated with filter_resp.

classmethod from_text_file(path, wavelength_unit=Unit('Angstrom'), **kwargs)

Load a :class:Filter from an input text file.

Parameters:

• path (str) – Path to the text file containing the filter information. The first and second columns must correspond to the wavelength and passband curve, respectively.

• wavelength_unit (:class:astropy.units.Unit) – Unit associated with wavelength values in the input file.

• **kwargs – Arguments to be passed to numpy.loadtxt()

Returns:

filter (Filter) – The Filter containing the input information.

classmethod from_svo(name, destination_dir=None)

Load a :class:Filter from the Spanish Vitural Observatory archive.

Parameters:

name (str) – SVO filter name. If the filter is not found locally, it will be downloaded from the archive.

Returns:

filter (Filter) – The Filter containing the input information.

Example

>>> from pst.observables import Filter
>>> panstarrs_r_filter = Filter.from_svo("PANSTARRS_PS1.r")

effective_wavelength()

Compute the effective wavelength of the filter.

Description

The effective wavelength is computed as

\[\lambda_{\rm eff} = \frac{\int{R(\lambda) \cdot \lambda d\lambda}}{\int{R(\lambda) d\lambda}}\]

returns:

eff_wl (astropy.units.Quantity) – The effective wavelength of the filter.

effective_bandwidth()

Compute the effective bandwidth of the filter.

Description

The effective bandwith is computed as

\[\Delta \lambda_{\rm BW} = \sqrt{8\log(2)} \left(\frac{\int{R(\lambda) \cdot \lambda^2 d\lambda}}{\int{R(\lambda) d\lambda}} - \lambda_{\rm eff}\right)^{1/2}\]

returns:

eff_bw (astropy.units.Quantity) – The effective bandwidth of the filter.

See also

effective_wavelength()

effective_transmission()

Compute the effective bandwidth of the filter.

Description

The effective transmission is computed as

\[R_{\rm eff} = \frac{\int{R(\lambda)^2 d\lambda}}{\int{R(\lambda) d\lambda}}\]

returns:

eff_tr (float) – The effective transmission of the filter.

interpolate(wavelength=None)

Interpolate and update the Filter response curve to an input wavelength.

Description

Interpolate linearly the Filter response curve to an input wavelength vector. The result will update the exising values of wavelength and response.

Parameters:

wavelength (numpy.ndarray or astropy.units.Quantity) – Wavelength vector to interpolate filt_resp. If type is numpy.ndarray, the value is converted in to an ``astropy.units.Quantity` in Angstrom.

returns:

response (np.ndarray) – Filter response curve interpolated to the input values of wavelength.

get_photons(spectra, spectra_err=None, mask_nan=True)

Compute the photon flux from an input spectra.

Description

The photon flux associated to the filter is computed by numerically integrating the input spectra with the filter response, using the trapezid method.

Parameters:

• spectra (np.ndarray or :class:astropy.units.Quantity) – Input spectra (flux density per wavelength unit) with same dimensions as the Filter wavelength.

• spectra_err (np.ndarray or :class:astropy.units.Quantity, optional) – Input spectra associated error.

• mask_nan (bool, optional) – If True, NaN values are masked.

returns:

• photon_flux (:class:astropy.units.Quantity) – Filter photon flux.

• photon_flux_err (:class:astropy.units.Quantity) – Filter photon flux associated error.

get_ab(spectra, spectra_err=None, mask_nan=True)

Compute the synthetic AB magnitude from an input spectra.

Description

The AB magnitude associated to the filter is computed by numerically integrating the input spectra with the filter response, using the trapezid method.

Parameters:

• spectra (np.ndarray or :class:astropy.units.Quantity) – Input spectra (flux density per wavelength unit) with same dimensions as the Filter wavelength.

• spectra_err (np.ndarray or :class:astropy.units.Quantity, optional) – Input spectra associated error.

• mask_nan (bool, optional) – If True, NaN values are masked.

returns:

• mag_ab (:class:astropy.units.Quantity) – AB magnitude.

• mag_ab_err (:class:astropy.units.Quantity) – AB magnitude associated error.

See also

get_photons()

get_fnu(spectra, spectra_err=None, mask_nan=True)

Compute synthetic flux density per frequency unit from a spectrum.

Parameters:

• spectra (np.ndarray or :class:astropy.units.Quantity) – Input spectra (flux density per wavelength unit) with same dimensions as the Filter wavelength.

• spectra_err (np.ndarray or :class:astropy.units.Quantity, optional) – Input spectra associated error.

• mask_nan (bool, optional) – If True, NaN values are masked.

Returns:

• f_nu (:class:astropy.units.Quantity) – Synthetic flux density in Jy.

• f_nu_err (:class:astropy.units.Quantity) – Associated uncertainty in Jy.

See also

get_photons()

get_flambda_vegamag(spectra, spectra_err=None, mask_nan=True)

Compute synthetic flux density per wavelength unit from a spectrum.

Parameters:

• spectra (np.ndarray or :class:astropy.units.Quantity) – Input spectra (flux density per wavelength unit) with same dimensions as the Filter wavelength.

• spectra_err (np.ndarray or :class:astropy.units.Quantity, optional) – Input spectra associated error.

• mask_nan (bool, optional) – If True, NaN values are masked.

Returns:

• f_lambda (:class:astropy.units.Quantity) – Synthetic flux density per wavelength.

• f_lambda_err (:class:astropy.units.Quantity) – Associated uncertainty estimate.

See also

get_photons()

plot(add_props=False, ax=None, show=False)

Plot the filter response curve.

Plot the original filter response curve together with the interpolated version computed using a new grid of wavelengths.

Parameters:

• add_props (bool) – If True, add vertical lines indicating the effective wavelength and bandwidth of the filter. Default is False.

• ax (matplotlib.axes.Axes, optional) – Matplotlib axis to plot on. If None, a new figure and axis are created.

• show (bool) – If True, display the plot by calling plt.show(). This requires an interactive matplotlib session. Default is False.

class pst.observables.TopHatFilter(central_wave, width, **kwargs)

Bases: Filter

Top hat photometric filter

See also

Filter

class pst.observables.FilterList(filters: List[Filter], wavelength: Quantity | None = None, response: ndarray | None = None, kernel_phot: Quantity | None = None, dlambda: Quantity | None = None, norm_phot: Quantity | None = None, names: List[str] | None = None)

Bases: object

A list of Filter instances for faster computations.

#TODO

filters: List[Filter]

wavelength: Quantity | None = None

response: ndarray | None = None

kernel_phot: Quantity | None = None

dlambda: Quantity | None = None

norm_phot: Quantity | None = None

names: List[str] | None = None

property n_bands: int

Number of filters contained in the list.

property effective_wavelength: Quantity

Array of effective wavelengths, one per filter.

wavelength_range(kappa_bw=2.0) → [<class 'astropy.units.quantity.Quantity'>, <class 'astropy.units.quantity.Quantity'>]

Get the net wavelength coverage by the filters.

interpolate(wavelength: ndarray | Quantity) → FilterList

Interpolate all filter responses to a common wavelength grid.

Parameters:

wavelength (array or Quantity) – Target wavelength grid. Must be 1D and monotonic increasing.

Returns:

self

get_photons(spectra: ndarray | Quantity, spectra_err: ndarray | Quantity | None = None, mask_nan: bool = True) → Tuple[Quantity, Quantity | None]

Vectorized photon flux through all bandpasses.

Parameters:

• spectra (array or Quantity) – Shape (…, n_wave). Flux density per wavelength.

• spectra_err (array or Quantity, optional) – Same shape as spectra.

• mask_nan (bool) – If True, NaNs in spectra are treated as zero contribution.

Returns:

• n_phot (Quantity) – Shape (…, n_bands)

• n_phot_err (Quantity or None) – Shape (…, n_bands)

get_fnu(spectra: ndarray | Quantity, spectra_err: ndarray | Quantity | None = None, mask_nan: bool = True) → Tuple[Quantity, Quantity | None]

Compute synthetic f_nu in AB system for all bands.

Returns:

• fnu (Quantity) – Shape (…, n_bands) in Jy

• fnu_err (Quantity or None) – Shape (…, n_bands) in Jy

abmag(spectra: ndarray | Quantity, spectra_err: ndarray | Quantity | None = None, mask_nan: bool = True) → Tuple[Quantity, Quantity | None]

Compute AB magnitudes for all bands.

plot(add_props=False, ax=None, show=False)

Plot the filters response curve.

Parameters:

show (bool) – If True

class pst.observables.EquivalentWidth(left_wl_range, central_wl_range, right_wl_range)

Bases: object

Equivalent width of an spectral region.

Description

Given a stellar spectra spectra \(F_\lambda\), the equivalent width is defined as the area of the spectral line, defined over a given spectral region (\(\lambda_{\rm C,\,min}, \lambda_{\rm C,\,max}\)), divided by the average flux of the continuum in a given spectral region. It is computed as:

\[EW = \int_{\lambda_{C,\,\rm min}}^{\lambda_{\rm C,\,max}} \left(1 - \frac{F_\lambda(\lambda)}{F_{\rm cont}(\lambda)}\right) d\lambda\]

where:

\[F_{cont}(\lambda) = \frac{F_{\rm B} \lambda_{\rm R} - F_{\rm R} \lambda_{\rm B}}{\lambda_{\rm R} - \lambda_{\rm B}} + \lambda \frac{F_{\rm R} - F_{\rm B}}{\lambda_{\rm R} - \lambda_{\rm B}}\]

and \(F_{\rm B}\) and \(F_{\rm R}\) are the average flux in the left and right spectral windows, respectively, defined as:

\[F_{\rm B} = \frac{1}{\lambda_{\rm B,\,max} - \lambda_{\rm B,\,min}} \int_{\lambda_{\rm B,\,min}}^{\lambda_{\rm B,\,max}} F_\lambda(\lambda) d\lambda\]

\[F_{\rm R} = \frac{1}{\lambda_{\rm R,\,max} - \lambda_{\rm R,\,min}} \int_{\lambda_{\rm R,\,min}}^{\lambda_{\rm R,\,max}} F_\lambda(\lambda) d\lambda\]

where \(\lambda_{\rm B,\,min}\) and \(\lambda_{\rm B,\,max}\) are the left spectral window boundaries, and \(\lambda_{\rm R,\,min}\) and \(\lambda_{\rm R,\,max}\) are the right spectral window boundaries.

Example

>>> from pst.observables import EquivalentWidth
>>> ew = EquivalentWidth(left_wl_range=(3700, 3900),
...                      central_wl_range=(4000, 4100),
...                      right_wl_range=(4200, 4400))
>>> wavelength = np.linspace(3600, 4600, 1000) * u.angstrom
>>> spectra = np.random.normal(1, 0.1, size=wavelength.size) * u.erg / u.s / u.cm**2 / u.angstrom
>>> ew_value, ew_err = ew.compute_ew(wavelength, spectra)

property left_wl_range: Quantity

Spectral range defining the left pseudocontinuum window \(\lambda_{\rm B,\,min}, \lambda_{\rm B,\,max}\).

property right_wl_range: Quantity

Spectral range defining the right pseudocontinuum window \(\lambda_{\rm R,\,min}, \lambda_{\rm R,\,max}\).

property central_wl_range: Quantity

Spectral range defining the equivalent width window \(\lambda_{\rm C,\,min}, \lambda_{\rm C,\,max}\).

compute_ew(wavelength, spectra, spectra_err=None)

Compute the equivalent width of a given input spectra.

Description

The equivalent width is computed using the definition given in the class description. Positive values of the equivalent width indicate an absorption line, while negative values indicate an emission line. The error on the equivalent width is computed using the error propagation formula, assuming null covariance between the spectral points.

Parameters:

• spectra (np.ndarray or :class:astropy.units.Quantity) – Input spectra. If the array is multidimensional, the first axis must correspond to the spectral direction.

• wavelength (np.ndarray or :class:astropy.units.Quantity) – Wavelength array associated to spectra.

• spectra_err (np.ndarray or :class:astropy.units.Quantity, optional) – If provided, computed the associated error of the equivalent width.

returns:

• ew (np.ndarray) – The equivalent width of the input spectra.

• ew_err (np.ndarray) – The associated error of the equivalent width.

classmethod from_json(path)

Load a EquivalentWidth from a JSON file.

Parameters:

path (str) – Path to the JSON file.

Returns:

ew (EquivvalentWidth)

classmethod from_name(name)

Load a EquivalentWidth from a JSON file.

Parameters:

name (str) – Name of the Lick index.

Returns:

ew (EquivvalentWidth)

Utilities

Module containing generic utility functions

pst.utils._const_interp_cumflux(edges_out, edges_in, f_centers_seg)

Interpolate the cumulative integrated flux assuming constant flux density within each input bin.

The function constructs the cumulative integral C(lambda) = integral f(lambda’) d(lambda’) evaluated at the input bin edges and linearly interpolates C onto the output edges. Differences of the interpolated cumulative at consecutive output edges yield the integrated flux per output bin. Because the interpolation operates on the cumulative, total flux is conserved by construction.

Parameters:

• edges_out (array_like of float, shape (M+1,)) – Target bin edges (must be strictly increasing).

• edges_in (array_like of float, shape (N+1,)) – Input bin edges for the current contiguous segment.

• f_centers_seg (array_like of float, shape (N,)) – Flux density sampled at input bin centers for this segment.

Returns:

cum_flux_out (ndarray of float, shape (M+1,)) – Cumulative integrated flux evaluated at edges_out.

Notes

This assumes f is constant within each input bin (zero order hold). It is fast and robust and is the default choice for most spectra.

pst.utils._f_edges_from_centers(f_centers)

Estimate flux density at bin edges from center samples. Intended for use with a piecewise linear model of f(lambda) inside each bin.

Interior edge values are computed as the average of adjacent centers. End edges are extrapolated linearly from the nearest two centers.

Parameters:

f_centers (array_like of float, shape (N,)) – Flux density at bin centers.

Returns:

f_edges (ndarray of float, shape (N+1,)) – Estimated flux density at bin edges.

pst.utils._lin_interp_cumflux(edges_out, edges_in, f_centers_seg)

Interpolate the cumulative integrated flux assuming linear flux density within each input bin.

The function estimates f(lambda) at the input edges from center values and integrates using the trapezoidal rule to obtain the cumulative C(lambda) at the edges. C is then linearly interpolated onto the output edges.

Parameters:

• edges_out (array_like of float, shape (M+1,)) – Target bin edges.

• edges_in (array_like of float, shape (N+1,)) – Input bin edges for the current segment.

• f_centers_seg (array_like of float, shape (N,)) – Flux density at input bin centers.

Returns:

cum_flux_out (ndarray of float, shape (M+1,)) – Cumulative integrated flux at edges_out.

pst.utils._spline_interp_cumflux(edges_out, edges_in, f_centers_seg)

Interpolate the cumulative integrated flux using a monotonic cubic (PCHIP) spline.

The cumulative C(lambda) is built assuming constant f within each input bin and is then evaluated on edges_out using a monotone cubic interpolant to avoid oscillations and overshoot.

Parameters:

• edges_out (array_like of float, shape (M+1,)) – Target bin edges.

• edges_in (array_like of float, shape (N+1,)) – Input bin edges for the current segment.

• f_centers_seg (array_like of float, shape (N,)) – Flux density at input bin centers.

Returns:

cum_flux_out (ndarray of float, shape (M+1,)) – Cumulative integrated flux at edges_out.

pst.utils._infer_wl_edges_from_centers(wave)

Infer bin edges from wavelength centers.

Midpoints between adjacent centers are used for interior edges, and the outer edges are extrapolated by half of the end spacing.

Parameters:

wave (array_like of float, shape (N,)) – Monotonic wavelength centers.

Returns:

edges (ndarray of float, shape (N+1,)) – Inferred wavelength bin edges.

pst.utils._segments_from_mask(mask)

Compute contiguous valid segments from a boolean mask.

Parameters:

mask (array_like of bool, shape (N,)) – True marks valid bins; False marks invalid or gap bins.

Returns:

segments (list of tuple) – List of (start, end) index pairs describing valid segments.

pst.utils.resample_via_cumulative_masked(w_in, f_in, *, w_out=None, edges_in=None, edges_out=None, mask_valid=None, fill_value=0.0, kind='const')

Flux conserving resampling using cumulative integrals. Robust to NaN or invalid values by treating them as gaps.

Each contiguous valid segment builds its own cumulative integral C(lambda) on input edges. C is evaluated on the output edges using one of the supported interpolation types, and the difference between consecutive output edges gives the integrated flux. Flux outside all valid segments can optionally be filled with a constant value.

Parameters:

• w_in (array_like of float) – Input bin centers (monotonic).

• f_in (array_like of float) – Flux density at input centers.

• w_out (array_like of float, optional) – Output bin centers (monotonic). Provide this or edges_out.

• edges_in (array_like of float, optional) – Input bin edges. If None, inferred from w_in.

• edges_out (array_like of float, optional) – Output bin edges. If None, inferred from w_out.

• mask_valid (array_like of bool, optional) – Mask of valid bins. If None, uses finite values of f_in.

• fill_value (float, default 0.0) – Flux density assumed outside valid coverage.

• kind ({“const”, “linear”, “cubic”}, default “const”) – Interpolation scheme for cumulative inside segments.

Returns:

f_out (ndarray of float) – Flux density per unit wavelength on the output grid.

pst.utils.resample_via_bin_frac(w_in, f_in, w_out=None, edges_in=None, edges_out=None, var_in=None, fill_value=0.0, return_matrix=False)

Flux-conserving resampling via bin-overlap fractions (matrix method).

Builds a sparse mapping matrix between input and output bins whose entries are the fractional overlaps (in wavelength) of each input bin with each output bin. Optionally propagates variances under the assumption of diagonal input covariance.

Parameters:

• w_in ((N,) array_like) – Input bin centers (strictly monotonic).

• f_in ((N,) array_like) – Input flux density per unit wavelength.

• w_out ((M,) array_like, optional) – Output bin centers. Provide this or edges_out.

• edges_in ((N+1,) array_like, optional) – Input edges. If None, inferred from w_in.

• edges_out ((M+1,) array_like, optional) – Output edges. If None, inferred from w_out.

• var_in ((N,) array_like, optional) – Input variances of flux density per input bin (uncorrelated).

• fill_value (float, default 0.0) – Flux density to add in output bins lying outside the input support (e.g., 0.0 for no extrapolation).

• return_matrix (bool, default False) – If True, also return the CSR sparse mapping matrix (N_out, N_in) of fractional overlaps.

Returns:

• f_out ((M,) ndarray) – Output flux density per unit wavelength.

• var_out ((M,) ndarray or None) – Propagated variance (None if var_in is None).

• sparse_matrix (scipy.sparse.csr_matrix, optional) – Mapping matrix such that F_out = M @ (f_in * delta_lambda_in).

Notes

• Exactly flux-conserving (up to FP tol.), independent of within-bin shape assumptions.

• Useful for variance propagation and covariance modeling.

• Complexity O(N_overlaps); efficient with CSR sparsity.

pst.utils.flux_conserving_interpolation(new_wave, wave, spectra, *, method='binfrac', spectra_err=None, **interp_args)

High-level wrapper for flux-conserving spectral resampling.

Provides two flux-conserving resampling paths: (i) cumulative piecewise method robust to gaps/outliers (method=”cumulative”)

via resample_via_cumulative_masked, and

2. exact bin-overlap matrix method (method=”binfrac”) via

resample_via_bin_frac with optional variance propagation.

Parameters:

• new_wave (array_like or Quantity) – Target wavelength centers. If Quantity, units are respected; if not, treated as unitless but must be consistent with wave.

• wave (array_like or Quantity) – Input wavelength centers.

• spectra (array_like or Quantity) – Flux density at wave. If Quantity, units are preserved.

• method ({“cumulative”,”binfrac”}, default “cumulative”) – Resampling engine. See functions referenced above for details and additional behavior (gaps, variance, matrix return).

• spectra_err (array_like or Quantity, optional) – 1-sigma uncertainty on spectra. If provided with method=”binfrac”, variances are propagated. For method=”cumulative”, spectra_err is resampled as a field (no covariance coupling).

• **interp_args – Passed through to the selected engine: - cumulative: w_out, edges_in, edges_out, mask_valid,

  fill_value, kind

  • binfrac: w_out, edges_in, edges_out, fill_value, return_matrix (ignored here)

Returns:

out (Quantity or tuple) – If spectra_err is None: resampled spectrum as Quantity with the same unit as spectra. If spectra_err is given: a tuple (spectrum, variance) where both are Quantity, and variance has squared units of spectra.

Notes

• For rigorous covariance propagation, prefer method=”binfrac”.

pst.utils.gaussian1d_conv(f, sigma, deltax)

Apply a gaussian convolution to a 1D array f(x).

params

fnp.array

1D array containing the data to be convolved with.

sigmanp.array

1D array containing the values of sigma at each value of x

deltaxfloat

Step size of x in “physical” units.

pst.utils.check_unit(quantity, default_unit=None, equivalence=None, **equiv_kwargs)

Check the units of an input quantity.

Parameters:

• quantity (np.ndarray or astropy.units.Quantity) – Input quantity.

• default_unit (astropy.units.Unit, default=None) – If quantity has not units, it corresponds to the unit assigned to it. Otherwise, it is used to check the equivalency with quantity.

pst.utils.check_parameter(x, default_unit=None, **param_kwargs)

Coerce input into a Parameter.

Parameters:

• x (number, astropy.units.Quantity, or Parameter) – Input value.

• default_unit (astropy.units.Unit, optional) – Unit to assume for bare numbers and conversion target.

• **param_kwargs – Forwarded to Parameter constructor when wrapping a non-Parameter.

Returns:

p (Parameter)

pst.utils.broadcast_to_axis(x: ndarray, target_ndim: int, axis: int) → ndarray

Expand 1D array x (wavelength axis) to match target ndim, placing wavelength on axis.

Legacy Compatibility

Warning

pst.models is retained as a compatibility layer and is deprecated. New code should import CEM classes from pst.cem instead.