diff --git a/examples/cosmo/luminosity_function_integration_test.py b/examples/cosmo/luminosity_function_integration_test.py new file mode 100644 index 0000000000..6439ac9dd7 --- /dev/null +++ b/examples/cosmo/luminosity_function_integration_test.py @@ -0,0 +1,97 @@ +""" +Camels luminosity function example +================================== + +Use test cosmological simulation data (from the [CAMELS simulations]( +https://www.camel-simulations.org/)) to generate SDSS photometry and +plot the luminosity functions. + +Note: this is an integration test, and requires the full camels +snapshot and subhalo catalogue. These can be downloaded from the +tests/data/ folder using the `download_camels.sh` script (an internet +connection is required). +""" + +import matplotlib.pyplot as plt +import numpy as np + +from synthesizer.conversions import lnu_to_absolute_mag +from synthesizer.emission_models import IncidentEmission +from synthesizer.filters import FilterCollection +from synthesizer.grid import Grid +from synthesizer.load_data.load_camels import load_CAMELS_IllustrisTNG + + +def calc_df(x, volume, binLimits): + hist, dummy = np.histogram(x, bins=binLimits) + hist = np.float64(hist) + phi = (hist / volume) / (binLimits[1] - binLimits[0]) + + phi_sigma = (np.sqrt(hist) / volume) / ( + binLimits[1] - binLimits[0] + ) # Poisson errors + + return phi, phi_sigma, hist + + +h = 0.6711 +grid_dir = "../../tests/test_grid" +grid_name = "test_grid" +grid = Grid(grid_name, grid_dir=grid_dir) +incident = IncidentEmission(grid) + +filter_codes = [f"SLOAN/SDSS.{f}" for f in ["g"]] +fc = FilterCollection(filter_codes=filter_codes, new_lam=grid.lam) + +gals = load_CAMELS_IllustrisTNG( + "../../tests/data/", + snap_name="CV_0_snap_086.hdf5", + group_name="CV_0_fof_subhalo_tab_086.hdf5", +) + +# Filter by stellar mass +mstar_mask = np.array([np.sum(g.stars.masses) > 1e9 for g in gals]) +gals = [g for g, _m in zip(gals, mstar_mask) if _m] + +# Calculate g-band magnitudes +specs = np.vstack([g.stars.get_spectra(incident).lnu for g in gals]) +phot = np.vstack( + [g.stars.get_photo_lnu(fc)["incident"]["SLOAN/SDSS.g"] for g in gals] +) +mags = lnu_to_absolute_mag(phot) + +# Calculate g-band luminosity function +binlims = np.linspace(-25, -16, 12) +bins = binlims[:-1] + (binlims[1] - binlims[0]) / 2 +phi, phi_sigma, _ = calc_df(mags, (25 / h) ** 3, binlims) + +# Plot luminosity function +fig, ax = plt.subplots(1, 1) + +ax.plot(bins, np.log10(phi), label="Current version") + +# Previous 'working' version LF for reference +# Commit hash: +# b683b81e5de1e8c6ab6938d047cab54e7c5a2fdf +phi_previous = [ + -np.inf, + -np.inf, + -4.63, + -3.63, + -3.23, + -2.66, + -2.59, + -2.85, + -2.73, + -2.89, + -3.33, +] + +ax.plot(bins, phi_previous, label="Previous version", ls="dashed") + +ax.legend() +ax.set_xlabel("$m_{g}^{AB}$") +ax.set_ylabel(r"$\phi \,/\, \mathrm{Mpc^{-3} \; dex^{-1}}$") + +plt.show() +# plt.savefig('test_lf.png'); plt.close() diff --git a/examples/cosmo/plot_parametric_young_stars.py b/examples/cosmo/plot_parametric_young_stars.py index e642b922d1..f5cc0e17b3 100644 --- a/examples/cosmo/plot_parametric_young_stars.py +++ b/examples/cosmo/plot_parametric_young_stars.py @@ -35,7 +35,7 @@ model = IncidentEmission(grid) # Select a single galaxy - gal = gals[1] + gal = gals[2] # Age limit at which we replace star particles age_lim = 500 * Myr diff --git a/tests/data/camels_snap.hdf5 b/tests/data/camels_snap.hdf5 index 0db0636027..c9ca760a1c 100644 Binary files a/tests/data/camels_snap.hdf5 and b/tests/data/camels_snap.hdf5 differ diff --git a/tests/data/camels_subhalo.hdf5 b/tests/data/camels_subhalo.hdf5 index c52e22df03..f9c2776592 100644 Binary files a/tests/data/camels_subhalo.hdf5 and b/tests/data/camels_subhalo.hdf5 differ diff --git a/tests/data/create_camels_small_snap.py b/tests/data/create_camels_small_snap.py index e2a5143c20..33cbb53e79 100644 --- a/tests/data/create_camels_small_snap.py +++ b/tests/data/create_camels_small_snap.py @@ -2,10 +2,10 @@ import numpy as np _dir = "." -snap_name = "LH_1_snap_031.hdf5" -group_name = "LH_1_fof_subhalo_tab_031.hdf5" +snap_name = "CV_0_snap_086.hdf5" +group_name = "CV_0_fof_subhalo_tab_086.hdf5" N = 10 # number of galaxies to extract -ignore_N = 1 # number of galaxies to ignore +ignore_N = 10 # number of galaxies to ignore with h5py.File(f"{_dir}/{snap_name}", "r") as hf: form_time = hf["PartType4/GFM_StellarFormationTime"][:] diff --git a/tests/data/download_camels.sh b/tests/data/download_camels.sh index ddc81ee9b8..b11343717e 100755 --- a/tests/data/download_camels.sh +++ b/tests/data/download_camels.sh @@ -1,11 +1,11 @@ #!/bin/bash dir=./ # $1 -snap=031 +snap=086 # $2 -lh=1 +cv=0 # $3 -wget "https://users.flatironinstitute.org/~camels/Sims/IllustrisTNG/LH/LH_${lh}/snap_$snap.hdf5" --output-document "${dir}/LH_${lh}_snap_${snap}.hdf5" +wget "https://users.flatironinstitute.org/~camels/Sims/IllustrisTNG/CV/CV_${cv}/snapshot_$snap.hdf5" --output-document "${dir}/CV_${cv}_snap_${snap}.hdf5" -wget "https://users.flatironinstitute.org/~camels/FOF_Subfind/IllustrisTNG/LH/LH_${lh}/fof_subhalo_tab_$snap.hdf5" --output-document "$dir/LH_${lh}_fof_subhalo_tab_$snap.hdf5" +wget "https://users.flatironinstitute.org/~camels/FOF_Subfind/IllustrisTNG/CV/CV_${cv}/groups_$snap.hdf5" --output-document "$dir/CV_${cv}_fof_subhalo_tab_$snap.hdf5"