In [1]:
import glob, copy
import argparse
import yaml, h5py
import numpy as NP
from astropy.io import ascii
from astropy import units as U
import astropy.cosmology as cosmology
import astropy.constants as FC
import matplotlib.pyplot as PLT
import matplotlib.patches as patches
import matplotlib.ticker as ticker
from matplotlib.lines import Line2D
from astroutils import cosmotile
from astroutils import mathops as OPS
from astroutils import DSP_modules as DSP
from astroutils import lookup_operations as LKP
from astroutils import constants as CNST
from prisim import delay_spectrum as DS
from IPython.core.debugger import set_trace
%matplotlib inline

In [2]:
cosmoPlanck15 = cosmology.Planck15 # Planck 2015 cosmology
cosmo100_Planck15 = cosmoPlanck15.clone(name='Modified Planck 2015 cosmology with h=1.0', H0=100.0) # Modified Planck 2015 cosmology with h=1.0, H= 100 km/s/Mpc

Define some useful functions


In [3]:
def read_raw_lightcone_cube(parms):

#     cube_source = parms['cube_source']
#     if not isinstance(cube_source, str):
#         raise TypeError('Input cube_source must be a string')
#     if cube_source.lower() not in ['21cmfast']:
#         raise ValueError('{0} cubes currently not supported'.format(cube_source))

    indir = parms['indir']
    infile_prefix = parms['infile_prefix']
    if infile_prefix is None:
        infile_prefix = ''
    infile_suffix = parms['infile_suffix']
    if infile_suffix is None:
        infile_suffix = ''

    fname_delimiter = parms['parseinfo']['delimiter']
    zstart_pos = parms['parseinfo']['zstart_pos']
    zend_pos = parms['parseinfo']['zend_pos']
    zstart_identifier = parms['parseinfo']['zstart_identifier']
    zend_identifier = parms['parseinfo']['zend_identifier']
    zstart_identifier_pos = parms['parseinfo']['zstart_identifier_pos']
    zend_identifier_pos = parms['parseinfo']['zend_identifier_pos']
    if zstart_identifier is not None:
        if zstart_identifier_pos.lower() not in ['before', 'after']:
            raise ValueError('zstart_identifier_pos must be set to "before" or "after"')
        elif zstart_identifier_pos.lower() == 'before':
            zstart_value_place = 1
        else:
            zstart_value_place = 0

    if zend_identifier is not None:
        if zend_identifier_pos.lower() not in ['before', 'after']:
            raise ValueError('zend_identifier_pos must be set to "before" or "after"')
        elif zend_identifier_pos.lower() == 'before':
            zend_value_place = 1
        else:
            zend_value_place = 0
            
    fullfnames = glob.glob(indir + infile_prefix + '*' + infile_suffix)
    fullfnames = NP.asarray(fullfnames)
    fnames = [fname.split('/')[-1] for fname in fullfnames]
    fnames = NP.asarray(fnames)
    if fnames[0].split(fname_delimiter)[-1] == 'lighttravel':
        dim = int(fnames[0].split(fname_delimiter)[-3])
        boxsize = float(fnames[0].split(fname_delimiter)[-2][:-3])
    else:
        dim = int(fnames[0].split(fname_delimiter)[-2])
        boxsize = float(fnames[0].split(fname_delimiter)[-1][:-3])
    boxres = boxsize / dim
    if zstart_identifier is not None:
        zstart_str = [fname.replace(fname_delimiter,' ').split()[zstart_pos].split(zstart_identifier)[zstart_value_place] for fname in fnames]
    else:
        zstart_str = [fname.replace(fname_delimiter,' ').split()[zstart_pos] for fname in fnames]
    if zend_identifier is not None:
        zend_str = [fname.replace(fname_delimiter,' ').split()[zend_pos].split(zend_identifier)[zend_value_place] for fname in fnames]
    else:
        zend_str = [fname.replace(fname_delimiter,' ').split()[zend_pos] for fname in fnames]
    
    zstart = NP.asarray(map(float, zstart_str))
    zend = NP.asarray(map(float, zend_str))
    sortind = NP.argsort(zstart)
    zstart = zstart[sortind]
    zend = zend[sortind]
    fnames = fnames[sortind]
    fullfnames = fullfnames[sortind]
    lightcone_cube = None
    for fi,fullfname in enumerate(fullfnames):
        lc_cube = cosmotile.fastread_21cmfast_cube(fullfname)
        if lightcone_cube is None:
            lightcone_cube = NP.copy(lc_cube)
        else:
            lightcone_cube = NP.concatenate((lightcone_cube, lc_cube), axis=0) # Line of sight seems to be the first axis
    return (lightcone_cube, boxres, zstart.min())

In [4]:
def monopole(lightcone_cube):
    """
    Assumed to have shape (nlos,...) where the average will be estimated across all except the first axis
    """
    if lightcone_cube.ndim <= 1:
        raise ValueError('Insufficient axes to average')
    avgaxes = NP.arange(1,lightcone_cube.ndim)
    return NP.mean(lightcone_cube, axis=tuple(avgaxes))

In [5]:
def write_raw_lightcone_cube(lightcone_cube, boxres, zmin, outfile, units='',
                             cosmo_in=cosmo100_Planck15):

    """
    Write a lightcone cube in 21cmFAST format to a file in HDF5 format
    """

    cosmo100 = cosmo.clone(name='Modified cosmology with h=1.0', H0=100.0) # Modified cosmology with h=1.0, H= 100 km/s/Mpc

    monopole_values = monopole(lightcone_cube)

    with h5py.File(outfile, 'w') as fileobj:
        hdr_grp = fileobj.create_group('header')
        hdr_grp['boxres'] = cosmo_in.h * boxres
        hdr_grp['boxres'].attrs['units'] = 'Mpc'
        hdr_grp['zmin'] = zmin
        cosmo_grp = fileobj.create_group('cosmo')
        cosmo_grp['H0'] = cosmo100.H0.value
        cosmo_grp['H0'].attrs['units'] = '{0}'.format(cosmo100.H0.unit)
        cosmo_grp['Om0'] = cosmo100.Om0
        cosmo_grp['Ode0'] = cosmo100.Ode0
        cosmo_grp['Ob0'] = cosmo100.Ob0
        cosmo_grp['w0'] = cosmo100.w(0.0)
        cosmo_grp['Tcmb0'] = cosmo100.Tcmb0.value
        cosmo_grp['Tcmb0'].attrs['units'] = '{0}'.format(cosmo100.Tcmb0.unit)
        cube_dset = fileobj.create_dataset('cube', lightcone_cube.shape, data=lightcone_cube)
        cube_dset.attrs['units'] = units
        monopole_dset = fileobj.create_dataset('monopole', monopole_values.shape, data=monopole_values)
        monopole_dset.attrs['units'] = units

In [6]:
def read_HDF5_lightcone_cube(infile):

    with h5py.File(infile, 'r') as fileobj:
        hdrinfo = {key: fileobj['header'][key].value for key in fileobj['header']}
        cosmoinfo = {key: fileobj['cosmo'][key].value for key in fileobj['cosmo']}
        lightcone_cube_units = fileobj['cube'].attrs['units']
        lightcone_cube = fileobj['cube'].value * U.Unit(lightcone_cube_units)
        if 'monopole_temperature' in fileobj:
            monopole_values = fileobj['monopole_temperature'] * U.Unit(lightcone_cube_units)
        elif 'monopole' in fileobj:
            monopole_values = fileobj['monopole'] * U.Unit(lightcone_cube_units)
        else:
            monopole_values = None

    boxres = hdrinfo['boxres'] * U.Unit('Mpc')
    zmin = hdrinfo['zmin']
    cosmo = cosmology.wCDM(cosmoinfo['H0'], cosmoinfo['Om0'], cosmoinfo['Ode0'], w0=cosmoinfo['w0'], Ob0=cosmoinfo['Ob0'], Tcmb0=cosmoinfo['Tcmb0']*U.K)
    
    return (lightcone_cube, monopole_values, boxres, zmin, cosmo)

In [7]:
def optical_depth_21cm(spin_T, x_H, overdensity, redshift, cosmo_in=cosmo100_Planck15, reference='Santos+2005'):

    """
    Compute 21cm optical depth based on the expression provided in the specified reference, T_spin, neutral fraction, 
    overdensity, and cosmological parameters
    """
    
    if reference.lower() not in ['santos+2005', 'pritchard_loeb_2010', 'zaldarriaga+2004', 'tozzi+2000', 'bharadwaj_ali_2005', 'cooray_2004', 'cooray_2005']:
        raise ValueError('Expression from reference {0} not supported'.format(reference))

    if not isinstance(spin_T, (int,float,NP.ndarray)):
        raise TypeError('Input spin_T must be scalar or numpy array')
    if not isinstance(x_H, (int,float,NP.ndarray)):
        raise TypeError('Input x_H must be scalar or numpy array')
    if not isinstance(overdensity, (int,float,NP.ndarray)):
        raise TypeError('Input overdensity must be scalar or numpy array')
    if not isinstance(redshift, (int,float,NP.ndarray)):
        raise TypeError('Input redshift must be scalar or numpy array')
        
    spin_T = spin_T.to('K')
    spin_T = NP.asarray(spin_T.value) * U.K
    x_H = NP.asarray(x_H)
    overdensity = NP.asarray(overdensity)
    redshift = NP.asarray(redshift)
    
    if spin_T.ndim == 0: 
        spin_T = spin_T.reshape(-1)
    if x_H.ndim == 0: 
        x_H = x_H.reshape(-1)
    if overdensity.ndim == 0: 
        overdensity = overdensity.reshape(-1)
    if redshift.ndim == 0: 
        redshift = redshift.reshape(-1)

    shape_Tspin = spin_T.shape
    shape_xHI = x_H.shape
    shape_overdensity = overdensity.shape
    shape_redshift = redshift.shape
    outshape_redshift = [1]*spin_T.ndim
    outshape_redshift[-1] = -1
    outshape_redshift = tuple(outshape_redshift)
    
    if (shape_Tspin != shape_xHI) or (shape_Tspin != shape_overdensity) or (shape_Tspin[-1] != redshift.size):
        raise ValueError('Dimension mismatch!')
        
    if reference.lower() == 'santos+2005':
        optical_depth = 8.6e-3 * x_H * (1+overdensity) * cosmo_in.Tcmb0 * (1+redshift.reshape(outshape_redshift)) / spin_T * (cosmo_in.Ob0 * cosmo_in.h**2 / 0.02) * NP.sqrt(0.15 / (cosmo_in.Om0 * cosmo_in.h**2)) * NP.sqrt((1+redshift.reshape(outshape_redshift))/10.0) * (0.7/cosmo_in.h)
        
    return optical_depth.decompose().value

In [8]:
def temperature_to_fluxdensity(temperature, beamSA, frequency, ref_freq=None, spindex=None):

    """
    Input temperature should be of shape (..., nlos) if spindex is not set or set to zero. 
    Input temperature should be of shape (..., 1) if spindex is set to a non-zero scalar.
    Input frequency must be of shape (nlos,)
    Input beamSA must be of shape (nbeam,nlos) or (nbeam,1) or (nbeam,) (The latter two cases imply that same beam size applies to all frequencies)
    Output fluxdensity will be of shape (nbeam,...,nlos)
    """

    temperature = U.Quantity(temperature.si.value, U.K, ndmin=1)
    beamSA = U.Quantity(beamSA.si.value, U.sr, ndmin=1)
    frequency = U.Quantity(frequency.si.value, U.Hz, ndmin=1)

    if spindex is not None:
        if not isinstance(spindex, (int,float)):
            raise TypeError('Input spindex must be a scalar')
        if (temperature.size != 1) and (NP.abs(spindex) > 1e-10):
            raise ValueError('When spindex is provided, temperature must be a scalar')
    else:
        spindex = 0.0

    inshape = temperature.shape
    if temperature.ndim == 1:
        if (temperature.shape[-1] != frequency.size) and (temperature.shape[-1] != 1): # Should be (1,) or (nlos,)
            raise ValueError('Input temperature and frequency have incompatible dimensions')
        inshape = (1,) + inshape
    temperature = temperature.reshape(-1,inshape[-1])
    
    if beamSA.ndim == 1:
        beamSA = NP.asarray(beamSA.value).reshape(-1,1) * U.sr # (nbeam,1)
    if beamSA.ndim == 2:
        if (beamSA.shape[-1] != frequency.size) and (beamSA.shape[-1] != 1): # should be (nbeam, nlos) or (nbeam,1)
            raise ValueError('Input beamSA has incompatible dimensions.')
    if beamSA.ndim > 2:
        raise ValueError('Input beamSA has incompatible dimensions')
        
    frequency = NP.asarray(frequency.value).reshape(-1) * U.Hz # (nlos,)
    freq_center = NP.mean(frequency)
    if ref_freq is None:
        ref_freq = NP.mean(frequency)
    
    wl = FC.c / frequency # (nlos,)

    spec_T = temperature * (frequency.reshape(1,-1)/ref_freq)**spindex # (npix,nlos)
    
    outshape = (beamSA.shape[0],) + inshape[:-1] + (frequency.shape[-1],)
    
    fluxdensity_spectrum = 2 * FC.k_B * spec_T[NP.newaxis,...] * beamSA[:,NP.newaxis,:].to('', equivalencies=U.dimensionless_angles()) / (wl.reshape(1,1,-1)**2) # S.(wl**2 / Omega) = 2.k.T
    return fluxdensity_spectrum.to('mJy').reshape(outshape)

In [9]:
def fluxdensity_to_temperature(fluxdensity, beamSA, frequency, ref_freq=None, spindex=None):

    """
    Input temperature should be of shape (..., nlos) if spindex is not set or set to zero. 
    Input temperature should be of shape (..., 1) if spindex is set to a non-zero scalar.
    Input frequency must be of shape (nlos,)
    Input beamSA must be of shape (nbeam,nlos) or (nbeam,1) or (nbeam,) (The latter two cases imply that same beam size applies to all frequencies)
    Output fluxdensity will be of shape (nbeam,...,nlos)
    """

    fluxdensity_unit = fluxdensity.si.unit
    fluxdensity = U.Quantity(fluxdensity.si.value, fluxdensity_unit, ndmin=1)
    beamSA = U.Quantity(beamSA.si.value, U.sr, ndmin=1)
    frequency = U.Quantity(frequency.si.value, U.Hz, ndmin=1)

    if spindex is not None:
        if not isinstance(spindex, (int,float)):
            raise TypeError('Input spindex must be a scalar')
        if (fluxdensity.size != 1) and (NP.abs(spindex) > 1e-10):
            raise ValueError('When spindex is provided, fluxdensity must be a scalar')
    else:
        spindex = 0.0

    inshape = fluxdensity.shape
    if fluxdensity.ndim == 1:
        if (fluxdensity.shape[-1] != frequency.size) and (fluxdensity.shape[-1] != 1):
            raise ValueError('Input fluxdensity and frequency have incompatible dimensions')
        inshape = (1,) + inshape
    fluxdensity = fluxdensity.reshape(-1,inshape[-1])
    
    if beamSA.ndim <= 1:
        beamSA = NP.asarray(beamSA.value).reshape(-1,1) * U.sr # (nbeam,1)
    if beamSA.ndim == 2:
        if (beamSA.shape[-1] != frequency.size) and (beamSA.shape[-1] != 1): # should be (nbeam, nlos) or (nbeam,1)
            raise ValueError('Input beamSA has incompatible dimensions.')
    if beamSA.ndim > 2:
        raise ValueError('Input beamSA has incompatible dimensions')

    frequency = NP.asarray(frequency.value).reshape(-1) * U.Hz # (nlos,)
    if ref_freq is None:
        ref_freq = NP.mean(frequency)
    
    wl = FC.c / frequency # (nlos,)

    fluxdensity_spectrum = fluxdensity * (frequency.reshape(1,-1)/ref_freq)**spindex # (npix,nlos)

    outshape = (beamSA.shape[0],) + inshape[:-1] + (frequency.shape[-1],)
    
    spec_T = fluxdensity_spectrum[NP.newaxis,...] * wl.reshape(1,1,-1)**2 / (2 * FC.k_B * beamSA[:,NP.newaxis,:].to('', equivalencies=U.dimensionless_angles())) # S.(wl**2 / Omega) = 2.k.T
    return spec_T.to('K').reshape(outshape)

In [10]:
def classical_source_confusion(theta_FWHM, frequency):

    """
    Reference: Condon et al. (2012)
    theta_FWHM: (ntheta,) with units of angle
    frequency: (nfreq,) with units of frequency
    
    Output: 1-sigma confusion assuming S/N=5 for source subtraction threshold
    """
    
    confusion_rms = 1.2 * U.uJy * (frequency/(3.02*U.GHz))**(-0.7) * (theta_FWHM/(8.0*U.arcsec))**(10/3.)
    return confusion_rms.to('uJy')

Start the code for computing and plotting

Read input parameters


In [11]:
parmsfile = '/lustre/aoc/users/nthyagar/codes/mine/python/projects/21cmforest/stats_analysis_theory_parms.yaml'
with open(parmsfile, 'r') as pfile:
    parms = yaml.safe_load(pfile)
    
projectdir = parms['dirstruct']['projectdir']
figdir = projectdir + parms['dirstruct']['figdir']
print(parms)


{'output': {'nzbins': 1, 'rest_frequency': 1420405751.77}, 'dirstruct': {'figdir': 'figures/', 'projectdir': '/lustre/aoc/users/nthyagar/projects/21cmforest/'}, 'sim': {'cosmo': {'Tcmb0': 2.7255, 'name': 'custom', 'h': 0.678, 'Ode0': None, 'Om0': 0.308, 'sigma8': 0.815, 'w0': -1, 'Y_He': 0.245, 'Ob0': 0.02226}, 'source': '21cmfast', 'units_in': {'boxdim': 'Mpc', 'boxval': 'mK'}}, 'actions': {'process_LC_HDF5': {'tabledir': 'tables/', 'cubefile': 'light_cone.hdf5', 'act': False, 'models': ['Faint_galaxies_fiducial_1024', 'Bright_galaxies_fiducial_1024', 'Faint_galaxies_coldIGM'], 'combined_dir': 'lightcone/combined/', 'modeldir': '/lustre/aoc/users/nthyagar/data/EoR_models/21cmFAST/Andrei_Mesinger/'}, 'process_raw_LC': {'input': {'qtyinfo': {'delta_T': {'infile_suffix': 'lighttravel', 'parseinfo': {'delimiter': '_', 'zstart_identifier_pos': 'before', 'zstart_identifier': 'zstart', 'zend_pos': 6, 'zstart_pos': 5, 'zend_identifier': 'zend', 'zend_identifier_pos': 'before'}, 'units': 'mK', 'indir': '/lustre/aoc/users/nthyagar/data/EoR_models/21cmFAST/Andrei_Mesinger/Bright_galaxies_fiducial_1024/lightcone/', 'infile_prefix': 'delta_T_v3_no_halos_'}, 'T_s': {'infile_suffix': 'lighttravel', 'parseinfo': {'delimiter': '_', 'zstart_identifier_pos': 'before', 'zstart_identifier': 'zstart', 'zend_pos': 2, 'zstart_pos': 1, 'zend_identifier': 'zend', 'zend_identifier_pos': 'before'}, 'units': 'K', 'indir': '/lustre/aoc/users/nthyagar/data/EoR_models/21cmFAST/Andrei_Mesinger/Bright_galaxies_fiducial_1024/lightcone/', 'infile_prefix': 'Ts_'}, 'x_H': {'infile_suffix': 'lighttravel', 'parseinfo': {'delimiter': '_', 'zstart_identifier_pos': 'before', 'zstart_identifier': 'zstart', 'zend_pos': 3, 'zstart_pos': 2, 'zend_identifier': 'zend', 'zend_identifier_pos': 'before'}, 'units': '', 'indir': '/lustre/aoc/users/nthyagar/data/EoR_models/21cmFAST/Andrei_Mesinger/Bright_galaxies_fiducial_1024/lightcone/', 'infile_prefix': 'xH_nohalos_'}, 'density': {'infile_suffix': 'lighttravel', 'parseinfo': {'delimiter': '_', 'zstart_identifier_pos': 'before', 'zstart_identifier': 'zstart', 'zend_pos': 4, 'zstart_pos': 3, 'zend_identifier': 'zend', 'zend_identifier_pos': 'before'}, 'units': '', 'indir': '/lustre/aoc/users/nthyagar/data/EoR_models/21cmFAST/Andrei_Mesinger/Bright_galaxies_fiducial_1024/lightcone/', 'infile_prefix': 'updated_smoothed_deltax_'}}}, 'cube_source': '21cmfast', 'qtys': ['delta_T', 'x_H', 'density', 'T_s'], 'act': False, 'output': {'outdir': '/lustre/aoc/users/nthyagar/data/EoR_models/21cmFAST/Andrei_Mesinger/Bright_galaxies_fiducial_1024/lightcone/combined/', 'outfile_prefix': None}}}, 'plot': {'1': {'models': ['Faint_galaxies_fiducial_1024', 'Bright_galaxies_fiducial_1024', 'Faint_galaxies_coldIGM'], 'Tsys': 500.0, 'frac_bw': 0.0667, 'seed': 100, 'z_qso': [10.0, 13.0, 21.0], 'data_loss': 0.0, 't_int': 60.0, 'spindex': -0.7, 'beam_majax': [25.0], 'rest_freq': 1420405751.77, 'n_antennas': 30, 'subband': {'fftpow': 2.0, 'shape': 'bhw', 'pad': 1.0}, 'flux_obs': [0.01, 0.1], 'snr_min': 10.0, 'action': False, 'z_pspec': [8.25, 11.5, 18.0], 'freq_obs': 150000000.0, 'nrand': 10000}, '1a': {'action': False}, '3': {'n_qso_ref': 100, 'A_over_T_lofar': 90.0, 'A_over_T_SKA2': 4000.0, 'A_over_T_SKA1': 800.0, 'Tsys_reffreq': 150000000, 'seed': 100, 'z_qso': [10.0, 12.0, 18.0], 'data_loss': 0.0, 'Tsys0': 384.0, 'A_over_T_ref': 800.0, 'snr_min': 5.0, 'rest_freq': 1420405751.77, 'A_over_T_EGMRT': 210.0, 'subband': {'fftpow': 1.0, 'shape': 'bhw', 'pad': 1.0}, 'flux_obs': [1.0, 10.0, 100.0], 'indir_pfx': '/lustre/aoc/users/nthyagar/data/EoR_models/21cmFAST/Andrei_Mesinger/', 'freq_obs': 150000000.0, 'indir_sfx': '/lightcone/combined/', 'models': ['Bright_galaxies_fiducial_1024', 'Faint_galaxies_fiducial_1024'], 'frac_bw': 0.0667, 'efficiency': 1.0, 'qtys': ['delta_T', 'x_H', 'density', 'T_s'], 't_int': 36000.0, 'spindex': -1.05, 'fov_fwhm': 5.0, 'beam_majax': [10.0, 10.0, 10.0], 'Tsys_spindex': -2.55, 'n_antennas': 30, 'action': True, 'z_pspec': [8.0, 9.5, 11.0], 'nrand': 10000, 'A_over_T_uGMRT': 70.0}, '2': {'indir_sfx': '/lightcone/combined/', 'models': ['Bright_galaxies_fiducial_1024', 'Faint_galaxies_fiducial_1024'], 'Tsys': 500.0, 'frac_bw': 0.0667, 'seed': 100, 'z_qso': [10.0, 13.0, 21.0], 'data_loss': 0.0, 'n_antennas': 30, 't_int': 60.0, 'spindex': -0.7, 'fov_fwhm': 180.0, 'snr_min': 10.0, 'rest_freq': 1420405751.77, 'qtys': ['delta_T', 'T_s'], 'subband': {'fftpow': 1.0, 'shape': 'bhw', 'pad': 1.0}, 'flux_obs': [0.01, 0.1], 'beam_majax': [25.0], 'action': False, 'z_pspec': [8.25, 11.5, 18.0], 'indir_pfx': '/lustre/aoc/users/nthyagar/data/EoR_models/21cmFAST/Andrei_Mesinger/', 'freq_obs': 150000000.0, 'nrand': 10000}, '2a': {'action': False}, '2b': {'action': False}, '2c': {'action': False}, '1b': {'action': False}}}

Read cosmological parameters and units


In [12]:
cosmoparms = parms['sim']['cosmo']
if cosmoparms['name'] is None:
    cosmo_in = None
elif cosmoparms['name'].lower() == 'custom':
    h = cosmoparms['h']
    H0 = 100.0 * h
    Om0 = cosmoparms['Om0']
    Ode0 = cosmoparms['Ode0']
    if Ode0 is None:
        Ode0 = 1.0 - Om0
    Ob0 = cosmoparms['Ob0'] / h**2
    w0 = cosmoparms['w0']
    Tcmb0 = cosmoparms['Tcmb0']
    cosmo_in = cosmology.wCDM(H0, Om0, Ode0, w0=w0, Ob0=Ob0, Tcmb0=Tcmb0*U.K)
elif cosmoparms['name'].lower() == 'wmap9':
    cosmo_in = cosmology.WMAP9
elif cosmoparms['name'].lower() == 'planck15':
    cosmo_in = cosmology.Planck15
else:
    raise ValueError('{0} preset not currently accepted for cosmology'.format(cosmoparms['name'].lower()))
# cosmo100 = cosmo.clone(name='Modified cosmology with h=1.0', H0=100.0) # Modified cosmology with h=1.0, H=100 km/s/Mpc

units_in = parms['sim']['units_in']
units_in['frequency'] = 'Hz'

In [13]:
print(cosmo_in)


wCDM(H0=67.8 km / (Mpc s), Om0=0.308, Ode0=0.692, w0=-1, Tcmb0=2.725 K, Neff=3.04, m_nu=[0. 0. 0.] eV, Ob0=0.0484)

Read info on cube operations to be performed


In [14]:
actions = [key for key in parms['actions'] if parms['actions'][key]['act']]

Convert raw lightcone to HDF5 format if required


In [11]:
if 'process_raw_LC' in actions:
    qtys = parms['actions']['process_raw_LC']['qtys']
    for qty in qtys:
        lightcone_cube, boxres, zmin = read_raw_lightcone_cube(parms['actions']['process_raw_LC']['input']['qtyinfo'][qty])
        outparms = parms['actions']['process_raw_LC']['output']
        outdir = outparms['outdir']
        outfile_prefix = outparms['outfile_prefix']
        if outfile_prefix is None:
            outfile = outdir + '{0}_light_cone'.format(qty)
        elif isinstance(outfile_prefix, str):
            outfile = outdir + outfile_prefix + '_light_cone'
        else:
            raise TypeError('Output filename prefix must be set to None or a string')
        write_raw_lightcone_cube(lightcone_cube, boxres, zmin, outfile+'.hdf5', units=parms['actions']['process_raw_LC']['input']['qtyinfo'][qty]['units'], cosmo_in=cosmo_in)

Work with the HDF5 lightcone


In [12]:
if 'process_LC_HDF5' in actions:
    models = parms['actions']['process_LC_HDF5']['models']
    modeldir = parms['actions']['process_LC_HDF5']['modeldir']
    combined_dir = parms['actions']['process_LC_HDF5']['combined_dir']
    cubedirs = [modeldir + model + '/' + combined_dir for model in models]
    cubefile = parms['actions']['process_LC_HDF5']['cubefile']
    tabledirs = [modeldir + model + '/' + parms['actions']['process_LC_HDF5']['tabledir'] for model in models]
    print(tabledirs)
    cube_delta_T = {}
    cube_xH = {}
    cube_density = {}
    cube_Tspins = {}
    lccubes = {}
    cube_monopole_values = {}
    boxresolutions = {}
    cube_zmins = {}
    cosmo_parms = {}
    cube_d_los_comov = {}
    cube_redshifts = {}
    cube_mean_Tspins = {}
    tab_redshifts = {}
    tab_Tspins = {}
    cube_solang = {}
    for modelind, model in enumerate(models):
        cubefilepath = cubedirs[modelind] + cubefile
        lightcone_cube, monopole_values, boxres, zmin, cosmo = read_HDF5_lightcone_cube(cubefilepath)
        d_los_comov_min = cosmo.comoving_distance(zmin)
        cube_d_los_comov[model] = d_los_comov_min + boxres * NP.arange(lightcone_cube.shape[0])
        redshifts = NP.asarray([cosmology.z_at_value(cosmo.comoving_distance, d_los) for d_los in cube_d_los_comov[model]])

        tabfile = glob.glob(tabledirs[modelind] + '*.txt')[0]
        tabinfo = ascii.read(tabfile)
        print(model+'\n-------------------------')
        print(tabinfo)
        print('\n')
        tab_redshifts[model] = tabinfo['col1'].data
        tab_Tspins[model] = tabinfo['col5'].data * U.K

        zind = NP.where(NP.logical_and(redshifts >= tab_redshifts[model].min(), redshifts <= tab_redshifts[model].max()))[0]
        cube_redshifts[model] = redshifts[zind]
        cube_d_los_comov[model] = cube_d_los_comov[model][zind]
        lccubes[model] = lightcone_cube[zind,:,:]
        cube_monopole_values[model] = monopole_values[zind]
        boxresolutions[model] = boxres
        pixang = boxres / cosmo.comoving_distance(cube_redshifts[model]) * U.rad
        cube_solang[model] = pixang**2
        # cube_solang[model] = (boxres * cosmo.arcsec_per_kpc_comoving(cube_redshifts[model]))**2
        cube_solang[model] = cube_solang[model].to('sr')
        cosmo_parms[model] = copy.deepcopy(cosmo)
        cube_zmins[model] = redshifts[zind].min()

        cube_mean_Tspins[model] = OPS.interpolate_array(tab_Tspins[model].value, tab_redshifts[model], cube_redshifts[model], kind='linear', axis=-1) * tab_Tspins[model].unit

Prepare plot information


In [16]:
plot_info = parms['plot']
plots = [key for key in plot_info if plot_info[key]['action']]

In [17]:
print(plots)
print(plot_info['3'])


['3']
{'n_qso_ref': 100, 'A_over_T_lofar': 90.0, 'A_over_T_SKA2': 4000.0, 'A_over_T_SKA1': 800.0, 'Tsys_reffreq': 150000000, 'seed': 100, 'z_qso': [10.0, 12.0, 18.0], 'data_loss': 0.0, 'Tsys0': 384.0, 'A_over_T_ref': 800.0, 'snr_min': 5.0, 'rest_freq': 1420405751.77, 'A_over_T_EGMRT': 210.0, 'subband': {'fftpow': 1.0, 'shape': 'bhw', 'pad': 1.0}, 'flux_obs': [1.0, 10.0, 100.0], 'indir_pfx': '/lustre/aoc/users/nthyagar/data/EoR_models/21cmFAST/Andrei_Mesinger/', 'freq_obs': 150000000.0, 'indir_sfx': '/lightcone/combined/', 'models': ['Bright_galaxies_fiducial_1024', 'Faint_galaxies_fiducial_1024'], 'frac_bw': 0.0667, 'efficiency': 1.0, 'qtys': ['delta_T', 'x_H', 'density', 'T_s'], 't_int': 36000.0, 'spindex': -1.05, 'fov_fwhm': 5.0, 'beam_majax': [10.0, 10.0, 10.0], 'Tsys_spindex': -2.55, 'n_antennas': 30, 'action': True, 'z_pspec': [8.0, 9.5, 11.0], 'nrand': 10000, 'A_over_T_uGMRT': 70.0}

** Go down to Plot 3, skip plots 1 and 2 ***

Plot 1: Read the information for processing the lightcone cubes for different tasks


In [324]:
if ('1' in plots) or ('1a' in plots) or ('1b' in plots) or ('1c' in plots):

    modelnames = plot_info['1']['models']
    beam_majax = NP.asarray(plot_info['1']['beam_majax']) * U.arcsec
    beam_sa_0 = (beam_majax**2).to('sr') # Beam solid angle at z=0
    flux_density_obs = NP.asarray(plot_info['1']['flux_obs']) * U.mJy
    beam_majax = beam_majax * NP.ones(flux_density_obs.size)
    beam_sa_0 = beam_sa_0 * NP.ones(flux_density_obs.size)
    flux_obs = (flux_density_obs / beam_sa_0).reshape(-1,1)
    freq_obs = plot_info['1']['freq_obs'] * U.Hz
    spindex = plot_info['1']['spindex']
    z_qso = NP.asarray(plot_info['1']['z_qso'])
    z_pspec = NP.asarray(plot_info['1']['z_pspec'])
    frac_bw = plot_info['1']['frac_bw']
    rest_freq = plot_info['1']['rest_freq'] * U.Hz
    if rest_freq is None:
        rest_freq = CNST.rest_freq_HI * U.Hz
    nrand = plot_info['1']['nrand']
    randseed = plot_info['1']['seed']
    subbandinfo = plot_info['1']['subband']
    winshape = subbandinfo['shape']
    fftpow = subbandinfo['fftpow']
    pad = subbandinfo['pad']
    Tsys = plot_info['1']['Tsys'] * U.K
    n_ant = plot_info['1']['n_antennas']
    t_int = plot_info['1']['t_int'] * U.s
    t_half = 0.5 * t_int
    detection_threshold = plot_info['1']['snr_min']
    data_loss_factor = plot_info['1']['data_loss']

In [325]:
print('Beam major axis: {0}'.format(beam_majax))
    print('Beam solid angle: {0}'.format(beam_sa_0))
    print('Flux Density observed at {0}: {1}'.format(freq_obs.to('MHz'), flux_density_obs))
    print('Flux observed at {0}: {1}'.format(freq_obs.to('MHz'), flux_obs.to('mJy arcmin-2')))
    print('Spectral Index: {0}'.format(spindex))
    print('System Temperature: {0}'.format(Tsys))
    print('Number of antennas: {0}'.format(n_ant))
    print('Integration time: {0}'.format(t_int))
    print('Detection threshold: {0}'.format(detection_threshold))
    print('Data loss factor: {0}'.format(data_loss_factor))


Beam major axis: [5. 5. 5. 5.] arcsec
Beam solid angle: [4.61508414e-10 4.61508414e-10 4.61508414e-10 4.61508414e-10] sr
Flux Density observed at 150.0 MHz: [  0.1   1.   10.  100. ] mJy
Flux observed at 150.0 MHz: [[   18.33464944]
 [  183.34649444]
 [ 1833.46494442]
 [18334.64944419]] mJy / arcmin2
Spectral Index: -1.05
System Temperature: 500.0 K
Number of antennas: 30
Integration time: 3600.0 s
Detection threshold: 10.0
Data loss factor: 0.0

Initialize variables


In [714]:
d_los_comov_center = cosmo.comoving_distance(z_pspec)
    zlo = z_pspec * (1 - 0.5 * frac_bw)
    zhi = z_pspec * (1 + 0.5 * frac_bw)
    d_los_comov_lo = cosmo.comoving_distance(zlo)
    d_los_comov_hi = cosmo.comoving_distance(zhi)
    d_los_comov_width = d_los_comov_hi - d_los_comov_lo
    
    print('Redshifts of background quasars = {0}'.format(z_qso))
    print('Redshifts where HI PS is estimated = {0}'.format(z_pspec))
    print('Center of LOS comoving distance to PS redshift: {0}/h'.format(d_los_comov_center))
    print('Width of LOS comoving distance at PS redshift: {0}/h'.format(d_los_comov_width))


Redshifts of background quasars = [10. 13. 21.]
Redshifts where HI PS is estimated = [ 8.28 11.3  18.  ]
Center of LOS comoving distance to PS redshift: [6256.07728224 6721.29807069 7322.42661988] Mpc/h
Width of LOS comoving distance at PS redshift: [105.37539567  94.30284466  78.24485193] Mpc/h

In [715]:
net_bg_temp_los = {}
    bg_temp_los = {}
    Tb_los_fluctuations = {}
    Tb_los_fluctuations_modified = {}
    modification_factor = {}
    pspec_kprll_onesided = {}
    noise_power_std_kprll_onesided = {}
    nmodes_onesided = {}
    nmodes_total_onesided = {}
    wlinv_kernel = {}
    kprll_onesided = {}
    nsamples_onesided = {}
    nsamples_required = {}
    nsamples_total_onesided = {}
    basic_sampstd_kprll_onesided = {}
    basic_uncertainty_kprll_onesided = {}

In [716]:
for model in modelnames:
        # flux_los = flux_obs * (rest_freq/freq_obs)**spindex / (1+cube_redshifts[model])**(1+spindex)
        # dtheta = (boxresolutions[model] * cosmo.arcsec_per_kpc_comoving(cube_redshifts[model]))
        # dOmega = dtheta**2
        # sa_filling_factor = (beam_majax / (boxresolutions[model]/))**2
        sa_filling_factor_pix_beam = cube_solang[model].reshape(1,-1) / beam_sa_0.reshape(-1,1)
        sa_filling_factor_pix_beam[sa_filling_factor_pix_beam.decompose() <= 1.0] = 1.0
        flux_los = flux_obs * (rest_freq/freq_obs)**spindex * (1+cube_redshifts[model].reshape(1,-1))**(3.0-spindex) / sa_filling_factor_pix_beam # Correction for filling factor between pixel and beam    
        bg_temp_los[model] = flux_los * (FC.c / rest_freq)**2 / (2 * FC.k_B)
        bg_temp_los[model] = bg_temp_los[model].si.value * U.K
#         sa_filling_factor = beam_sa_0.reshape(-1,1) / cube_solang[model].reshape(1,-1)
#         sa_filling_factor[sa_filling_factor.decompose() >= 1.0] = 1.0
#         flux_los = flux_obs * (rest_freq/freq_obs)**spindex * (1+cube_redshifts[model].reshape(1,-1))**(3.0-spindex) * sa_filling_factor 
#         bg_temp_los[model] = flux_los * (FC.c / rest_freq)**2 / (2 * FC.k_B)

        los_ind_to_modify = cube_redshifts[model].reshape(1,-1) <= z_qso.reshape(-1,1) # nz_qso x nlos
        net_bg_temp_los[model] = bg_temp_los[model][:,NP.newaxis,:] * los_ind_to_modify[NP.newaxis,:,:].astype(NP.float) + cosmo_parms[model].Tcmb0 * (1 + cube_redshifts[model].reshape(1,1,-1)) # nsa x nz_qso x nlos

        rstate = NP.random.RandomState(seed=randseed)
        randpix = rstate.randint(lccubes[model].shape[1]*lccubes[model].shape[2], size=nrand)
        Tb_los_fluctuations[model] = lccubes[model].reshape(cube_redshifts[model].size, -1)[:,randpix] - cube_monopole_values[model].reshape(-1,1) # nlos x nxy
        Tb_los_fluctuations[model] = Tb_los_fluctuations[model].to('K').T # nxy x nlos
        modification_factor[model] = (cube_mean_Tspins[model].reshape(1,1,-1) - net_bg_temp_los[model][:,:,:]) / (cube_mean_Tspins[model].reshape(1,1,-1) - cosmo_parms[model].Tcmb0 * (1+cube_redshifts[model].reshape(1,1,-1))) # nsa x nz_qso x nlos
        Tb_los_fluctuations_modified[model] = Tb_los_fluctuations[model][NP.newaxis,NP.newaxis,:,:] * modification_factor[model][:,:,NP.newaxis,:] # nsa x nz_qso x nxy x nlos
        Tb_los_fluctuations_modified[model] = Tb_los_fluctuations_modified[model].to('K')
        
        print(model+'\n-----------------------')
        print('Pix/Beam filling fraction: {0}\n'.format(sa_filling_factor_pix_beam.squeeze()))
        print('Specific Intensity in a cube pixel along LOS: {0}\n'.format(flux_los.to('mJy arcmin-2').squeeze()))
        print('Brightness temperature from background quasar as seen by each pixel at LOS:\n{0}\n'.format(bg_temp_los[model].to('K').squeeze()))
        print('Modification factor in background radiation temperature: {0}\n'.format(modification_factor[model].decompose()))
        print('Maximum modification factor in background radiation temperature: {0}\n'.format(NP.amax(modification_factor[model].decompose(), axis=-1).squeeze()))


Faint_galaxies_fiducial_1024
-----------------------
Pix/Beam filling fraction: [[2.62376559 2.62273566 2.62170633 ... 1.19718551 1.19686804 1.1965507 ]
 [2.62376559 2.62273566 2.62170633 ... 1.19718551 1.19686804 1.1965507 ]]

Specific Intensity in a cube pixel along LOS: [[3.45150551e+00 3.45903996e+00 3.46659201e+00 ... 5.46859578e+03
  5.49383491e+03 5.51920445e+03]
 [3.45150551e+01 3.45903996e+01 3.46659201e+01 ... 5.46859578e+04
  5.49383491e+04 5.51920445e+04]] mJy / arcmin2

Brightness temperature from background quasar as seen by each pixel at LOS:
[[6.58051329e-01 6.59487820e-01 6.60927666e-01 ... 1.04262233e+03
  1.04743433e+03 1.05227119e+03]
 [6.58051329e+00 6.59487820e+00 6.60927666e+00 ... 1.04262233e+04
  1.04743433e+04 1.05227119e+04]] K

Modification factor in background radiation temperature: [[[0.99953985 0.9995392  0.99953855 ... 1.         1.         1.        ]
  [0.99953985 0.9995392  0.99953855 ... 1.         1.         1.        ]
  [0.99953985 0.9995392  0.99953855 ... 1.         1.         1.        ]]

 [[0.9953985  0.99539199 0.99538547 ... 1.         1.         1.        ]
  [0.9953985  0.99539199 0.99538547 ... 1.         1.         1.        ]
  [0.9953985  0.99539199 0.99538547 ... 1.         1.         1.        ]]]

Maximum modification factor in background radiation temperature: [[1.00000000e+00 1.11004819e+02 1.11004819e+02]
 [1.00000000e+00 1.10104819e+03 1.10104819e+03]]

Bright_galaxies_fiducial_1024
-----------------------
Pix/Beam filling fraction: [[2.62376559 2.62273566 2.62170633 ... 1.19718551 1.19686804 1.1965507 ]
 [2.62376559 2.62273566 2.62170633 ... 1.19718551 1.19686804 1.1965507 ]]

Specific Intensity in a cube pixel along LOS: [[3.45150551e+00 3.45903996e+00 3.46659201e+00 ... 5.46859578e+03
  5.49383491e+03 5.51920445e+03]
 [3.45150551e+01 3.45903996e+01 3.46659201e+01 ... 5.46859578e+04
  5.49383491e+04 5.51920445e+04]] mJy / arcmin2

Brightness temperature from background quasar as seen by each pixel at LOS:
[[6.58051329e-01 6.59487820e-01 6.60927666e-01 ... 1.04262233e+03
  1.04743433e+03 1.05227119e+03]
 [6.58051329e+00 6.59487820e+00 6.60927666e+00 ... 1.04262233e+04
  1.04743433e+04 1.05227119e+04]] K

Modification factor in background radiation temperature: [[[0.99903752 0.99903526 0.99903299 ... 1.         1.         1.        ]
  [0.99903752 0.99903526 0.99903299 ... 1.         1.         1.        ]
  [0.99903752 0.99903526 0.99903299 ... 1.         1.         1.        ]]

 [[0.9903752  0.99035258 0.9903299  ... 1.         1.         1.        ]
  [0.9903752  0.99035258 0.9903299  ... 1.         1.         1.        ]
  [0.9903752  0.99035258 0.9903299  ... 1.         1.         1.        ]]]

Maximum modification factor in background radiation temperature: [[ 480.08721277  480.08721277  480.08721277]
 [4791.8721277  4791.8721277  4791.8721277 ]]

Faint_galaxies_coldIGM
-----------------------
Pix/Beam filling fraction: [[3.64587014 3.64418321 3.64249744 ... 1.9870684  1.98638957 1.9857111 ]
 [3.64587014 3.64418321 3.64249744 ... 1.9870684  1.98638957 1.9857111 ]]

Specific Intensity in a cube pixel along LOS: [[4.39648368e+00 4.40913152e+00 4.42181869e+00 ... 1.32088121e+03
  1.32767941e+03 1.33451663e+03]
 [4.39648368e+01 4.40913152e+01 4.42181869e+01 ... 1.32088121e+04
  1.32767941e+04 1.33451663e+04]] mJy / arcmin2

Brightness temperature from background quasar as seen by each pixel at LOS:
[[8.38217387e-01 8.40628778e-01 8.43047666e-01 ... 2.51834347e+02
  2.53130466e+02 2.54434026e+02]
 [8.38217387e+00 8.40628778e+00 8.43047666e+00 ... 2.51834347e+03
  2.53130466e+03 2.54434026e+03]] K

Modification factor in background radiation temperature: [[[0.99948356 0.99948185 0.99948013 ... 1.         1.         1.        ]
  [0.99948356 0.99948185 0.99948013 ... 1.         1.         1.        ]
  [0.99948356 0.99948185 0.99948013 ... 1.         1.         1.        ]]

 [[0.99483559 0.99481846 0.99480126 ... 1.         1.         1.        ]
  [0.99483559 0.99481846 0.99480126 ... 1.         1.         1.        ]
  [0.99483559 0.99481846 0.99480126 ... 1.         1.         1.        ]]]

Maximum modification factor in background radiation temperature: [[1.00000000e+00 5.18429009e+02 5.18429009e+02]
 [1.00000000e+00 5.17529009e+03 5.17529009e+03]]

Plot 1a: Mean $T_\mathrm{spin}$ and background radiation temperatures, $T_\mathrm{\gamma}$ from the QSO and the CMB as seen by the observer at some measured flux density when extrapolated back in redshift


In [735]:
if '1a' in plots:
        bgrad_colrs = ['blue', '0.8', 'black']
        zqso_ls = ['-.', '--', ':']
        ncols = beam_majax.size
        for mdlind, mdl in enumerate(modelnames):
            fig, axs = PLT.subplots(ncols=ncols, sharey=True, figsize=(ncols*3.5, 4.5), squeeze=False)
            if axs.size == 1:
                axs = axs.reshape(-1)
            for beamind, beamsize in enumerate(beam_majax):
                line_spin, = axs[0,beamind].plot(1+cube_redshifts[mdl], cube_mean_Tspins[mdl], ls='-', lw=2, color='cyan', label=r'$T_\mathrm{spin}$')
                line_cmb, = axs[0,beamind].plot(1+cube_redshifts[mdl], cosmo_parms[mdl].Tcmb0 * (1+cube_redshifts[mdl]), ls='-', lw=3, color=bgrad_colrs[0], label=r'$T_\mathrm{bg}=T_\mathrm{cmb}$')
                line_qso, = axs[0,beamind].plot(1+cube_redshifts[mdl], bg_temp_los[mdl][beamind,:], ls='-', lw=4, color=bgrad_colrs[1], label=r'$T_\mathrm{bg}=T_\mathrm{qso}$')
                line_cmb_qso = []
                for zqi, zq in enumerate(z_qso):
                    line_cmb_qso += axs[0,beamind].plot(1+cube_redshifts[mdl], net_bg_temp_los[mdl][beamind,zqi,:], ls=zqso_ls[zqi], lw=1.5, color=bgrad_colrs[2], label=r'$T_\mathrm{bg}=T_\mathrm{qso}(z_\mathrm{qso}=$'+'{0:.1f})'.format(zq) + r'$+T_\mathrm{cmb}$')
                line_cmb_qso = tuple(line_cmb_qso)
                if beamind == 0:
                    if mdl == 'Faint_galaxies_coldIGM':
                        first_legend = axs[0,beamind].legend(handles=(line_spin, line_qso, line_cmb), loc='upper left', shadow=False, fontsize=7, labelspacing=0)
                    else:
                        first_legend = axs[0,beamind].legend(handles=(line_spin, line_qso, line_cmb), loc='lower left', shadow=False, fontsize=7, labelspacing=0)
                    second_legend = axs[0,beamind].legend(handles=line_cmb_qso, loc='upper right', shadow=False, fontsize=7, labelspacing=0)
                    axs[0,beamind].add_artist(first_legend)
                    axs[0,beamind].add_artist(second_legend)

                axs[0,beamind].text(0.5, 0.75, r'$S_\mathrm{qso}=$'+'{0:.2g} mJy'.format((flux_obs[beamind]*beam_sa_0[beamind]).to('mJy').value[0]), transform=axs[0,beamind].transAxes, ha='center', va='top', fontsize=12)
                axs[0,beamind].set_ylim(ymin=1, ymax=5e4)
                axs[0,beamind].set_xscale('log')
                axs[0,beamind].set_xlabel(r'$1+z$', fontsize=12, weight='medium')
                axs[0,beamind].set_yscale('log')
#                 axs[0,beamind].set_xlim(32.5, 6.5)
                tickformatter = ticker.FuncFormatter(lambda y, _: '{:.16g}'.format(y))
                axs[0,beamind].xaxis.set_major_formatter(tickformatter)
                axs[0,beamind].xaxis.set_minor_formatter(tickformatter)
                axs[0,beamind].tick_params(axis='both', which='both', labelsize=10)
                axs[0,beamind].invert_xaxis()
 
            fig.subplots_adjust(hspace=0, wspace=0)
            fig.subplots_adjust(top=0.98, bottom=0.1, left=0.1, right=0.98)

            big_ax = fig.add_subplot(111)
            big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
            # big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
            big_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
            big_ax.set_xticks([])
            big_ax.set_yticks([])
#             big_ax.set_xlabel(r'$1+z$', fontsize=12, weight='medium', labelpad=20)
            big_ax.set_ylabel(r'$T$ [K]', fontsize=12, weight='medium', labelpad=25)

            PLT.savefig(figdir+'{0}_background_brightness_temperature.pdf'.format(mdl), bbox_inches=0)


Plot 1b: Modification factors to brightness temperature fluctuations


In [271]:
if '1b' in plots:        
        zqso_ls = ['-', '--', ':']
        zqso_modification_colrs = ['0.85', '0.7', '0.55']
        shade_maxyval_zqso = [2.5e3, 5e3, 1e4]
        ncols = beam_majax.size
        mdl_simplified_names = ['Bright Gal.', 'Faint Gal.', 'Cold IGM']
        for mdlind, mdl in enumerate(modelnames):
            fig, axs = PLT.subplots(nrows=2, ncols=ncols, sharex=True, figsize=(ncols*3.5, 2*2.5), squeeze=False)
#             if axs.size == 1:
#                 axs = axs.reshape(-1)
            for beamind, beamsize in enumerate(beam_majax):
                increase_modification_ind_list = []
                for zqi, zq in enumerate(z_qso):
                    ind_increase_modification = NP.where(NP.abs(modification_factor[mdl][beamind,zqi,:]) > 1.0)[0]
                    increase_modification_ind_list += [ind_increase_modification]

                axs[0,beamind].plot(1+cube_redshifts[mdl], cube_mean_Tspins[mdl] - cosmo_parms[mdl].Tcmb0 * (1+cube_redshifts[mdl]), ls='-', lw=4, color='cyan', label=r'$T_\mathrm{s} - T_\mathrm{cmb}$')
                for zqi, zq in reversed(list(enumerate(z_qso))):
                    axs[0,beamind].vlines(1+cube_redshifts[mdl][increase_modification_ind_list[zqi]], -shade_maxyval_zqso[zqi], shade_maxyval_zqso[zqi], color=zqso_modification_colrs[zqi], alpha=0.5)
                    axs[0,beamind].plot(1+cube_redshifts[mdl], cube_mean_Tspins[mdl] - net_bg_temp_los[mdl][beamind,zqi,:], ls=zqso_ls[zqi], lw=2, color='black', label=r'$T_\mathrm{s} - [T_\mathrm{qso}(z=$'+'{0:.1f}'.format(zq)+r'$)+T_\mathrm{cmb}]$')
                    
                if beamind == 0:
                    axs[0,beamind].legend(loc='upper left', shadow=False, fontsize=7, labelspacing=0)
                                    
                axs[0,beamind].set_yscale('symlog', linthreshy=10)
                axs[0,beamind].set_xscale('log')
                axs[0,beamind].set_ylabel(r'$T_\mathrm{s} - T_\mathrm{bg}$ [K]', fontsize=12, weight='medium')
                
                for zqi, zq in reversed(list(enumerate(z_qso))):    
                    axs[1,beamind].vlines(1+cube_redshifts[mdl][increase_modification_ind_list[zqi]], -shade_maxyval_zqso[zqi], shade_maxyval_zqso[zqi], color=zqso_modification_colrs[zqi], alpha=0.5)
                    axs[1,beamind].plot(1+cube_redshifts[mdl], modification_factor[mdl][beamind,zqi,:], ls=zqso_ls[zqi], lw=2, color='black', label=r'CMB + $z_\mathrm{qso}=$'+'{0:.1f}'.format(zq))
                axs[1,beamind].set_yscale('symlog', linthreshy=1)
                axs[1,beamind].set_xscale('log')
                axs[1,beamind].set_ylabel('Modification factor', fontsize=12, weight='medium')
                
                tickformatter = ticker.FuncFormatter(lambda y, _: '{:.16g}'.format(y))
                axs[1,beamind].xaxis.set_major_formatter(tickformatter)
                axs[1,beamind].xaxis.set_minor_formatter(tickformatter)
                
                axs[0,beamind].set_xlim(40,6)
                axs[1,beamind].set_xlim(40,6)

                fig.subplots_adjust(hspace=0, wspace=0)
                
                fig.subplots_adjust(top=0.98)
                fig.subplots_adjust(bottom=0.11)
                fig.subplots_adjust(left=0.22)
                fig.subplots_adjust(right=0.98)

                big_ax = fig.add_subplot(111)
                big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
                # big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
                big_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
                big_ax.set_xticks([])
                big_ax.set_yticks([])
                big_ax.set_xlabel(r'$1+z$', fontsize=12, weight='medium', labelpad=20)

            PLT.savefig(figdir+'{0}_background_brightness_temperature_modification.pdf'.format(mdl), bbox_inches=0)



In [18]:
for model in modelnames:
        d_los_comov_wts = NP.empty((z_pspec.size, cube_d_los_comov[model].size), dtype=NP.float_) # nz x nlos
        frac_width = DSP.window_N2width(n_window=None, shape=winshape, fftpow=fftpow, area_normalize=False, power_normalize=True)
        window_loss_factor = 1 / frac_width
        n_window = NP.round(window_loss_factor * (d_los_comov_width/boxresolutions[model]).decompose().value).astype(NP.int)
        d_los_comov_center_ind, ind_los_pixels, d_los_comov_offset = LKP.find_1NN(cube_d_los_comov[model].to('Mpc').value.reshape(-1,1), d_los_comov_center.to('Mpc').value.reshape(-1,1), distance_ULIM=0.51*boxresolutions[model].to('Mpc').value, remove_oob=True)
        sortind = NP.argsort(ind_los_pixels)
        d_los_comov_center_ind = d_los_comov_center_ind[sortind]
        ind_los_pixels = ind_los_pixels[sortind]
        d_los_comov_offset = d_los_comov_offset[sortind]
        n_window = n_window[sortind]

        for i,ind_los_pix in enumerate(ind_los_pixels):
            window = NP.sqrt(frac_width * n_window[i]) * DSP.window_fftpow(n_window[i], shape=winshape, fftpow=fftpow, centering=True, peak=None, area_normalize=False, power_normalize=True)
            window_lospix = cube_d_los_comov[model].to('Mpc').value[ind_los_pix] + boxresolutions[model].to('Mpc').value * (NP.arange(n_window[i]) - int(n_window[i]/2))
            ind_window_lospix, inds_lospix, dlos = LKP.find_1NN(cube_d_los_comov[model].to('Mpc').value.reshape(-1,1), window_lospix.reshape(-1,1), distance_ULIM=0.51*boxresolutions[model].to('Mpc').value, remove_oob=True)
            sind = NP.argsort(ind_window_lospix)
            ind_window_lospix = ind_window_lospix[sind]
            inds_lospix = inds_lospix[sind]
            dlos = dlos[sind]
            window = window[ind_window_lospix]
            window = NP.pad(window, ((inds_lospix.min(), cube_d_los_comov[model].to('Mpc').value.size-1-inds_lospix.max())), mode='constant', constant_values=((0.0,0.0)))
            d_los_comov_wts[i,:] = window # nz x nlos
        npad = int(cube_d_los_comov[model].size * pad)
        kprll_twosided = 2 * NP.pi * DSP.spectral_axis(cube_d_los_comov[model].size + npad, delx=boxresolutions[model].to('Mpc').value, use_real=False, shift=True) / U.Mpc 
        ndim_padtuple = [(0,0), (0,0), (0,0), (0,npad)] 
        wlinv_krnl = DSP.FT1D(NP.pad(d_los_comov_wts, [(0,0), (0,npad)], mode='constant'), ax=-1, inverse=False, use_real=False, shift=True) * boxresolutions[model]
        power_kprll_twosided = DSP.FT1D(NP.pad(d_los_comov_wts[NP.newaxis,:,NP.newaxis,:] * Tb_los_fluctuations_modified[model].value, ndim_padtuple, mode='constant'), ax=-1, inverse=False, use_real=False, shift=True) * boxresolutions[model] * Tb_los_fluctuations_modified[model].unit # nsa x nz x nxy x nlos
        power_kprll_twosided = NP.abs(power_kprll_twosided) ** 2 / d_los_comov_width.reshape(1,-1,1,1)
        downsample_factor = NP.min((cube_d_los_comov[model].size + npad) * boxresolutions[model] / d_los_comov_width).decompose().value

        # kprll_twosided = DSP.downsampler(kprll_twosided.value, downsample_factor, axis=-1, method='FFT') * kprll_twosided.unit
        # wlinv_kernel[model] = DSP.downsampler(wlinv_krnl.value, downsample_factor, axis=-1, method='FFT') * wlinv_krnl.unit
        # power_kprll_twosided = DSP.downsampler(power_kprll_twosided.value, downsample_factor, axis=-1, method='FFT') * power_kprll_twosided.unit

        negind = NP.arange(kprll_twosided.size//2 + 1)
        posind = kprll_twosided.size//2 + NP.arange((kprll_twosided.size+1)//2)
        nsamples_onesided[model] = NP.ones(kprll_twosided.size//2 + 1, dtype=NP.float)
        nsamples_onesided[model][1:(kprll_twosided.size + 1)//2] = 2

        kprll_onesided[model] = NP.abs(kprll_twosided[:kprll_twosided.size//2 + 1][::-1])
        power_folded = OPS.reverse(power_kprll_twosided[:,:,:,negind].value, axis=3) * power_kprll_twosided.unit
        power_folded[:,:,:,:(kprll_twosided.size+1)//2] = power_folded[:,:,:,:(kprll_twosided.size+1)//2] + power_kprll_twosided[:,:,:,posind]
        pspec_kprll_onesided[model] = 0.5 * NP.mean(power_folded, axis=2)
        pspec_kprll_onesided[model] = DSP.downsampler(pspec_kprll_onesided[model].value, downsample_factor, axis=-1, method='interp', kind='cubic') * pspec_kprll_onesided[model].unit
        kprll_onesided[model] = DSP.downsampler(kprll_onesided[model].value, downsample_factor, axis=-1, method='interp', kind='linear') * kprll_onesided[model].unit
        nsamples_onesided[model] = DSP.downsampler(nsamples_onesided[model], downsample_factor, axis=-1, method='interp', kind='linear')

        nmodes_onesided[model] = d_los_comov_width.reshape(-1,1) * kprll_onesided[model].reshape(1,-1) / (2 * NP.pi)
        nsamples_total_onesided[model] = nsamples_onesided[model].reshape(1,-1) * nmodes_onesided[model]
        freq_resolution = boxresolutions[model] * (cosmo_parms[model].H0 * rest_freq * cosmo_parms[model].efunc(z_pspec) ) / (FC.c * (1+z_pspec)**2)
        bw = d_los_comov_width * (cosmo_parms[model].H0 * rest_freq * cosmo_parms[model].efunc(z_pspec) ) / (FC.c * (1+z_pspec)**2)
        noise_power_std_kprll_onesided[model] = ((Tsys / NP.sqrt(2 * n_ant * (n_ant-1) * freq_resolution.reshape(1,-1,1) * t_half * (1.0-data_loss_factor))) * NP.sqrt((bw/freq_resolution).decompose().value.reshape(1,-1,1)) * freq_resolution.reshape(1,-1,1))**2 * d_los_comov_width.reshape(1,-1,1) / bw.reshape(1,-1,1)**2 / NP.sqrt(nsamples_onesided[model].reshape(1,1,-1)) # Factor 2 is for dual-pol, See SIRA II Eq.(9-23), d_los_comov_width/boxresolutions[model] or bw/freq_resolution is for number of cells going into the Fourier transform, Factor of d_los_comov_width comes due to the integral of the window function along LOS, NP.sqrt(nsamples_onesided[model]) is the reduction in noise in power spectrum units
        # The expression below is equivalent to the above but done in comoving space whereas the one above is in frequency space
        # noise_power_std_kprll_onesided[model] = ((Tsys / NP.sqrt(2 * n_ant * (n_ant-1) * freq_resolution.reshape(1,-1,1) * t_half)) * NP.sqrt((bw/freq_resolution).decompose().value.reshape(1,-1,1)) * boxresolutions[model])**2 / d_los_comov_width.reshape(1,-1,1) / NP.sqrt(nsamples_onesided[model].reshape(1,1,-1)) # Factor 2 is for dual-pol, See SIRA II Eq.(9-23), d_los_comov_width/boxresolutions[model] or bw/freq_resolution is for number of cells going into the Fourier transform, Factor of d_los_comov_width comes due to the integral of the window function along LOS, NP.sqrt(nsamples_onesided[model]) is the reduction in noise in power spectrum units
        basic_sampstd_kprll_onesided[model] = pspec_kprll_onesided[model] / NP.sqrt(nsamples_onesided[model].reshape(1,1,-1))
        basic_uncertainty_kprll_onesided[model] = NP.sqrt(basic_sampstd_kprll_onesided[model]**2 + noise_power_std_kprll_onesided[model]**2)
        nsamples_required[model] = (detection_threshold / (pspec_kprll_onesided[model] / basic_uncertainty_kprll_onesided[model]))**2


	Renormalized the shaping window to unit power.
	Renormalized the shaping window to unit power.
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
	Renormalized the shaping window to unit power.
	Renormalized the shaping window to unit power.
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
	Renormalized the shaping window to unit power.
	Renormalized the shaping window to unit power.
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.

In [920]:
print(pspec_kprll_onesided[model].unit)
    print(noise_power_std_kprll_onesided[model].to('K2 Mpc').unit)


K2 Mpc
K2 Mpc

Plot 2: Based on re-estimating the $\delta T_\mathrm{b}$ fluctuations from $T_\mathrm{spin}$ cubes instead of $\langle T_\mathrm{spin}\rangle$


In [18]:
if ('2' in plots) or ('2a' in plots) or ('2b' in plots) or ('2c' in plots):

    modelnames = plot_info['2']['models']
    beam_majax = NP.asarray(plot_info['2']['beam_majax']) * U.arcsec
    beam_sa_0 = (beam_majax**2).to('sr') # Beam solid angle at z=0
    fov_fwhm_rad = NP.radians(plot_info['2']['fov_fwhm']) # Field of view FWHM (in rad)
    fov_hwhm_dircos = NP.sin(0.5*fov_fwhm_rad) # Field of HWHM in direction-cosine units
    flux_density_obs = NP.asarray(plot_info['2']['flux_obs']) * U.mJy
    beam_majax = beam_majax * NP.ones(flux_density_obs.size)
    beam_sa_0 = beam_sa_0 * NP.ones(flux_density_obs.size)
    flux_obs = (flux_density_obs / beam_sa_0).reshape(-1,1)
    freq_obs = plot_info['2']['freq_obs'] * U.Hz
    bl_max = 1.22 * (FC.c / freq_obs) / (beam_majax.to('rad').value) # Max baseline corresponding to beam major axis
    spindex = plot_info['2']['spindex']
    z_qso = NP.asarray(plot_info['2']['z_qso'])
    z_pspec = NP.asarray(plot_info['2']['z_pspec'])
    frac_bw = plot_info['2']['frac_bw']
    rest_freq = plot_info['2']['rest_freq'] * U.Hz
    if rest_freq is None:
        rest_freq = CNST.rest_freq_HI * U.Hz
    nrand = plot_info['2']['nrand']
    randseed = plot_info['2']['seed']
    subbandinfo = plot_info['2']['subband']
    winshape = subbandinfo['shape']
    fftpow = subbandinfo['fftpow']
    pad = subbandinfo['pad']
    Tsys = plot_info['2']['Tsys'] * U.K
    n_ant = plot_info['2']['n_antennas']
    t_int = plot_info['2']['t_int'] * U.s
    t_half = 0.5 * t_int
    detection_threshold = plot_info['2']['snr_min']
    data_loss_factor = plot_info['2']['data_loss']    
    
    indir_pfx = plot_info['2']['indir_pfx']
    indir_sfx = plot_info['2']['indir_sfx']
    
    nrand = plot_info['2']['nrand']
    randseed = plot_info['2']['seed']
    rstate = NP.random.RandomState(randseed)
    
    qtys = plot_info['2']['qtys']

In [19]:
print(beam_majax.to('arcsec'))
    print(fov_fwhm_rad)
    print(fov_hwhm_dircos)
    print(bl_max.to('m'))
    print(flux_obs)


[25. 25.] arcsec
3.141592653589793
1.0
[20117.51802179 20117.51802179] m
[[ 680722.72473844]
 [6807227.24738435]] mJy / sr

In [20]:
if ('2' in plots) or ('2a' in plots) or ('2b' in plots) or ('2c' in plots):

    sampled_los_qtys = {}    
    for modelind,model in enumerate(modelnames):
        sampled_los_qtys[model] = {}
        indir = indir_pfx + model + indir_sfx
        for qty in qtys:
            sampled_los_qtys[model][qty] = {}
            cubefile = indir + '{0}_light_cone.hdf5'.format(qty)
            cube, los_monopole, boxres, zmin, cosmo = read_HDF5_lightcone_cube(cubefile)
            rstate = NP.random.RandomState(randseed)
            pixind = rstate.randint(0, high=cube.shape[1]*cube.shape[2], size=nrand)
            sampled_los_qtys[model][qty]['cube'] = cube.reshape(cube.shape[0],-1)[:,pixind].T # (npix, nlos)
            sampled_los_qtys[model][qty]['monopole'] = los_monopole
        d_los_comov_min = cosmo.comoving_distance(zmin)
        sampled_los_qtys[model]['d_los_comov'] = d_los_comov_min + boxres * NP.arange(cube.shape[0]) # (nlos,)
        zgrid = NP.logspace(NP.log10(zmin-0.5), NP.log10(100.0), num=50000)
        d_los_comov_grid = cosmo.comoving_distance(zgrid)
        sampled_los_qtys[model]['redshifts'] = NP.interp(sampled_los_qtys[model]['d_los_comov'].to('Mpc').value, d_los_comov_grid.to('Mpc').value, zgrid)
#         sampled_los_qtys[model]['redshifts'] = NP.asarray([cosmology.z_at_value(cosmo.comoving_distance, d_los) for d_los in sampled_los_qtys[model]['d_los_comov']]) # (nlos,)
        sampled_los_qtys[model]['d_los_comov_center'] = cosmo.comoving_distance(z_pspec)
        zlo = z_pspec * (1 - 0.5 * frac_bw)
        zhi = z_pspec * (1 + 0.5 * frac_bw)
        d_los_comov_lo = cosmo.comoving_distance(zlo)
        d_los_comov_hi = cosmo.comoving_distance(zhi)
        sampled_los_qtys[model]['d_los_comov_width'] = d_los_comov_hi - d_los_comov_lo

        sampled_los_qtys[model]['boxres'] = boxres
        sampled_los_qtys[model]['pixang'] = boxres / cosmo.comoving_distance(sampled_los_qtys[model]['redshifts']) * U.rad # (nlos,)
        sampled_los_qtys[model]['pixsa'] = sampled_los_qtys[model]['pixang']**2
        sampled_los_qtys[model]['pixsa'] = sampled_los_qtys[model]['pixsa'].to('sr')
        sampled_los_qtys[model]['cosmo'] = copy.deepcopy(cosmo)
        sampled_los_qtys[model]['zmin'] = zmin
        
        sa_filling_factor_pix_beam = sampled_los_qtys[model]['pixsa'].reshape(1,-1) / beam_sa_0.reshape(-1,1) # (nbeam, nlos)
        sa_filling_factor_pix_beam[sa_filling_factor_pix_beam.decompose() <= 1.0] = 1.0 
        flux_los = flux_obs * (rest_freq/freq_obs)**spindex * (1+sampled_los_qtys[model]['redshifts'].reshape(1,-1))**(3.0-spindex) / sa_filling_factor_pix_beam # Correction for filling factor between pixel and beam, (nbeam, nlos)  
        sampled_los_qtys[model]['T_bg'] = {'cube': flux_los * (FC.c / rest_freq)**2 / (2 * FC.k_B)}  # (nbeam, nlos)
        sampled_los_qtys[model]['T_bg']['cube'] = sampled_los_qtys[model]['T_bg']['cube'].si.value * U.K  # (nbeam, nlos)
        
        los_ind_to_modify = sampled_los_qtys[model]['redshifts'].reshape(1,-1) <= z_qso.reshape(-1,1) # (nz_qso, nlos)
        sampled_los_qtys[model]['T_bg_net'] = {'cube': sampled_los_qtys[model]['T_bg']['cube'][:,NP.newaxis,:] * los_ind_to_modify[NP.newaxis,:,:].astype(NP.float) + sampled_los_qtys[model]['cosmo'].Tcmb0 * (1 + sampled_los_qtys[model]['redshifts'].reshape(1,1,-1))} # (nbeam, nz_qso, nlos)

        sampled_los_qtys[model]['delta_T_nomonopole'] = {}
        sampled_los_qtys[model]['delta_T_nomonopole']['cube'] = sampled_los_qtys[model]['delta_T']['cube'] - sampled_los_qtys[model]['delta_T']['monopole'].reshape(1,-1) # (npix, nlos)
        sampled_los_qtys[model]['delta_T_nomonopole']['cube'] = sampled_los_qtys[model]['delta_T_nomonopole']['cube'].to('K') # (npix, nlos)
 
        sampled_los_qtys[model]['modification'] = {}
        sampled_los_qtys[model]['modification']['cube'] = (sampled_los_qtys[model]['T_s']['cube'][NP.newaxis,NP.newaxis,:,:] - sampled_los_qtys[model]['T_bg_net']['cube'][:,:,NP.newaxis,:]) / (sampled_los_qtys[model]['T_s']['cube'][NP.newaxis,NP.newaxis,:,:] - sampled_los_qtys[model]['cosmo'].Tcmb0 * (1+sampled_los_qtys[model]['redshifts'].reshape(1,1,1,-1))) # (nbeam, nz_qso, npix, nlos)

        sampled_los_qtys[model]['delta_T_modified'] = {}
        sampled_los_qtys[model]['delta_T_modified_nomonopole'] = {}
        sampled_los_qtys[model]['delta_T_modified']['cube'] = sampled_los_qtys[model]['delta_T']['cube'][NP.newaxis,NP.newaxis,:,:] * sampled_los_qtys[model]['modification']['cube'] # (nbeam, nz_qs0, npix, nlos)
        sampled_los_qtys[model]['delta_T_modified']['cube'] = sampled_los_qtys[model]['delta_T_modified']['cube'].to('K') # (nbeam, nz_qso, npix, nlos)
        sampled_los_qtys[model]['delta_T_modified_nomonopole']['cube'] = sampled_los_qtys[model]['delta_T_modified']['cube'] - sampled_los_qtys[model]['delta_T']['monopole'].reshape(1,1,1,-1)

In [21]:
print(subbandinfo)
            print('Redshifts of background quasars = {0}'.format(z_qso))
            print('Redshifts where HI PS is estimated = {0}'.format(z_pspec))
            print('Center of LOS comoving distance to PS redshift: {0}/h'.format(sampled_los_qtys[model]['d_los_comov_center']))
            print('Width of LOS comoving distance at PS redshift: {0}/h'.format(sampled_los_qtys[model]['d_los_comov_width']))


{'fftpow': 1.0, 'shape': 'bhw', 'pad': 1.0}
Redshifts of background quasars = [10. 13. 21.]
Redshifts where HI PS is estimated = [ 8.25 11.5  18.  ]
Center of LOS comoving distance to PS redshift: [6250.34249052 6746.00634124 7322.42661988] Mpc/h
Width of LOS comoving distance at PS redshift: [105.50351902  93.67938598  78.24485193] Mpc/h

In [23]:
randpixind = NP.arange(10000)
#         losind = NP.arange(1783,1856)
        losind = NP.arange(764,862)
        mltdimind_0 = NP.ix_(randpixind,losind)
        mltdimind_1 = NP.ix_([0],[0],randpixind,losind)
        print(mltdimind_0)
        print(randpixind)
        print(sampled_los_qtys[model]['redshifts'][losind])
        print(sampled_los_qtys[model]['modification']['cube'].shape)
        print(NP.abs(sampled_los_qtys[model]['modification']['cube'][mltdimind_1]).min())
        print(NP.abs(sampled_los_qtys[model]['modification']['cube'][mltdimind_1]).max())
#         print(sampled_los_qtys[model]['modification']['cube'][mltdimind_1])
        print(sampled_los_qtys[model]['delta_T']['cube'][mltdimind_0].min())
        print(sampled_los_qtys[model]['delta_T']['cube'][mltdimind_0].max())
        print(sampled_los_qtys[model]['delta_T_nomonopole']['cube'][mltdimind_0].to('mK').min())
        print(sampled_los_qtys[model]['delta_T_nomonopole']['cube'][mltdimind_0].to('mK').max())
        print(sampled_los_qtys[model]['delta_T_modified']['cube'][mltdimind_1].to('mK').min())
        print(sampled_los_qtys[model]['delta_T_modified']['cube'][mltdimind_1].to('mK').max())
        print(sampled_los_qtys[model]['delta_T_modified_nomonopole']['cube'][mltdimind_1].to('mK').min())
        print(sampled_los_qtys[model]['delta_T_modified_nomonopole']['cube'][mltdimind_1].to('mK').max())
#         print(sampled_los_qtys[model]['modification']['cube'].squeeze())
#         print(sampled_los_qtys[model]['T_s']['cube'])
#         print(sampled_los_qtys[model]['delta_T_modified']['cube'].shape)
#         print(sampled_los_qtys[model]['delta_T_modified_nomonopole']['cube'].shape)


(array([[   0],
       [   1],
       [   2],
       ...,
       [9997],
       [9998],
       [9999]]), array([[764, 765, 766, 767, 768, 769, 770, 771, 772, 773, 774, 775, 776,
        777, 778, 779, 780, 781, 782, 783, 784, 785, 786, 787, 788, 789,
        790, 791, 792, 793, 794, 795, 796, 797, 798, 799, 800, 801, 802,
        803, 804, 805, 806, 807, 808, 809, 810, 811, 812, 813, 814, 815,
        816, 817, 818, 819, 820, 821, 822, 823, 824, 825, 826, 827, 828,
        829, 830, 831, 832, 833, 834, 835, 836, 837, 838, 839, 840, 841,
        842, 843, 844, 845, 846, 847, 848, 849, 850, 851, 852, 853, 854,
        855, 856, 857, 858, 859, 860, 861]]))
[   0    1    2 ... 9997 9998 9999]
[8.01435414 8.01967552 8.0250016  8.0303324  8.0356679  8.04100811
 8.04635306 8.05170273 8.05705714 8.06241628 8.06778018 8.07314882
 8.07852222 8.08390038 8.08928331 8.09467102 8.1000635  8.10546077
 8.11086283 8.11626969 8.12168134 8.12709781 8.1325191  8.1379452
 8.14337612 8.14881188 8.15425247 8.15969791 8.16514819 8.17060333
 8.17606332 8.18152819 8.18699792 8.19247253 8.19795203 8.20343641
 8.20892569 8.21441986 8.21991895 8.22542294 8.23093186 8.2364457
 8.24196446 8.24748817 8.25301681 8.25855041 8.26408895 8.26963246
 8.27518093 8.28073437 8.28629279 8.2918562  8.29742459 8.30299797
 8.30857636 8.31415976 8.31974817 8.3253416  8.33094005 8.33654353
 8.34215205 8.34776562 8.35338424 8.35900791 8.36463664 8.37027044
 8.37590931 8.38155326 8.38720231 8.39285644 8.39851567 8.40418
 8.40984945 8.41552401 8.42120369 8.42688851 8.43257846 8.43827355
 8.44397379 8.44967918 8.45538973 8.46110545 8.46682634 8.47255242
 8.47828367 8.48402012 8.48976177 8.49550862 8.50126068 8.50701796
 8.51278046 8.51854819 8.52432116 8.53009936 8.53588282 8.54167153
 8.54746551 8.55326475]
(2, 3, 10000, 3072)
0.956509145717
0.999070055599
0.0 mK
34.0754165649 mK
-11.7094955444 mK
22.6278018951 mK
0.0 mK
33.4708734294 mK
-11.7094954476 mK
22.0232584393 mK

In [26]:
if '2a' in plots:        
        zqso_ls = ['-', '--', ':']
        zqso_modification_colrs = ['lightblue', 'lightgreen', 'lightsalmon']
        zqso_mean_modification_colrs = ['blue', 'green', 'red']
        shade_maxyval_zqso = [0.5e3, 1e3, 2.5e3]
        zqso_modification_shade_colrs = ['0.85', '0.7', '0.55']
        ncols = beam_majax.size
        mdl_simplified_names = ['Bright Gal.', 'Faint Gal.', 'Cold IGM']
        for mdlind, mdl in enumerate(modelnames):
            print('\n{0}\n============================='.format(mdl))
            temperature_contrast_factor_orig = 1.0 - sampled_los_qtys[model]['cosmo'].Tcmb0 * (1 + sampled_los_qtys[model]['redshifts'].reshape(1,-1)) / sampled_los_qtys[mdl]['T_s']['cube'] # (npix, nlos)
            temperature_contrast_factor_orig = temperature_contrast_factor_orig.decompose().value # (npix, nlos)
            temperature_contrast_factor_modified = 1.0 - sampled_los_qtys[model]['T_bg_net']['cube'][:,:,NP.newaxis,:] / sampled_los_qtys[mdl]['T_s']['cube'][NP.newaxis,NP.newaxis,:,:] # (nbeam, nz_qso, npix, nlos)
            temperature_contrast_factor_modified = temperature_contrast_factor_modified.decompose().value # (nbeam, nz_qso, npix, nlos)

#             mean_modification_factor = 10**NP.mean(NP.log10(sampled_los_qtys[model]['modification']['cube']), axis=2) # (nbeam, nz_qso, nlos)
            mean_modification_factor = NP.std(sampled_los_qtys[mdl]['modification']['cube'], axis=2) # (nbeam, nz_qso, nlos)
            std_modification_factor = NP.std(sampled_los_qtys[mdl]['modification']['cube'], axis=2) # (nbeam, nz_qso, nlos)

            mean_modification_factor_v2 = NP.mean(temperature_contrast_factor_modified, axis=-2) / NP.mean(temperature_contrast_factor_orig, axis=-2) # (nbeam, nz_qso, nlos)

            fig, axs = PLT.subplots(ncols=ncols, sharey=True, figsize=(ncols*2, 3.5), squeeze=False)
            if axs.size == 1:
                axs = axs.reshape(-1)

            for beamind, beamsize in enumerate(beam_majax):
                increase_modification_ind_list = []
                for zqi, zq in enumerate(z_qso):
                    ind_increase_modification = NP.where(NP.abs(mean_modification_factor_v2[beamind,zqi,:]) > 1.0)[0]
                    increase_modification_ind_list += [ind_increase_modification]
                for zqi, zq in reversed(list(enumerate(z_qso))):
                    axs[0,beamind].vlines(1+sampled_los_qtys[mdl]['redshifts'][increase_modification_ind_list[zqi]], -shade_maxyval_zqso[zqi]/20, shade_maxyval_zqso[zqi]/20, color=zqso_modification_shade_colrs[zqi])
                axs[0,beamind].plot(1+sampled_los_qtys[mdl]['redshifts'], NP.mean(temperature_contrast_factor_orig, axis=-2), ls='-', lw=4, color='cyan', label=r'$T_\mathrm{bg}=T_\mathrm{cmb}$')
                for zqi, zq in reversed(list(enumerate(z_qso))):    
#                     axs[0,beamind].fill_between(1+sampled_los_qtys[mdl]['redshifts'], mean_modification_factor[beamind,zqi,:]-std_modification_factor[beamind,zqi,:], mean_modification_factor[beamind,zqi,:]+std_modification_factor[beamind,zqi,:], color=zqso_modification_colrs[zqi])
#                     axs[0,beamind].plot(1+sampled_los_qtys[mdl]['redshifts'], mean_modification_factor[beamind,zqi,:], lw=2, color=zqso_mean_modification_colrs[zqi], label=r'$z_\mathrm{qso}=$'+'{0:.1f}'.format(zq))
                    axs[0,beamind].plot(1+sampled_los_qtys[mdl]['redshifts'], NP.mean(temperature_contrast_factor_modified[beamind,zqi,:,:], axis=-2), ls=zqso_ls[zqi], lw=2, color='black', label=r'$T_\mathrm{bg} = [T_\mathrm{qso}(z=$'+'{0:.1f}'.format(zq)+r'$)+T_\mathrm{cmb}]$')
                        
                axs[0,beamind].set_yscale('symlog', linthreshy=1.0)
                axs[0,beamind].set_xscale('log')
                axs[0,beamind].set_xlim(32.5,6.5)
                axs[0,beamind].set_xlabel(r'$1+z$', fontsize=12, weight='medium')
                
                if beamind == 0:
                    axs[0,beamind].set_ylabel(r'$1-\frac{T_\mathrm{bg}}{T_\mathrm{s}}$', fontsize=12, weight='medium', labelpad=-10)
                    legend = axs[0,beamind].legend(loc='upper left', shadow=False, fontsize=7)

                axs[0,beamind].text(0.98, 0.02, r'$S_\mathrm{qso}=$'+'{0:.2f} mJy'.format((flux_obs[beamind]*beam_sa_0[beamind]).to('mJy').value[0]), transform=axs[0,beamind].transAxes, ha='right', va='bottom', fontsize=10)
                
                tickformatter = ticker.FuncFormatter(lambda x, _: '{:.6g}'.format(x))
                axs[0,beamind].xaxis.set_major_formatter(tickformatter)
                axs[0,beamind].xaxis.set_minor_formatter(tickformatter)
                axs[0,beamind].yaxis.set_major_formatter(tickformatter)
#                 axs[0,beamind].yaxis.set_minor_formatter(tickformatter)   
                
                fig.subplots_adjust(hspace=0, wspace=0.0)
                
                fig.subplots_adjust(top=0.98)
                fig.subplots_adjust(bottom=0.13)
                fig.subplots_adjust(left=0.15)
                fig.subplots_adjust(right=0.98)

                big_ax = fig.add_subplot(111)
                big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
                # big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
                big_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
                big_ax.set_xticks([])
                big_ax.set_yticks([])
#                 big_ax.set_xlabel(r'$1+z$', fontsize=12, weight='medium', labelpad=20)
            PLT.savefig(figdir+'{0}_mean_temperature_contrast_factor.pdf'.format(mdl), bbox_inches=0)

            ###############################

            fig, axs = PLT.subplots(ncols=ncols, sharey=True, figsize=(ncols*2, 3.5), squeeze=False)
            if axs.size == 1:
                axs = axs.reshape(-1)

            for beamind, beamsize in enumerate(beam_majax):                
                for zqi, zq in reversed(list(enumerate(z_qso))):
                    axs[0,beamind].vlines(1+sampled_los_qtys[mdl]['redshifts'][increase_modification_ind_list[zqi]], -shade_maxyval_zqso[zqi], shade_maxyval_zqso[zqi], color=zqso_modification_shade_colrs[zqi], alpha=0.5)
                for zqi, zq in enumerate(z_qso):
#                     axs[0,beamind].plot(1+sampled_los_qtys[mdl]['redshifts'], mean_modification_factor_v2[beamind,zqi,:], ls=zqso_ls[zqi], lw=2, color='black', label=r'CMB + $z_\mathrm{qso}=$'+'{0:.1f}'.format(zq))
                    axs[0,beamind].plot(1+sampled_los_qtys[mdl]['redshifts'], mean_modification_factor_v2[beamind,zqi,:], ls=zqso_ls[zqi], lw=2, color='black', label=r'$z_\mathrm{qso}=$'+'{0:.1f}'.format(zq))

                axs[0,beamind].set_yscale('symlog', linthreshy=1)
#                 axs[0,beamind].set_yscale('log')
                axs[0,beamind].set_xscale('log')
                axs[0,beamind].set_xlim(32.5,6.5)
                axs[0,beamind].set_xlabel(r'$1+z$', fontsize=12, weight='medium')
                if beamind == 0:
                    axs[0,beamind].set_ylabel('Contrast Modification factor', fontsize=12, weight='medium')
                    legend = axs[0,beamind].legend(loc='lower left', shadow=False, fontsize=7.5)

                axs[0,beamind].text(0.98, 0.98, r'$S_\mathrm{qso}=$'+'{0:.2f} mJy'.format((flux_obs[beamind]*beam_sa_0[beamind]).to('mJy').value[0]), transform=axs[0,beamind].transAxes, ha='right', va='top', fontsize=11)
                
#                 tickformatter = ticker.FuncFormatter(lambda y, _: '{:.16g}'.format(y))
                axs[0,beamind].xaxis.set_major_formatter(tickformatter)
                axs[0,beamind].xaxis.set_minor_formatter(tickformatter)
                axs[0,beamind].yaxis.set_major_formatter(tickformatter)
#                 axs[0,beamind].yaxis.set_minor_formatter(tickformatter)   
                
                fig.subplots_adjust(hspace=0, wspace=0)
                
                fig.subplots_adjust(top=0.98)
                fig.subplots_adjust(bottom=0.13)
                fig.subplots_adjust(left=0.19)
                fig.subplots_adjust(right=0.98)

#                 big_ax = fig.add_subplot(111)
#                 big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
#                 # big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
#                 big_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
#                 big_ax.set_xticks([])
#                 big_ax.set_yticks([])
#                 big_ax.set_xlabel(r'$1+z$', fontsize=12, weight='medium', labelpad=20)
            PLT.savefig(figdir+'{0}_mean_temperature_contrast_modification_factor.pdf'.format(mdl), bbox_inches=0)


Bright_galaxies_fiducial_1024
=============================

Faint_galaxies_fiducial_1024
=============================

In [28]:
if '2b' in plots:
        ncols = z_qso.size
#         nrows = beam_majax.size
        nrows = 2
        for mdlind, mdl in enumerate(modelnames):
            for zqi, zq in enumerate(z_qso):
                los_ind = NP.where(NP.abs(sampled_los_qtys[mdl]['d_los_comov'] - sampled_los_qtys[mdl]['cosmo'].comoving_distance(z_pspec[zqi])) <= 0.5 * sampled_los_qtys[mdl]['d_los_comov_width'][zqi])[0]

                fig, axs = PLT.subplots(nrows=nrows, ncols=beam_majax.size, sharex=True, figsize=(ncols*1.75, 3.5), squeeze=False)
                if axs.size == 1:
                    axs = axs.reshape(-1)
                for beamind, beamsize in enumerate(beam_majax):
                    
                    axs[0,beamind].plot(sampled_los_qtys[mdl]['delta_T_modified_nomonopole']['cube'][beamind,zqi,:,los_ind].to('mK').value.reshape(-1), (NP.zeros(nrand).reshape(-1,1)+sampled_los_qtys[mdl]['d_los_comov'][los_ind].reshape(1,-1)).ravel(), ls='none', marker='.', ms=0.5, color='black')
                    axs[0,beamind].plot(sampled_los_qtys[mdl]['delta_T_nomonopole']['cube'][:,los_ind].to('mK').value.reshape(-1), (NP.zeros(nrand).reshape(-1,1)+sampled_los_qtys[mdl]['d_los_comov'][los_ind].reshape(1,-1)).ravel(), ls='none', marker='.', ms=1, color='gray', alpha=0.5)
                    dnsty_orig, be_orig, patches = axs[1,beamind].hist(sampled_los_qtys[mdl]['delta_T_nomonopole']['cube'][:,los_ind].to('mK').value.reshape(-1), bins='auto', density=False, log=True, histtype='step', color='gray')
                    dnsty, be, patches = axs[1,beamind].hist(sampled_los_qtys[mdl]['delta_T_modified_nomonopole']['cube'][beamind,zqi,:,los_ind].to('mK').value.reshape(-1), bins='auto', density=False, log=True, histtype='step', color='black')
#                     print(be)
#                     print(dnsty)
#                     print(los_ind)
#                     print(sampled_los_qtys[mdl]['d_los_comov'][los_ind])

                    axs[0,beamind].set_xlim(-1e4,5e2)  
                    axs[0,beamind].set_ylim(sampled_los_qtys[mdl]['d_los_comov'][los_ind].to('Mpc').value.min(), sampled_los_qtys[mdl]['d_los_comov'][los_ind].to('Mpc').value.max()) 
                    axs[0,beamind].set_xscale('symlog', linthreshx=10.0)
                    axs[1,beamind].set_ylim(0.9,9e4)
                    if beamind == 0:
                        axs[0,beamind].set_ylabel(r'$D_\mathrm{C}\,[h^{-1}$ Mpc]', fontsize=12, weight='medium')
                        axs[1,beamind].set_ylabel(r'$N$', fontsize=12, weight='medium')
                    else:
#                         empty_tick_labels = ['']*len(axs[0,0].get_yticklabels())
#                         axs[0,beamind].set_yticklabels(empty_tick_labels)

                        axs[0,beamind].yaxis.set_major_formatter(ticker.NullFormatter())
                        axs[1,beamind].yaxis.set_major_formatter(ticker.NullFormatter())
                    axs[1,beamind].text(0.1, 0.95, r'$z=$'+'{0:.1f}'.format(z_pspec[zqi]), transform=axs[1,beamind].transAxes, ha='left', va='top')
                    axs[1,beamind].text(0.99, 0.95, r'$S_\mathrm{qso}=$'+'{0:.2g} mJy'.format((flux_obs[beamind]*beam_sa_0[beamind]).to('mJy').value[0]), transform=axs[1,beamind].transAxes, ha='right', va='top', fontsize=8.5)

                fig.subplots_adjust(hspace=0, wspace=0)
                fig.subplots_adjust(top=0.98)
                fig.subplots_adjust(bottom=0.13)
                fig.subplots_adjust(left=0.15)
                fig.subplots_adjust(right=0.98)

                big_ax = fig.add_subplot(111)
                big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
                # big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
                big_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
                big_ax.set_xticks([])
                big_ax.set_yticks([])
                big_ax.set_xlabel(r'$\delta T_\mathrm{b}$ [mK]', fontsize=12, weight='medium', labelpad=20)
#             PLT.savefig(figdir+'{0}_Tb_pixel_scatter_distribution.pdf'.format(mdl), bbox_inches=0)
            print(figdir+'{0}_Tb_pixel_scatter_distribution.pdf'.format(mdl))


/lustre/aoc/users/nthyagar/projects/21cmforest/figures/Bright_galaxies_fiducial_1024_Tb_pixel_scatter_distribution.pdf
/lustre/aoc/users/nthyagar/projects/21cmforest/figures/Faint_galaxies_fiducial_1024_Tb_pixel_scatter_distribution.pdf

In [29]:
if '2b' in plots:
        delta_T_range = [-8e3, 5e2]
        ncols = 2
#         nrows = beam_majax.size
        nrows = 1
        bgrad_colrs = ['darkturquoise']
        bgrad_colrs += map(str, NP.linspace(0.6, 0.3, beam_majax.size).tolist())
        for mdlind, mdl in enumerate(modelnames):
            for zqi, zq in enumerate(z_qso):
                los_ind = NP.where(NP.abs(sampled_los_qtys[mdl]['d_los_comov'] - sampled_los_qtys[mdl]['cosmo'].comoving_distance(z_pspec[zqi])) <= 0.5 * sampled_los_qtys[mdl]['d_los_comov_width'][zqi])[0]
                fig, axs = PLT.subplots(nrows=nrows, ncols=ncols, sharey=True, figsize=(ncols*2.5, 3.5), squeeze=False)
                if axs.size == 1:
                    axs = axs.reshape(-1)
                axs[0,0].plot(sampled_los_qtys[mdl]['d_los_comov'][los_ind], sampled_los_qtys[mdl]['delta_T_nomonopole']['cube'][nrand//500,los_ind].to('mK').value, ls='-', lw=1, color=bgrad_colrs[0], )
                dnsty_orig, be_orig, patches = axs[0,1].hist(sampled_los_qtys[mdl]['delta_T_nomonopole']['cube'][:,los_ind].to('mK').value.reshape(-1), bins='sqrt', range=delta_T_range, density=True, log=True, orientation='horizontal', histtype='step', color=bgrad_colrs[0])
                for beamind, beamsize in enumerate(beam_majax):
                    axs[0,0].plot(sampled_los_qtys[mdl]['d_los_comov'][los_ind], sampled_los_qtys[mdl]['delta_T_modified_nomonopole']['cube'][beamind,zqi,nrand//500,los_ind].to('mK').value, ls='-', lw=1, color=bgrad_colrs[beamind+1])
                    dnsty, be, patches = axs[0,1].hist(sampled_los_qtys[mdl]['delta_T_modified_nomonopole']['cube'][beamind,zqi,:,los_ind].to('mK').value.reshape(-1), bins='auto', range=delta_T_range, density=True, log=True, orientation='horizontal', histtype='step', color=bgrad_colrs[beamind+1])
                    axs[0,0].set_ylim(delta_T_range)
                    axs[0,0].set_yscale('symlog', linthreshy=10.0)
                    axs[0,0].text(0.98, 0.98, r'$z=$'+'{0:.1f}'.format(z_pspec[zqi]), transform=axs[0,0].transAxes, ha='right', va='top')
                    axs[0,0].set_ylabel(r'$\delta T_\mathrm{b}$ [mK]', fontsize=12, weight='medium')
                    axs[0,0].set_xlabel(r'$D_\mathrm{C}\,[h^{-1}$ Mpc]', fontsize=12, weight='medium')
                    axs[0,1].set_xlim(4e-6,None)
                    axs[0,1].set_xlabel('Density Function', fontsize=12, weight='medium')
                fig.subplots_adjust(hspace=0, wspace=0)
                fig.subplots_adjust(top=0.98)
                fig.subplots_adjust(bottom=0.15)
                fig.subplots_adjust(left=0.15)
                fig.subplots_adjust(right=0.98)

#                 big_ax = fig.add_subplot(111)
#                 big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
#                 # big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
#                 big_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
#                 big_ax.set_xticks([])
#                 big_ax.set_yticks([])
#                 big_ax.set_xlabel(r'$\delta T_\mathrm{b}$ [mK]', fontsize=12, weight='medium', labelpad=20)

            PLT.savefig(figdir+'delta_T_temp_pre_post_modification_z_{0:.1f}_{1}_.pdf'.format(z_pspec[zqi],mdl), bbox_inches=0)
            print(figdir+'delta_T_temp_pre_post_modification_z_{0:.1f}_{1}.pdf'.format(z_pspec[zqi],mdl))


/lustre/aoc/users/nthyagar/projects/21cmforest/figures/delta_T_temp_pre_post_modification_z_18.0_Bright_galaxies_fiducial_1024.pdf
/lustre/aoc/users/nthyagar/projects/21cmforest/figures/delta_T_temp_pre_post_modification_z_18.0_Faint_galaxies_fiducial_1024.pdf