Working with data from simulations

In cosmological simulations it is very typical to have massive particles that represent a single stellar population with a given age and metallicity. In order to compare predictions from these simulations with observational data it is ofter required to convert from these quantities to flux, magnitudes or surface brightness. In this example, we explore how to do this with milespy with an example simulated galaxy.

from astropy import units as u
from astropy.visualization import quantity_support
quantity_support()
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import numpy as np

We load the data from a simple data file part of pynbody test suite. It contains a simulation of a galaxy performed with the GASOLINE code.

import  requests
import os
import  tarfile

url = "https://zenodo.org/records/12687409/files/gasoline.tar.gz"

if not os.path.exists("gasoline_ahf"):
    with requests.get(url, stream=True) as rx, tarfile.open(fileobj=rx.raw, mode="r:gz") as tarobj:
        tarobj.extractall()

The data is loaded with pynbody and converted to astropy Quantity arrays.

import pynbody
s = pynbody.load('gasoline_ahf/g15784.lr.01024')
h = s.halos()
s.physical_units()
pynbody.analysis.center(h[0])
pynbody.analysis.sideon(h[0])

x = np.array(s.s['x'].in_units('kpc')) << u.kpc
y = np.array(s.s['y'].in_units('kpc')) << u.kpc
z = np.array(s.s['z'].in_units('kpc')) << u.kpc
vx = np.array(s.s['vx'].in_units('km s**-1')) << (u.km/u.s)
vy = np.array(s.s['vy'].in_units('km s**-1')) << (u.km/u.s)
vz = np.array(s.s['vz'].in_units('km s**-1')) << (u.km/u.s)
mass = np.array(s.s['mass'].in_units('Msol')) << u.Msun
met = np.log10(np.array(s.s['metals'])/0.0134) << u.dex
age = np.array(s.s['age'].in_units('Gyr')) << u.Gyr

pynbody.halo : Unable to load AHF substructure file; continuing without. To expose the underlying problem as an exception, pass ignore_missing_substructure=False to the AHFCatalogue constructor
/tmp/ipykernel_891727/1505409435.py:15: RuntimeWarning: divide by zero encountered in log10
  met = np.log10(np.array(s.s['metals'])/0.0134) << u.dex

Then, we mask the data that is not of interest and we clip the age and metallicities, as the library does not contain SSP with ages < 0.6 Gyr and thus will throw errors.

m = np.sqrt(x**2+y**2+z**2) < 20*u.kpc
x = x[m]
y = y[m]
z = z[m]
vz = vz[m]
mass = mass[m]
met = met[m]
age = age[m]
np.clip(met, a_min=-2.32*u.dex, a_max=0.22*u.dex, out=met)
np.clip(age, a_min=0.0631*u.Gyr, a_max=14*u.Gyr, out=age)

plt.plot(x,y, 'k,', alpha=0.1)
plt.gca().set_aspect('equal')

For this exercise we will use the EMILES SSP library

from milespy import SSPLibrary
lib = SSPLibrary(source='EMILES_SSP')
plt.plot(lib.models.spectral_axis, lib.models[10].flux);

However, we do not require all the spectral range for this analysis. Thus, we trim the models to the region of interest. Note that not doing this could cause a very large memory allocation, as there are ~50000 spectral points and ~150000 particles, resulting in ~60 GB of required RAM! If you really need the full spectral range, you could perform the following analysis in batches of particles to make it more manageable for your system.

lib.trim(3000.*u.AA, 9000*u.AA)
plt.plot(lib.models.spectral_axis, lib.models[10].flux);

If this is still not enough or you do not need the spectral resolution, you can also resample the whole library.

lib.resample(np.arange(3010, 9000, 3)<<u.AA)
plt.plot(lib.models.spectral_axis, lib.models[10].flux);

We decide to use closest as it is way faster than interpolate and having 150000 particles we can assume that the mean deviations will be small… but that is something you must consider carefully!

%time specs = lib.closest(age=age, met=met, imf_slope=np.full(met.shape, 1.3), mass=mass)

100%|██████████████████████████████████████████████████████████████████████████████████████████████████| 227645/227645 [00:08<00:00, 27788.12it/s]

CPU times: user 8.76 s, sys: 777 ms, total: 9.54 s
Wall time: 9.41 s

As a side note, you can save and restore the generated spectra as shown in the cell below if you do not want to repeat the costly interpolation all the time.

#np.savez('all', flux=specs.flux.value, wave=specs.spectral_axis.value)
#loaded = np.loadz('all')
#flux = u.Quantity(loaded['flux'], unit=u.L_sun / u.AA)
#wave = u.Quantity(loaded['wave'], unit=u.AA)
#specs = spectra(flux=flux, spectral_axis=wave)

Taking into account line-of-sight velocities

In simulations, the velocity information of each SSP is also accesible, thus we can take this into account in the generated spectra. We can compute, for example, the spectra of the galaxy with and without the velocity shifts.

In the current projection, the velocity in the Z-axis would be the line-of-sight velocity for the shifts,

v_los = vz

and we use it to generate new spectra with this velocity shift.

%time shifted = specs.velocity_shift(v_los)

100%|███████████████████████████████████████████████████████████████████████████████████████████████████| 227645/227645 [00:26<00:00, 8568.38it/s]

CPU times: user 26.4 s, sys: 440 ms, total: 26.9 s
Wall time: 26.8 s

We can then compute the spectra of the galaxy as a whole, and compare it to the spectra obtained without velocity information.

galaxy_spectrum = specs.sum(axis=0)
galaxy_spectrum_los = shifted.sum(axis=0)

plt.plot(galaxy_spectrum.spectral_axis, galaxy_spectrum.flux)
plt.plot(galaxy_spectrum_los.spectral_axis, galaxy_spectrum_los.flux, alpha=0.9)

x1, x2, y1, y2 = 6520, 6600, 2.7e6, 4.5e6  # subregion of the original image
axins = plt.gca().inset_axes(
    [0.4, 0.1, 0.4, 0.4],
    xlim=(x1, x2), ylim=(y1, y2), xticklabels=[], yticklabels=[])
axins.plot(galaxy_spectrum.spectral_axis, galaxy_spectrum.flux)
axins.plot(galaxy_spectrum_los.spectral_axis, galaxy_spectrum_los.flux, alpha=0.9)

plt.gca().indicate_inset_zoom(axins, edgecolor="black");

In this simple example, the velocity distribution is nearly Guassian, and thus we obtain a simple smoothing of the spectra. However, this procedure can be the basis for more complex analysis.

Building a spectral cube

As we have a spectrum per particle, it could be of interest to build a spectral cube, where we sum the contribution of all particles in a given bin (weighted by their mass). For this, we first define the binning configuration and we used binned_statistic_2d to get the bin to which each particle belongs.

from scipy.stats import binned_statistic_2d

nbins = 100
bins = np.linspace(-20,20,nbins)
_, xed, yed, nn = binned_statistic_2d(x, y, mass, statistic='sum', bins=bins, expand_binnumbers=True)

Then we create an empty spectral cube and we start summing the contribution to each bin

from milespy.spectra import Spectra

cube = Spectra(spectral_axis=specs.spectral_axis, flux=u.Quantity(np.zeros((nbins+1,nbins+1,specs.npix)), unit=specs.flux.unit))
for i,(xx,yy) in enumerate(nn.T):
    cube.data[yy,xx,:] += specs.data[i, :]

We compute the zero moment (i.e., intensity) of the whole cube

from specutils.analysis import moment

mom = moment(cube, order=0)
plt.imshow(mom.value, norm=LogNorm());

Spectral regions

Using spectral_slab we can take a look at specific part of the spectra for its analysis

from specutils.manipulation import spectral_slab

subspec = spectral_slab(cube, 6530.*u.AA, 6590*u.AA)
ha_flux = np.sum(subspec.flux, axis=(0,1))
plt.plot(subspec.spectral_axis, ha_flux);

ha_mom = moment(subspec, order=0)
plt.imshow(ha_mom.value, norm=LogNorm());

Magnitudes

We are going to compute the i, g and r magnitudes of the galaxy, assuming that it is at a distance of 10 Mpc

import milespy.filter as flib

fnames = flib.search("OAJ_JPAS.(iSDSS|gSDSS|rSDSS)")
filts = flib.get(fnames)
abs_mags = cube.magnitudes(filters=filts)

d = 10 * u.Mpc
distance_modulus = 5 * np.log10( d.to(u.pc)/(10*u.pc) ) - 5
mags = {k: (abs_mags[k] + distance_modulus) for k in abs_mags.keys()}

plt.imshow(mags['OAJ_JPAS.rSDSS'], cmap='Greys', vmin=9, vmax=16, extent=(-0.2, 0.2, -0.2, 0.2))
plt.gca().set_facecolor("k")
plt.xlabel('x [kpc]')
plt.ylabel('y [kpc]')
cbar = plt.colorbar().set_label(r'm$_\mathrm{r}$', fontsize=20)

A color map can also be obtained

plt.imshow(mags['OAJ_JPAS.gSDSS']-mags['OAJ_JPAS.rSDSS'], cmap='bone')
plt.colorbar();

False-color image

Lastly, with some tinkering, the magnitude maps can be combined to create a false-color image. However, there is no fixed way to approach this. In our case we renormalize each filter independently.

def norm(f):
    mask = np.isfinite(f)
    minf = f[mask].min()
    maxf = f[mask].max()
    f[mask] = (f[mask] - minf)/(maxf-minf)
    return 1-1.4*f

r = norm(mags['OAJ_JPAS.iSDSS'])
g = norm(mags['OAJ_JPAS.rSDSS'])
b = norm(mags['OAJ_JPAS.gSDSS'])
rgb = np.dstack((r,g,b))
plt.imshow(rgb);

WARNING:matplotlib.image:Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

But there are other stablished tools in astropy, for example. However, these also required some parameter tuning to get a nice-looking image…

from astropy.visualization import make_lupton_rgb
rgb_default = make_lupton_rgb(10**(-mags['OAJ_JPAS.iSDSS']), 10**(-mags['OAJ_JPAS.rSDSS']), 10**(-mags['OAJ_JPAS.gSDSS']), stretch=0.8, Q=2.5)
plt.imshow(rgb_default);