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
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')
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)
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)
In [14]:
actions = [key for key in parms['actions'] if parms['actions'][key]['act']]
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)
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
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'])
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))
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))
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()))
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)
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
In [920]:
print(pspec_kprll_onesided[model].unit)
print(noise_power_std_kprll_onesided[model].to('K2 Mpc').unit)
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)
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']))
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)
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)
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))
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))