PRO dustem_radiance

loadct,39
win=0

dirtmp='~/tmp/'
readcol,dirtmp+'/data/ISRF.DAT',lambda,Inu,skip=7,/silent
readcol,dirtmp+'out/EXT.RES',lam,PAH0a,PAH1a,VSGa,BEa,aSila,PAH0s,PAH1s,VSGs,BEs,aSils,skip=8,/silent
nua = 3e10/(lam*1e-4)
readcol,dirtmp+'out/SED.RES',lame,PAH0e,PAH1e,VSGe,BEe,aSile,/silent
nue = 3e10/(lame*1e-4)

totalabs=PAH0a+PAH1a+VSGa+BEa+aSila
totalsca=PAH0s+PAH1s+VSGs+BEs+aSils
totalext=totalabs+totalsca
totalemis=PAH0e+PAH1e+VSGe+BEe+aSile

mH = 1.67e-24

n=n_elements(lam)
isrf = fltarr(n)
for i=0, n-1 do isrf(i)=interpol(Inu,lambda,lam(i))
nua1 = 3e10 / 0.1e-4
nua2 = 3e10 / 20e-4
; Absorption
AaSil 	= integral(nua,aSila*isrf,nua1,nua2)*mH
ABE 	= integral(nua,BEa*isrf,nua1,nua2)*mH
AVSG 	= integral(nua,VSGa*isrf,nua1,nua2)*mH
APAH 	= integral(nua,(PAH1a+PAH0a)*isrf,nua1,nua2)*mH
Abs 	= integral(nua,totalabs*isrf,nua1,nua2)*mH
Ext 	= integral(nua,totalext*isrf,nua1,nua2)*mH

;Scattering
SaSil 	= integral(nua,aSils*isrf,nua1,nua2)*mH
SBE 	= integral(nua,BEs*isrf,nua1,nua2)*mH

; Emission
nue1 = 3e10/0.1e-4
nue2 = 100.d9
EPAH  	= integral(nue,(PAH0e+PAH1e)/nue,nue1,nue2)
EVSG 	= integral(nue,VSGe/nue,nue1,nue2)
EBE 	= integral(nue,BEe/nue,nue1,nue2)
EaSil 	= integral(nue,aSile/nue,nue1,nue2)
Eemis 	= integral(nue,totalemis/nue,nue1,nue2)

print
print,"DUSTEM MC11"
print,"==========="
;print,"Energy absorbed (W/m2/1d20H over 4PI steradians)"
print,"Energy reemitted by dust (W/m2/sr/1d20H)"
coef = 4*!pi / 1d17
print
print,"aSil	"
print,"	ABS	= ",Aasil/coef
print,"	EMI 	= ",Easil/coef
print,"BE	"
print,"	ABS	= ",ABE/coef
print,"	EMI	= ",EBE/coef
print,"VSG	  "
print,"	ABS	= ",AVSG/coef
print,"	EMI	= ",EVSG/coef
print,"PAH	"
print,"	ABS	= ",APAH/coef
print,"	EMI	= ",EPAH/coef
print 
print,"Total BG"
print,"	ABS	= ",(Aasil+ABE)/coef
print,"	EMIS	= ",(EaSil+EBE)/coef
print 
print,"Total "
print,"	ABS	= ",(Aasil+ABE+AVSG+APAH)/coef
print,"	EMI	= ",Eemis/coef
print
print,"ABS/EXT ",Abs/Ext
print,"ABS/EXT BG only ",(AaSil+ABE)/(AaSil+SaSil+ABE+SBE)
;print,"Total ext (for check with Mathis only)",Ext/coef
print

win+=1
window,win
plot,lam,totalext*isrf*1d16,/xlog,xr=[0.1,5],/ylog,yr=[1e-1,1e2]
oplot,lam,aSila*isrf*1d16,col=250
oplot,lam,BEa*isrf*1d16,col=150
oplot,lam,VSGa*isrf*1d16,col=100
oplot,lam,(PAH0a+PAH1a)*isrf*1d16,col=30



;P06b Eq. 11
;***********
; PIP82
Tobs = 19.9
beta = 1.65
tauperH = 7.1e-27
sigmaS=5.67e-8
nu0=353d9
k=1.38e-23
h=6.63e-34
Rdata=tauperH*sigmaS/!pi*Tobs^4*(k*Tobs/h/nu0)^beta*factorial(3.+beta)/factorial(3.)*riemannzeta(4+beta)/riemannzeta(4.)*1d20
; Color correction ?
cc = 1.11
Rdata *= cc

print
print,"P06b"
print,"===="
print,"Radiance from P06b :",Rdata," W/m2/sr/1d20H"

dir='~/Software/DUSTEM'
ext=read_xcat(dir+'/Data/EXTs/Mathis1990_DISM.xcat',/silent)
ext.spec *= 1.d-21
ext.error *= 1.d-21
isrf=read_xcat(dir+'/Data/ISRFs/ISRF.xcat',/silent)

;wav=isrf.lambda
;n=n_elements(wav)
;ns = fltarr(n)
;for i=0, n-1 do ns(i)=interpol(ext.spec,ext.wave,wav(i))
;
;E = integral(wav*1e-4,ns*isrf.Inu,0.1e-4,20.e-4)
;print
;print,"Mathis Extinction curve"
;print,"======================="
;print,"Total absorbed energy = ",E*1.d17*1.d20," W/m2/1d20H"
;print,"Total energy reemitted = ",E*1.d17*1.d20/(4*!pi),"W/m2/sr/1d20H"

;win+=1
;window,win
;plot,wav,1.d37*ns*isrf.inu,xr=[0.1,20],/xlog

stop

END