[1]:
import numpy as np
import matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
import py21cmfast as p21c
from py21cmfast import global_params
from py21cmfast import plotting
random_seed = 1605
EoR_colour = matplotlib.colors.LinearSegmentedColormap.from_list('mycmap',\
[(0, 'white'),(0.33, 'yellow'),(0.5, 'orange'),(0.68, 'red'),\
(0.83333, 'black'),(0.9, 'blue'),(1, 'cyan')])
plt.register_cmap(cmap=EoR_colour)
This result was obtained using 21cmFAST at commit 2bb4807c7ef1a41649188a3efc462084f2e3b2e0
This notebook shows how to include the effect of the DM-baryon relative velocities, and the new EOS2021 parameters.
Based on Muñoz+21 (https://arxiv.org/abs/2110.13919). See https://drive.google.com/drive/folders/1-50AO-i3arCnfHc22YWXJacs4u-xsPL6?usp=sharing for the large (1.5Gpc) AllGalaxies simulation with the same parameters.
It is recommended to do the other tutorials first
Fiducial and lightcones¶
Let’s fix the initial condition for this tutorial.
[2]:
output_dir = '/Users/julian/Dropbox/Research/EOS_21/'
HII_DIM = 64
BOX_LEN = 200 #cell size of ~3 Mpc or below for relative velocities
# USE_FFTW_WISDOM make FFT faster AND use relative velocities. , 'USE_INTERPOLATION_TABLES': True or code is too slow
user_params = {"HII_DIM":HII_DIM, "BOX_LEN": BOX_LEN, "USE_FFTW_WISDOM": True, 'USE_INTERPOLATION_TABLES': True,
"FAST_FCOLL_TABLES": True,
"USE_RELATIVE_VELOCITIES": True, "POWER_SPECTRUM": 5}
#set FAST_FCOLL_TABLES to TRUE if using minihaloes, it speeds up the generation of tables by ~x30 (see Appendix of 2110.13919)
#USE_RELATIVE_VELOCITIES is important for minihaloes. If True, POWER_SPECTRUM has to be set to 5 (CLASS) to get the transfers.
initial_conditions = p21c.initial_conditions(user_params=user_params,
random_seed=random_seed,
direc=output_dir
#, regenerate=True
)
[3]:
plotting.coeval_sliceplot(initial_conditions, "lowres_vcb");
plotting.coeval_sliceplot(initial_conditions, "lowres_density");
Let’s run a ‘fiducial’ model and see its lightcones
Note that the reference model has
F_STAR7_MINI ~ F_STAR10
and
F_ESC7_MINI ~ 1%, as low, but conservative fiducial
Also we take L_X_MINI=L_X out of simplicity (and ignorance)
[4]:
# the lightcones we want to plot later together with their color maps and min/max
lightcone_quantities = ('brightness_temp','Ts_box','xH_box',"dNrec_box",'z_re_box','Gamma12_box','J_21_LW_box',"density")
cmaps = [EoR_colour,'Reds','magma','magma','magma','cubehelix','cubehelix','viridis']
vmins = [-150, 1e1, 0, 0, 5, 0, 0, -1]
vmaxs = [ 30, 1e3, 1, 2, 9, 1,10, 1]
astro_params_vcb = {"ALPHA_ESC": -0.3, "F_ESC10": -1.35,
"ALPHA_STAR": 0.5, "F_STAR10": -1.25, "t_STAR" :0.5,
"F_STAR7_MINI": -2.5, "ALPHA_STAR_MINI": 0, "F_ESC7_MINI" : -1.35,
"L_X": 40.5, "L_X_MINI": 40.5, "NU_X_THRESH": 500.0,
"A_VCB": 1.0, "A_LW": 2.0}
astro_params_novcb=astro_params_vcb
astro_params_novcb.update({'A_VCB': 0.0})
#setting 'A_VCB': 0 sets to zero the effect of relative velocities (fiducial value is 1.0)
#the parameter 'A_LW' (with fid value of 2.0) does the same for LW feedback.
flag_options_fid = {"INHOMO_RECO":True, 'USE_MASS_DEPENDENT_ZETA':True, 'USE_TS_FLUCT':True,
'USE_MINI_HALOS':True, 'FIX_VCB_AVG':False}
flag_options_fid_vavg = flag_options_fid
flag_options_fid_vavg.update({'FIX_VCB_AVG': True})
#the flag FIX_VCB_AVG side-steps the relative-velocity ICs, and instead fixes all velocities to some average value.
#It gets the background right but it's missing VAOs and 21cm power at large scales
[ ]:
ZMIN=5.
lightcone_fid_vcb = p21c.run_lightcone(
redshift = ZMIN,
init_box = initial_conditions,
flag_options = flag_options_fid,
astro_params = astro_params_vcb,
lightcone_quantities=lightcone_quantities,
global_quantities=lightcone_quantities,
random_seed = random_seed,
direc = output_dir,
write=True#, regenerate=True
)
fig, axs = plt.subplots(len(lightcone_quantities),1,
figsize=(20,10))#(getattr(lightcone_fid_vcb, lightcone_quantities[0]).shape[2]*0.01,
#getattr(lightcone_fid_vcb, lightcone_quantities[0]).shape[1]*0.01*len(lightcone_quantities)))
for ii, lightcone_quantity in enumerate(lightcone_quantities):
axs[ii].imshow(getattr(lightcone_fid_vcb, lightcone_quantity)[1],
vmin=vmins[ii], vmax=vmaxs[ii],cmap=cmaps[ii])
axs[ii].text(1, 0.05, lightcone_quantity,horizontalalignment='right',verticalalignment='bottom',
transform=axs[ii].transAxes,color = 'red',backgroundcolor='white',fontsize = 15)
axs[ii].xaxis.set_tick_params(labelsize=10)
axs[ii].yaxis.set_tick_params(labelsize=0)
plt.tight_layout()
fig.subplots_adjust(hspace = 0.01)
[ ]:
#also run one without velocities and with fixed vcb=vavg (for comparison)
lightcone_fid_novcb = p21c.run_lightcone(
redshift = ZMIN,
init_box = initial_conditions,
flag_options = flag_options_fid,
astro_params = astro_params_novcb,
lightcone_quantities=lightcone_quantities,
global_quantities=lightcone_quantities,
random_seed = random_seed,
direc = output_dir,
write=True#, regenerate=True
)
lightcone_fid_vcbavg = p21c.run_lightcone(
redshift = ZMIN,
init_box = initial_conditions,
flag_options = flag_options_fid_vavg,
astro_params = astro_params_vcb,
lightcone_quantities=lightcone_quantities,
global_quantities=lightcone_quantities,
random_seed = random_seed,
direc = output_dir,
write=True#, regenerate=True
)
[ ]:
[ ]:
#plus run one with only atomic-cooling galaxies but same otherwise
flag_options_NOMINI=flag_options_fid
flag_options_NOMINI.update({'USE_MINI_HALOS': False})
lightcone_fid_NOMINI = p21c.run_lightcone(
redshift = ZMIN,
init_box = initial_conditions,
flag_options = flag_options_NOMINI,
astro_params = astro_params_vcb,
lightcone_quantities=lightcone_quantities,
global_quantities=lightcone_quantities,
random_seed = random_seed,
direc = output_dir,
write=True#, regenerate=True
)
[ ]:
#compare vcb and novcb
fig, axs = plt.subplots(2 ,1,
figsize=(20,6))
axs[0].imshow(getattr(lightcone_fid_vcb, 'brightness_temp')[1],
vmin=vmins[0], vmax=vmaxs[0],cmap=cmaps[0])
axs[1].imshow(getattr(lightcone_fid_novcb, 'brightness_temp')[1],
vmin=vmins[0], vmax=vmaxs[0],cmap=cmaps[0])
axs[0].text(1, 0.05, 'vcb' ,horizontalalignment='right',verticalalignment='bottom',
transform=axs[0].transAxes,color = 'red',backgroundcolor='white',fontsize = 15)
axs[1].text(1, 0.05, 'novcb' ,horizontalalignment='right',verticalalignment='bottom',
transform=axs[1].transAxes,color = 'red',backgroundcolor='white',fontsize = 15)
# axs[0].xaxis.set_tick_params(labelsize=10)
# axs[1].yaxis.set_tick_params(labelsize=0)
plt.tight_layout()
fig.subplots_adjust(hspace = 0.01)
[ ]:
#plot tau
tau_vcb=tau_novcb=tau_NOMINI=np.array([])
for il,lightcone in enumerate([lightcone_fid_vcb,lightcone_fid_novcb,lightcone_fid_NOMINI]):
z_e=np.array([]);
tau_e=np.array([]);
for i in range(len(lightcone.node_redshifts)-1):
tauz=p21c.compute_tau(redshifts=lightcone.node_redshifts[-1:-2-i:-1],
global_xHI=lightcone.global_xHI[-1:-2-i:-1])
tau_e=np.append(tau_e,tauz)
z_e=np.append(z_e,lightcone.node_redshifts[-2-i])
#add lower zs where we manually set xH=1
zlow=np.linspace(lightcone.node_redshifts[-1]-0.1, 0.1, 14)
for zl in zlow:
tauz=p21c.compute_tau(redshifts=np.array([zl]), global_xHI=np.array([lightcone.global_xHI[-1]]))
tau_e=np.append([tauz],tau_e)
z_e=np.append([zl],z_e)
if(il==0):
tau_vcb=tau_e
elif (il==1):
tau_novcb=tau_e
else:
tau_NOMINI=tau_e
linestyles = ['-', '-.',':']
colors = ['black','gray','#377eb8']
lws = [3,1,2]
fig, axs = plt.subplots(1, 1, sharex=True, figsize=(8,4))
kk=0
axs.plot(z_e, tau_vcb, label = 'vcb',
color=colors[kk],linestyle=linestyles[kk], lw=lws[kk])
kk=1
axs.plot(z_e, tau_novcb, label = 'no vcb',
color=colors[kk],linestyle=linestyles[kk], lw=lws[kk])
kk=2
axs.plot(z_e, tau_NOMINI, label = 'no MINI',
color=colors[kk],linestyle=linestyles[kk],lw=lws[kk])
axs.set_ylim(0., 0.1)
axs.set_xlabel('redshift',fontsize=15)
axs.xaxis.set_tick_params(labelsize=15)
axs.set_xlim(0.,20.)
axs.set_ylabel('$\\tau$',fontsize=15)
axs.yaxis.set_tick_params(labelsize=15)
plt.tight_layout()
fig.subplots_adjust(hspace = 0.0,wspace=0.0)
tauPmin=0.0561-0.0071
tauPmax=0.0561+0.0071
axs.axhspan(tauPmin, tauPmax, alpha=0.34, color='black')
axs.grid()
#Planck2020: tau=0.0561±0.0071
[ ]:
#check that the tau z=15-30 is below 0.02 as Planck requires
print(z_e[-1],z_e[55])
tau_vcb[-1]-tau_vcb[55]
[ ]:
linestyles = ['-', '-.',':']
colors = ['black','gray','#377eb8']
lws = [3,1,2]
labels = ['vcb', 'no vcb', 'no MINI']
fig, axs = plt.subplots(1, 1, sharex=True, figsize=(8,4))
for kk,lightcone in enumerate([lightcone_fid_vcb,lightcone_fid_novcb,lightcone_fid_NOMINI]):
axs.plot(lightcone.node_redshifts, lightcone.global_xHI, label = labels[kk],
color=colors[kk],linestyle=linestyles[kk], lw=lws[kk])
axs.set_ylim(0., 1.)
axs.set_xlabel('redshift',fontsize=15)
axs.xaxis.set_tick_params(labelsize=15)
axs.set_xlim(5.,20.)
axs.set_ylabel('$x_{HI}$',fontsize=15)
axs.yaxis.set_tick_params(labelsize=15)
plt.tight_layout()
fig.subplots_adjust(hspace = 0.0,wspace=0.0)
axs.grid()
[ ]:
[ ]:
#choose a redshift to print coeval slices and see if there are VAOs. Usually best then T21~T21min/2
zz=zlist21[40]
print(zz)
[ ]:
#We plot a coeval box, but we compare the vcb case against the vcb=vavg, since the no velocity (vcb=0) case has a background evolution that is too different.
coeval_fid_vcb = p21c.run_coeval(
redshift = zz,
init_box = initial_conditions,
flag_options = flag_options_fid,
astro_params = astro_params_vcb,
random_seed = random_seed,
direc = output_dir,
write=True#, regenerate=True
)
coeval_fid_vcbavg = p21c.run_coeval(
redshift = zz,
init_box = initial_conditions,
flag_options = flag_options_fid_vavg,
astro_params = astro_params_vcb,
random_seed = random_seed,
direc = output_dir,
write=True#, regenerate=True
)
[ ]:
[ ]:
T21slice_vcb=coeval_fid_vcb.brightness_temp
T21avg_vcb = np.mean(T21slice_vcb)
dT21slice_vcb = T21slice_vcb - T21avg_vcb
T21slice_novcb=coeval_fid_vcbavg.brightness_temp
T21avg_novcb = np.mean(T21slice_novcb)
dT21slice_novcb = T21slice_novcb - T21avg_novcb
sigma21=np.sqrt(np.var(dT21slice_vcb))
T21maxplot = 3.0*sigma21
T21minplot = -2.0 * sigma21
origin = 'lower'
extend = 'both'
origin = None
extend = 'neither'
xx = np.linspace(0, BOX_LEN, HII_DIM, endpoint=False)
yy = xx
indexv=0
fig, ax = plt.subplots(2, 2, constrained_layout=True, figsize=(10,8),
sharex='col', sharey='row',
gridspec_kw={'hspace': 0, 'wspace': 0})
cs0=ax[0,0].contourf(xx, yy, dT21slice_novcb[indexv], extend=extend, origin=origin,
vmin=T21minplot, vmax=T21maxplot,cmap='bwr')
fig.colorbar(cs0, ax=ax[0,0], shrink=0.9, location='left')
cs1=ax[0,1].contourf(xx, yy, dT21slice_vcb[indexv], extend=extend, origin=origin,
vmin=T21minplot, vmax=T21maxplot,cmap='bwr')
fig.colorbar(cs1, ax=ax[0,1], shrink=0.9)
deltaslice=initial_conditions.lowres_density
deltaavg = np.mean(deltaslice)
ddeltaslice = deltaslice - deltaavg
vcbslice=initial_conditions.lowres_vcb
vcbavg = np.mean(vcbslice)
dvcbslice = vcbslice
print(vcbavg)
csd=ax[1,0].contourf(xx, yy, ddeltaslice[indexv])
fig.colorbar(csd, ax=ax[1,0], shrink=0.9, location='left')
csv=ax[1,1].contourf(xx, yy, dvcbslice[indexv])
fig.colorbar(csv, ax=ax[1,1], shrink=0.9, extend=extend)
plt.show()
plt.tight_layout()
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
global_quantities = ('brightness_temp','Ts_box','xH_box',"dNrec_box",'z_re_box','Gamma12_box','J_21_LW_box',"density")
#choose some to plot...
plot_quantities = ('brightness_temp','Ts_box','xH_box',"dNrec_box",'Gamma12_box','J_21_LW_box')
ymins = [-120, 1e1, 0, 0, 0, 0]
ymaxs = [ 30, 1e3, 1, 1, 1,5]
linestyles = ['-', '-',':','-.','-.',':']
colors = ['gray','black','#e41a1c','#377eb8','#e41a1c','#377eb8']
lws = [2,2,2,2]
textss = ['vcb','MCGs']
factorss = [[0, 1],] * len(textss)
labelss = [['NO', 'reference'],] * len(textss)
fig, axss = plt.subplots(len(plot_quantities), len(labelss),
sharex=True, figsize=(4*len(labelss),2*len(plot_quantities)))
for pp, texts in enumerate(textss):
labels = labelss[pp]
factors = factorss[pp]
axs = axss[:,pp]
for kk, label in enumerate(labels):
factor = factors[kk]
if kk==0:
if pp == 0:
lightcone = lightcone_fid_NOMINI
else:
lightcone = lightcone_fid_novcb
else:
lightcone = lightcone_fid_vcb
freqs = 1420.4 / (np.array(lightcone.node_redshifts) + 1.)
for jj, global_quantity in enumerate(plot_quantities):
axs[jj].plot(freqs, getattr(lightcone, 'global_%s'%global_quantity.replace('_box','')),
color=colors[kk],linestyle=linestyles[kk], label = labels[kk],lw=lws[kk])
axs[0].text(0.01, 0.99, texts,horizontalalignment='right',verticalalignment='bottom',
transform=axs[0].transAxes,fontsize = 15)
for jj, global_quantity in enumerate(plot_quantities):
axs[jj].set_ylim(ymins[jj], ymaxs[jj])
axs[-1].set_xlabel('Frequency/MHz',fontsize=15)
axs[-1].xaxis.set_tick_params(labelsize=15)
axs[0].set_xlim(1420.4 / (35 + 1.), 1420.4 / (5.5 + 1.))
zlabels = np.array([ 6, 7, 8, 10, 13, 18, 25, 35])
ax2 = axs[0].twiny()
ax2.set_xlim(axs[0].get_xlim())
ax2.set_xticks(1420.4 / (zlabels + 1.))
ax2.set_xticklabels(zlabels.astype(np.str))
ax2.set_xlabel("redshift",fontsize=15)
ax2.xaxis.set_tick_params(labelsize=15)
ax2.grid(False)
if pp == 0:
axs[0].legend(loc='lower right', ncol=2,fontsize=13,fancybox=True,frameon=True)
for jj, global_quantity in enumerate(plot_quantities):
axs[jj].set_ylabel('global_%s'%global_quantity.replace('_box',''),fontsize=15)
axs[jj].yaxis.set_tick_params(labelsize=15)
else:
for jj, global_quantity in enumerate(plot_quantities):
axs[jj].set_ylabel('global_%s'%global_quantity.replace('_box',''),fontsize=0)
axs[jj].yaxis.set_tick_params(labelsize=0)
plt.tight_layout()
fig.subplots_adjust(hspace = 0.0,wspace=0.0)
# varying parameters
let’s vary the parameters that describe mini-halos and see the impact to the global signal. Warning: It may take a while to run all these boxes!
We keep other parameters fixed and vary one of following by a factor of 1/3 and 3:
F_STAR7_MINI
F_ESC7_MINI
L_X_MINI
A_LW
We also have a NOmini model where mini-halos are not included
[ ]:
#defining those color, linstyle, blabla
linestyles = ['-', '-',':','-.','-.',':']
colors = ['gray','black','#e41a1c','#377eb8','#e41a1c','#377eb8']
lws = [1,3,2,2,2,2]
textss = ['varying '+r'$f_{*,7}^{\rm mol}$',\
'varying '+r'$f_{\rm esc}^{\rm mol}$',\
'varying '+r'$L_{\rm x}^{\rm mol}$',\
'varying '+r'$A_{\rm LW}$']
factorss = [[0, 1, 0.33, 3.0],] * len(textss)
labelss = [['No Velocity', 'Fiducial', '/3', 'x3'],] * len(textss)
[ ]:
global_quantities = ('brightness_temp','Ts_box','xH_box',"dNrec_box",'z_re_box','Gamma12_box','J_21_LW_box',"density")
#choose some to plot...
plot_quantities = ('brightness_temp','Ts_box','xH_box',"dNrec_box",'Gamma12_box','J_21_LW_box')
ymins = [-120, 1e1, 0, 0, 0, 0]
ymaxs = [ 30, 1e3, 1, 1, 1,10]
fig, axss = plt.subplots(len(plot_quantities), len(labelss),
sharex=True, figsize=(4*len(labelss),2*len(global_quantities)))
for pp, texts in enumerate(textss):
labels = labelss[pp]
factors = factorss[pp]
axs = axss[:,pp]
for kk, label in enumerate(labels):
flag_options = flag_options_fid.copy()
astro_params = astro_params_vcb.copy()
factor = factors[kk]
if label == 'No Velocity':
lightcone = lightcone_fid_novcb
elif label == 'Fiducial':
lightcone = lightcone_fid_vcb
else:
if pp == 0:
astro_params.update({'F_STAR7_MINI': astro_params_vcb['F_STAR7_MINI']+np.log10(factor)})
elif pp == 1:
astro_params.update({'F_ESC7_MINI': astro_params_vcb['F_ESC7_MINI']+np.log10(factor)})
elif pp == 2:
astro_params.update({'L_X_MINI': astro_params_vcb['L_X_MINI']+np.log10(factor)})
elif pp == 3:
astro_params.update({'A_LW': astro_params_vcb['A_LW']*factor})
else:
print('Make a choice!')
lightcone = p21c.run_lightcone(
redshift = ZMIN,
init_box = initial_conditions,
flag_options = flag_options_fid,
astro_params = astro_params,
global_quantities=global_quantities,
random_seed = random_seed,
direc = output_dir
)
freqs = 1420.4 / (np.array(lightcone.node_redshifts) + 1.)
for jj, global_quantity in enumerate(plot_quantities):
axs[jj].plot(freqs, getattr(lightcone, 'global_%s'%global_quantity.replace('_box','')),
color=colors[kk],linestyle=linestyles[kk], label = labels[kk],lw=lws[kk])
axs[0].text(0.01, 0.99, texts,horizontalalignment='left',verticalalignment='top',
transform=axs[0].transAxes,fontsize = 15)
for jj, global_quantity in enumerate(plot_quantities):
axs[jj].set_ylim(ymins[jj], ymaxs[jj])
axs[-1].set_xlabel('Frequency/MHz',fontsize=15)
axs[-1].xaxis.set_tick_params(labelsize=15)
axs[0].set_xlim(1420.4 / (35 + 1.), 1420.4 / (5.5 + 1.))
zlabels = np.array([ 6, 7, 8, 10, 13, 18, 25, 35])
ax2 = axs[0].twiny()
ax2.set_xlim(axs[0].get_xlim())
ax2.set_xticks(1420.4 / (zlabels + 1.))
ax2.set_xticklabels(zlabels.astype(np.str))
ax2.set_xlabel("redshift",fontsize=15)
ax2.xaxis.set_tick_params(labelsize=15)
ax2.grid(False)
if pp == 0:
axs[0].legend(loc='lower right', ncol=2,fontsize=13,fancybox=True,frameon=True)
for jj, global_quantity in enumerate(plot_quantities):
axs[jj].set_ylabel('global_%s'%global_quantity.replace('_box',''),fontsize=15)
axs[jj].yaxis.set_tick_params(labelsize=15)
else:
for jj, global_quantity in enumerate(plot_quantities):
axs[jj].set_ylabel('global_%s'%global_quantity.replace('_box',''),fontsize=0)
axs[jj].yaxis.set_tick_params(labelsize=0)
plt.tight_layout()
fig.subplots_adjust(hspace = 0.0,wspace=0.0)
[ ]:
[ ]:
# define functions to calculate PS, following py21cmmc
from powerbox.tools import get_power
def compute_power(
box,
length,
n_psbins,
log_bins=True,
ignore_kperp_zero=True,
ignore_kpar_zero=False,
ignore_k_zero=False,
):
# Determine the weighting function required from ignoring k's.
k_weights = np.ones(box.shape, dtype=int)
n0 = k_weights.shape[0]
n1 = k_weights.shape[-1]
if ignore_kperp_zero:
k_weights[n0 // 2, n0 // 2, :] = 0
if ignore_kpar_zero:
k_weights[:, :, n1 // 2] = 0
if ignore_k_zero:
k_weights[n0 // 2, n0 // 2, n1 // 2] = 0
res = get_power(
box,
boxlength=length,
bins=n_psbins,
bin_ave=False,
get_variance=False,
log_bins=log_bins,
k_weights=k_weights,
)
res = list(res)
k = res[1]
if log_bins:
k = np.exp((np.log(k[1:]) + np.log(k[:-1])) / 2)
else:
k = (k[1:] + k[:-1]) / 2
res[1] = k
return res
def powerspectra(brightness_temp, n_psbins=50, nchunks=10, min_k=0.1, max_k=1.0, logk=True):
data = []
chunk_indices = list(range(0,brightness_temp.n_slices,round(brightness_temp.n_slices / nchunks),))
if len(chunk_indices) > nchunks:
chunk_indices = chunk_indices[:-1]
chunk_indices.append(brightness_temp.n_slices)
for i in range(nchunks):
start = chunk_indices[i]
end = chunk_indices[i + 1]
chunklen = (end - start) * brightness_temp.cell_size
power, k = compute_power(
brightness_temp.brightness_temp[:, :, start:end],
(BOX_LEN, BOX_LEN, chunklen),
n_psbins,
log_bins=logk,
)
data.append({"k": k, "delta": power * k ** 3 / (2 * np.pi ** 2)})
return data
[ ]:
# do 5 chunks but only plot 1 - 4, the 0th has no power for minihalo models where xH=0
nchunks = 4
k_fundamental = 2*np.pi/BOX_LEN
k_max = k_fundamental * HII_DIM
Nk=np.floor(HII_DIM/1).astype(int)
fig, axss = plt.subplots(nchunks, len(labelss), sharex=True,sharey=True,figsize=(4*len(labelss),3*(nchunks)),subplot_kw={"xscale":'log', "yscale":'log'})
for pp, texts in enumerate(textss):
labels = labelss[pp]
factors = factorss[pp]
axs = axss[:,pp]
for kk, label in enumerate(labels):
flag_options = flag_options_fid.copy()
astro_params = astro_params_vcb.copy()
factor = factors[kk]
if label == 'No Velocity':
lightcone = lightcone_fid_novcb
elif label == 'Fiducial':
lightcone = lightcone_fid_vcb
else:
if pp == 0:
astro_params.update({'F_STAR7_MINI': astro_params_vcb['F_STAR7_MINI']+np.log10(factor)})
elif pp == 1:
astro_params.update({'F_ESC7_MINI': astro_params_vcb['F_ESC7_MINI']+np.log10(factor)})
elif pp == 2:
astro_params.update({'L_X_MINI': astro_params_vcb['L_X_MINI']+np.log10(factor)})
elif pp == 3:
astro_params.update({'A_LW': astro_params_vcb['A_LW']+np.log10(factor)})
else:
print('Make a choice!')
lightcone = p21c.run_lightcone(
redshift = ZMIN,
init_box = initial_conditions,
flag_options = flag_options_fid,
astro_params = astro_params,
global_quantities=global_quantities,
random_seed = random_seed,
direc = output_dir
)
PS = powerspectra(lightcone, min_k = k_fundamental, max_k = k_max)
for ii in range(nchunks):
axs[ii].plot(PS[ii+1]['k'], PS[ii+1]['delta'], color=colors[kk],linestyle=linestyles[kk], label = labels[kk],lw=lws[kk])
if pp == len(textss)-1 and kk == 0:
axs[ii].text(0.99, 0.01, 'Chunk-%02d'%(ii+1),horizontalalignment='right',verticalalignment='bottom',
transform=axs[ii].transAxes,fontsize = 15)
axs[0].text(0.01, 0.99, texts,horizontalalignment='left',verticalalignment='top',
transform=axs[0].transAxes,fontsize = 15)
axs[-1].set_xlabel("$k$ [Mpc$^{-3}$]",fontsize=15)
axs[-1].xaxis.set_tick_params(labelsize=15)
if pp == 0:
for ii in range(nchunks):
axs[ii].set_ylim(2e-1, 2e2)
axs[ii].set_ylabel("$k^3 P$", fontsize=15)
axs[ii].yaxis.set_tick_params(labelsize=15)
else:
for ii in range(nchunks-1):
axs[ii].set_ylim(2e-1, 2e2)
axs[ii].set_ylabel("$k^3 P$", fontsize=0)
axs[ii].yaxis.set_tick_params(labelsize=0)
axss[0,0].legend(loc='lower left', ncol=2,fontsize=13,fancybox=True,frameon=True)
plt.tight_layout()
fig.subplots_adjust(hspace = 0.0,wspace=0.0)
Note that I’ve run these simulations in parallel before this tutorial. With these setup, each took ~6h to finish. Here, running means read the cached outputs.
global properties – optical depth¶
[ ]:
#defining those color, linstyle, blabla
linestyles = ['-', '-',':','-.','-.',':']
colors = ['gray','black','#e41a1c','#377eb8','#e41a1c','#377eb8']
lws = [1,3,2,2,2,2]
textss_tau = ['varying '+r'$f_{*,7}^{\rm mol}$',\
'varying '+r'$f_{\rm esc}^{\rm mol}$',\
'varying '+r'$A_{\rm LW}$']
factorss_tau = [[0, 1, 0.33, 3.0],] * len(textss_tau)
labelss_tau = [['No Velocity', 'Fiducial', '/3', 'x3'],] * len(textss_tau)
[ ]:
plot_quantities = ['tau_e']
ymins = [0]
ymaxs = [0.2]
fig, axss_tau = plt.subplots(len(plot_quantities), len(labelss_tau),
sharex=True, figsize=(4*len(labelss_tau),3*len(plot_quantities)))
for pp, texts in enumerate(textss_tau):
labels = labelss_tau[pp]
factors = factorss_tau[pp]
axs = axss_tau[pp]
for kk, label in enumerate(labels):
flag_options = flag_options_fid.copy()
astro_params = astro_params_vcb.copy()
factor = factors[kk]
if label == 'No Velocity':
lightcone = lightcone_fid_novcb
elif label == 'Fiducial':
lightcone = lightcone_fid_vcb
else:
if pp == 0:
astro_params.update({'F_STAR7_MINI': astro_params_vcb['F_STAR7_MINI']+np.log10(factor)})
elif pp == 1:
astro_params.update({'F_ESC7_MINI': astro_params_vcb['F_ESC7_MINI']+np.log10(factor)})
elif pp == 2:
astro_params.update({'A_LW': astro_params_vcb['A_LW']*factor})
else:
print('Make a choice!')
lightcone = p21c.run_lightcone(
redshift = ZMIN,
init_box = initial_conditions,
flag_options = flag_options_fid,
astro_params = astro_params,
global_quantities=global_quantities,
random_seed = random_seed,
direc = output_dir
)
z_e=np.array([]);
tau_e=np.array([]);
for i in range(len(lightcone.node_redshifts)-1):
tauz=p21c.compute_tau(redshifts=lightcone.node_redshifts[-1:-2-i:-1],
global_xHI=lightcone.global_xHI[-1:-2-i:-1])
tau_e=np.append(tau_e,tauz)
z_e=np.append(z_e,lightcone.node_redshifts[-2-i])
#print(i,lightcone.node_redshifts[i],tauz)
#add lower zs where we manually set xH=1
zlow=np.linspace(lightcone.node_redshifts[-1]-0.1, 0.1, 14)
for zl in zlow:
tauz=p21c.compute_tau(redshifts=np.array([zl]), global_xHI=np.array([lightcone.global_xHI[-1]]))
tau_e=np.append([tauz],tau_e)
z_e=np.append([zl],z_e)
# freqs = 1420.4 / (np.array(lightcone.node_redshifts) + 1.)
for jj, global_quantity in enumerate(plot_quantities):
axs.plot(z_e, tau_e,
color=colors[kk],linestyle=linestyles[kk], label = labels[kk],lw=lws[kk])
axs.text(0.01, 0.99, texts,horizontalalignment='left',verticalalignment='top',
transform=axs.transAxes,fontsize = 15)
axs.set_ylim(ymins[0], ymaxs[0])
axs.set_xlabel('redshift',fontsize=15)
axs.xaxis.set_tick_params(labelsize=15)
axs.set_xlim(0.,20.)
if pp == 0:
for ii in range(nchunks):
axs.set_ylabel('$\\tau$',fontsize=15)
axs.yaxis.set_tick_params(labelsize=15)
else:
for ii in range(nchunks-1):
axs.yaxis.set_tick_params(labelsize=0)
plt.tight_layout()
fig.subplots_adjust(hspace = 0.0,wspace=0.0)
[ ]:
[ ]:
21-cm power spectra¶
[ ]:
# do 5 chunks but only plot 1 - 4, the 0th has no power for minihalo models where xH=0
nchunks = 4
fig, axss = plt.subplots(nchunks, len(labelss), sharex=True,sharey=True,figsize=(4*len(labelss),3*(nchunks)),subplot_kw={"xscale":'log', "yscale":'log'})
for pp, texts in enumerate(textss):
labels = labelss[pp]
factors = factorss[pp]
axs = axss[:,pp]
for kk, label in enumerate(labels):
factor = factors[kk]
if kk==0:
lightcone = lightcone_fid_vcbavg
else:
lightcone = lightcone_fid_vcb
PS = powerspectra(lightcone, min_k = k_fundamental, max_k = k_max)
for ii in range(nchunks):
axs[ii].plot(PS[ii+1]['k'], PS[ii+1]['delta'], color=colors[kk],linestyle=linestyles[kk], label = labels[kk],lw=lws[kk])
if pp == len(textss)-1 and kk == 0:
axs[ii].text(0.99, 0.01, 'Chunk-%02d'%(ii+1),horizontalalignment='right',verticalalignment='bottom',
transform=axs[ii].transAxes,fontsize = 15)
axs[0].text(0.01, 0.99, texts,horizontalalignment='left',verticalalignment='top',
transform=axs[0].transAxes,fontsize = 15)
axs[-1].set_xlabel("$k$ [Mpc$^{-3}$]",fontsize=15)
axs[-1].xaxis.set_tick_params(labelsize=15)
if pp == 0:
for ii in range(nchunks):
axs[ii].set_ylim(2e-1, 2e2)
axs[ii].set_ylabel("$k^3 P$", fontsize=15)
axs[ii].yaxis.set_tick_params(labelsize=15)
else:
for ii in range(nchunks-1):
axs[ii].set_ylim(2e-1, 2e2)
axs[ii].set_ylabel("$k^3 P$", fontsize=0)
axs[ii].yaxis.set_tick_params(labelsize=0)
axss[0,0].legend(loc='lower left', ncol=2,fontsize=13,fancybox=True,frameon=True)
plt.tight_layout()
fig.subplots_adjust(hspace = 0.0,wspace=0.0)
[ ]:
nchunks=5
textss = ['vcb','vcb']
factorss = [[0, 1],] * len(textss)
labelss = [['Regular', 'Avg'],] * len(textss)
k_fundamental = 2*np.pi/BOX_LEN
k_max = k_fundamental * HII_DIM
Nk=np.floor(HII_DIM/1).astype(int)
PSv= powerspectra(lightcone_fid_vcb, min_k = k_fundamental, max_k = k_max)
PSvavg= powerspectra(lightcone_fid_vcbavg, min_k = k_fundamental, max_k = k_max)
[ ]:
klist= PSv[0]['k']
P21diff = [(PSv[i]['delta']-PSvavg[i]['delta'])/PSvavg[i]['delta'] for i in range(nchunks)]
import matplotlib.pyplot as plt
fig, axss = plt.subplots(nchunks, 1, sharex=True,sharey=True,figsize=(2*len(labelss),3*(nchunks)),subplot_kw={"xscale":'linear', "yscale":'linear'})
for ii in range(nchunks):
axss[ii].plot(klist, P21diff[ii])
plt.xscale('log')
axss[0].legend(loc='lower left', ncol=2,fontsize=13,fancybox=True,frameon=True)
plt.tight_layout()
fig.subplots_adjust(hspace = 0.0,wspace=0.0)
[ ]: