Quickstart Guide to DSPS

This demo notebook begins by downloading the DSPS default option for the SSP spectral library. These data are stored at this URL in a flat hdf5 file with column names as expected by the dsps.load_ssp_templates function, which we will demonstrate below.

When downloading and storing SSP libraries, you can optionally use the DSPS_DRN environment variable to specify the default location where DSPS will look for SSP libraries. But here we’ll just save the downloaded data to tempdata.h5, directly pass the filename to the data loader. The load_ssp_templates that we’ll use to load these SSPs is just a convenience function - all of the DSPS functions that we’ll demonstrate in this notebook accept plain arrays and floats as inputs, and so you can store your SSP data on disk in whatever format you like.

! curl https://portal.nersc.gov/project/hacc/aphearin/DSPS_data/ssp_data_fsps_v3.2_lgmet_age.h5 > tempdata.h5

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100 58.7M  100 58.7M    0     0  30.4M      0  0:00:01  0:00:01 --:--:-- 30.4M

Inspect the SSP data

from dsps import load_ssp_templates
ssp_data = load_ssp_templates(fn='tempdata.h5')

print(ssp_data._fields)

print('\nssp_lgmet.shape = {}'.format(ssp_data.ssp_lgmet.shape))
print('ssp_lg_age_gyr.shape = {}'.format(ssp_data.ssp_lg_age_gyr.shape))
print('ssp_wave.shape = {}'.format(ssp_data.ssp_wave.shape))
print('ssp_flux.shape = {}'.format(ssp_data.ssp_flux.shape))

('ssp_lgmet', 'ssp_lg_age_gyr', 'ssp_wave', 'ssp_flux', 'ssp_emline_wave', 'ssp_emline_luminosity')

ssp_lgmet.shape = (12,)
ssp_lg_age_gyr.shape = (107,)
ssp_wave.shape = (5994,)
ssp_flux.shape = (12, 107, 5994)

The returned ssp_data is a namedtuple storing 4 ndarrays for the age-metallicity grid of the SSP spectra. Galaxy SEDs are calculated via probability-weighted sums of these spectral templates. For a galaxy observed at some \(t_{\rm obs}\), we’ll calculate the restframe SED of two different models in the cells below:

1. a galaxy with a tabulated star formation history (SFH), and metallicity Z distributed as a lognormal about some median Z, using the calc_rest_sed_sfh_table_lognormal_mdf function.

2. a galaxy with SFH table and also tabulated history of metallicity (ZH), using the calc_rest_sed_sfh_table_met_table function.

In the cells below, we’ll randomly generate an SFH and ZH for a galaxy, and then plot the results.

import numpy as np

gal_t_table = np.linspace(0.05, 13.8, 100) # age of the universe in Gyr
gal_sfr_table = np.random.uniform(0, 10, gal_t_table.size) # SFR in Msun/yr

gal_lgmet = -2.0 # log10(Z)
gal_lgmet_scatter = 0.2 # lognormal scatter in the metallicity distribution function

The SED calculating functions require you specify the time of the observation, t_obs, rather than the redshift, z_obs. We’ll use the age_at_z function in dsps.cosmology to calculate the relationship between these two quantities, assuming the default redshift of DSPS. You could also use this same function to compute gal_t_table in case your input SFH is tabulated as a function of redshift.

from dsps.cosmology import age_at_z, DEFAULT_COSMOLOGY

print(DEFAULT_COSMOLOGY)

z_obs = 0.5
t_obs = age_at_z(z_obs, *DEFAULT_COSMOLOGY) # age of the universe in Gyr at z_obs
t_obs = t_obs[0] # age_at_z function returns an array, but SED functions accept a float for this argument

CosmoParams(Om0=0.3075, w0=-1.0, wa=0.0, h=0.6774)

from dsps import calc_rest_sed_sfh_table_lognormal_mdf
from dsps import calc_rest_sed_sfh_table_met_table

sed_info = calc_rest_sed_sfh_table_lognormal_mdf(
    gal_t_table, gal_sfr_table, gal_lgmet, gal_lgmet_scatter,
    ssp_data.ssp_lgmet, ssp_data.ssp_lg_age_gyr, ssp_data.ssp_flux, t_obs)


gal_lgmet_table = np.linspace(-3, -2, gal_t_table.size)

sed_info2 = calc_rest_sed_sfh_table_met_table(
    gal_t_table, gal_sfr_table, gal_lgmet_table, gal_lgmet_scatter,
    ssp_data.ssp_lgmet, ssp_data.ssp_lg_age_gyr, ssp_data.ssp_flux, t_obs)

from matplotlib import pyplot as plt

fig, ax = plt.subplots(1, 1)
__=ax.loglog()
__=ax.plot(ssp_data.ssp_wave, sed_info.rest_sed)
__=ax.plot(ssp_data.ssp_wave, sed_info2.rest_sed)

Calculating photometry

Now we’ll use dsps.photometry to calculate the apparent and absolute magnitudes of an SED observed through a broadband filter. One additional ingredient we need for this calculation is the transmission curve of some filter. For this, we’ll download one from the same public URL where we previously downloaded the SSP spectra, and just write the result to a temporary file. But the load_transmission_curve function also supports the use of the DSPS_DRN environment variable in case you want DSPS to remember your default data location. The only difference is that you should store your transmission curves in the filters subdirectory of DSPS_DRN. And as above, you can also ignore this data-loading convenience function and store your transmission curves wherever you like and in whatever format you prefer.

! curl https://portal.nersc.gov/project/hacc/aphearin/DSPS_data/filters/lsst_g_transmission.h5 > tempfilter.h5

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0

/home/docs/checkouts/readthedocs.org/user_builds/dsps/conda/latest/lib/python3.11/pty.py:89: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.
  pid, fd = os.forkpty()

100  134k  100  134k    0     0   307k      0 --:--:-- --:--:-- --:--:--  307k

from dsps import load_transmission_curve
lsst_g = load_transmission_curve(fn="tempfilter.h5")

print(lsst_g._fields)

print('\nlsst_g.wave.shape = {}'.format(lsst_g.wave.shape))
print('lsst_g.transmission.shape = {}'.format(lsst_g.transmission.shape))

('wave', 'transmission')

lsst_g.wave.shape = (8500,)
lsst_g.transmission.shape = (8500,)

Calculating absolute magnitude

Since we have already calculated above the restframe SED at the time the galaxy is observed, then we can directly integrate the SED against the filter transmission curve to compute the absolute magnitude of the galaxy.

from dsps import calc_rest_mag
rest_mag = calc_rest_mag(ssp_data.ssp_wave, sed_info.rest_sed, lsst_g.wave, lsst_g.transmission)

Calculating apparent magnitude

To calculate the apparent magnitude, we need to redshift the SED and also apply the appropriate cosmological dimming factor to the restframe flux. These calculations are done under-the-hood with the flat_wcdm.py module in dsps.cosmology, and so we just need to pass in the same cosmological parameters used above to calculate t_obs from z_obs.

from dsps import calc_obs_mag

obs_mag = calc_obs_mag(ssp_data.ssp_wave, sed_info.rest_sed, lsst_g.wave, lsst_g.transmission,
                      z_obs, *DEFAULT_COSMOLOGY)

! rm tempdata.h5

! rm tempfilter.h5