Dust extinction

# NBVAL_SKIP
import os
os.environ['SPS_HOME'] = '/home/annalena/sps_fsps'

Dust extinction models in Rubix

This notebook shows the basics of the dust extinction models implemented in Rubix. We have closely followed the implementation by the dust extinction package. Currently we only support a subset of all available models, namely the Cardelli, Clayton, & Mathis (1989) Milky Way R(V) dependent model, the Gordon et al. (2023) Milky Way R(V) dependent model and the Fitzpatrick & Massa (1990) 6 parameter ultraviolet shape model.

We will demonstrate how to use these models to calculate and visualize the effects of dust extinction on stellar spectra. Additionally, we will show how to integrate these models into a Rubix pipeline to simulate the impact of dust on galaxy observations.

First, we import the dust models from Rubix.

# NBVAL_SKIP
from rubix.spectra.dust.extinction_models import Cardelli89, Gordon23

# NBVAL_SKIP
import matplotlib.pyplot as plt
import numpy as np
import jax.numpy as jnp

We visulaize some of the aspects of the models, i.e. their A(x)/Av as a function of wavelength.

# NBVAL_SKIP
fig, ax = plt.subplots()

# generate the curves and plot them
x = np.arange(0.5,10.0,0.1) # in 1/microns
Rvs = [2.0,3.0,4.0,5.0,6.0]
for cur_Rv in Rvs:
    ext_model = Cardelli89(Rv=cur_Rv)
    ax.plot(x,ext_model(x),label='R(V) = ' + str(cur_Rv))

ax.set_xlabel(r'$x$ [$\mu m^{-1}$]')
ax.set_ylabel(r'$A(x)/A(V)$')

# for 2nd x-axis with lambda values
axis_xs = np.array([0.1, 0.12, 0.15, 0.2, 0.3, 0.5, 1.0])
new_ticks = 1 / axis_xs
new_ticks_labels = ["%.2f" % z for z in axis_xs]
tax = ax.twiny()
tax.set_xlim(ax.get_xlim())
tax.set_xticks(new_ticks)
tax.set_xticklabels(new_ticks_labels)
tax.set_xlabel(r"$\lambda$ [$\mu$m]")

ax.legend(loc='best')

<matplotlib.legend.Legend at 0x722a582340b0>

We can now also use those models and show their effects on a black body spectrum. For that, we instantiate the Cardelli model, create a black body spectrum with astropy and apply the dust extinction with a fiducial Rv of 3.1 to the spectrum for a range of Av parameters.

# NBVAL_SKIP
# Let's import some packages
from astropy.modeling.models import BlackBody
import astropy.units as u
from matplotlib.ticker import ScalarFormatter

# NBVAL_SKIP
# initialize cardelli model with Rv=3.1
ext = Cardelli89(Rv=3.1)

# NBVAL_SKIP
# generate wavelengths between 3 and 10 microns
#    within the valid range for the Cardelli R(V) dependent model
lam = np.logspace(np.log10(3), np.log10(10.0), num=1000)

# setup the inputs for the blackbody function
wavelengths = lam*1e4 # Angstroem
temperature = 10000 # Kelvin

# NBVAL_SKIP
# get the blackbody flux
bb_lam = BlackBody(10000*u.K, scale=1.0 * u.erg / (u.cm ** 2 * u.AA * u.s * u.sr))
flux = bb_lam(wavelengths)

# NBVAL_SKIP
# get the extinguished blackbody flux for different amounts of dust
flux_ext_av05 = flux*ext.extinguish(lam, Av=0.5)
flux_ext_av15 = flux*ext.extinguish(lam, Av=1.5)
flux_ext_ebv10 = flux*ext.extinguish(lam, Ebv=1.0)

# NBVAL_SKIP
# plot the intrinsic and extinguished fluxes
fig, ax = plt.subplots()

ax.plot(wavelengths, flux, label='Intrinsic')
ax.plot(wavelengths, flux_ext_av05, label='$A(V) = 0.5$')
ax.plot(wavelengths, flux_ext_av15, label='$A(V) = 1.5$')
ax.plot(wavelengths, flux_ext_ebv10, label='$E(B-V) = 1.0$')

ax.set_xlabel('$\lambda$ [$\AA$]')
ax.set_ylabel('$Flux$')

ax.set_xscale('log')
ax.xaxis.set_major_formatter(ScalarFormatter())
ax.set_yscale('log')

ax.set_title('Example extinguishing a blackbody')

ax.legend(loc='best')
plt.tight_layout()

<>:10: SyntaxWarning: invalid escape sequence '\l'
<>:10: SyntaxWarning: invalid escape sequence '\l'
/tmp/ipykernel_521197/1777885020.py:10: SyntaxWarning: invalid escape sequence '\l'
  ax.set_xlabel('$\lambda$ [$\AA$]')

We see that the Cardelli model has some limited range in wavelength. Now let’s try the same for the Gordon et al. model which has a broader wavelength support.

# NBVAL_SKIP
# generate wavelengths between 0.092 and 31 microns
#    within the valid range for the Gordon23 R(V) dependent relationship
lam = jnp.logspace(np.log10(0.092), np.log10(31.0), num=1000)

# setup the inputs for the blackbody function
wavelengths = lam*1e4 # Angstroem
temperature = 10000 # Kelvin

# get the blackbody flux
bb_lam = BlackBody(10000*u.K, scale=1.0 * u.erg / (u.cm ** 2 * u.AA * u.s * u.sr))
flux = bb_lam(wavelengths)

# initialize the model
ext = Gordon23(Rv=3.1)

# get the extinguished blackbody flux for different amounts of dust
flux_ext_av05 = flux*ext.extinguish(lam, Av=0.5)
flux_ext_av15 = flux*ext.extinguish(lam, Av=1.5)
flux_ext_ebv10 = flux*ext.extinguish(lam, Ebv=1.0)

# plot the intrinsic and extinguished fluxes
fig, ax = plt.subplots()

ax.plot(wavelengths, flux, label='Intrinsic')
ax.plot(wavelengths, flux_ext_av05, label='$A(V) = 0.5$')
ax.plot(wavelengths, flux_ext_av15, label='$A(V) = 1.5$')
ax.plot(wavelengths, flux_ext_ebv10, label='$E(B-V) = 1.0$')

ax.set_xlabel(r'$\lambda$ [$\AA$]')
ax.set_ylabel('$Flux$')

ax.set_xscale('log')
ax.xaxis.set_major_formatter(ScalarFormatter())
ax.set_yscale('log')

ax.set_title('Example extinguishing a blackbody')

ax.legend(loc='best')
plt.tight_layout()
plt.show()

We see, as expected, the impact of dust is most important for short wavelength, i.e. the blue part of the spectrum.

Run the RUBIX pipeline with dust

We now turn to running the RUBIX pipeline with dust included. For this, we first need to setup the config accordingly. That is as easy as replacing "pipeline":{"name": "calc_ifu"} with "pipeline":{"name": "calc_dusty_ifu"} in the config.

In order to comapre a dusty and non dusty IFU cube, we first run a normal RUBIX pipeline.

#NBVAL_SKIP

import matplotlib.pyplot as plt
from rubix.core.pipeline import RubixPipeline 
import os
config = {
    "pipeline":{"name": "calc_ifu"},
    
    "logger": {
        "log_level": "DEBUG",
        "log_file_path": None,
        "format": "%(asctime)s - %(name)s - %(levelname)s - %(message)s",
    },
    "data": {
        "name": "IllustrisAPI",
        "args": {
            "api_key": os.environ.get("ILLUSTRIS_API_KEY"),
            "particle_type": ["stars", "gas"],
            "simulation": "TNG50-1",
            "snapshot": 99,
            "save_data_path": "data",
        },
        
        "load_galaxy_args": {
        "id": 11,
        "reuse": False,
        },
        
        "subset": {
            "use_subset": True,
            "subset_size": 50000,
        },
    },
    "simulation": {
        "name": "IllustrisTNG",
        "args": {
            "path": "data/galaxy-id-11.hdf5",
        },
    
    },
    "output_path": "output",

    "telescope":
        {"name": "MUSE",
         "psf": {"name": "gaussian", "size": 5, "sigma": 0.6},
         "lsf": {"sigma": 0.5},
         "noise": {"signal_to_noise": 1,"noise_distribution": "normal"},},
    "cosmology":
        {"name": "PLANCK15"},
        
    "galaxy":
        {"dist_z": 0.1,
         "rotation": {"type": "edge-on"},
        },
        
    "ssp": {
        "template": {
            "name": "BruzualCharlot2003"
        },
        "dust": {
            "extinction_model": "Cardelli89", #"Gordon23", 
            "dust_to_gas_ratio": 0.01, # need to check Remyer's paper
            "dust_to_metals_ratio": 0.4, # do we need this ratio if we set the dust_to_gas_ratio?
            "dust_grain_density": 3.5, # g/cm^3 #check this value
            "Rv": 3.1,
        },
    },        
}

2025-11-11 10:18:11,641 - rubix - INFO - 
   ___  __  _____  _____  __
  / _ \/ / / / _ )/  _/ |/_/
 / , _/ /_/ / _  |/ /_>  <
/_/|_|\____/____/___/_/|_|


2025-11-11 10:18:11,642 - rubix - INFO - Rubix version: 0.0.post626+g42b4b7505.d20251110
2025-11-11 10:18:11,643 - rubix - INFO - JAX version: 0.7.2
2025-11-11 10:18:11,643 - rubix - INFO - Running on [CpuDevice(id=0)] devices
2025-11-11 10:18:11,644 - rubix - WARNING - python-fsps is not installed. Please install it to use this function. Install using pip install fsps and check the installation page: https://dfm.io/python-fsps/current/installation/ for more details. Especially, make sure to set all necessary environment variables.

#NBVAL_SKIP
pipe = RubixPipeline(config)

inputdata = pipe.prepare_data()
rubixdata = pipe.run_sharded(inputdata)

/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml
  warnings.warn(
2025-11-11 10:18:11,911 - rubix - INFO - Getting rubix data...
2025-11-11 10:18:11,912 - rubix - INFO - Rubix galaxy file already exists, skipping conversion
/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/jax/_src/numpy/scalar_types.py:50: UserWarning: Explicitly requested dtype float64 requested in asarray is not available, and will be truncated to dtype float32. To enable more dtypes, set the jax_enable_x64 configuration option or the JAX_ENABLE_X64 shell environment variable. See https://github.com/jax-ml/jax#current-gotchas for more.
  return asarray(x, dtype=self.dtype)
/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/core/data.py:491: UserWarning: Explicitly requested dtype <class 'jax.numpy.float64'> requested in array is not available, and will be truncated to dtype float32. To enable more dtypes, set the jax_enable_x64 configuration option or the JAX_ENABLE_X64 shell environment variable. See https://github.com/jax-ml/jax#current-gotchas for more.
  rubixdata.galaxy.center = jnp.array(data["subhalo_center"], dtype=jnp.float64)
2025-11-11 10:18:11,995 - rubix - INFO - Centering stars particles
2025-11-11 10:18:12,927 - rubix - WARNING - The Subset value is set in config. Using only subset of size 50000 for stars
2025-11-11 10:18:12,981 - rubix - INFO - Centering gas particles
2025-11-11 10:18:13,515 - rubix - WARNING - The Subset value is set in config. Using only subset of size 50000 for gas
2025-11-11 10:18:13,518 - rubix - INFO - Data loaded with 50000 star particles and 50000 gas particles.
2025-11-11 10:18:13,518 - rubix - INFO - Data preparation completed in 1.61 seconds.
2025-11-11 10:18:13,519 - rubix - INFO - Setting up the pipeline...
2025-11-11 10:18:13,520 - rubix - DEBUG - Pipeline Configuration: {'Transformers': {'rotate_galaxy': {'name': 'rotate_galaxy', 'depends_on': None, 'args': [], 'kwargs': {}}, 'filter_particles': {'name': 'filter_particles', 'depends_on': 'rotate_galaxy', 'args': [], 'kwargs': {}}, 'spaxel_assignment': {'name': 'spaxel_assignment', 'depends_on': 'filter_particles', 'args': [], 'kwargs': {}}, 'calculate_datacube_particlewise': {'name': 'calculate_datacube_particlewise', 'depends_on': 'spaxel_assignment', 'args': [], 'kwargs': {}}, 'convolve_psf': {'name': 'convolve_psf', 'depends_on': 'calculate_datacube_particlewise', 'args': [], 'kwargs': {}}, 'convolve_lsf': {'name': 'convolve_lsf', 'depends_on': 'convolve_psf', 'args': [], 'kwargs': {}}, 'apply_noise': {'name': 'apply_noise', 'depends_on': 'convolve_lsf', 'args': [], 'kwargs': {}}}}
2025-11-11 10:18:13,521 - rubix - DEBUG - Rotation Type found: edge-on
2025-11-11 10:18:13,522 - rubix - INFO - Calculating spatial bin edges...
/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml
  warnings.warn(
2025-11-11 10:18:13,535 - rubix - INFO - Getting cosmology...
2025-11-11 10:18:13,765 - rubix - INFO - Calculating spatial bin edges...
2025-11-11 10:18:13,777 - rubix - INFO - Getting cosmology...
2025-11-11 10:18:13,789 - rubix - INFO - Getting cosmology...
2025-11-11 10:18:13,831 - rubix - DEBUG - Method not defined, using default method: cubic
/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml
  warnings.warn(
2025-11-11 10:18:13,898 - rubix - DEBUG - Method not defined, using default method: cubic
/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml
  warnings.warn(
2025-11-11 10:18:14,162 - rubix - INFO - Assembling the pipeline...
2025-11-11 10:18:14,163 - rubix - INFO - Compiling the expressions...
2025-11-11 10:18:14,164 - rubix - INFO - Number of devices: 1
2025-11-11 10:18:14,343 - rubix - INFO - Rotating galaxy with alpha=90.0, beta=0.0, gamma=0.0
2025-11-11 10:18:14,344 - rubix - INFO - Rotating galaxy for simulation: IllustrisTNG
2025-11-11 10:18:14,344 - rubix - INFO - Rotating gas
2025-11-11 10:18:14,559 - rubix - INFO - Filtering particles outside the aperture...
2025-11-11 10:18:14,569 - rubix - INFO - Assigning particles to spaxels...
2025-11-11 10:18:14,609 - rubix - INFO - Calculating Data Cube (combined per‐particle)…
2025-11-11 10:18:14,787 - rubix - DEBUG - Datacube shape: (25, 25, 3721)
2025-11-11 10:18:14,787 - rubix - INFO - Convolving with PSF...
2025-11-11 10:18:14,792 - rubix - INFO - Convolving with LSF...
2025-11-11 10:18:14,798 - rubix - INFO - Applying noise to datacube with signal to noise ratio: 1 and noise distribution: normal
2025-11-11 10:18:16,985 - rubix - INFO - Total time for sharded pipeline run: 3.47 seconds.

Now we run the pipeline including the effects of dust.

Next to setting "pipeline":{"name": "calc_ifu"} there are some more nobs under the section ssp for dust that we can tweek if needed.

Options to consider are as follows:

• the exact “extinction_model” to use. Currently Rubix supports “Cardelli89” or “Gordon23”

• the “dust_to_gas_model” to use. This currently refers to the fitting formula used by Remy-Ruyer et al. 2014. See their Table 1 for more info.

• the “Xco” model used by Remy-Ruyer et al 2014. Either “Z” or “MW”

• the “dust_grain_density” which depends on the type of dust at hand, see e.g. the NIST tables.

• the “Rv” value in case one uses an Rv dependent dust model.

#NBVAL_SKIP

import matplotlib.pyplot as plt
from rubix.core.pipeline import RubixPipeline 
import os
config = {
    "pipeline":{"name": "calc_dusty_ifu"},
    
    "logger": {
        "log_level": "DEBUG",
        "log_file_path": None,
        "format": "%(asctime)s - %(name)s - %(levelname)s - %(message)s",
    },
    "data": {
        "name": "IllustrisAPI",
        "args": {
            "api_key": os.environ.get("ILLUSTRIS_API_KEY"),
            "particle_type": ["stars", "gas"],
            "simulation": "TNG50-1",
            "snapshot": 99,
            "save_data_path": "data",
        },
        
        "load_galaxy_args": {
        "id": 11,
        "reuse": True,
        },
        
        "subset": {
            "use_subset": True,
            "subset_size": 50000,
        },
    },
    "simulation": {
        "name": "IllustrisTNG",
        "args": {
            "path": "data/galaxy-id-11.hdf5",
        },
    
    },
    "output_path": "output",

    "telescope":
        {"name": "MUSE",
         "psf": {"name": "gaussian", "size": 5, "sigma": 0.6},
         "lsf": {"sigma": 0.5},
         "noise": {"signal_to_noise": 1,"noise_distribution": "normal"},},
    "cosmology":
        {"name": "PLANCK15"},
        
    "galaxy":
        {"dist_z": 0.1,
         "rotation": {"type": "edge-on"},
        },
        
    "ssp": {
        "template": {
            "name": "BruzualCharlot2003"
        },
        "dust": {
            "extinction_model": "Cardelli89", #"Gordon23", 
            "dust_to_gas_model": "broken power law fit", # from Remyer's paper see their Table 1
            "Xco": "Z", # from Remyer's paper, see their Table 1
            "dust_grain_density": 3.0, # #check this value, reverse engeneered from Ibarrra-Medel 2018
            "Rv": 3.1,
        },
    },        
}

#NBVAL_SKIP
pipe = RubixPipeline(config)

inputdata = pipe.prepare_data()
rubixdata_dust = pipe.run_sharded(inputdata)

/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml
  warnings.warn(
2025-11-11 10:18:23,421 - rubix - INFO - Getting rubix data...
2025-11-11 10:18:23,424 - rubix - INFO - Rubix galaxy file already exists, skipping conversion
/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/jax/_src/numpy/scalar_types.py:50: UserWarning: Explicitly requested dtype float64 requested in asarray is not available, and will be truncated to dtype float32. To enable more dtypes, set the jax_enable_x64 configuration option or the JAX_ENABLE_X64 shell environment variable. See https://github.com/jax-ml/jax#current-gotchas for more.
  return asarray(x, dtype=self.dtype)
/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/core/data.py:491: UserWarning: Explicitly requested dtype <class 'jax.numpy.float64'> requested in array is not available, and will be truncated to dtype float32. To enable more dtypes, set the jax_enable_x64 configuration option or the JAX_ENABLE_X64 shell environment variable. See https://github.com/jax-ml/jax#current-gotchas for more.
  rubixdata.galaxy.center = jnp.array(data["subhalo_center"], dtype=jnp.float64)
2025-11-11 10:18:23,452 - rubix - INFO - Centering stars particles
2025-11-11 10:18:23,777 - rubix - WARNING - The Subset value is set in config. Using only subset of size 50000 for stars
2025-11-11 10:18:23,783 - rubix - INFO - Centering gas particles
2025-11-11 10:18:23,834 - rubix - WARNING - The Subset value is set in config. Using only subset of size 50000 for gas
2025-11-11 10:18:23,835 - rubix - INFO - Data loaded with 50000 star particles and 50000 gas particles.
2025-11-11 10:18:23,836 - rubix - INFO - Data preparation completed in 0.42 seconds.
2025-11-11 10:18:23,837 - rubix - INFO - Setting up the pipeline...
2025-11-11 10:18:23,838 - rubix - DEBUG - Pipeline Configuration: {'Transformers': {'rotate_galaxy': {'name': 'rotate_galaxy', 'depends_on': None, 'args': [], 'kwargs': {}}, 'filter_particles': {'name': 'filter_particles', 'depends_on': 'rotate_galaxy', 'args': [], 'kwargs': {}}, 'spaxel_assignment': {'name': 'spaxel_assignment', 'depends_on': 'filter_particles', 'args': [], 'kwargs': {}}, 'calculate_extinction': {'name': 'calculate_extinction', 'depends_on': 'spaxel_assignment', 'args': [], 'kwargs': {}}, 'calculate_dusty_datacube_particlewise': {'name': 'calculate_dusty_datacube_particlewise', 'depends_on': 'calculate_extinction', 'args': [], 'kwargs': {}}, 'convolve_psf': {'name': 'convolve_psf', 'depends_on': 'calculate_dusty_datacube_particlewise', 'args': [], 'kwargs': {}}, 'convolve_lsf': {'name': 'convolve_lsf', 'depends_on': 'convolve_psf', 'args': [], 'kwargs': {}}, 'apply_noise': {'name': 'apply_noise', 'depends_on': 'convolve_lsf', 'args': [], 'kwargs': {}}}}
2025-11-11 10:18:23,838 - rubix - DEBUG - Rotation Type found: edge-on
2025-11-11 10:18:23,842 - rubix - INFO - Calculating spatial bin edges...
/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml
  warnings.warn(
2025-11-11 10:18:23,866 - rubix - INFO - Getting cosmology...
2025-11-11 10:18:23,886 - rubix - INFO - Calculating spatial bin edges...
2025-11-11 10:18:23,899 - rubix - INFO - Getting cosmology...
2025-11-11 10:18:23,912 - rubix - INFO - Getting cosmology...
2025-11-11 10:18:23,934 - rubix - DEBUG - Method not defined, using default method: cubic
/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml
  warnings.warn(
2025-11-11 10:18:23,969 - rubix - DEBUG - Method not defined, using default method: cubic
/home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/factory.py:26: UserWarning: No telescope config provided, using default stored in /home/annalena/.conda/envs/rubix/lib/python3.12/site-packages/rubix/telescope/telescopes.yaml
  warnings.warn(
2025-11-11 10:18:23,995 - rubix - INFO - Assembling the pipeline...
2025-11-11 10:18:23,996 - rubix - INFO - Compiling the expressions...
2025-11-11 10:18:23,997 - rubix - INFO - Number of devices: 1
2025-11-11 10:18:24,193 - rubix - INFO - Rotating galaxy with alpha=90.0, beta=0.0, gamma=0.0
2025-11-11 10:18:24,194 - rubix - INFO - Rotating galaxy for simulation: IllustrisTNG
2025-11-11 10:18:24,195 - rubix - INFO - Rotating gas
2025-11-11 10:18:24,364 - rubix - INFO - Filtering particles outside the aperture...
2025-11-11 10:18:24,370 - rubix - INFO - Assigning particles to spaxels...
2025-11-11 10:18:24,374 - rubix - INFO - Applying dust extinction to the spaxel data...
2025-11-11 10:18:24,375 - rubix - INFO - Applying dust extinction to the spaxel data using vmap...
2025-11-11 10:18:24,481 - rubix - INFO - Calculating Data Cube (combined per‐particle)…
2025-11-11 10:18:24,508 - rubix - DEBUG - Datacube shape: (25, 25, 3721)
2025-11-11 10:18:24,508 - rubix - INFO - Convolving with PSF...
2025-11-11 10:18:24,511 - rubix - INFO - Convolving with LSF...
2025-11-11 10:18:24,514 - rubix - INFO - Applying noise to datacube with signal to noise ratio: 1 and noise distribution: normal
2025-11-11 10:18:26,659 - rubix - INFO - Total time for sharded pipeline run: 2.82 seconds.

Let’s compare one example spaxel spectrum with and without dust.

#NBVAL_SKIP
wave = pipe.telescope.wave_seq

spectra = rubixdata # Spectra of all stars
dusty_spectra = rubixdata_dust # Spectra of all stars
print(spectra.shape)
print(dusty_spectra.shape)

plt.plot(wave, spectra[12,12,:])
plt.plot(wave, dusty_spectra[12,12,:])

(25, 25, 3721)
(25, 25, 3721)

[<matplotlib.lines.Line2D at 0x7229802724e0>]

Now let’s visualize a nice edge-on galaxy in SDSS broad-band images with some nice dust lanes…

# NBVAL_SKIP
from rubix.telescope.filters import load_filter, convolve_filter_with_spectra

# NBVAL_SKIP
# load all fliter curves for SLOAN
curves = load_filter("SLOAN")

# NBVAL_SKIP
wave = pipe.telescope.wave_seq
filters,images = curves.apply_filter_curves(rubixdata_dust, wave).values()

for i_dust,name in zip(images, filters):
    plt.figure()
    plt.imshow(i_dust)
    plt.colorbar()
    plt.title(name)

Sanity check: overlay gas column density map over the dusty emission image

# NBVAL_SKIP
# The input data are rotated in the same way as the particles are rotaterd to calculate the IFU. 
# This step is necessary, because we only have the raw input data and the pipeline only returns the datacube and not the per particle information.
from rubix.core.rotation import get_galaxy_rotation
rotate = get_galaxy_rotation(config)

inputdatadata = rotate(inputdata)

2025-11-11 10:18:53,867 - rubix - DEBUG - Rotation Type found: edge-on
2025-11-11 10:18:53,869 - rubix - INFO - Rotating galaxy with alpha=90.0, beta=0.0, gamma=0.0
2025-11-11 10:18:53,870 - rubix - INFO - Rotating galaxy for simulation: IllustrisTNG
2025-11-11 10:18:53,870 - rubix - INFO - Rotating gas

# NBVAL_SKIP
idx = np.where(inputdata.gas.mass != 0)
gas_map = np.histogram2d(inputdata.gas.coords[:,0][idx], inputdata.gas.coords[:,1][idx], bins=(25,25), weights=np.squeeze(inputdata.gas.mass)[idx])

# NBVAL_SKIP
plt.figure()
plt.imshow(gas_map[0].T, cmap='inferno')
plt.colorbar()
plt.title("gas map")

Text(0.5, 1.0, 'gas map')

# NBVAL_SKIP
plt.figure()
plt.imshow(i_dust)
plt.imshow(gas_map[0].T, cmap='inferno', alpha=0.6)
plt.colorbar()
plt.title("emission and gas map overlayed")

Text(0.5, 1.0, 'emission and gas map overlayed')

And in comparison to this, the same galaxy without dust…

# NBVAL_SKIP
wave = pipe.telescope.wave_seq
filters,images = curves.apply_filter_curves(rubixdata, wave).values()

for i,name in zip(images, filters):
    plt.figure()
    plt.imshow(i)
    plt.colorbar()
    plt.title(name)