Compare SPS grid assignment methods

This example compares a the cloud in cell (CIC) and nearest grid point (NGP) grid assignment methods. These methods dictate how mass is assigned to spectra in the SPS grids.

/home/runner/work/synthesizer/synthesizer/src/synthesizer/particle/galaxy.py:251: RuntimeWarning: Current mass of stars not provided, setting stellar_mass_weighted_age to `None`
  self.calculate_integrated_stellar_properties()
/home/runner/work/synthesizer/synthesizer/src/synthesizer/particle/galaxy.py:111: RuntimeWarning: In `load_gas`: one of either `masses` or `metallicities` is not provided, setting `gas` object to `None`
  self.load_gas(gas=gas)
/home/runner/work/synthesizer/synthesizer/src/synthesizer/particle/galaxy.py:118: RuntimeWarning: Current mass of stars not provided, setting stellar_mass_weighted_age to `None`
  self.calculate_integrated_stellar_properties()
/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/unyt/array.py:1949: RuntimeWarning: invalid value encountered in divide
  out_arr = func(

import matplotlib.pyplot as plt
from synthesizer.grid import Grid
from synthesizer.parametric import SFH, ZDist
from synthesizer.parametric import Stars as ParametricStars
from synthesizer.particle.galaxy import Galaxy as ParticleGalaxy
from synthesizer.particle.stars import sample_sfhz
from unyt import Myr

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

# Define the SFH and metallicity distribution
Z_p = {"metallicity": 0.01}
metal_dist = ZDist.DeltaConstant(**Z_p)
sfh_p = {"duration": 100 * Myr}
sfh = SFH.Constant(**sfh_p)

# Define the parametric stars
sfzh = ParametricStars(
    grid.log10age,
    grid.metallicity,
    sf_hist=sfh,
    metal_dist=metal_dist,
    initial_mass=10**9,
)

# How many particles?
nstar = 10**5

# Get the stars object
stars = sample_sfhz(
    sfzh.sfzh,
    sfzh.log10ages,
    sfzh.log10metallicities,
    nstar,
    initial_mass=10**9 / nstar,
)

# Create galaxy object
particle_galaxy = ParticleGalaxy(stars=stars)

# Calculate the stars SEDs using both grid assignment schemes
cic_sed = particle_galaxy.stars.get_spectra_incident(
    grid,
    grid_assignment_method="cic",
)
ngp_sed = particle_galaxy.stars.get_spectra_incident(
    grid,
    grid_assignment_method="ngp",
)

# Setup the plot
fig = plt.figure()
ax = fig.add_subplot(111)
resi_ax = ax.twinx()
ax.grid(True)

resi_ax.semilogx(
    ngp_sed.lam,
    ngp_sed.lnu / cic_sed.lnu,
    color="m",
    linestyle="dotted",
    label="Residual",
    alpha=0.6,
)
ax.loglog(
    cic_sed.lam,
    cic_sed.lnu,
    label="CIC",
)
ax.loglog(
    ngp_sed.lam,
    ngp_sed.lnu,
    label="NGP",
)


resi_ax.set_ylabel("NGP / CIC")
x_units = str(cic_sed.lam.units)
y_units = str(cic_sed.lnu.units)
ax.set_xlabel(r"$\lambda/[\mathrm{" + x_units + r"}]$")
ax.set_ylabel(r"$L_{\nu}/[\mathrm{" + y_units + r"}]$")

ax.legend()
resi_ax.legend(loc="upper left")
ax.set_xlim(10**2.0, 10**5.2)
ax.set_ylim(10**25.0, 10**30.0)
plt.show()