py: read grids and global_xH using dragons

grab_reion_data.py

import numpy as np
from dragons import meraxes

fname = "/lustre/projects/p021_astro/smutch/meraxes/tiamat/paper_runs/uvb/output/meraxes.hdf5"

# grab the snapshot list
snaplist, _, _ = meraxes.read_snaplist(fname)

# read all of the global_xH values
global_xH = meraxes.read_global_xH(fname, snaplist, quiet=True)

# use the global_xH values to filter out the snapshots that have no grids
snaplist = snaplist[~np.isnan(global_xH)]
global_xH = global_xH[~np.isnan(global_xH)]

# now snaplist only contains the snapshots that have grids
snap = snaplist[6]
print "Reading grid for snapshot %d..." % snap
xH = meraxes.read_grid(fname, snap, "xH")

print "Neutral fraction = %.3f" % xH.mean()
print "Reported global neutral fraction = %.3f" % global_xH[6]

# mock up full global neutral fraction array by assuming xH = 1.0 before
# reionization and xH = 0.0 after...
snaplist, _, _ = meraxes.read_snaplist(fname)
global_xH = meraxes.read_global_xH(fname, snaplist, quiet=True)
reion_start, reion_end = np.where(~np.isnan(global_xH))[0][[0, -1]]
global_xH[:reion_start] = 1.0
global_xH[reion_end+1:] = 0.0

print "Global_xH history:"
print global_xH