In [1]:
import numpy as np
import matplotlib.pyplot as plt

font = {'size'   : 15}
matplotlib.rc('font', **font)

In [2]:
cp /Users/lia/Academic/CygX-3/6601_paper1/best_fits.py /Users/lia/Dropbox/notebooks/cygx3_redo/best_fits.py

Calculate scattering opacity for different scattering models and distributions


In [3]:
import sigma_scat as ss
import dust
import constants as c

In [4]:
NH, d2g   = 1.e21, 0.009
dust_mass = NH * c.mp() * d2g 
ERANGE    = np.power( 10.0, np.arange(-0.6,1.0,0.005) )

Set up grain size distributions and materials


In [5]:
RHO_SIL, RHO_GRA, RHO_AVG = 3.8, 2.2, 3.0 # g cm^-3; see Draine book

In [6]:
MRN_RANGE = np.arange(0.005,0.25001,0.05)
MRN_sil   = dust.Dustdist( rad=MRN_RANGE, p=3.5, rho=RHO_SIL )
MRN_gra   = dust.Dustdist( rad=MRN_RANGE, p=3.5, rho=RHO_GRA )
MRN_avg   = dust.Dustdist( rad=MRN_RANGE, p=3.5, rho=RHO_AVG )

In [7]:
# choices motivated by crazy fit
BIG_RANGE = np.arange(0.005, 1.5, 0.05)
BIG_sil   = dust.Dustdist( rad=BIG_RANGE, p=3.5, rho=RHO_SIL )
BIG_gra   = dust.Dustdist( rad=BIG_RANGE, p=3.5, rho=RHO_GRA )
BIG_avg   = dust.Dustdist( rad=BIG_RANGE, p=3.5, rho=RHO_AVG )

In [8]:
## Rayleigh-Gans plus Drude approximation
RGD_mrn = ss.Kappascat( E=ERANGE, dist=dust.Dustspectrum(rad=MRN_avg, md=dust_mass), scatm=ss.makeScatmodel('RG','Drude') )
RGD_big = ss.Kappascat( E=ERANGE, dist=dust.Dustspectrum(rad=BIG_avg, md=dust_mass), scatm=ss.makeScatmodel('RG','Drude') )

In [9]:
## Mie scattering for the small grain MRN distribution
Mie_mrn_sil = ss.Kappascat( E=ERANGE, dist=dust.Dustspectrum(rad=MRN_sil, md=dust_mass), scatm=ss.makeScatmodel('Mie','Silicate') )
Mie_mrn_gra = ss.Kappascat( E=ERANGE, dist=dust.Dustspectrum(rad=MRN_gra, md=dust_mass), scatm=ss.makeScatmodel('Mie','Graphite') )

In [10]:
## Mie scattering for the slightly bigger grain distribution
Mie_big_sil = ss.Kappascat( E=ERANGE, dist=dust.Dustspectrum(rad=BIG_sil, md=dust_mass), scatm=ss.makeScatmodel('Mie','Silicate') )
Mie_big_gra = ss.Kappascat( E=ERANGE, dist=dust.Dustspectrum(rad=BIG_gra, md=dust_mass), scatm=ss.makeScatmodel('Mie','Graphite') )

Plot it


In [15]:
MD = dust_mass * 10.0 # for N_H = 1.e22

plt.plot( RGD_mrn.E, RGD_mrn.kappa * MD, '0.4', lw=2, label='RG-Drude' )
plt.plot( Mie_mrn_sil.E, Mie_mrn_sil.kappa * MD, 'g', lw=2, label='Mie-Silicate' )
plt.plot( Mie_mrn_gra.E, Mie_mrn_gra.kappa * MD, 'b', lw=2, label='Mie-Graphite' )

plt.legend( loc='upper right', fontsize=12 )

plt.plot( RGD_big.E, RGD_big.kappa * MD, '0.4', lw=1, ls='-' )
plt.plot( Mie_big_sil.E, Mie_big_sil.kappa * MD, 'g', lw=1, ls='-' )
plt.plot( Mie_big_gra.E, Mie_big_gra.kappa * MD, 'b', lw=1, ls='-' )

plt.loglog()
plt.xlim(0.3,10)
plt.ylim(1.e-2,10.0)
plt.xlabel( "Energy [keV]", size=15 )
plt.ylabel( r"Scattering Opacity [$\tau$ per $N_{\rm H}/10^{22}$]", size=15 )

plt.text( 0.5, 0.1, '$0.25\ \mu m$\ncut-off', size=12)
plt.text( 1.2, 3.0, '$1.5\ \mu m$\ncut-off', size=12)
#plt.savefig('/Users/lia/Dropbox/Writing/CygX3-paper1/figure5.pdf', format='pdf')


What is $\kappa_{\rm keV}$ for MRN RG-Drude?


In [12]:
from scipy.interpolate import interp1d

In [13]:
kappa_interp = interp1d( RGD_mrn.E, RGD_mrn.kappa )
kappa_1keV   = kappa_interp(1.0)
print "Tau (1 keV) for RG-Drude with MRN:", kappa_1keV


Tau (1 keV) for RG-Drude with MRN: 3326.17984598

In [14]:
plt.plot( RGD_mrn.E, RGD_mrn.kappa )
plt.plot( 1.0, kappa_1keV, 'ko', markersize=10 )
plt.loglog()


Out[14]:
[]

In [14]: