In [1]:
import numpy as np
import matplotlib.pyplot as plt
import cosmology as cosm
# Some plot settings
rcParams['font.size'] = 12
#rcParams['text.usetex'] = True
In [2]:
## UNITS
C = 3.e10 # cm/s
GPC = 1.e9 * 3.e18 # cm
ARCSEC = 5.e-6 # radians
ARCMIN = 3.e-4 # radians
YEAR = 3.e7 # seconds
In [3]:
## Values to use an example
ZS = 2.0
COSMO = cosm.Cosmology()
In [4]:
print "Omega_Lambda:", COSMO.l
print "Omega_Matter:", COSMO.m
source_distance = cosm.DChi(ZS, 0.0) # Gpc
print "Distance to a z=1 source [Gpc]:", cosm.DChi(1.0, 0.0)
print "Distance to a z=" + np.str(ZS)[0:2] + " source [Gpc]:", source_distance
print "Distance to a z=4 source [Gpc]:", cosm.DChi(4.0, 0.0)
Corrales & Paerels (2012): $$ \Delta \chi_{\rm max} \approx 4.5\ {\rm kpc}\ \chi_{Gpc} E^{-2}_{keV} a_{0.1}^{-2} $$
And in general (Corrales & Paerels 2012, Eq 1 of Miralde-Escude 1999): $$ \Delta \chi = \chi_{OS} \alpha^2 (1-x) $$
Keep in mind that the minimal x value is: $$ x_{min} \approx \frac{\alpha}{2 \sigma} $$ where $$ \sigma = \frac{10.4'}{a(0.1\ \mu{\rm m})\ E({\rm keV})} $$
In [5]:
def charsig(a, E):
""" returns the characteristic scattering angle, in arcseconds
a: [um] Grain size
E: [keV] Photon energy"""
charsig0 = 1.04 * 60 # arcsec per 1 um grain size per keV
return charsig0 / a / E
def delta_time(alpha, x, zs=ZS):
# alpha: [arcsec] Observation angle
time_os = cosm.DChi(zs, 0.0) * (GPC/C) / YEAR # years
return time_os * (alpha*ARCSEC)**2 * (1-x)
Set x=0 to get the coefficient (T):
$$ \Delta t = T (1-x) $$
In [6]:
# Time delay coefficient for 30 arcseconds
print delta_time(30, 0.0) # years
This agrees with my pen and paper calculation, good
In [7]:
alpha = np.logspace(-1,3,300)
plt.plot(alpha, delta_time(alpha, 0.0), 'k-')
plt.loglog()
plt.xlabel("Observation angle [alpha]")
plt.ylabel("Time delay coefficient [years]")
# Overplot the characteristic halo size
plt.axvline(charsig(0.1,1.0), ls='--', color='r')
Out[7]:
In [8]:
def tfade(a=0.1, E=1.0, zs=ZS):
time_os = cosm.DChi(zs, 0.0) * (GPC/C) / YEAR # years
sig_rad = charsig(a, E) * ARCSEC # radiansg
return time_os * sig_rad**2 / 2.0
In [9]:
print "Tfade for 1 keV halo from 0.1 um grains [years]:", tfade(a=0.1, E=1.0)
print "Tfade for 1 keV halo from 1.0 um grains [years]:", tfade(a=1.0, E=1.0)
In [10]:
ztest = np.linspace(0.5,4.0,20)
tfade_max = [tfade(a=0.1, zs=z) for z in ztest]
tfade_min = [tfade(a=1.0, zs=z) for z in ztest]
plt.plot(ztest, tfade_max, 'k-', lw=2, label='0.1 um grains')
plt.plot(ztest, tfade_min, 'k-', lw=2, alpha=0.5, label='1.0 um grains')
plt.semilogy()
plt.legend(loc='upper left', frameon=False)
plt.xlabel("Redshift of source ($z$)")
plt.ylabel("Halo lifetime [years]")
Out[10]:
In [11]:
def tmax(alpha, a=0.1, E=1.0, zs=ZS):
""" Calculate the maximum time delay associated with
observation angle and typical grain population
alpha: [arcmin] Observation angle
a: [micron] Grain size
E: [keV] Photon energy
zs: Redshift of source"""
time_os = cosm.DChi(zs, 0.0) * (GPC/C) / YEAR
alpha_rad = alpha * ARCSEC
xterm = 1 - alpha/(2*charsig(a,E))
return time_os * alpha_rad**2 * xterm
In [12]:
plt.plot(alpha, tmax(alpha, a=0.1), 'k-', lw=2, label='0.1 um grains')
plt.plot(alpha, tmax(alpha, a=1.0), 'k-', lw=2, alpha=0.5, label='1.0 um grains')
plt.loglog()
plt.ylim(1.0,1.e5)
plt.xlim(1.0,1.e3)
plt.legend(loc='upper left', frameon=False)
plt.xlabel("Observation angle [arcsec]")
plt.ylabel("Lookback time [years]")
# Overplot location of characteristic scattering angle
plt.axvline(charsig(0.1,1.0), ls='--', color='k')
plt.axvline(charsig(1.0,1.0), ls='--', color='0.5')
plt.axhline(1.e3, ls=':', color='r')
Out[12]:
In [13]:
## Pretty plot settings
#from astroML.plotting import setup_text_plots
#setup_text_plots(fontsize=8, usetex=True)
#import matplotlib as mpl
#mpl.rc('axes', labelsize='large')
#mpl.rc('xtick', labelsize='large')
#mpl.rc('ytick', labelsize='large')
In [14]:
fig = plt.figure(figsize=(6,2))
plt.plot(alpha, tmax(alpha, a=0.1), 'k-', lw=2, label='0.1 um grains')
plt.plot(alpha, tmax(alpha, a=1.0), 'k-', lw=2, alpha=0.5, label='1.0 um grains')
plt.loglog()
plt.ylim(1.0,1.e5)
plt.xlim(0.1,1.e3)
plt.legend(loc='upper left', frameon=False)#, fontsize=12)
plt.xlabel("Observation angle [arcsec]", size=15)
plt.ylabel("$\delta t$ [years]", size=15)
#plt.savefig('/Users/lia/Dropbox/HEAD2014/time_delay.pdf', format='pdf')
Out[14]:
In [15]:
print "Time delay at 10'':", tmax(30.0, a=0.1, zs=2)
In [ ]: