In [1]:
cd ~/Dropbox/notebooks/ghost_halos
In [2]:
import numpy as np
import matplotlib.pyplot as plt
#rcParams['font.size'] = 12
font = {'size' : 15}
matplotlib.rc('font', **font)
Infrastructure needed to work with saved halo calculations
See mie_scattering_halos.ipynb
In [3]:
from astropy.io import fits
def load_halo(fname):
try:
hdulist = fits.open(fname)
except:
print "Cannot open "+fname
return
tdata = hdulist[1].data
hdulist.close()
return tdata['Obsangle'], tdata['Intensity'] # arcsec, arcsec^-2
Halos from uniformly distributed, diffuse IGM dust
In [4]:
# 0.1 - 1.0 um "grey" dust grain size distributions
alpha, uIGM_grey_gra = load_halo('halo_uIGM_grey_gra.fits')
alpha, uIGM_grey_sil = load_halo('halo_uIGM_grey_sil.fits')
# 1.0 um grains only
alpha, uIGM_big_gra = load_halo('halo_uIGM_gra.fits')
alpha, uIGM_big_sil = load_halo('halo_uIGM_sil.fits')
# MRN silicate (approximates SMC-like extinction curve)
alpha, uIGM_mrn_sil = load_halo('halo_uIGM_mrn_sil.fits')
Halos from a foreground galaxy
In [5]:
alpha, sIGM_mrn_sil = load_halo('halo_sIGM_mrn.fits')
alpha, sIGM_grey_gra = load_halo('halo_sIGM_grey_gra.fits')
alpha, sIGM_grey_sil = load_halo("halo_sIGM_grey_sil.fits")
Infrastructure needed for PSF template
See plot_radprofile.ipynb
Fit parameters saved as plot_radprofile_fitparams.txt
In [6]:
from astropy.modeling import models
XTRANS = 0.5 # pixels, core-wing transition radius, min x val for wing portion
def wing_model(x, amp, p, xtrans=XTRANS):
wmod = models.PowerLaw1D(amplitude=amp, x_0=1.0, alpha=p)
expmod = np.exp(-xtrans/x)
return wmod(x) * expmod
def core_model(x, amp, p, g):
cmod = models.Beta1D(amplitude=amp, x_0=0.0, alpha=p, gamma=g)
return cmod(x)
def psf_model(params, x):
a1, p1, g1, a2, p2, c0 = params
result = core_model(x, a1, p1, g1) + wing_model(x, a2, p2)
return result # arcsec^-2
In [7]:
from astropy.io import ascii
pfilename = 'plot_radprofile_fitparams.txt'
ptable = ascii.read(pfilename)
pfit = np.zeros(6)
for i in range(6):
pfit[i] = ptable[i][0]
pnames = ['Beta1d Amp', 'Beta1d power', 'Beta1d gamma', 'Powerlaw Amp', 'Powerlaw power', 'Background']
pdict = dict(zip(pnames,pfit))
for key in pdict: print key, ':', pdict[key]
Plot up the PSF template
In [8]:
def plot_IGM_dust_halos(ax):
ax.plot(alpha, uIGM_grey_gra, 'k-', lw=2, label='Grey IGM')
ax.plot(alpha, uIGM_big_gra, 'k-', label='$1\ \mu m$ IGM')
ax.plot(alpha, uIGM_mrn_sil, 'r-', label='SMC-like IGM')
plt.legend(loc='upper right', fontsize=12)
plt.text(1.5, 0.015, 'PSF', color='0.5', size=15)
ax.plot(alpha, psf_model(pfit,alpha), '-', lw=2, color='0.5', alpha=0.8)
plt.loglog()
plt.xlabel('Observation Angle [arcsec]', size=15)
plt.ylabel(r'Normalized Intensity [arcsec$^{-2}$]', size=15)
plt.title('Uniformly distributed IGM')
return
ax = plt.subplot(111)
plot_IGM_dust_halos(ax)
# Plot dust from foreground galaxy CGM
#plt.plot(alpha, sIGM_mrn_sil, 'k--', lw=1, alpha=0.3)
#plt.plot(alpha, sIGM_grey_gra, 'k-', lw=2, alpha=0.3)
In [9]:
## From dust near Her X-1
alpha, fg_herx1_int = load_halo('fg_herx1.fits')
## From dust in a MW CGM
alpha, fg_sIGM_mrn = load_halo('fg_halo_IGMscreen.fits')
In [10]:
def plot_CGM_dust_halos(ax):
ax.plot(alpha, fg_sIGM_mrn, 'r:', label='MW CGM (MRN)')
ax.plot(alpha, sIGM_mrn_sil, ls='dashdot', color='r', label='$z=0.5$ CGM (MRN)')
#ax.plot(alpha, sIGM_grey_gra, 'k-', lw=2, label='$z=0.5$ CGM (grey)')
#ax.plot(alpha, fg_herx1_int, 'r-', label='Her X-1 Foreground')
plt.legend(loc='upper right')
#plt.text(1.5, 0.015, 'PSF', color='0.5', size=15)
ax.plot(alpha, psf_model(pfit,alpha), '-', lw=2, color='0.5', alpha=0.8)
plt.loglog()
plt.xlabel('Observation Angle [arcsec]', size=15)
plt.ylabel('Normalized intensity [arcsec$^{-2}$]', size=15)
plt.title('Possible sources of foreground scattering')
# 1. Milky Way foreground -- use MRN silicate
# 2. CGM from a z=0.5 galaxy -- use MRN silicate
# 3. Contamination to PSF (Galactic ISM)
return
ax = plt.subplot(111)
plot_CGM_dust_halos(ax)
In [11]:
fig = plt.figure(figsize=(12,4))
ax1 = plt.subplot(121)
plot_IGM_dust_halos(ax1)
ax2 = plt.subplot(122)
plot_CGM_dust_halos(ax2)
plt.ylabel('')
ax2.yaxis.set_ticklabels('')
fig.subplots_adjust(wspace=0.0)
In [12]:
def plot_IGM_ratios(ax):
ratio_grey = uIGM_grey_gra / psf_model(pfit, alpha)
ratio_1um = uIGM_big_gra / psf_model(pfit, alpha)
ratio_red = uIGM_mrn_sil / psf_model(pfit, alpha)
ax.plot( alpha, ratio_grey, color='k', lw=2 )
ax.plot( alpha, ratio_1um, color='k' )
ax.plot( alpha, ratio_red, color='r' )
return
ax = plt.subplot(111)
plot_IGM_ratios(ax)
plt.semilogx()
plt.title('Uniform IGM halo ratios')
plt.xlabel('Observation angle [arcsec]')
plt.ylabel(r'$I_h/I_{\rm PSF}$')
Out[12]:
In [13]:
def plot_CGM_ratios(ax):
ratio_CGM = sIGM_mrn_sil / psf_model(pfit, alpha)
ratio_MW = fg_sIGM_mrn / psf_model(pfit, alpha)
ax.plot( alpha, ratio_CGM, color='r', ls='--' )
ax.plot( alpha, ratio_MW, color='r', ls='dashdot' )
return
ax = plt.subplot(111)
plot_CGM_ratios(ax)
plt.semilogx()
plt.title('CGM halo ratios')
plt.xlabel('Observation angle [arcsec]')
plt.ylabel(r'$I_h/I_{\rm PSF}$')
Out[13]:
In [14]:
import matplotlib.gridspec as gridspec
In [20]:
fig = plt.figure( figsize=(8,8) )
gs = gridspec.GridSpec(2, 1, height_ratios=[2,1])
ax1 = plt.subplot(gs[0])
plot_IGM_dust_halos(ax1)
ax1.plot(alpha, sIGM_mrn_sil, ls='--', color='r', label='$z_g=0.5$ screen')
ax1.plot(alpha, fg_sIGM_mrn, 'r', ls='dashdot', label='MW foreground')
ax1.set_title('')
ax1.xaxis.set_ticklabels('')
ax1.set_xlabel('')
ax1.set_ylim(3.e-12, 1.0)
plt.legend(loc='upper right', fontsize=12)
ax2 = plt.subplot(gs[1])
plot_IGM_ratios(ax2)
plot_CGM_ratios(ax2)
ax2.set_xlabel( 'Observation angle [arcsec]' )
ax2.set_ylabel(r'$I_h / I_{\rm PSF}$')
ax2.yaxis.set_ticks(np.arange(0.05,0.4,0.1))
ax2.set_ylim(0, 0.4)
plt.semilogx()
gs.update(hspace=0.0)
#plt.savefig('figure02.pdf', format='pdf')
In [15]:
In [15]: