{ "cells": [ { "cell_type": "markdown", "id": "21609608-8da1-40d1-bf4d-12d8ca84f267", "metadata": {}, "source": [ "# Zeus21 with Population II and III Stars" ] }, { "cell_type": "markdown", "id": "5f2e29a3", "metadata": {}, "source": [ "This quick tutorial covers all the new features in the updated public version of Zeus21 (see our paper [Cruz et al. 2024](https://arxiv.org/abs/2407.18294)) for more information.\n", "\n", "In this new version, we have included the influence of Population III stars as well as the suppressive feedback effects from Lyman-Werner radiation and the relative velocity between CDM and baryons." ] }, { "cell_type": "code", "execution_count": 1, "id": "2c9e8ec0", "metadata": {}, "outputs": [], "source": [ "import math\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "\n", "#import Zeus\n", "import zeus21\n", "\n", "#set up the CLASS cosmology\n", "from classy import Class\n", "\n", "#and the user parameters\n", "UserParams = zeus21.User_Parameters(precisionboost=1.2)\n" ] }, { "cell_type": "markdown", "id": "f3646152-80fa-4166-bcd8-1645ffbebe5b", "metadata": {}, "source": [ "## Step 1: Set up Cosmology" ] }, { "cell_type": "markdown", "id": "e56fce7d", "metadata": {}, "source": [ "We begin by running CLASS, where the associated parameters can be altered below. Then we save the cosmo parameters, the correlation functions, and the halo mass function at all desired z and M." ] }, { "cell_type": "code", "execution_count": 2, "id": "01adaa46-dcdd-4405-a6a1-5b800584dc4d", "metadata": {}, "outputs": [], "source": [ "# cosmo inputs from table 2 last column of 1807.06209, as 21cmFAST\n", "ombh2 = 0.02242 \n", "omch2 = 0.11933\n", "tau_re = 0.0544\n", "hLittle = 0.6766\n", "ns = 0.9665\n", "As = np.exp(3.047)*10**(-10.)\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "b5c07181", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CLASS has run, we store the cosmology.\n", "HMF interpolator built. This ends the cosmology part -- moving to astrophysics.\n" ] } ], "source": [ "#set up user parameters\n", "UserParams = zeus21.User_Parameters(FLAG_FORCE_LINEAR_CF=False,zmin_T21=10.)\n", "\n", "#define all cosmology (including derived) parameters, and save them to the CosmoParams structure\n", "CosmoParams = zeus21.Cosmo_Parameters(UserParams=UserParams, omegab = ombh2, omegac = omch2, h_fid = hLittle, As = As, ns = ns, tau_fid = tau_re, USE_RELATIVE_VELOCITIES = True, Flag_emulate_21cmfast=False)\n", "print('CLASS has run, we store the cosmology.')\n", "\n", "# Compute the HMF structure that stores associated quantities\n", "HMFintclass = zeus21.HMF_interpolator(User_Parameters=UserParams,Cosmo_Parameters=CosmoParams)\n", "print('HMF interpolator built. This ends the cosmology part -- moving to astrophysics.')" ] }, { "cell_type": "markdown", "id": "8eeda1be-f47c-4cb4-b9f2-6998e06c6641", "metadata": {}, "source": [ "## Step 2: Astrophysics" ] }, { "cell_type": "markdown", "id": "faf27bcd", "metadata": {}, "source": [ "After that we set up the astro parameters, calculate the SFRD and with it all the global signal and related quantities. In principle one can re-run the astrophysics part only if you're certain of the cosmology parameters. \n", "\n", "With the current implementation, the computation of global quantities below takes less than a second on a laptop.\n", "\n", "For more details on parameters, please refer to those used in Table 1 of our paper linked above.\n", "\n", "NOTE: Pop III stars and LW feedback are turned on using the USE_POP_III and USE_LW_FEEDBACK flags below." ] }, { "cell_type": "code", "execution_count": null, "id": "fbcff876-3535-475d-9cb0-c124b9de2fff", "metadata": { "scrolled": true, "tags": [] }, "outputs": [], "source": [ "### POP II quantities first\n", "\n", "################################\n", "### Model Parameters\n", "accretion_model = \"exp\" # Accretion model. \"exp\" for exponential, \"EPS\" for EPS. \"RP16\" for the dynamically averaged fitting function in Rodríguez-Puebla+16. Default is \"exp\"\n", "\n", "################################\n", "### SFR(Mh) Parameteres\n", "alphastar = 0.69 # alphastar powerlaw index for low masses, default 0.5\n", "betastar = -1.68 # betastar powerlaw index for high masses, default -0.5\n", "epsstar = 10**-1.11 # epsilonstar = fstar at Mc\n", "Mc = 10**11.93 # Pivot mass at which the power law cuts for model 0, default Mc = 3e11\n", "dlog10epsstardz = -0.08 # dlog10epsilonstar/dz, default 0\n", "\n", "################################\n", "### Escape fraction parameters\n", "fesc10 = 0.1 # fesc(M) parameter. Power law normalized (fesc10) at M=1e10 Msun with index alphaesc\n", "alphaesc = 0.0\n", "L40_xray = 10**0.5 # L40_xray: soft-band (E<2 keV) lum/SFR in Xrays in units of 10^40 erg/s/(Msun/yr)\n", "E0_xray = 500. # E0_xray: minimum energy in eV\n", "alpha_xray = -1.0 # Xray SED power-law index\n", "Emax_xray_norm=2000 # max energy in eV to normalize SED. Keep at 2000 eV normally\n", "\n", "################################\n", "### LyA parameters\n", "N_alpha_perbaryon_II = 9690 # number of Pop II photons between LyA and Ly Cont. per baryon (from BL05)\n", "N_alpha_perbaryon_III = 17900 # number of Pop III photons between LyA and Ly Cont. per baryon value of 17900 is from Klessen & Glover 2023 (2303.12500), table A2\n", "\n", "################################\n", "### MTURN Parameters: \n", "Mturn_fixed = None # Mturn_fixed: None if use Matom(z) at each z, Some value if fixed Mturn\n", "FLAG_MTURN_SHARP= False # Mturn_sharp: False if regular exponential cutoff, True if sharp cutoff, active only if Mturn_fixed is on\n", "\n", "\n", "\n", "\n", "################################\n", "# Pop III Quantities\n", "alphastar_III = 0 \n", "betastar_III = 0\n", "epsstar_III = 10**-2.5709708032788794/3 #AV sugmaUV = 1.3\n", "Mc_III = 1e7\n", "dlog10epsstardz_III = 0.0\n", "\n", "fesc7_III = 10**(-1.35)\n", "alphaesc_III = -0.3\n", "L40_xray_III = 10**0.5\n", "alpha_xray_III = -1.0\n", "\n", "\n", "USE_POPIII = True\n", "USE_LW_FEEDBACK = True\n", "\n", "A_LW = 2.0\n", "beta_LW = 0.6\n", "\n", "A_vcb = 1.0\n", "beta_vcb = 1.8\n", "\n", "\n", "\n", "AstroParams = zeus21.Astro_Parameters(CosmoParams=CosmoParams, \n", " quadratic_SFRD_lognormal=False,\n", " FLAG_USE_PSD=False, \n", " \n", " \n", " accretion_model = accretion_model,\n", " \n", " alphastar = alphastar, \n", " betastar = betastar, \n", " epsstar = epsstar, \n", " Mc = Mc, \n", " dlog10epsstardz = dlog10epsstardz,\n", " \n", " fesc10 = fesc10, \n", " alphaesc = alphaesc,\n", " L40_xray = L40_xray, \n", " E0_xray = E0_xray, \n", " alpha_xray = alpha_xray, \n", " Emax_xray_norm = Emax_xray_norm, \n", " \n", " N_alpha_perbaryon_II = N_alpha_perbaryon_II, \n", " N_alpha_perbaryon_III = N_alpha_perbaryon_III,\n", " \n", " Mturn_fixed = Mturn_fixed, \n", " FLAG_MTURN_SHARP = FLAG_MTURN_SHARP,\n", " \n", " \n", " USE_POPIII = USE_POPIII, \n", " USE_LW_FEEDBACK = USE_LW_FEEDBACK,\n", "\n", " alphastar_III = alphastar_III, \n", " betastar_III = betastar_III,\n", " epsstar_III = epsstar_III,\n", " Mc_III = Mc_III,\n", " dlog10epsstardz_III = dlog10epsstardz_III,\n", "\n", " fesc7_III = fesc7_III,\n", " alphaesc_III = alphaesc_III,\n", " L40_xray_III = L40_xray_III,\n", " alpha_xray_III = alpha_xray_III,\n", " \n", " A_LW = A_LW,\n", " beta_LW = beta_LW,\n", " \n", " A_vcb = A_vcb,\n", " beta_vcb = beta_vcb,\n", " )\n", "\n", "\n", "CoeffStructure = zeus21.get_T21_coefficients(UserParams=UserParams, CosmoParams=CosmoParams ,AstroParams=AstroParams, HMFinterp=HMFintclass)\n", "SFRD_class = zeus21.sfrd.SFRD_class(UserParams, CosmoParams, AstroParams, HMFintclass)\n", "zlist= CoeffStructure.zintegral\n", "print('SFRD and coefficients stored. Move ahead.')\n" ] }, { "cell_type": "code", "execution_count": null, "id": "4cb4b017", "metadata": {}, "outputs": [], "source": [ "### POP II quantities first\n", "\n", "################################\n", "### Model Parameters\n", "accretion_model = \"exp\" # ACCRETION MODEL: 0 for exponential, 1 for EPS, default EXP\n", "\n", "################################\n", "### SFR(Mh) Parameteres\n", "alphastar = 0.5 # alphastar powerlaw index for low masses, default 0.5\n", "betastar = -0.5 # betastar powerlaw index for high masses, default -0.5\n", "epsstar = 10**-1. # epsilonstar = fstar at Mc\n", "Mc = 3e11 # Pivot mass at which the power law cuts for model 0, default Mc = 3e11\n", "dlog10epsstardz = 0.0 # dlog10epsilonstar/dz, default 0\n", "\n", "################################\n", "### Escape fraction parameters\n", "fesc10 = 0.1 # fesc(M) parameter. Power law normalized (fesc10) at M=1e10 Msun with index alphaesc\n", "alphaesc = 0.0\n", "L40_xray = 10**0.5 # L40_xray: soft-band (E<2 keV) lum/SFR in Xrays in units of 10^40 erg/s/(Msun/yr)\n", "E0_xray = 500. # E0_xray: minimum energy in eV\n", "alpha_xray = -1.0 # Xray SED power-law index\n", "Emax_xray_norm=2000 # max energy in eV to normalize SED. Keep at 2000 eV normally\n", "\n", "################################\n", "### LyA parameters\n", "N_alpha_perbaryon_II = 9690 # number of Pop II photons between LyA and Ly Cont. per baryon (from BL05)\n", "N_alpha_perbaryon_III = 17900 # number of Pop III photons between LyA and Ly Cont. per baryon value of 17900 is from Klessen & Glover 2023 (2303.12500), table A2\n", "\n", "################################\n", "### MTURN Parameters: \n", "Mturn_fixed = None # Mturn_fixed: None if use Matom(z) at each z, Some value if fixed Mturn\n", "FLAG_MTURN_SHARP= False # Mturn_sharp: False if regular exponential cutoff, True if sharp cutoff, active only if Mturn_fixed is on\n", "\n", "################################\n", "ZMIN = 10.0 # down to which z we compute the evolution\n", "\n", "\n", "\n", "################################\n", "# Pop III Quantities\n", "alphastar_III = 0 \n", "betastar_III = 0\n", "epsstar_III = 10**(-3.0)\n", "Mc_III = 1e7\n", "dlog10epsstardz_III = 0.0\n", "\n", "fesc7_III = 10**(-1.35)\n", "alphaesc_III = -0.3\n", "L40_xray_III = 10**0.5\n", "alpha_xray_III = -1.0\n", "\n", "\n", "USE_POPIII = True\n", "USE_LW_FEEDBACK = True\n", "\n", "A_LW = 2.0\n", "beta_LW = 0.6\n", "\n", "A_vcb = 1.0\n", "beta_vcb = 1.8\n", "\n", "\n", "#set up your astro parameters too, here the peak of f*(Mh) as an example\n", "\n", "AstroParams = zeus21.Astro_Parameters(CosmoParams=CosmoParams, \n", " quadratic_SFRD_lognormal=False,\n", " FLAG_USE_PSD=False, \n", " \n", " accretion_model = accretion_model,\n", " \n", " alphastar = alphastar, \n", " betastar = betastar, \n", " epsstar = epsstar, \n", " Mc = Mc, \n", " dlog10epsstardz = dlog10epsstardz,\n", " \n", " fesc10 = fesc10, \n", " alphaesc = alphaesc,\n", " L40_xray = L40_xray, \n", " E0_xray = E0_xray, \n", " alpha_xray = alpha_xray, \n", " Emax_xray_norm = Emax_xray_norm, \n", " \n", " N_alpha_perbaryon_II = N_alpha_perbaryon_II, \n", " N_alpha_perbaryon_III = N_alpha_perbaryon_III,\n", " \n", " Mturn_fixed = Mturn_fixed, \n", " FLAG_MTURN_SHARP = FLAG_MTURN_SHARP,\n", " \n", " USE_POPIII = USE_POPIII, \n", " USE_LW_FEEDBACK = USE_LW_FEEDBACK,\n", "\n", " alphastar_III = alphastar_III, \n", " betastar_III = betastar_III,\n", " epsstar_III = epsstar_III,\n", " Mc_III = Mc_III,\n", " dlog10epsstardz_III = dlog10epsstardz_III,\n", "\n", " fesc7_III = fesc7_III,\n", " alphaesc_III = alphaesc_III,\n", " L40_xray_III = L40_xray_III,\n", " alpha_xray_III = alpha_xray_III,\n", " \n", " A_LW = A_LW,\n", " beta_LW = beta_LW,\n", " \n", " A_vcb = A_vcb,\n", " beta_vcb = beta_vcb)\n", "\n", "CoeffStructure = zeus21.get_T21_coefficients(UserParams=UserParams, CosmoParams=CosmoParams ,AstroParams=AstroParams, HMFinterp=HMFintclass)\n", "SFRD_class = zeus21.sfrd.SFRD_class(UserParams, CosmoParams, AstroParams, HMFintclass)\n", "zlist= CoeffStructure.zintegral\n", "print('SFRD and coefficients stored. Move ahead.')\n", "\n" ] }, { "cell_type": "markdown", "id": "fa72898f", "metadata": {}, "source": [ "## Step 3: Plotting globally-averaged quantities" ] }, { "cell_type": "code", "execution_count": null, "id": "507bde44", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize = (12, 6.75))\n", "\n", "plt.plot(zlist, CoeffStructure.T21avg, color = \"#E64D4E\", linewidth = 3)\n", "\n", "plt.xlim([10, 35])\n", "plt.ylim([-120, 30])\n", "\n", "plt.xlabel(r'$z$', fontsize = 30)\n", "plt.ylabel(r'$T_{21}$ [mK]', fontsize = 30)\n", "plt.xticks(fontsize=30)\n", "plt.yticks(fontsize=30)\n", "plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(axis=\"y\", labelsize=30, pad = 10)\n", "plt.tick_params(axis=\"x\", labelsize=30, pad = 10)\n", "\n", "plt.tight_layout()\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "3a57e7ae", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize = (12, 6.75))\n", "\n", "plt.plot(zlist,CoeffStructure.T_CMB , color = \"#EE9063\", linewidth = 3, linestyle = 'dotted', label = r'$T_\\mathrm{CMB}$')\n", "\n", "\n", "plt.plot(zlist,CoeffStructure.Tk_ad , color = \"#0B92B1\", linewidth = 3, label = r'$T_\\mathrm{ad}$')\n", "plt.plot(zlist,CoeffStructure.Tk_xray , color = \"#E64D4E\", linewidth = 3, label = r'$T_x$')\n", "plt.plot(zlist,CoeffStructure.Tk_avg, color = \"#665191\", linewidth = 3, label = r'$T_k$')\n", "\n", "plt.plot(zlist,CoeffStructure._invTs_avg**-1 , color = \"#77AC54\", linewidth = 3, label = r'$T_s$')\n", "\n", "\n", "plt.xlim([10, 35])\n", "plt.ylim([0, 100])\n", "\n", "plt.xlabel(r'$z$', fontsize = 30)\n", "plt.ylabel(r'$T$ [K]', fontsize = 30)\n", "plt.xticks(fontsize=30)\n", "plt.yticks(fontsize=30)\n", "plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(axis=\"y\", labelsize=30, pad = 10)\n", "plt.tick_params(axis=\"x\", labelsize=30, pad = 10)\n", "plt.legend(fontsize=20, frameon = False)\n", "plt.tight_layout()\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "6815b8c5", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize = (12, 6.75))\n", "\n", "plt.semilogy(CoeffStructure.zintegral, CoeffStructure.SFRD_II_avg, color=\"#E64D4E\", linewidth=3.0, label = 'Pop II')\n", "plt.semilogy(CoeffStructure.zintegral, CoeffStructure.SFRD_III_avg, color=\"#0B92B1\", linewidth=3.0, label = 'Pop III')\n", "\n", "plt.semilogy(CoeffStructure.zintegral, CoeffStructure.SFRD_II_avg + CoeffStructure.SFRD_III_avg, color=\"#665191\", linewidth=3.0, label = 'Pop II + III')\n", "\n", "plt.xlim([10, 35])\n", "plt.ylim(1e-7, 1e-1)\n", "\n", "plt.xlabel(r'$z$', fontsize = 30)\n", "plt.ylabel(r'$\\dot{\\overline{\\rho}}_*(z)$', fontsize = 30)\n", "\n", "plt.xticks(fontsize=30)\n", "plt.yticks(fontsize=30)\n", "plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(axis=\"y\", labelsize=30, pad = 10)\n", "plt.tick_params(axis=\"x\", labelsize=30, pad = 10)\n", "plt.legend(fontsize=20, frameon = False)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "b1a55a1e", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize = (12, 6.75))\n", "\n", "plt.semilogy(CoeffStructure.zintegral, CoeffStructure.J_21_LW_II, color=\"#E64D4E\", linewidth=3.0, label = 'Pop II')\n", "plt.semilogy(CoeffStructure.zintegral, CoeffStructure.J_21_LW_III, color=\"#0B92B1\", linewidth=3.0, label = 'Pop III')\n", "\n", "plt.semilogy(CoeffStructure.zintegral, CoeffStructure.J_21_LW_II + CoeffStructure.J_21_LW_III, color=\"#665191\", linewidth=3.0, label = 'Pop II + III')\n", "\n", "plt.xlim([10, 35])\n", "plt.ylim(1e-3, 1e2)\n", "\n", "plt.xlabel(r'$z$', fontsize = 30)\n", "plt.ylabel(r'$J_{21}$', fontsize = 30)\n", "\n", "plt.xticks(fontsize=30)\n", "plt.yticks(fontsize=30)\n", "plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(axis=\"y\", labelsize=30, pad = 10)\n", "plt.tick_params(axis=\"x\", labelsize=30, pad = 10)\n", "plt.legend(fontsize=20, frameon = False)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "80eb8f5b", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize = (12, 6.75))\n", "\n", "plt.semilogy(CoeffStructure.zintegral, CoeffStructure.xa_avg, color=\"#665191\", linewidth=3.0, label = 'Pop II + III')\n", "\n", "plt.xlim([10, 35])\n", "plt.ylim(1e-3, 1e2)\n", "\n", "plt.xlabel(r'$z$', fontsize = 30)\n", "plt.ylabel(r'$x_\\alpha$', fontsize = 30)\n", "\n", "plt.xticks(fontsize=30)\n", "plt.yticks(fontsize=30)\n", "plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(axis=\"y\", labelsize=30, pad = 10)\n", "plt.tick_params(axis=\"x\", labelsize=30, pad = 10)\n", "plt.legend(fontsize=20, frameon = False)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "a85d3ff1", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize = (12, 6.75))\n", "\n", "plt.semilogy(CoeffStructure.zintegral,CoeffStructure.xe_avg, color=\"#E64D4E\", linewidth=3.0)\n", "\n", "plt.xlim([10, 35])\n", "plt.ylim([1e-4, 1e-3])\n", "\n", "plt.xlabel(r'$z$', fontsize = 30)\n", "plt.ylabel(r'$x_e$', fontsize = 30)\n", "plt.xticks(fontsize=30)\n", "plt.yticks(fontsize=30)\n", "plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(axis=\"y\", labelsize=30, pad = 10)\n", "plt.tick_params(axis=\"x\", labelsize=30, pad = 10)\n", "# plt.legend(fontsize=20, frameon = False)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "77eb7769", "metadata": {}, "source": [ "## Step 4: Fluctuations" ] }, { "cell_type": "markdown", "id": "36a3cd23", "metadata": {}, "source": [ "The last structure to run is Power_Spectra, which uses our lognormal and log-chisquare prescription to compute the statistics of IGM fluctuations in the form of two-point functions.\n", "\n", "Note that the new version of Zeus includes a modified prescription for computing power spectra. The form for the 21-cm brightness temperature used in the old version of Zeus is\n", "\n", "\\begin{equation}\n", "T_{21}=T_0(z)\\left(1+\\delta-\\delta_v\\right) x_{\\mathrm{HI}}\\left(\\frac{x_\\alpha}{1+x_\\alpha}\\right)\\left(1-\\frac{T_{\\mathrm{CMB}}}{T_c}\\right)\n", "\\end{equation}\n", "\n", "where the leftmost $\\delta$ term is the total matter overdensity. Because $T_{21}$ tracks neutral hydrogen fluctuations, we make the substitution $\\delta \\rightarrow \\delta_b$. Thus, the large-scale structure (LSS) terms that fold in to computing $\\Delta^2_{21}(k,z)$ now depend on the baryon power spectrum $P_\\mathrm{b}(k,z)$ and the baryon-dark matter cross power spectrum $P_\\mathrm{b,cdm}(k,z)$.\n", "\n", "This new feature can be toggled with UserParams.USE_BARYON_FLAG = 1 for the new $\\delta_b$ prescription or UserParams.USE_BARYON_FLAG = 0 for the old $\\delta$ form." ] }, { "cell_type": "code", "execution_count": null, "id": "f80f6e74", "metadata": {}, "outputs": [], "source": [ "RSDMODE = 1\n", "\n", "UserParams.USE_BARYON_FLAG = 1\n", "PS21 = zeus21.Power_Spectra(UserParams, CosmoParams, AstroParams, CoeffStructure, RSD_MODE = RSDMODE)\n", "klist = PS21.klist_PS" ] }, { "cell_type": "markdown", "id": "ff154eb4", "metadata": {}, "source": [ "## Plotting Results" ] }, { "cell_type": "code", "execution_count": null, "id": "62acc2ca", "metadata": {}, "outputs": [], "source": [ "from scipy.interpolate import RegularGridInterpolator\n", "\n", "interp = RegularGridInterpolator((zlist, klist), PS21.Deltasq_T21, method = 'cubic', bounds_error=False, fill_value=0)\n", "interpLin = RegularGridInterpolator((zlist, klist), PS21.Deltasq_T21_lin, method = 'cubic', bounds_error=False, fill_value=0)" ] }, { "cell_type": "code", "execution_count": null, "id": "0a03d011", "metadata": {}, "outputs": [], "source": [ "#choose a z to plot\n", "kchoose=0.3 \n", "zlistHighRes = np.geomspace(zlist[0], zlist[-1], 1000)\n", "\n", "powerSpectrum = interp((zlistHighRes, kchoose))\n", "powerSpectrumLin = interpLin((zlistHighRes, kchoose))\n", "\n", "plt.figure(figsize = (12, 6.75))\n", "\n", "plt.semilogy(zlistHighRes, powerSpectrum, color=\"#E64D4E\", linewidth=3.0, label = 'Nonlinear')\n", "plt.semilogy(zlistHighRes, powerSpectrumLin, color=\"#0B92B1\", linewidth=3.0, linestyle = 'dashed', label = 'Linear')\n", "\n", "plt.xlim([10, 25])\n", "plt.ylim([1,200])\n", "\n", "plt.xlabel(r'$z$', fontsize = 30)\n", "plt.ylabel(r'$\\Delta^2_{21}\\,\\rm[mK^2]$', fontsize = 30)\n", "\n", "plt.xticks(fontsize=30)\n", "plt.yticks(fontsize=30)\n", "plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(axis=\"y\", labelsize=30, pad = 10)\n", "plt.tick_params(axis=\"x\", labelsize=30, pad = 10)\n", "plt.legend(fontsize=20, frameon = False)\n", "plt.tight_layout()\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "730f80cd", "metadata": {}, "outputs": [], "source": [ "#choose a z to plot\n", "zchoose=16\n", "klistHighRes = np.geomspace(1e-4, 1e1, 1000)\n", "\n", "powerSpectrum = interp((zchoose, klistHighRes))\n", "powerSpectrumLin = interpLin((zchoose, klistHighRes))\n", "\n", "plt.figure(figsize = (12, 6.75))\n", "\n", "plt.loglog(klistHighRes, powerSpectrum, color=\"#E64D4E\", linewidth=3.0, label = 'Nonlinear')\n", "plt.loglog(klistHighRes, powerSpectrumLin, color=\"#0B92B1\", linewidth=3.0, linestyle = 'dashed', label = 'Linear')\n", "\n", "plt.xlim([1e-3,1e0])\n", "plt.ylim([1e-1,1e3])\n", "\n", "plt.xlabel(r'$k\\,\\rm [Mpc^{-1}]$', fontsize = 30)\n", "plt.ylabel(r'$\\Delta^2_{21}\\,\\rm[mK^2]$', fontsize = 30)\n", "plt.xticks(fontsize=30)\n", "plt.yticks(fontsize=30)\n", "plt.tick_params(which='major', length=12, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(which='minor', length=5, width=2, direction='in', top = True, bottom = True, left = True, right = True)\n", "plt.tick_params(axis=\"y\", labelsize=30, pad = 10)\n", "plt.tick_params(axis=\"x\", labelsize=30, pad = 10)\n", "plt.legend(fontsize=20, frameon = False)\n", "plt.tight_layout()\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "4c7cc3ff", "metadata": {}, "source": [ "For questions, don't hesitate to reach out to hcruz2@jhu.edu!" ] }, { "cell_type": "code", "execution_count": null, "id": "7c58d071", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "21zeus_hack", "language": "python", "name": "21zeus_hack" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.13" } }, "nbformat": 4, "nbformat_minor": 5 }