244 Mpc/h

Here we download and analyse the simulation data produced in for a cosmological volume of size 244/h Mpc in each direction.

[1]:
import numpy as np
import matplotlib.pyplot as plt
import pickle, os
from glob import glob
from StoReS import *
import tools21cm as t2c
from tqdm import tqdm

N-body simulation with 244 Mpc/h in each direction

Cosmology of the simulation:

\(\Omega_\mathrm{m} = 0.27\), \(\Omega_\mathrm{b} = 0.044\), \(h_\mathrm{0} = 0.7\), \(\sigma_\mathrm{8} = 0.8\), \(n_\mathrm{s} = 0.96\)

We define the box length and the path to the folder for downloading the data.

[2]:
box_len  = 244/0.7
save_dir = './work/'
eor_hist = {}

Density fields - 244 Mpc/h

N-Body (CUBEP3M) snapshots will be downloaded.

[3]:
dn_url = 'https://ttt.astro.su.se/~gmell/244Mpc/densities/nc250/coarser_densities/'
url_dict = {'dens': dn_url}
c2r = C2RAY(work_dir={'dens': save_dir+'density/'}, verbose=False)
c2r.set_links(url_dict)
zs = c2r.zs_dict['dens']
for i,z in tqdm(enumerate(zs)):
    dn, ff = c2r.get_density_data_z(z)
    # os.remove(ff)
print('Data saved in {}'.format(c2r.work_dir['dens']))
77it [00:07, 10.79it/s]
Data saved in ./work/density/

244Mpc/h Reionisation Simulation Suite

LB1

This simulation models the reionization process assuming spin saturation (\(T_{S}\gg T_{CMB}\)). is described in See Dixon et al. (2016) the detailed description.

[4]:
xf_url = 'https://ttt.astro.su.se/~gmell/244Mpc/LB1/'
dn_url = 'https://ttt.astro.su.se/~gmell/244Mpc/densities/nc250/coarser_densities/'
url_dict = {'xfrac': xf_url, 'dens': dn_url, }
c2r = C2RAY(work_dir={'dens': save_dir+'density/', 'xfrac': save_dir+'LB1/'}, verbose=False)
c2r.set_links(url_dict)
zs = np.intersect1d(c2r.zs_dict['xfrac'], c2r.zs_dict['dens'])
xs = []
for i,z in tqdm(enumerate(zs)):
    data = c2r.get_data_z(z)
    xf = data['xfrac']
    dn = data['dens']
    xs.append(xf.mean())
    # os.remove(data['xfrac_filename'])
    # os.remove(data['dens_filename'])
xs = np.array(xs)
eor_hist['LB1'] = {'z': zs, 'x': xs}
print('Data saved in {}'.format(c2r.work_dir['xfrac']))

LB2

This simulation models the reionization process assuming spin saturation (\(T_{S}\gg T_{CMB}\)). is described in See Dixon et al. (2016) the detailed description.

[5]:
xf_url = 'https://ttt.astro.su.se/~gmell/244Mpc/LB2/'
dn_url = 'https://ttt.astro.su.se/~gmell/244Mpc/densities/nc250/coarser_densities/'
url_dict = {'xfrac': xf_url, 'dens': dn_url, }
c2r = C2RAY(work_dir={'dens': save_dir+'density/', 'xfrac': save_dir+'LB2/'}, verbose=False)
c2r.set_links(url_dict)
zs = np.intersect1d(c2r.zs_dict['xfrac'], c2r.zs_dict['dens'])
xs = []
for i,z in tqdm(enumerate(zs)):
    data = c2r.get_data_z(z)
    xf = data['xfrac']
    dn = data['dens']
    xs.append(xf.mean())
    # os.remove(data['xfrac_filename'])
    # os.remove(data['dens_filename'])
xs = np.array(xs)
eor_hist['LB2'] = {'z': zs, 'x': xs}
print('Data saved in {}'.format(c2r.work_dir['xfrac']))

LB3

This simulation models the reionization process assuming spin saturation (\(T_{S}\gg T_{CMB}\)). is described in See Dixon et al. (2016) the detailed description.

[6]:
xf_url = 'https://ttt.astro.su.se/~gmell/244Mpc/LB3/'
dn_url = 'https://ttt.astro.su.se/~gmell/244Mpc/densities/nc250/coarser_densities/'
url_dict = {'xfrac': xf_url, 'dens': dn_url, }
c2r = C2RAY(work_dir={'dens': save_dir+'density/', 'xfrac': save_dir+'LB3/'}, verbose=False)
c2r.set_links(url_dict)
zs = np.intersect1d(c2r.zs_dict['xfrac'], c2r.zs_dict['dens'])
xs = []
for i,z in tqdm(enumerate(zs)):
    data = c2r.get_data_z(z)
    xf = data['xfrac']
    dn = data['dens']
    xs.append(xf.mean())
    # os.remove(data['xfrac_filename'])
    # os.remove(data['dens_filename'])
xs = np.array(xs)
eor_hist['LB3'] = {'z': zs, 'x': xs}
print('Data saved in {}'.format(c2r.work_dir['xfrac']))

LB4

This simulation models the reionization process assuming spin saturation (\(T_{S}\gg T_{CMB}\)). is described in See Dixon et al. (2016) the detailed description.

[7]:
xf_url = 'https://ttt.astro.su.se/~gmell/244Mpc/LB4/'
dn_url = 'https://ttt.astro.su.se/~gmell/244Mpc/densities/nc250/coarser_densities/'
url_dict = {'xfrac': xf_url, 'dens': dn_url, }
c2r = C2RAY(work_dir={'dens': save_dir+'density/', 'xfrac': save_dir+'LB4/'}, verbose=False)
c2r.set_links(url_dict)
zs = np.intersect1d(c2r.zs_dict['xfrac'], c2r.zs_dict['dens'])
xs = []
for i,z in tqdm(enumerate(zs)):
    data = c2r.get_data_z(z)
    xf = data['xfrac']
    dn = data['dens']
    xs.append(xf.mean())
    # os.remove(data['xfrac_filename'])
    # os.remove(data['dens_filename'])
xs = np.array(xs)
eor_hist['LB4'] = {'z': zs, 'x': xs}
print('Data saved in {}'.format(c2r.work_dir['xfrac']))

Reionization history of the above simulations

[8]:
fig, ax = plt.subplots(1,1,figsize=(7,6))
# fig.suptitle('Dixon et al. (2016)')
ax.plot(eor_hist['LB1']['z'], 1-eor_hist['LB1']['x'], lw=3, ls='-',  label='LB1')
ax.plot(eor_hist['LB2']['z'], 1-eor_hist['LB2']['x'], lw=3, ls='--', label='LB2')
ax.plot(eor_hist['LB3']['z'], 1-eor_hist['LB3']['x'], lw=3, ls='-.', label='LB3')
ax.plot(eor_hist['LB4']['z'], 1-eor_hist['LB4']['x'], lw=3, ls=':',  label='LB4')
ax.set_xlabel('z', fontsize=14)
ax.set_ylabel('$x_\mathrm{HI}$', fontsize=14)
ax.legend()
ax.axis([6,14,-0.01,1.01])
plt.tight_layout()
plt.show()
../_images/examples_244Mpc_18_0.png

Power spectrum at z=6.9 from LB1

[9]:
xf_url = 'https://ttt.astro.su.se/~gmell/244Mpc/LB1/'
dn_url = 'https://ttt.astro.su.se/~gmell/244Mpc/densities/nc250/coarser_densities/'
url_dict = {'dens': dn_url, 'xfrac': xf_url}

c2r = C2RAY(work_dir={'dens': save_dir+'density/', 'xfrac': save_dir+'LB1/'})
c2r.set_links(url_dict)
z = 6.905
data = c2r.get_data_z(z)
xf = data['xfrac']
dn = data['dens']

Downloading simulation...
100% [..................................................] 125000028 / 125000028...done
The file already exists.
[11]:
dt = t2c.calc_dt(xf, dn, z=z)
ps, ks = t2c.power_spectrum_1d(dt, kbins=20, box_dims=box_len)
[12]:
fig, ax = plt.subplots(1,1,figsize=(7,6))
fig.suptitle('LB1')
ax.loglog(ks, ps*ks**3/2/np.pi**2, lw=3, ls='-',  label='z={:.1f}'.format(z))
ax.set_xlabel('k [1/Mpc]', fontsize=14)
ax.set_ylabel('$\Delta^2_\mathrm{21}$ [mK$^2$]', fontsize=14)
ax.legend()
# ax.axis([6,14,-0.01,1.01])
plt.tight_layout()
plt.show()
../_images/examples_244Mpc_22_0.png