Generate active galaxy

Finall, in this example we’re going to demonstrate how to make a composite galaxy, including with imaging. For more information on defining parametric morphology see the Imaging examples.

from unyt import K, Msun, Myr, deg, kelvin, kpc, yr

from synthesizer.emission_models import (
    AttenuatedEmission,
    DustEmission,
    GalaxyEmissionModel,
    PacmanEmission,
    UnifiedAGN,
)
from synthesizer.emission_models.attenuation import PowerLaw
from synthesizer.emission_models.dust.emission import Greybody
from synthesizer.filters import UVJ
from synthesizer.grid import Grid
from synthesizer.parametric import SFH, BlackHole, Stars, ZDist
from synthesizer.parametric.galaxy import Galaxy
from synthesizer.parametric.morphology import Sersic2D

Let’s begin by defining the geometry of the images:

# Define geometry of the images
resolution = 0.01 * kpc  # resolution in kpc
npix = 301
fov = resolution * npix

Define the grid_dir:

grid_dir = "../../../tests/test_grid"

Diffuse ISM dust

Later on we’re going to apply “diffuse” (ISM) dust attenuation to both the stars and blackholes so let’s define these now. Note, these could be different for each component.

tau_v = 0.5
# dust curve slope
alpha = -1.0
dust_curve = PowerLaw(slope=alpha)
dust_emission_model = Greybody(30 * K, 1.2)

The Stars

Let’s first build the stellar component of our galaxy, including setting the morphology so we can make an image later.

First define the grid. We need this to create our star formation and metal enrichment history.

grid_name = "test_grid"
grid = Grid(grid_name, grid_dir=grid_dir)

Let’s define our star formation and metal enrichment history:

stellar_mass = 10**10 * Msun
sfh = SFH.Constant(max_age=100 * Myr)
metal_dist = ZDist.Normal(mean=0.01, sigma=0.05)

Next, let’s define the morphology:

morph = Sersic2D(
    r_eff=1.0 * kpc, sersic_index=1.0, ellipticity=0.5, theta=35.0
)

Now let’s initialise our stars object:

# Get the 2D star formation and metal enrichment history for the
# given SPS grid. This is (age, Z).
stars = Stars(
    grid.log10age,
    grid.metallicity,
    sf_hist=sfh,
    metal_dist=metal_dist,
    initial_mass=stellar_mass,
    morphology=morph,
)

The blackholes

Let’s define the properties of our blackholes. We don’t need to define the morphology since blackholes automatically assume a PointSource geometry.

black_hole = BlackHole(
    mass=10**6.5 * Msun,
    inclination=60 * deg,
    accretion_rate=0.5 * Msun / yr,
    metallicity=0.01,
)

The Emission Model

Define the emission model including contributions from the stars and blackholes.

# Get the NLR and BLR grids
nlr_grid = Grid("test_grid_agn-nlr", grid_dir="../../../tests/test_grid")
blr_grid = Grid("test_grid_agn-blr", grid_dir="../../../tests/test_grid")

# Get the model
pacman = PacmanEmission(
    label="stellar_total",
    grid=grid,
    tau_v=tau_v,
    dust_curve=dust_curve,
    dust_emission=dust_emission_model,
)

# Get the model
uni_agn = UnifiedAGN(
    nlr_grid,
    blr_grid,
    covering_fraction_nlr=0.1,
    covering_fraction_blr=0.1,
    torus_emission_model=Greybody(1000 * kelvin, 1.5),
    label="agn_intrinsic",
)

# Include attenuation and dust emission
att_model = AttenuatedEmission(
    dust_curve=dust_curve,
    apply_dust_to=uni_agn,
    tau_v=tau_v,
    emitter="blackhole",
    label="agn_attenuated",
)
dust_model = DustEmission(
    label="agn_dust_emission",
    dust_emission_model=dust_emission_model,
    dust_lum_intrinsic=uni_agn,
    dust_lum_attenuated=att_model,
    emitter="blackhole",
)

total = GalaxyEmissionModel(
    label="total",
    combine=(dust_model, pacman),
)

The Galaxy

Now that we have our components we can pass them to create an instance of a Galaxy.

# Initialise Galaxy object
galaxy = Galaxy(stars=stars, black_holes=black_hole)
print(galaxy)

+----------------------------------------------------------------------------------------+
|                                         GALAXY                                         |
+---------------+------------------------------------------------------------------------+
| Attribute     | Value                                                                  |
+---------------+------------------------------------------------------------------------+
| galaxy_type   | 'Parametric'                                                           |
+---------------+------------------------------------------------------------------------+
| stars         | <synthesizer.parametric.stars.Stars object at 0x7fe898baf5e0>          |
+---------------+------------------------------------------------------------------------+
| black_holes   | <synthesizer.parametric.blackholes.BlackHole object at 0x7fe84abe8910> |
+---------------+------------------------------------------------------------------------+
| name          | 'parametric galaxy'                                                    |
+---------------+------------------------------------------------------------------------+
| sfzh (51, 13) | 0.00e+00 -> 5.01e+08 (Mean: 1.51e+07)                                  |
+---------------+------------------------------------------------------------------------+

And now we can get the spectra by passing the EmissionModel to the get_spectra method.

galaxy.get_spectra(total)

<synthesizer.sed.Sed at 0x7fe84abeb0a0>

Now we have a galaxy with spectra isolated to the components. To get the combined “total”, “emergent”, and “intrinsic” spectra from all components we can call the galaxy.get_spectra_combined.

galaxy.get_spectra_combined()

Spectra

We can use the plot_spectra function to make a quick spectrum of the galaxy.

fig, ax = galaxy.plot_spectra(
    combined_spectra=True,
    stellar_spectra=True,
    black_hole_spectra=True,
    quantity_to_plot="luminosity",
    figsize=(10, 6),
)

Images

To make images we need filters to define the photometry. For details see the filter and imaging docs.

filters = UVJ(new_lam=grid.lam)

With the spectra and filters made we can make images of each component and combine them.

# Get photometry
galaxy.get_photo_lnu(filters)

# Make images
img = galaxy.get_images_luminosity(
    resolution=resolution,
    fov=fov,
    emission_model=total,
    img_type="smoothed",
)

# Make and plot an rgb image
img.make_rgb_image(rgb_filters={"R": "J", "G": "V", "B": "U"})
fig, ax, _ = img.plot_rgb_image(show=True)