Spectral Data Cubes from Galaxy Particle distributions

In this example we show how to create a spectral data cube from particle data.

We first load some demo CAMELS data and a grid, as demonstrated elsewhere in the documentation.

import time

import numpy as np
from astropy.cosmology import Planck18 as cosmo
from unyt import kpc

from synthesizer.emission_models import IntrinsicEmission
from synthesizer.grid import Grid
from synthesizer.imaging import SpectralCube
from synthesizer.kernel_functions import Kernel
from synthesizer.load_data.load_camels import load_CAMELS_IllustrisTNG

# Define the grid
grid_name = "test_grid"
grid_dir = "../../../tests/test_grid/"
grid = Grid(grid_name, grid_dir=grid_dir)

# Create galaxy object
gal = load_CAMELS_IllustrisTNG(
    "../../../tests/data/",
    snap_name="camels_snap.hdf5",
    group_name="camels_subhalo.hdf5",
    physical=True,
)[0]

# Calculate the stellar rest frame SEDs for all particles in erg / s / Hz
model = IntrinsicEmission(grid, fesc=0.1, per_particle=True)
sed = gal.stars.get_spectra(model)

# Calculate the observed SED in nJy
sed.get_fnu(cosmo, gal.redshift)

unyt_array([[0.        , 0.        , 0.        , ..., 0.02796633, 0.0279594 ,
        0.02795347],
       [0.        , 0.        , 0.        , ..., 0.02788938, 0.02787349,
        0.02786148],
       [0.        , 0.        , 0.        , ..., 0.02595419, 0.02592956,
        0.0259072 ],
       ...,
       [0.        , 0.        , 0.        , ..., 0.00839355, 0.00839358,
        0.00839135],
       [0.        , 0.        , 0.        , ..., 0.02674974, 0.02674306,
        0.02673475],
       [0.        , 0.        , 0.        , ..., 0.0143164 , 0.01431551,
        0.0143167 ]], 'nJy')

Spectral Data Cube Creation

We now have most of the ingredients we need to generate a spectral data cube from our galaxy. The only parts that are missing are the wavelength array of our spectral data cube, its resolution and FOV, and a kernel for smoothing particles over. We’ll define these below and move on to making the spectral data cube.

# Define the width of the image
width = 30 * kpc

# Define image resolution (here we arbitrarily set it to 100
# pixels along an axis)
resolution = width / 200

# Define the wavelength array
lam = np.linspace(10**3.5, 10**4.5, 1000)

print(
    "Data cube spatial width is %.2f kpc with a %.2f kpc spaxel resolution"
    % (width.value, resolution.value)
)

# Get the SPH kernel
kernel = Kernel()
kernel_data = kernel.get_kernel()

Data cube spatial width is 30.00 kpc with a 0.15 kpc spaxel resolution

Synthesizer allows you to make either a “2D histogram” data cube, where particles are sorted into individual pixels, or data cubes where particles are smoothed over their SPH kernels. We’ll focus on the latter, but to do the former you simply call get_data_cube_hist with an sed, coordinates and the Sed quantity you want to populate the data cube with.

To make a smoothed data cube we first instantiate the SpectralCube, and then call get_data_cube_smoothed, which takes a sed, coordinates, smoothing_lengths, a kernel, the threshold of the kernel, and the Sed quantity you want to populate the data cube with. The possible quantities are "lnu", "luminosity" or "llam" for rest frame luminosities, or "fnu", "flam" or "flux" for fluxes (the latter 3 require get_fnu or get_fnu0 to have been called). We will make a cube populated with "fnu", i.e. the spectral flux density.

cube_start = time.time()

# Get the data cube
cube = SpectralCube(resolution=resolution, lam=lam, fov=width)

# And get the cube itself
cube.get_data_cube_smoothed(
    sed,
    coordinates=gal.stars.coordinates,
    smoothing_lengths=gal.stars.smoothing_lengths,
    kernel=kernel_data,
    kernel_threshold=1,
    quantity="fnu",
)

print("Spectral data cube created, took:", time.time() - cube_start)

Spectral data cube created, took: 0.5700433254241943

And that’s it. We now have a spectral data cube to analyse. We can visualise the data cube by making an animation.

# Animate the data cube
ani = cube.animate_data_cube(fps=240, show=True)

Galaxy helper method

If you don’t want to use the low level SpectralCube object we also include a helper method on a galaxy.

cube = gal.get_data_cube(
    resolution,
    width,
    lam,
    cube_type="hist",
    stellar_spectra="intrinsic",
    kernel=kernel,
    kernel_threshold=1,
    quantity="flux",
)