py: Sobacchi global model params fitting.

sobacchi_params_from_run.py

#!/usr/bin/env python

"""Calculate the Sobacchi & Mesinger global model parameters from a run.

Usage: sobacchi_params_from_run.py <fname>
"""

import numpy as np
import matplotlib.pyplot as plt
from dragons import meraxes
from docopt import docopt
from astropy import cosmology
import astropy.units as U

__author__ = "Simon Mutch"
__date__ = "2015-05-22"

def grab_xHII(fname):
    """Grab the snapshot, xHI and redshift lists."""
    snap, redshift, _ = meraxes.read_snaplist(fname)
    xHI = meraxes.read_global_xH(fname, snap, quiet=True) 
    sel = np.isfinite(xHI)
    snap, xHI, redshift = snap[sel][::-1], xHI[sel][::-1], redshift[sel][::-1]
    return snap, xHI, redshift

def find_midpoint(redshift, xHI):
    """Interpolate the redshift of xHI=0.5."""
    midpoint = np.interp(0.5, xHI, redshift)
    return midpoint

def find_duration(redshift, xHI):
    """Find the redshift interval spanned between xHI=0.4 & 0.6."""
    end = np.interp(0.4, xHI, redshift)
    start = np.interp(0.6, xHI, redshift)
    duration = start - end
    return duration

def delta_crit(z, cosmo):
    """Critical overdensity at a given redshift."""
    d = cosmo.Om(z) - 1.0
    return 18.*(np.pi**2) + 82.0*d - 39.0*(d**2)

def Tvir_to_Mvir(T, z, cosmo, mu=0.6):
    """Convert Tvir to Mvir."""
    z_term = ((1.+z)/10.)**(-3./2.)
    T_term = (T / 1.98e4)**(3./2.)
    cosmo_term = (cosmo.Om0/cosmo.Om(z) * delta_crit(z, cosmo)/18./(np.pi**2))**(-1./2.)
    mol_term = (mu/0.6)**(-3./2.)
    factor = 1.0e8*cosmo.h
    return factor * mol_term * cosmo_term * T_term * z_term

def find_sound_crossing_timescale(midpoint, cosmo):
    """Work out the redshift interval of a sound crossing time for a Tvir=5e4K
    halo at z=z_{XHI=0.5}."""
    T0 = 5.0e4  # Kelvin
    M0 = Tvir_to_Mvir(T0, midpoint, cosmo)
    sc_time = 200.0 * (M0/1.0e8)**(1./3.) * ((1+midpoint)/10.)**(-1.) * (cosmo.Om0*(cosmo.h**2)/0.15)**(-1./3.) * U.Myr
    mp_lbt = cosmo.lookback_time(midpoint)
    delta_z = midpoint - cosmology.z_at_value(cosmo.lookback_time, mp_lbt - sc_time)
    return delta_z

if __name__ == '__main__':
    args = docopt(__doc__)
    fname = args['<fname>']
    little_h = 0.7
    meraxes.set_little_h(little_h)
    params = meraxes.read_input_params(fname)
    cosmo = cosmology.core.FlatLambdaCDM(H0=little_h * 100., Om0=params['OmegaM'], name="Tiamat")

    snap, xHI, redshift = grab_xHII(fname)
    midpoint = find_midpoint(redshift, xHI)
    duration = find_duration(redshift, xHI)
    delta_z_sc = find_sound_crossing_timescale(midpoint, cosmo)

    print "Found:"
    print "z_re = %g" % midpoint
    print "Dz_re = %g" % duration
    print "Dz_sc = %g" % delta_z_sc
    print