In [2]:
cd ~/Dropbox/notebooks/ghost_halos
In [3]:
import numpy as np
import matplotlib.pyplot as plt
rcParams['font.size'] = 12
from halo import *
Set up dust distributions
In [4]:
RHO_GRA, RHO_SIL, RHO_DRUDE = 2.2, 3.8, 3.0 # g cm^-3
big_gra = dust.Grain(1.0, rho=RHO_GRA)
big_sil = dust.Grain(1.0, rho=RHO_SIL)
grey_gra = dust.Dustdist(rad=np.linspace(0.1,1.0,20), p=3.5, rho=RHO_GRA)
grey_sil = dust.Dustdist(rad=np.linspace(0.1,1.0,20), p=3.5, rho=RHO_SIL)
mrn_rgd = dust.Dustdist(rad=np.linspace(0.005,0.3,20), p=3.5, rho=RHO_DRUDE)
mrn_gra = dust.Dustdist(rad=np.linspace(0.005,0.3,20), p=3.5, rho=RHO_GRA)
mrn_sil = dust.Dustdist(rad=np.linspace(0.005,0.3,20), p=3.5, rho=RHO_SIL)
Set up scattering models
In [5]:
rgdrude = ss.makeScatmodel('RG','Drude')
mie_gra = ss.makeScatmodel('Mie','Graphite')
mie_sil = ss.makeScatmodel('Mie','Silicate')
In [6]:
## What is the optical depth of MRN dust at 1 keV?
md_nh21 = 1.e21 * c.mp() * 0.009
kappa_Gal = ss.Kappascat(E=1.0, \
dist=dust.Dustspectrum(rad=mrn_rgd, md=md_nh21), \
scatm=rgdrude)
print 'Taux (1 keV) for Galactic (MRN) grains, per 10^21 H:', \
kappa_Gal.kappa[0] * kappa_Gal.dist.md
Initialize cosmology and other relevant parameters for halo objects
In [7]:
E0 = 1.0 # keV
ALPHA = np.logspace(-1.0,3.0,100) # arcsec (observation angle)
ZS = 2.0 # Redshift of source quasar
ZG = 0.5 # Redshift of a foreground galaxy
lcdm = cosmo.Cosmology(d=3.e-6) # Omega_d^IGM for z~0.5 from Menard & Fukugita (2012)
print "LCDM cosmology:"
print "H0 =", lcdm.h0
print "Omega_M =", lcdm.m
print "Omega_Lambda =", lcdm.l
print "Omega_d(IGM) =", lcdm.d
Set up pyfits to save the outputs as fits files (saves time for plotting later)
In [8]:
from astropy.io import fits
In [9]:
def save_halo_obj(halo, fname, galtype=False):
# Set up the table
col1 = fits.Column(name='Obsangle', format='E', array=halo.alpha) # arcsec
col2 = fits.Column(name='Intensity', format='E', array=halo.intensity ) # arcsec**-2
cols = fits.ColDefs([col1, col2])
tbhdu = fits.BinTableHDU.from_columns(cols)
# Set up the fits file with header
prihdr = fits.Header()
prihdr['COMMENT'] = "lia@space.mit.edu, mie_scattering_halos.ipynb"
if galtype:
prihdr['HALOTYPE'] = "Galactic"
prihdr['MD'] = np.str(halo.dist.md)
prihdr['ISMTYPE'] = halo.htype.ismtype
else:
prihdr['HALOTYPE'] = "Cosmological"
prihdr['IGMTYPE'] = halo.htype.igmtype
prihdr['ZS'] = np.str(halo.htype.zs)
if halo.htype.igmtype == 'Screen':
prihdr['ZG'] = np.str(halo.htype.zg)
prihdr['MD'] = np.str(halo.dist.md)
if halo.htype.igmtype == 'Uniform':
prihdr['OMEGA_D'] = np.str(halo.htype.cosm.d)
prihdr['SCATM'] = halo.scatm.stype
if type(halo.taux) == numpy.float64: prihdr['TAUX'] = np.str(halo.taux)
else: prihdr['TAUX'] = np.str(halo.taux[0])
prihdu = fits.PrimaryHDU(header=prihdr)
# Write the fits
thdulist = fits.HDUList([prihdu, tbhdu])
thdulist.writeto(fname)
return
In [10]:
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
In [11]:
def load_halo_taux(fname):
try:
hdulist = fits.open(fname)
return hdulist[0].header['TAUX']
except:
print "Cannot open "+fname
return
In [12]:
#halo_uIGM_gra = Halo(E0=E0, alpha=ALPHA, rad=big_gra, scatm=mie_gra)
#UniformIGM(halo_uIGM_gra, zs=ZS, cosm=lcdm)
#save_halo_obj(halo_uIGM_gra, "halo_uIGM_gra.fits")
In [13]:
#halo_uIGM_sil = Halo(E0=E0, alpha=ALPHA, rad=big_sil, scatm=mie_sil)
#UniformIGM(halo_uIGM_sil, zs=ZS, cosm=lcdm)
#save_halo_obj(halo_uIGM_sil, "halo_uIGM_sil.fits")
In [14]:
#i_g, i_s = halo_uIGM_gra.intensity, halo_uIGM_sil.intensity
alpha, i_g = load_halo('halo_uIGM_gra.fits')
alpha, i_s = load_halo('halo_uIGM_sil.fits')
plt.plot(alpha, i_g, 'b-', label='Graphite')
plt.plot(alpha, i_s, 'g-', label='Silicate')
plt.loglog()
plt.legend(loc='upper right')
plt.xlabel('Observation angle [arcsec]')
plt.ylabel('Normalized intensity [arcsec$^{-2}$]')
plt.title('1 um IGM dust')
Out[14]:
In [15]:
#halo_uIGM_grey_gra = Halo(E0=E0, alpha=ALPHA, rad=grey_gra, scatm=mie_gra)
#UniformIGM(halo_uIGM_grey_gra, zs=ZS, cosm=lcdm)
#save_halo_obj(halo_uIGM_grey_gra, "halo_uIGM_grey_gra.fits")
In [16]:
#halo_uIGM_grey_sil = Halo(E0=E0, alpha=ALPHA, rad=grey_sil, scatm=mie_sil)
#UniformIGM(halo_uIGM_grey_sil, zs=ZS, cosm=lcdm)
#save_halo_obj(halo_uIGM_grey_sil, "halo_uIGM_grey_sil.fits")
In [17]:
#ig_g, ig_s = halo_uIGM_grey_gra.intensity, halo_uIGM_grey_sil.intensity
alpha, ig_g = load_halo('halo_uIGM_grey_gra.fits')
alpha, ig_s = load_halo('halo_uIGM_grey_sil.fits')
plt.plot(alpha, i_g, 'b-', label='1 um Graphite')
plt.plot(alpha, ig_g, 'b--', label='')
plt.plot(alpha, i_s, 'g-', label='1 um Silicate')
plt.plot(alpha, ig_s, 'g--', label='0.1-1 um dist')
plt.legend(loc='lower left', frameon=False)
plt.loglog()
plt.xlabel('Observation angle [arcsec]')
plt.ylabel('Normalized intensity [arcsec$^{-2}$]')
plt.title('Uniformly distributed IGM dust')
Out[17]:
There is not much difference between the 1 um only and 0.1-1 um smooth distribution of dust grains
There is also not much difference between Graphite and Silicate dust grains
In [18]:
#print 'Taux (1 keV), 1 um dust model:', halo_uIGM_gra.taux, halo_uIGM_sil.taux
#print 'Taux (1 keV), 0.1-1 um dust model:', halo_uIGM_grey_gra.taux, halo_uIGM_grey_sil.taux
print 'Taux (1 keV), 1 um dust model:', \
load_halo_taux("halo_uIGM_gra.fits"), load_halo_taux("halo_uIGM_sil.fits")
print 'Taux (1 keV), 0.1-1 um dust model:', \
load_halo_taux("halo_uIGM_grey_gra.fits"), load_halo_taux("halo_uIGM_grey_sil.fits")
In [19]:
#halo_uIGM_mrn_sil = Halo(E0=E0, alpha=ALPHA, rad=mrn_sil, scatm=rgdrude)
#UniformIGM(halo_uIGM_mrn_sil, zs=ZS, cosm=lcdm)
#save_halo_obj(halo_uIGM_mrn_sil, "halo_uIGM_mrn_sil.fits")
In [20]:
alpha, i_mrn_sil = load_halo("halo_uIGM_mrn_sil.fits")
In [21]:
plt.plot(alpha, ig_g, 'b--', label='0.1-1 um Graphite')
plt.plot(alpha, ig_s, 'g--', label='0.1-1 um Silicate')
plt.plot(alpha, i_mrn_sil, 'r-', label='MRN silicate')
plt.legend(loc='lower left', frameon=False)
plt.loglog()
plt.xlabel('Observation angle [arcsec]')
plt.ylabel('Normalized intensity [arcsec$^{-2}$]')
plt.title('Uniformly Distributed IGM dust')
Out[21]:
From Menard 2010:
$$ M_d = \frac{A_V}{1.086\ \kappa_V} \approx 4 \times 10^{-3} \left( \frac{r_p}{100\ h^{-1} \ {\rm kpc}} \right)^{-0.86} \kappa_V^{-1} \ {\rm g} \ {\rm cm}^2 $$They use extinction opacity for SMC-like dust from WD01:
$$ \kappa_V \approx 1.5 \times 10^{4} \ {\rm cm}^2 \ {\rm g}^{-1}$$What is the opacity of my simple MRN-like dust model?
In [22]:
V_um = 0.55 # microns, V-band center
V_keV = c.kev2lam() / (0.55*1.e-4) #cm keV / cm
kappa_gra = ss.Kappascat(E=V_keV, scatm=mie_gra, dist=dust.Dustspectrum(rad=mrn_gra)).kappa[0]
kappa_sil = ss.Kappascat(E=V_keV, scatm=mie_sil, dist=dust.Dustspectrum(rad=mrn_sil)).kappa[0]
print "Kappa_V (MRN, graphite):", kappa_gra
print "Kappa_V (MRN, silicate):", kappa_sil
In [23]:
def menard_mass_column(rp, kappa):
# rp -- impact parameter [kpc]
# kappa -- V-band opacity [cm^2 g^-1]
return 4.e-3 * np.power(rp/100.0, -0.86) / kappa
In [24]:
kappa_WD01 = 1.5e4 # g cm^-2
r_impact = np.logspace(1.0,3.0,100)
md_M10 = menard_mass_column(r_impact, kappa_WD01)
md_mrn_gra = menard_mass_column(r_impact, kappa_gra)
md_mrn_sil = menard_mass_column(r_impact, kappa_sil)
plt.plot(r_impact, md_M10, 'k-', lw=2, alpha=0.5)
plt.plot(r_impact, md_mrn_gra, 'b--')
plt.plot(r_impact, md_mrn_sil, 'g--')
plt.loglog()
plt.xlabel('Impact parameter [$h^{-1}$ kpc]')
plt.ylabel('Dust mass column [g cm$^{-2}$]')
Out[24]:
This makes sense, SMC-like grains are mostly Silicate (no strong 2175 Angs bump)
Dust mass column at $r_p = 100 h^{-1}$ kpc
In [25]:
kappa_grey_gra = ss.Kappascat(E=V_keV, scatm=mie_gra,
dist=dust.Dustspectrum(rad=grey_gra)).kappa[0]
kappa_grey_sil = ss.Kappascat(E=V_keV, scatm=mie_sil,
dist=dust.Dustspectrum(rad=grey_sil)).kappa[0]
MD100_menard = menard_mass_column(100.0, kappa_WD01)
MD100_grey_gra = menard_mass_column(100.0, kappa_grey_gra)
MD100_grey_sil = menard_mass_column(100.0, kappa_grey_sil)
print "M_d(100 kpc), Menard 2010:", MD100_menard
print "M_d(100 kpc), Grey graphite:", MD100_grey_gra
print "M_d(100 kpc), Grey silicate:", MD100_grey_sil
We get similar dust mass columns, regardless of grey dust model
Note that the MRN silicate grain kappa_V is similar to the value used by MSFR10.
Therefore use MRN distribution of silicate grains as an approximate for WD01 SMC dust (used by MSFR10 to fit extinction curve and estimate dust mass).
In [26]:
#halo_sIGM_mrn = Halo(E0, alpha=ALPHA, rad=mrn_sil, scatm=mie_sil)
#ScreenIGM(halo_sIGM_mrn, zs=ZS, zg=ZG, md=MD100_menard)
#save_halo_obj(halo_sIGM_mrn, 'halo_sIGM_mrn.fits')
In [27]:
#halo_sIGM_grey_gra = Halo(E0, alpha=ALPHA, rad=grey_gra, scatm=mie_gra)
#ScreenIGM(halo_sIGM_grey_gra, zs=ZS, zg=ZG, md=MD100_grey_gra)
#save_halo_obj(halo_sIGM_grey_gra, 'halo_sIGM_grey_gra.fits')
In [28]:
#halo_sIGM_grey_sil = Halo(E0, alpha=ALPHA, rad=grey_sil, scatm=mie_sil)
#ScreenIGM(halo_sIGM_grey_sil, zs=ZS, zg=ZG, md=MD100_grey_sil)
#save_halo_obj(halo_sIGM_grey_sil, 'halo_sIGM_grey_sil.fits')
In [29]:
alpha, sMRN_int = load_halo("halo_sIGM_mrn.fits")
alpha, sg_gra = load_halo("halo_sIGM_grey_gra.fits")
alpha, sg_sil = load_halo("halo_sIGM_grey_sil.fits")
plt.plot(ALPHA, sMRN_int, 'r-', label='MRN CGM')
plt.plot(ALPHA, sg_gra, 'b:', lw=2, label='0.1-1 um Graphite CGM')
plt.plot(ALPHA, sg_sil, 'g:', lw=2, label='0.1-1 um Silicate CGM')
plt.loglog()
plt.legend(loc='lower left', frameon=False)
plt.xlabel('Observation angle [arcsec]')
plt.ylabel('Normalized intensity [arcsec$^{-2}$]')
plt.title('Scattering from dust in $z=' + np.str(ZG) + '$ CGM')
Out[29]:
Again, not much difference between graphite and silicate case
In [30]:
from galhalo import *
In [31]:
fg_halo = Halo(E0=E0, alpha=ALPHA, rad=mrn_rgd, scatm=rgdrude)
UniformISM(fg_halo, NH=1.e19, nx=5000)
print "Mass column [g cm^-2]:", fg_halo.dist.md
print "Optical depth:", fg_halo.taux[0]
In [32]:
fg_halo_screen = Halo(E0=E0, alpha=ALPHA, rad=mrn_rgd, scatm=rgdrude)
DiscreteISM(fg_halo_screen, NH=1.e19, xg=0.999)
print fg_halo_screen.taux[0]
In [33]:
fg_halo_IGMscreen = Halo(E0=E0, alpha=ALPHA, rad=mrn_rgd, scatm=rgdrude)
ScreenIGM(fg_halo_IGMscreen, zs=ZS, zg=0.0, md=1.e19*c.mp()*0.009)
#save_halo_obj(fg_halo_IGMscreen, "fg_halo_IGMscreen.fits")
In [34]:
plt.plot(ALPHA, fg_halo.intensity, 'k-', label='Uniform Galactic')
plt.plot(ALPHA, fg_halo_screen.intensity, 'k--', label='Galactic Screen')
plt.plot(ALPHA, fg_halo_IGMscreen.intensity, 'k:', label='Cosmological Screen')
plt.legend(loc='upper right', frameon=False)
plt.loglog()
plt.xlabel('Observation angle [arcsec]')
plt.ylabel('Normalized intensity [arcsec$^{-2}$]')
plt.title('Possible foreground MW contamination, $N_H=10^{19}$, MRN dust')
Out[34]:
In [35]:
#fg_halo_grey_gra = Halo(E0=E0, alpha=ALPHA, rad=grey_gra, scatm=mie_gra)
#ScreenIGM(fg_halo_grey_gra, zs=ZS, zg=0.0, md=1.e19*c.mp()*0.009)
#save_halo_obj(fg_halo_grey_gra, 'fg_halo_grey_gra.fits')
In [36]:
#fg_halo_grey_sil = Halo(E0=E0, alpha=ALPHA, rad=grey_sil, scatm=mie_sil)
#ScreenIGM(fg_halo_grey_sil, zs=ZS, zg=0.0, md=1.e19*c.mp()*0.009)
#save_halo_obj(fg_halo_grey_sil, 'fg_halo_grey_sil.fits')
In [37]:
alpha, fg_ggra = load_halo('fg_halo_grey_gra.fits')
alpha, fg_gsil = load_halo('fg_halo_grey_sil.fits')
plt.plot(ALPHA, fg_halo_IGMscreen.intensity, 'k-', label='MRN dust')
plt.plot(ALPHA, fg_ggra, 'b--', label='Grey Graphite' )
plt.plot(ALPHA, fg_gsil, 'g--', label='Grey Silicate')
plt.legend(loc='lower left', frameon=False, fontsize=12)
plt.loglog()
plt.title('Foreground dust, Milky Way CGM, $N_H = 10^{19}$')
plt.xlabel('Observation angle [arcsec]')
plt.ylabel('Normalized intensity [arcsec$^{-2}$]')
Out[37]:
Note that here, Graphite vs Silicate makes a difference
In [38]:
NH_convert = 5.8e21 # cm^2 mag^-1
NH_herx1 = 0.014 * NH_convert
print NH_herx1
In [39]:
#fg_herx1 = Halo(E0=E0, alpha=ALPHA, rad=mrn_rgd, scatm=rgdrude)
#UniformISM(fg_herx1, NH=NH_herx1, nx=5000)
#save_halo_obj(fg_herx1, 'fg_herx1.fits', galtype=True)
#kappa_herx1 = ss.Kappascat(E=E0, scatm=rgdrude, dist=dust.Dustspectrum(rad=mrn_rgd))
#print fg_herx1.taux[0], kappa_herx1.kappa[0] * fg_herx1.dist.md
In [40]:
alpha, fg_herx1_int = load_halo('fg_herx1.fits')
plt.plot(alpha, fg_herx1_int, 'k-')
plt.loglog()
plt.xlabel('Observation angle [arcsec]')
plt.ylabel('Normalized intensity [arcsec$^{-2}$]')
plt.title('Predicted scattering halo for Her X-1, MRN dust')
Out[40]:
In [41]:
from astropy.modeling import models
Only difference between this code and plot_radprofile is inclusion of background, and removal of SFLUX term. (Want to model normalized surface brightness).
In [42]:
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 [43]:
from astropy.io import ascii
In [44]:
pfilename = 'plot_radprofile_fitparams.txt'
ptable = ascii.read(pfilename)
pfit = np.zeros(6)
for i in range(6):
pfit[i] = ptable[i][0]
print ptable
In [45]:
fig = plt.figure(figsize=(8,5.3))
# Plot IGM dust halos
plt.plot(alpha, ig_s, 'k-', lw=2, label='Grey IGM')
plt.plot(alpha, i_s, 'k-', lw=1, label='1 um IGM')
## Don't show this stuff for now, because it is superfluous
# (It's not really visible anyways, and could be confused with foreground dust)
# Plot Galhalo dust halos
#plt.plot(ALPHA, sg_sil, 'k--', lw=2, label='Grey CGM')
#plt.plot(ALPHA, sMRN_int, 'r--', lw=2, label='MRN CGM')
# Plot Foreground dust halo
#plt.plot(ALPHA, fg_halo_IGMscreen.intensity, 'r-', alpha=0.8, label='MRN Foreground')
plt.plot(ALPHA, fg_gsil, '--', color='0.3', label='Grey Foreground')
#plt.legend(loc='upper right', frameon=False, fontsize=15)
plt.text(1.5, 0.015, 'PSF\ntemplate', color='b', size=15)
# Plot PSF model
plt.plot(alpha, psf_model(pfit,alpha), '-', lw=1, color='b', alpha=0.8)
plt.loglog()
plt.xlabel('Observation Angle [arcsec]', size=15)
plt.ylabel('Normalized intensity [arcsec$^{-2}$]', size=15)
#plt.savefig('/Users/lia/Dropbox/HEAD2014/mie_scattering_halos.pdf', format='pdf')
Out[45]:
In [46]:
import radprofile as rp
In [47]:
## Open function to load the profile object
def open_profile(fname):
result = rp.Profile()
ffile = fits.open(fname)
pdata = ffile[1].data
result.rleft = pdata['R'][:,0]
result.rright = pdata['R'][:,1]
try:
result.surbri = pdata['SUR_FLUX'] # photons/cm**2/pixel**2/s
result.surbri_err = pdata['SUR_FLUX_ERR']
except:
print "Could not find value:", value
print "Your options are:"
print pdata.columns.names
return result
In [48]:
FNORM = 0.0007349181010116529 # 0.8 - 1.2 keV bin
fprofile = open_profile('herx1_c67_1keV_radprofile_flux.fits')
# Normalize by the flux in this bin
fprofile.divide(FNORM)
In [53]:
fig = plt.figure(figsize=(8,5.3))
# Plot IGM dust halos
plt.plot(alpha, ig_s, 'k-', lw=3, label='Grey IGM')
# Plot Galhalo dust halos
plt.plot(ALPHA, sg_sil, 'k--', lw=3, label='Grey CGM')
plt.plot(ALPHA, sMRN_int, 'r--', lw=3, label='MRN CGM')
# Plot Foreground dust halo
plt.plot(ALPHA, fg_halo_IGMscreen.intensity, 'r-', alpha=0.5, label='MRN Foreground')
plt.plot(ALPHA, fg_gsil, 'k-', alpha=0.5, label='Grey Foreground')
# Her X-1 Scattering halo
plt.plot(alpha, fg_herx1_int, 'm-', lw=2, label='Her X-1 halo')
plt.legend(loc='upper right', frameon=False)
plt.plot(alpha, psf_model(pfit,alpha), 'k-', lw=1, alpha=0.5)
plt.legend(loc='upper right', frameon=False)
arcsec_per_pix = 0.5 # arcsec/pix
plt.plot(alpha, psf_model(pfit,alpha), 'k-', lw=2, alpha=0.5)
plt.errorbar(fprofile.rmid*arcsec_per_pix, fprofile.surbri/arcsec_per_pix**2, \
yerr = fprofile.surbri_err/arcsec_per_pix**2, \
capsize=0, lw=2, color='k', marker='', ls='')
plt.plot(alpha, psf_model(pfit, alpha)+pfit[-1]/FNORM, 'b--', lw=2)
plt.loglog()
plt.xlabel('Observation Angle [arcsec]')
plt.ylabel('Normalized intensity [arcsec$^{-2}$]')
Out[53]:
In [50]:
pfit[-1]
Out[50]:
Minimal flux to get into RASS Bright Sources Catalog: $5 \times 10^{-13}$ cgs
From IGM Dust Calculation (evernote): Chandra flux 0.04 ct/s = $3 \times 10^{-13}$ cgs around 1 keV
Estimate of Chandra count rate: $$ R_{Chandra} = F \ \frac{0.04\ {\rm ct/s}}{3 \times 10^{-13}\ {\rm cgs}} $$
The typical Chandra 0.5-2 keV background is 0.1 ct/s/chip = $4 \times 10^{-7}$ ct/s/arcsec$^2$
In [51]:
RASSBSC_flux = 5.e-13 # cgs
RChandra = 0.04 / 3.e-13 # ct/s per flux (cgs)
SXRB = 4.e-7 # ct/s/arcsec^2
In [52]:
src_rate = RASSBSC_flux * 10 * RChandra # A tiny tail of bright quasars
fig = plt.figure(figsize=(8,5.3))
# Plot the typical background
plt.axhline(SXRB, color='0.8', ls='-', label='')
plt.text(50.0, 6.e-7, '0.5-2 keV background', color='0.5', size=12)
plt.xlim(0.1,1.e3)
# Plot IGM dust halos
plt.plot(alpha, ig_s * src_rate, 'k-', lw=2, label='Grey IGM')
plt.plot(alpha, i_s * src_rate, 'k-', lw=1, label='1 um IGM')
# Plot Foreground dust halo
plt.plot(ALPHA, fg_gsil * src_rate, '--', color='0.3', label='')
plt.text(0.3, 1.e-8, 'Grey foreground dust', color='0.3', size=12)
plt.legend(loc='upper right', frameon=False, fontsize=15)
# Plot PSF model
#plt.plot(alpha, psf_model(pfit,alpha)*src_rate/1.e3, '-', lw=1, color='b', alpha=0.8)
#plt.text(10, 1.e-10, 'PSF\n(dimmed)', color='b', size=15)
plt.loglog()
plt.ylim(1.e-11, 1.e-3)
plt.xlabel('Observation Angle [arcsec]', size=15)
plt.ylabel('Surface Brightness [ct s$^{-1}$ arcsec$^{-2}$]', size=15)
#plt.savefig('/Users/lia/Dropbox/HEAD2014/mie_scattering_halos_SB.pdf', format='pdf')
Out[52]:
In [52]: