Generating photometry from a `Galaxy`¶

If you are working with a Galaxy (or many Galaxy objects) containing multiple spectra, you can use the galaxy level get_photo_lnu and get_photo_fluxes methods to generate photometry for all spectra in the galaxy.

Before we can demonstrate this we firs need to make a parametric galaxy with a Stars and Blackhole component, define an emission model for a whole galaxy and generate their rest–frame and observed spectra.

[1]:

import numpy as np
from astropy.cosmology import Planck18 as cosmo
from unyt import K, Mpc, Msun, Myr, cm, yr

from synthesizer.emission_models import (
    BimodalPacmanEmission,
    Blackbody,
    DustEmission,
    EmissionModel,
    Greybody,
    UnifiedAGN,
)
from synthesizer.emission_models.attenuation import PowerLaw
from synthesizer.filters import FilterCollection
from synthesizer.grid import Grid
from synthesizer.parametric import SFH, ZDist
from synthesizer.parametric import Stars as ParametricStars
from synthesizer.particle import BlackHoles, Galaxy
from synthesizer.particle.stars import sample_sfzh

# Get the grids which we'll need for extraction
grid_name = "test_grid"
grid = Grid(grid_name)
nlr_grid = Grid("test_grid_agn-nlr")
blr_grid = Grid("test_grid_agn-blr")

# Define the metallicity history
zh = ZDist.DeltaConstant(metallicity=0.01)

# Define the star formation history
sfh_p = {"max_age": 100 * Myr}
sfh = SFH.Constant(**sfh_p)

# Initialise the parametric Stars object
param_stars = ParametricStars(
    grid.log10age,
    grid.metallicity,
    sf_hist=sfh,
    metal_dist=zh,
    initial_mass=10**9 * Msun,
)

# Define the number of stellar particles we want
n = 500

# Sample the parametric SFZH, producing a particle Stars object
# we will also pass some keyword arguments for some example attributes
part_stars = sample_sfzh(
    sfzh=param_stars.sfzh,
    log10ages=param_stars.log10ages,
    log10metallicities=param_stars.log10metallicities,
    nstar=n,
    current_masses=np.full(n, 10**8.7 / n) * Msun,
    redshift=1,
    coordinates=np.random.normal(0, 0.01, (n, 3)) * Mpc,
    centre=np.zeros(3) * Mpc,
    tau_v=np.random.rand(n),
)

# Make fake properties
n = 4
masses = 10 ** np.random.uniform(low=7, high=9, size=n) * Msun
coordinates = np.random.normal(0, 0.01, (n, 3)) * Mpc
accretion_rates = 10 ** np.random.uniform(low=-2, high=1, size=n) * Msun / yr
metallicities = np.full(n, 0.01)
tau_v = 0.5 * np.ones(n)

# And get the black holes object
blackholes = BlackHoles(
    masses=masses,
    coordinates=coordinates,
    accretion_rates=accretion_rates,
    metallicities=metallicities,
    ionisation_parameter_nlr=0.01,
    hydrogen_density_nlr=1e4 * cm**-3,
    ionisation_parameter_blr=0.1,
    hydrogen_density_blr=1e10 * cm**-3,
    tau_v=tau_v,
)

# And create the galaxy
gal = Galaxy(
    stars=part_stars,
    black_holes=blackholes,
    redshift=1,
)

# Get the stellar pacman model
pc_model = BimodalPacmanEmission(
    grid=grid,
    tau_v_ism=1.0,
    tau_v_birth=0.7,
    dust_curve_ism=PowerLaw(slope=-1.3),
    dust_curve_birth=PowerLaw(slope=-0.7),
    dust_emission_ism=Blackbody(temperature=100 * K),
    dust_emission_birth=Blackbody(temperature=30 * K),
    fesc=0.2,
    fesc_ly_alpha=0.9,
    label="stellar_total",
)
pc_model.save_spectra("stellar_total", "intrinsic")

# Initialise the UnifiedAGN model
uniagn = UnifiedAGN(
    nlr_grid,
    blr_grid,
    ionisation_parameter_nlr=0.01,
    hydrogen_density_nlr=1e4 * cm**-3,
    ionisation_parameter_blr=0.1,
    hydrogen_density_blr=1e10 * cm**-3,
    torus_emission_model=Greybody(1000 * K, 1.5),
    diffuse_dust_curve=PowerLaw(slope=-1.0),
    per_particle=True,
)
# relabel intrinsic and attenuated to avoid clashes
uniagn.relabel("intrinsic", "agn_intrinsic")
uniagn.relabel("attenuated", "agn_attenuated")


gal_intrinsic = EmissionModel(
    label="total_intrinsic",
    combine=(uniagn["agn_intrinsic"], pc_model["intrinsic"]),
    emitter="galaxy",
)

gal_attenuated = EmissionModel(
    label="total_attenuated",
    combine=(uniagn["agn_attenuated"], pc_model["attenuated"]),
    related_models=(gal_intrinsic,),
    emitter="galaxy",
)

# And now include the dust emission
gal_dust = DustEmission(
    dust_emission_model=Greybody(30 * K, 1.2),
    dust_lum_intrinsic=gal_intrinsic,
    dust_lum_attenuated=gal_attenuated,
    emitter="galaxy",
    label="dust_emission",
)

gal_total = EmissionModel(
    label="total",
    combine=(gal_attenuated, gal_dust),
    related_models=(gal_intrinsic,),
    emitter="galaxy",
)
# Get the spectra
sed = gal.get_spectra(gal_total)

# Get fluxes
gal.get_observed_spectra(cosmo)

/tmp/ipykernel_7934/2561947077.py:14: FutureWarning: The filters module has been moved to the instruments module. Please update your imports synthesizer.filters -> synthesizer.instruments
  from synthesizer.filters import FilterCollection
/opt/hostedtoolcache/Python/3.10.19/x64/lib/python3.10/site-packages/unyt/array.py:1832: RuntimeWarning: divide by zero encountered in log10
  out_arr = func(np.asarray(inp), out=out_func, **kwargs)
/opt/hostedtoolcache/Python/3.10.19/x64/lib/python3.10/site-packages/unyt/array.py:1832: RuntimeWarning: overflow encountered in exp
  out_arr = func(np.asarray(inp), out=out_func, **kwargs)
/opt/hostedtoolcache/Python/3.10.19/x64/lib/python3.10/site-packages/unyt/array.py:1972: RuntimeWarning: overflow encountered in multiply
  out_arr = func(

We can then combine with a FilterCollection to generate the galaxy-level photometry.

[2]:

# Get the filter collection
fs = [
    f"JWST/NIRCam.{f}"
    for f in ["F090W", "F150W", "F200W", "F277W", "F356W", "F444W"]
]

filters = FilterCollection(
    filter_codes=fs,
    new_lam=grid.lam,
)

# Get the photometry
gal.get_photo_lnu(filters)
gal.get_photo_fnu(filters)

The photometry produced by these methods are stored in the photo_lnu and photo_fluxes attributes, either at the base galaxy level or the individual components. These attributes are dictionaries containing the photometry for each spectra key.

For example, on the Galaxy.Stars object:

[3]:

print(gal.stars.photo_lnu)
print(gal.stars.photo_fnu)

{'intrinsic': <synthesizer.photometry.PhotometryCollection object at 0x7fa7acf26080>}
{'intrinsic': <synthesizer.photometry.PhotometryCollection object at 0x7fa7acf25ed0>}

Or on the galaxy level

[4]:

print(gal.photo_lnu)

{'total': <synthesizer.photometry.PhotometryCollection object at 0x7fa743257460>, 'total_attenuated': <synthesizer.photometry.PhotometryCollection object at 0x7fa743257610>, 'dust_emission': <synthesizer.photometry.PhotometryCollection object at 0x7fa7432577c0>, 'total_intrinsic': <synthesizer.photometry.PhotometryCollection object at 0x7fa743257970>}

As before we can print the photometry.

[5]:

print(gal.photo_fnu["total_attenuated"])

-----------------------------------------------------
|                 PHOTOMETRY (FLUX)                 |
|------------------------------------|--------------|
| JWST/NIRCam.F090W (λ = 9.02e+03 Å) | 8.67e+04 nJy |
|------------------------------------|--------------|
| JWST/NIRCam.F150W (λ = 1.50e+04 Å) | 7.65e+04 nJy |
|------------------------------------|--------------|
| JWST/NIRCam.F200W (λ = 1.99e+04 Å) | 8.22e+04 nJy |
|------------------------------------|--------------|
| JWST/NIRCam.F277W (λ = 2.76e+04 Å) | 1.64e+05 nJy |
|------------------------------------|--------------|
| JWST/NIRCam.F356W (λ = 3.57e+04 Å) | 4.04e+05 nJy |
|------------------------------------|--------------|
| JWST/NIRCam.F444W (λ = 4.40e+04 Å) | 6.77e+05 nJy |
-----------------------------------------------------

Or plot them.

[6]:

gal.photo_fnu["total_attenuated"].plot_photometry(show=True)

[6]:

(<Figure size 350x500 with 2 Axes>,
 <Axes: xlabel='$\\lambda_\\mathrm{obs}/[\\mathrm{\\AA}]$', ylabel='$F/[\\mathrm{\\rm{nJy}}]$'>)

../../_images/observables_photometry_galaxy_phot_11_1.png

[7]:

gal.photo_fnu["total_intrinsic"].plot_photometry(show=True)

[7]:

(<Figure size 350x500 with 2 Axes>,
 <Axes: xlabel='$\\lambda_\\mathrm{obs}/[\\mathrm{\\AA}]$', ylabel='$F/[\\mathrm{\\rm{nJy}}]$'>)

../../_images/observables_photometry_galaxy_phot_12_1.png

[8]:

gal.photo_fnu["total"].plot_photometry(show=True)

[8]:

(<Figure size 350x500 with 2 Axes>,
 <Axes: xlabel='$\\lambda_\\mathrm{obs}/[\\mathrm{\\AA}]$', ylabel='$F/[\\mathrm{\\rm{nJy}}]$'>)

../../_images/observables_photometry_galaxy_phot_13_1.png

Calculating light radii¶

Once we have photometry we can calculate the radius enclosing a given fraction of the light for a component. Here we’ll calculate the half light radius for both the intrinsic emission and the total emission in “F444W” in terms of luminosity, but before we can do that we need to get the particle spectra and call get_particle_photo_lnu to first calculate the per particle photometry (above we used the galaxy level methods to calculate integrated spectra).

[9]:

# Get the particle spectra
pc_model.set_per_particle(True)
gal.stars.get_spectra(pc_model)

# Get the particle photometry
gal.stars.get_particle_photo_lnu(filters)

int_r50 = gal.stars.get_half_luminosity_radius(
    "intrinsic", "JWST/NIRCam.F444W"
)
tot_r50 = gal.stars.get_half_luminosity_radius(
    "stellar_total", "JWST/NIRCam.F444W"
)
print(int_r50, tot_r50)

0.014129268905544463 Mpc 0.014129481002467097 Mpc

Similarly to the “attr” radii we can compute for any particle component, we can also compute the radius enclosing any fraction of the light for any particle component.

[10]:

int_r20 = gal.stars.get_luminosity_radius(
    "intrinsic", "JWST/NIRCam.F444W", frac=0.2
)
tot_r20 = gal.stars.get_luminosity_radius(
    "stellar_total", "JWST/NIRCam.F444W", frac=0.2
)
int_r80 = gal.stars.get_luminosity_radius(
    "intrinsic", "JWST/NIRCam.F444W", frac=0.8
)
tot_r80 = gal.stars.get_luminosity_radius(
    "stellar_total", "JWST/NIRCam.F444W", frac=0.8
)
print(int_r20, int_r80, tot_r20, tot_r80)

0.009441265688010491 Mpc 0.020139019050080406 Mpc 0.00961959083042451 Mpc 0.020147561890801007 Mpc

Calculating weighted attributes¶

Similarly to the weighted attributes shown in the stellar component docs, we can also calculate luminosity and flux weighted attributes for the galaxy’s components. Again these exist in generic and specific attribute forms. First we demonstrate the generic form, this requires the attribute to weight, the spectra type, and the filter code to weight by as arguments.

[11]:

print(
    gal.stars.get_lum_weighted_attr(
        "log10ages", "stellar_total", "JWST/NIRCam.F444W"
    )
)

6.866358186061508

We can also do this for fluxes.

[12]:

# We haven't computed fluxes yet so lets do that first
gal.stars.get_spectra(pc_model)
gal.get_observed_spectra(cosmo)
gal.stars.get_particle_photo_fnu(filters)

print(
    gal.stars.get_flux_weighted_attr(
        "log10ages", "stellar_total", "JWST/NIRCam.F444W"
    )
)

7.18137629894157

We also provide some specific methods for commonly weighted attributes. These methods still require the spectra type and filter code as arguments, but the attribute to weight is fixed.

[13]:

print(gal.stars.get_lum_weighted_age("stellar_total", "JWST/NIRCam.F444W"))
print(gal.stars.get_flux_weighted_age("stellar_total", "JWST/NIRCam.F444W"))
print(
    gal.stars.get_lum_weighted_metallicity(
        "stellar_total", "JWST/NIRCam.F444W"
    )
)
print(
    gal.stars.get_flux_weighted_metallicity(
        "stellar_total", "JWST/NIRCam.F444W"
    )
)
print(
    gal.stars.get_lum_weighted_optical_depth(
        "stellar_total", "JWST/NIRCam.F444W"
    )
)
print(
    gal.stars.get_flux_weighted_optical_depth(
        "stellar_total", "JWST/NIRCam.F444W"
    )
)

18316785.39617686 yr
27314145.97131986 yr
0.01
0.010000000000000002
0.4719336724806078
0.45958470074320845

Generating photometry from a Galaxy¶

Calculating light radii¶

Calculating weighted attributes¶

Generating photometry from a `Galaxy`¶