py:mass resolution

mass_resolution.py

# ---
# jupyter:
#   jupytext:
#     formats: ipynb,py:percent
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.11.3
#   kernelspec:
#     display_name: Python [conda env:std]
#     language: python
#     name: conda-env-std-py
# ---

# %%
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from astropy import constants as C
from astropy import cosmology
from astropy import units as U

# %%
cosmo = cosmology.Planck18
sns.set_theme("talk", "ticks", font_scale=1.1)


# %%
def delta_crit(z: float = 0):
    d = cosmo.Om(z) - 1.0
    return 18.0 * (np.pi ** 2) + 82.0 * d - 39.0 * (d ** 2)


def virial_mass(redshift: float, temp: U.Quantity = 1.0e4 * U.Kelvin):
    mass = (2.0 * temp * C.k_B / (C.G * 0.59 * C.m_p)) ** 1.5 * np.sqrt(
        3.0
        / (
            delta_crit(redshift)
            * 4.0
            * np.pi
            * cosmo.critical_density(redshift)
            * cosmo.Om(redshift)
        )
    )
    return mass.to(U.Msun)


# %%
TIAMAT_PARTMASS = 3.899487e6 * U.Msun
GENESIS_L210_N4320_PARTMASS = 1.47436e7 * U.Msun
HALO_NP = 30
ZCOOL_TARGET = 20.0

# %%
redshift = np.linspace(4.0, 21.0, 100)
mass = np.log10(virial_mass(redshift).value)

fig, ax = plt.subplots(1, 1, figsize=(12, 6))
ax.plot(redshift, mass, lw=4)
ax.fill_between(redshift, mass, 7.0, alpha=0.2)

tiamat_res = np.log10(TIAMAT_PARTMASS.to(U.Msun).value * HALO_NP)
ax.axhline(tiamat_res, ls="--", lw=3, color=sns.color_palette()[1])
ax.annotate(
    r"Tiamat $\left(h^{-1}67.8\, {\rm Mpc} \right)^3, N_p = 2160^3$",
    (19.5, tiamat_res + 0.05),
    ha="right",
    size="small",
)

genesis_L210_N4320_res = np.log10(
    GENESIS_L210_N4320_PARTMASS.to(U.Msun).value * HALO_NP
)
ax.axhline(genesis_L210_N4320_res, ls="--", lw=3, color=sns.color_palette()[2])
ax.annotate(
    r"Genesis $\left( h^{-1}210\, {\rm Mpc} \right)^3, N_p = 4320^3$",
    (19.5, genesis_L210_N4320_res + 0.05),
    ha="right",
    size="small",
)

extended_res = np.interp(ZCOOL_TARGET, redshift, mass)
npart = (10.0 ** (genesis_L210_N4320_res - extended_res)) * 4320 ** 3

ax.axhline(extended_res, ls="--", lw=3, color=sns.color_palette()[3])
ax.annotate(
    r"Genesis - extended $\mathbf{\left( h^{-1}210\, {\rm \bf Mpc} \right)^3, N_p^{\rm \bf eff} = "
    + f"{np.cbrt(npart):.0f}"
    + r"^3}$",
    (5.5, extended_res + 0.05),
    weight="bold",
    ha="left",
    size="small",
)

ax.set(
    ylabel=r"$\log_{10}\left(M / {\rm M_\odot}\right)$",
    xlabel="redshift",
    xlim=(5.0, 20.0),
    ylim=(7.25, 9.0),
)
ax.xaxis.set_major_locator(plt.MultipleLocator(2.0))
ax.set_title(r"$M_{\rm cool}\left(T_{\rm vir} \sim 10^4\, {\rm K}\right)$", loc="left")
sns.despine(ax=ax)

fname = "../plots/mass_resolution.pdf"
plt.savefig(fname, bbox_inches="tight", pad_inches=0.05)

# %%
# !realpath $fname