In [1]:
import numpy as np
import matplotlib.pyplot as plt
rcParams['font.size'] = 12
In [2]:
from halo import *
Cosmological parameters
In [3]:
IGM_upperlim_cosmology = cosmo.Cosmology(d=1.e-5)
print "H0 =", IGM_upperlim_cosmology.h0
print "Omega_Lambda =", IGM_upperlim_cosmology.l
print "Omega_M =", IGM_upperlim_cosmology.m
print "Omega_d_IGM =", IGM_upperlim_cosmology.d
Scattering Models
In [4]:
RHO_SIL, RHO_GRA = 3.8, 2.2 # g cm^-3
RGDRUDE = ss.makeScatmodel('RG','Drude')
GRAPHITE = ss.makeScatmodel('Mie','Graphite')
SILICATE = ss.makeScatmodel('Mie','Silicate')
In [70]:
## View dust distribution defaults
test = dust.Dustdist()
print 'Default rho =', test.rho
print 'Default p =', test.p
Raw data last run on: January 28, 2015
Calculate 0.5 keV optical depth as a function of grain size
In [5]:
ENERGY = 0.5
Z_S = 2.0
grainsize = np.linspace(0.1,1.0,20)
'''
taux_rg = [cosmo.CosmTauX(Z_S, E=ENERGY, dist=dust.Grain(a), scatm=RGDRUDE, cosm=IGM_upperlim_cosmology)
for a in grainsize]
taux_gra = [cosmo.CosmTauX(Z_S, E=ENERGY, dist=dust.Grain(a,rho=RHO_GRA), scatm=GRAPHITE, cosm=IGM_upperlim_cosmology)
for a in grainsize]
taux_sil = [cosmo.CosmTauX(Z_S, E=ENERGY, dist=dust.Grain(a,rho=RHO_SIL), scatm=SILICATE, cosm=IGM_upperlim_cosmology)
for a in grainsize]
'''
Calculate optical depth as a function of energy for two grain sizes: 0.1 and 1.0 um
In [6]:
small_a, big_a = 0.1, 1.0 # um
energy_range = np.linspace(0.5,1.5,20)
'''
taux_small_rg = [cosmo.CosmTauX(Z_S, E=energy, dist=dust.Grain(small_a), scatm=RGDRUDE, cosm=IGM_upperlim_cosmology)
for energy in energy_range]
taux_small_gra = [cosmo.CosmTauX(Z_S, E=energy, dist=dust.Grain(small_a,rho=RHO_GRA), scatm=GRAPHITE, cosm=IGM_upperlim_cosmology)
for energy in energy_range]
taux_small_sil = [cosmo.CosmTauX(Z_S, E=energy, dist=dust.Grain(small_a,rho=RHO_SIL), scatm=SILICATE, cosm=IGM_upperlim_cosmology)
for energy in energy_range]
taux_big_rg = [cosmo.CosmTauX(Z_S, E=energy, dist=dust.Grain(big_a), scatm=RGDRUDE, cosm=IGM_upperlim_cosmology)
for energy in energy_range]
taux_big_gra = [cosmo.CosmTauX(Z_S, E=energy, dist=dust.Grain(big_a,rho=RHO_GRA), scatm=GRAPHITE, cosm=IGM_upperlim_cosmology)
for energy in energy_range]
taux_big_sil = [cosmo.CosmTauX(Z_S, E=energy, dist=dust.Grain(big_a,rho=RHO_SIL), scatm=SILICATE, cosm=IGM_upperlim_cosmology)
for energy in energy_range]
'''
Calculate optical depth for the power law distribution of Silicate and Graphite grains
p=3.5 power law distribution between 0.1 um and 1.0 um
In [68]:
radii = np.linspace(0.1, 1.0, 50)
'''
taux_gra_dist = [cosmo.CosmTauX(Z_S, E=energy, dist=dust.Dustdist(rad=radii,rho=RHO_GRA), scatm=GRAPHITE, cosm=IGM_upperlim_cosmology) \
for energy in energy_range]
taux_sil_dist = [cosmo.CosmTauX(Z_S, E=energy, dist=dust.Dustdist(rad=radii,rho=RHO_SIL), scatm=SILICATE, cosm=IGM_upperlim_cosmology) \
for energy in energy_range]
'''
Out[68]:
January 29, 2015 - There was something wrong with one of the silicate calculations. Fill in that missing data point.
In [48]:
print energy_range[-3], taux_sil_dist[-3]
In [49]:
replacement = cosmo.CosmTauX(Z_S, E=energy_range[-3], dist=dust.Dustdist(rad=radii,rho=RHO_SIL), \
scatm=SILICATE, cosm=IGM_upperlim_cosmology)
In [50]:
print replacement
In [51]:
taux_sil_dist[-3] = replacement
Write everything to a file for later reference
Files last saved on: January 29, 2015
In [16]:
from astropy.io import ascii
In [22]:
'''
outfile = 'figure01_left.txt'
ascii.write([grainsize, taux_rg, taux_gra, taux_sil], outfile, \
names=['Grainsize','RG-Drude','Mie-Graphite','Mie-Silicate'])
'''
In [20]:
'''
outfile = 'figure01_right.txt'
ascii.write([energy_range, taux_small_rg, taux_small_gra, taux_small_sil, taux_big_rg, taux_big_gra, taux_big_sil], outfile, \
names=['Energy','Small-RG-Drude','Small-Mie-Graphite','Small-Mie-Silicate','Big-RG-Drude','Big-Mie-Graphite','Big-Mie-Silicate'])
'''
In [52]:
'''
outfile = 'figure01_right_color.txt'
ascii.write([energy_range, taux_gra_dist, taux_sil_dist], outfile, \
names=['Energy','DDist-Mie-Graphite','DDist-Mie-Silicate'])
'''
Load data
In [24]:
infile = 'figure01_left.txt'
data = ascii.read(infile)
taux_rg = data['RG-Drude']
taux_gra = data['Mie-Graphite']
taux_sil = data['Mie-Silicate']
In [26]:
infile = 'figure01_right.txt'
data = ascii.read(infile)
taux_small_rg = data['Small-RG-Drude']
taux_small_gra = data['Small-Mie-Graphite']
taux_small_sil = data['Small-Mie-Silicate']
taux_big_rg = data['Big-RG-Drude']
taux_big_gra = data['Big-Mie-Graphite']
taux_big_sil = data['Big-Mie-Silicate']
In [53]:
infile = 'figure01_right_color.txt'
data = ascii.read(infile)
taux_gra_dist = data['DDist-Mie-Graphite']
taux_sil_dist = data['DDist-Mie-Silicate']
Plot it up
In [54]:
def plot_tau_vs_grainsize(ax, ylim=(0.0,0.15)):
ax.plot(grainsize, taux_rg, 'k-', lw=2, label='RG-Drude')
plt.plot(grainsize, taux_gra, '--', color='k', lw=1, label='Mie-Graphite')
plt.plot(grainsize, taux_sil, color='k', ls='dashdot', lw=1, label='Mie-Silicate')
plt.legend(loc='upper right', fontsize=12)
plt.ylim(ylim)
plt.text(0.15, 0.13, r'$0.5\ {\rm keV}$', size=15)
return
In [55]:
def plot_tau_vs_energy(ax, ylim=(0.0,0.05)):
gcolor = {'small':'0.0', 'big':'0.3'}
ax.plot(energy_range, taux_small_rg, color=gcolor['small'], lw=2, label='RG-Drude')
ax.plot(energy_range, taux_small_gra, ls='--', color=gcolor['small'], lw=1, label='RG-Drude')
ax.plot(energy_range, taux_small_sil, ls='dashdot', color=gcolor['small'], lw=1, label='RG-Drude')
#plt.legend(loc='upper right')
ax.plot(energy_range, taux_big_rg, color=gcolor['big'], lw=2)
ax.plot(energy_range, taux_big_gra, ls='--', color=gcolor['big'], lw=1)
ax.plot(energy_range, taux_big_sil, ls='dashdot', color=gcolor['big'], lw=1)
plt.ylim(ylim)
plt.xlim(0.5,1.5)
plt.text(0.6, 0.007, r'$0.1\ \mu m$', color=gcolor['small'], size=15)
plt.text(1.0, 0.04, r'$1.0\ \mu m$', color=gcolor['big'], size=15)
return
In [64]:
def plot_tau_vs_energy_ddist(ax, ylim=(0.0, 0.05)):
gcolor = {'sil':'g', 'gra':'b'}
ax.plot(energy_range, taux_gra_dist, color=gcolor['gra'], ls='-', lw=2, alpha=0.8)
ax.plot(energy_range, taux_sil_dist, color=gcolor['sil'], ls='-', lw=2, alpha=0.8)
plt.text(0.77, 0.035, 'gra', color='b', size=12)
plt.text(0.60, 0.028, 'sil', color='g', size=12)
return
In [66]:
fig = plt.figure(figsize=(12,4))
left = plt.subplot(121)
plot_tau_vs_grainsize(left)
plt.xlabel(r'Grain size [$\mu$m]')
plt.ylabel(r'$\tau_{\rm sca}$', size=20)
right=plt.subplot(122)
plot_tau_vs_energy(right)
plt.xlabel(r'Energy [keV]')
plot_tau_vs_energy_ddist(right)
fig.subplots_adjust(wspace=0.1)
#plt.savefig('figure01.pdf', format='pdf')
In [9]: