Skip to content

Bug fixes to mini-halos.ipynb and relative-velocities.ipynb #525

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
raxquax123 wants to merge 2 commits into from
+540 −390
Closed

Bug fixes to mini-halos.ipynb and relative-velocities.ipynb #525

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
dbf629d
bug fixes to relative-velocities.ipynb and mini-halos.iypnb
raxquax123 Jun 12, 2025
0f12d9f
Merge remote-tracking branch 'origin/main'
raxquax123 Jun 12, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
138 changes: 77 additions & 61 deletions docs/tutorials/mini-halos.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,10 @@
" (1, \"cyan\"),\n",
" ],\n",
")\n",
"plt.register_cmap(cmap=EoR_colour)"
"try:\n",
" plt.colormaps.register(cmap=EoR_colour)\n",
"except ValueError:\n",
" print(\"EoR colormap already registered. Continuing without re-registering.\")"
]
},
{
Expand Down Expand Up @@ -70,13 +73,13 @@
"cache = p21c.OutputCache(output_dir)\n",
"\n",
"# USE_FFTW_WISDOM make FFT faster\n",
"inputs = p21c.InputParameters.from_template('Qin20').evolve_input_structs(\n",
"inputs = p21c.InputParameters.from_template('Qin20', random_seed=random_seed,).evolve_input_structs(\n",
" HII_DIM = 128,\n",
" BOX_LEN = 250,\n",
")\n",
"\n",
"initial_conditions = p21c.initial_conditions(\n",
" inputs=inputs, random_seed=random_seed, direc=output_dir\n",
" inputs=inputs, cache=cache,\n",
")"
]
},
Expand Down Expand Up @@ -149,7 +152,7 @@
" resolution=inputs.simulation_options.cell_size,\n",
")\n",
"\n",
"lightcone_fid = p21c.run_lightcone(\n",
"_, _, _, lightcone_fid = p21c.run_lightcone(\n",
" lightconer=lcn,\n",
" inputs=inputs,\n",
" lightcone_quantities=lightcone_quantities,\n",
Expand All @@ -160,15 +163,15 @@
" len(lightcone_quantities),\n",
" 1,\n",
" figsize=(\n",
" getattr(lightcone_fid, lightcone_quantities[0]).shape[2] * 0.01,\n",
" getattr(lightcone_fid, lightcone_quantities[0]).shape[1]\n",
" lightcone_fid.lightcones[lightcone_quantities[0]].shape[2] * 0.01,\n",
" lightcone_fid.lightcones[lightcone_quantities[0]].shape[1]\n",
" * 0.01\n",
" * len(lightcone_quantities),\n",
" ),\n",
")\n",
"for ii, lightcone_quantity in enumerate(lightcone_quantities):\n",
" axs[ii].imshow(\n",
" getattr(lightcone_fid, lightcone_quantity)[1],\n",
" lightcone_fid.lightcones[lightcone_quantity][1],\n",
" vmin=vmins[ii],\n",
" vmax=vmaxs[ii],\n",
" cmap=cmaps[ii],\n",
Expand Down Expand Up @@ -201,7 +204,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"let's vary paremeters that describe mini-halos and see the impact to the global signal"
"Let's vary parameters that describe mini-halos and see the impact to the global signal"
]
},
{
Expand Down Expand Up @@ -235,14 +238,14 @@
"lws = [1, 3, 2, 2, 2, 2]\n",
"\n",
"texts = [\n",
" \"varying \" + r\"$f_{*,7}^{\\rm mol}$\",\n",
" \"varying \" + r\"$f_{\\rm esc}^{\\rm mol}$\",\n",
" \"varying \" + r\"$L_{\\rm x}^{\\rm mol}$\",\n",
" \"varying \" + r\"$1-f_{\\rm H_2}^{\\rm shield}$\",\n",
" \"varying \" + \"$f_{*,7}$\",\n",
" \"varying \" + \"$f_{esc}$\",\n",
" \"varying \" + \"$L_x$\",\n",
" \"varying \" + \"$f_{H_2}$\",\n",
"]\n",
"var_params = [\"F_STAR7_MINI\",\"F_ESC7_MINI\",\"L_X_MINI\",\"F_H2_SHIELD\"]\n",
"factors = [0, 1, 0.1, 0.5, 2, 10]\n",
"labels = [\"NOmini\", \"reference\", \"x0.1\", \"x0.5\", \"x2\", \"x10\"]"
"factors = [0.0000001, 1, 0.1, 0.5, 2, 10]\n",
"labelss = [\"NOmini\", \"reference\", \"x0.1\", \"x0.5\", \"x2\", \"x10\"]"
]
},
{
Expand Down Expand Up @@ -289,7 +292,7 @@
" \"cumulative_recombinations\",\n",
" \"z_reion\",\n",
" \"ionisation_rate_G12\",\n",
" \"J_21_LW\",\n",
" #\"J_21_LW\",\n",
" \"density\",\n",
")\n",
"# choose some to plot...\n",
Expand All @@ -299,62 +302,71 @@
" \"neutral_fraction\",\n",
" \"cumulative_recombinations\",\n",
" \"ionisation_rate_G12\",\n",
" \"J_21_LW\",\n",
" #\"J_21_LW\",\n",
")\n",
"ymins = [-120, 1e1, 0, 0, 0, 0]\n",
"ymaxs = [30, 1e3, 1, 1, 1, 10]\n",
"\n",
"lcn = p21c.RectilinearLightconer.with_equal_cdist_slices(\n",
" min_redshift=min(inputs.node_redshifts),\n",
" max_redshift=max(inputs.node_redshifts),\n",
" quantities=global_quantities,\n",
" resolution=inputs.simulation_options.cell_size,\n",
")\n",
"\n",
"fig, axss = plt.subplots(\n",
" len(plot_quantities),\n",
" len(labelss),\n",
" len(var_params),\n",
" sharex=True,\n",
" figsize=(4 * len(labelss), 2 * len(global_quantities)),\n",
" figsize=(4 * len(var_params), 2 * len(global_quantities)),\n",
")\n",
"\n",
"for pp, text in enumerate(texts):\n",
" axs = axss[:, pp]\n",
" for kk, label in enumerate(labels):\n",
" for kk, label in enumerate(labelss):\n",
" factor = factors[kk]\n",
"\n",
" if var_params[pp] == \"F_H2_SHIELD\":\n",
" if factor > 1:\n",
" continue\n",
" factor = min(1, factor)\n",
" inputs_var = inputs.evolve_input_structs(\n",
" USE_MINI_HALOS=(label!=\"NOmini\"),\n",
" F_H2_SHIELD=1. - (1. - inputs.astro_params.F_H2_SHIELD)*factor,\n",
" )\n",
" else:\n",
" if pp == 2 or pp == 3:\n",
" factor = min(1, factor)\n",
" inputs_var = inputs.evolve_input_structs(\n",
" USE_MINI_HALOS=(label!=\"NOmini\"),\n",
" **{var_params[pp] : inputs.astro_params.get(var_params[pp]) * factor}\n",
" **{var_params[pp] : getattr(inputs.astro_params, var_params[pp]) * factor}\n",
" )\n",
" if label == \"reference\":\n",
" lightcone = lightcone_fid\n",
" else:\n",
" lightcone = p21c.run_lightcone(\n",
" initial_conditions=initial_conditions,\n",
" _, _, _, lightcone = p21c.run_lightcone(\n",
" lightconer=lcn,\n",
" inputs=inputs_var,\n",
" global_quantities=global_quantities,\n",
" direc=cache,\n",
" initial_conditions = initial_conditions,\n",
" cache=cache,\n",
" )\n",
"\n",
" freqs = 1420.4 / (np.array(lightcone.node_redshifts) + 1.0)\n",
" freqs = 1420.4 / (np.array(lightcone.inputs.node_redshifts) + 1.0)\n",
" for jj, global_quantity in enumerate(plot_quantities):\n",
" axs[jj].plot(\n",
" freqs,\n",
" getattr(\n",
" lightcone, \"global_{}\".format(global_quantity.replace(\"_box\", \"\"))\n",
" ),\n",
" lightcone.global_quantities[global_quantity],\n",
" color=colors[kk],\n",
" linestyle=linestyles[kk],\n",
" label=labels[kk],\n",
" label=labelss[kk],\n",
" lw=lws[kk],\n",
" )\n",
"\n",
" axs[0].text(\n",
" 0.01,\n",
" 0.99,\n",
" texts,\n",
" text,\n",
" horizontalalignment=\"left\",\n",
" verticalalignment=\"top\",\n",
" transform=axs[0].transAxes,\n",
Expand Down Expand Up @@ -391,7 +403,7 @@
" )\n",
" axs[jj].yaxis.set_tick_params(labelsize=0)\n",
"\n",
"plt.tight_layout()\n",
"#plt.tight_layout()\n",
"fig.subplots_adjust(hspace=0.0, wspace=0.0)"
]
},
Expand All @@ -416,6 +428,9 @@
"# define functions to calculate PS, following py21cmmc\n",
"from powerbox.tools import get_power\n",
"\n",
"BOX_LEN = 250\n",
"HII_DIM = 64\n",
"\n",
"\n",
"def compute_power(\n",
" box,\n",
Expand Down Expand Up @@ -460,33 +475,37 @@
"\n",
"\n",
"def powerspectra(\n",
" brightness_temp, n_psbins=50, nchunks=10, min_k=0.1, max_k=1.0, logk=True\n",
"):\n",
" lightcone, n_psbins=50, nchunks=10, logk=True\n",
"): \n",
" brightness_temp = lightcone.lightcones[\"brightness_temp\"]\n",
" data = []\n",
" chunk_indices = list(\n",
" range(\n",
" 0,\n",
" brightness_temp.n_slices,\n",
" round(brightness_temp.n_slices / nchunks),\n",
" np.shape(brightness_temp)[2],\n",
" round(np.shape(brightness_temp)[2] / nchunks),\n",
" )\n",
" )\n",
"\n",
" if len(chunk_indices) > nchunks:\n",
" chunk_indices = chunk_indices[:-1]\n",
" chunk_indices.append(brightness_temp.n_slices)\n",
" chunk_indices.append(np.shape(brightness_temp)[2])\n",
"\n",
" for i in range(nchunks):\n",
" start = chunk_indices[i]\n",
" end = chunk_indices[i + 1]\n",
" chunklen = (end - start) * brightness_temp.cell_size\n",
" chunklen = (end - start) * lightcone.inputs.simulation_options.cell_size.value\n",
"\n",
" power, k = compute_power(\n",
" brightness_temp.brightness_temp[:, :, start:end],\n",
" brightness_temp[:, :, start:end],\n",
" (BOX_LEN, BOX_LEN, chunklen),\n",
" n_psbins,\n",
" log_bins=logk,\n",
" )\n",
" k = k[~np.isnan(power)]\n",
" power = power[~np.isnan(power)]\n",
" data.append({\"k\": k, \"delta\": power * k**3 / (2 * np.pi**2)})\n",
"\n",
" return data"
]
},
Expand Down Expand Up @@ -514,54 +533,58 @@
],
"source": [
"# do 5 chunks but only plot 1 - 4, the 0th has no power for minihalo models where xH=0\n",
"nchunks = 4\n",
"nchunks = 5\n",
"\n",
"fig, axss = plt.subplots(\n",
" nchunks,\n",
" len(labelss),\n",
" (nchunks-1),\n",
" len(var_params),\n",
" sharex=True,\n",
" sharey=True,\n",
" figsize=(4 * len(labelss), 3 * (nchunks)),\n",
" figsize=(4 * len(var_params), 3 * (nchunks-1)),\n",
" subplot_kw={\"xscale\": \"log\", \"yscale\": \"log\"},\n",
")\n",
"\n",
"for pp, text in enumerate(texts):\n",
" labels = labels[pp]\n",
" factors = factors[pp]\n",
" labels = labelss[pp]\n",
" #factor = factors[pp]\n",
" axs = axss[:, pp]\n",
" for kk, label in enumerate(labels):\n",
" for kk, label in enumerate(labelss):\n",
" factor = factors[kk]\n",
"\n",
" if var_params[pp] == \"F_H2_SHIELD\":\n",
" if factor > 1:\n",
" continue\n",
" factor = min(1, factor)\n",
" inputs_var = inputs.evolve_input_structs(\n",
" USE_MINI_HALOS=(label!=\"NOmini\"),\n",
" F_H2_SHIELD=1. - (1. - inputs.astro_params.F_H2_SHIELD)*factor,\n",
" )\n",
" else:\n",
" if pp == 2 or pp == 3:\n",
" factor = min(1, factor)\n",
" inputs_var = inputs.evolve_input_structs(\n",
" USE_MINI_HALOS=(label!=\"NOmini\"),\n",
" **{var_params[pp] : inputs.astro_params.get(var_params[pp]) * factor}\n",
" **{var_params[pp] : getattr(inputs.astro_params, var_params[pp]) * factor}\n",
" )\n",
" if label == \"reference\":\n",
" lightcone = lightcone_fid\n",
" else:\n",
" lightcone = p21c.run_lightcone(\n",
" initial_conditions=initial_conditions,\n",
" _, _, _, lightcone = p21c.run_lightcone(\n",
" lightconer=lcn,\n",
" inputs=inputs_var,\n",
" global_quantities=global_quantities,\n",
" direc=cache,\n",
" initial_conditions = initial_conditions,\n",
" cache=cache,\n",
" )\n",
"\n",
" PS = powerspectra(lightcone)\n",
" for ii in range(nchunks):\n",
" PS = powerspectra(lightcone, nchunks=nchunks)\n",
" for ii in range(nchunks - 1):\n",
" axs[ii].plot(\n",
" PS[ii + 1][\"k\"],\n",
" PS[ii + 1][\"delta\"],\n",
" color=colors[kk],\n",
" linestyle=linestyles[kk],\n",
" label=labels[kk],\n",
" label=labelss[kk],\n",
" lw=lws[kk],\n",
" )\n",
"\n",
Expand All @@ -579,7 +602,7 @@
" axs[0].text(\n",
" 0.01,\n",
" 0.99,\n",
" texts,\n",
" text,\n",
" horizontalalignment=\"left\",\n",
" verticalalignment=\"top\",\n",
" transform=axs[0].transAxes,\n",
Expand All @@ -590,7 +613,7 @@
" axs[-1].xaxis.set_tick_params(labelsize=15)\n",
"\n",
" if pp == 0:\n",
" for ii in range(nchunks):\n",
" for ii in range(nchunks - 1):\n",
" axs[ii].set_ylim(2e-1, 2e2)\n",
" axs[ii].set_ylabel(\"$k^3 P$\", fontsize=15)\n",
" axs[ii].yaxis.set_tick_params(labelsize=15)\n",
Expand All @@ -609,15 +632,8 @@
"cell_type": "markdown",
"metadata": {},
"source": [
" Now you know how minihalo can shape the 21-cm signal!"
" Now you know how minihalos can shape the 21-cm signal!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand Down
Loading
Loading