Viewing contents of file '../idllib/astron/contrib/varosi/vlibm/allpro/dust_emission.pro'
;+
; NAME:
; dust_emission
; PURPOSE:
; Compute dust temperature distributions and infrared emission spectrum
; from the absorbed Luminosities as computed by MGEP_RAD_TRANS.
; All input luminosities and output spectra are in solar units.
;
; CALLING:
; dust_IR_struct = dust_emission( Lumins, KAPPA_DUST= , ... )
;
; INPUT:
; Lumins = structure of all Luminosities in model,
; as computed and returned by procedure MGEP_Rad_Trans.
;
; required tags are: Lumins.Lum_dust: fltarr( ncomp, nphase )
; Lumins.dust_comp: strarr( ncomp )
; Lumins.dust_phase: strarr( nphase )
;
; KEYWORD INPUTS:
;
; KAPPA_DUST = structure variable with dust absorption/scattering coeffs.:
; KAPPA_ICM_DUST = optional, structure with ICM dust abs. & scat. coeffs.
;
; Required Structure Tags:
; Dust_Kappa.wavelen (Nwave) (wavelengths, in microns)
; Dust_Kappa.Kabs (Nwave,Ncomp) (absorption per gram)
; Dust_Kappa.Tmax (Ncomp) (sublimation temperatures)
;
; DUST_MASS = matrix (Ncomp,Nphase) of total dust mass
; for each component & phase of medium, in solar mass units.
;
; /PLOT : plot intermediate results.
; /VERBOSE : print informational messages about results.
; /FULL_STRUCT : return all infrared emission spectra in structure
; (default is to only return dust temperatures).
;
; POWER_LAW_TDIST = integer array of Nphase elements, giving power law
; flag for ICM and clumps in the case of a two-phase medium.
; Values of 0 : assume single dust temperature for that phase.
; Values of 1 : assume that the dust temperature follows
; an inverse power law distribution, like around a point source,
; and in such case the following indices will be used:
;
; DINDEX = radial inverse power law index of dust density variation,
; (e.g. DINDEX=2 implies density goes as 1/r^2, default = 0).
; Must be an array of Nphase elements, giving density index for
; ICM and clumps in the case of a two-phase medium.
;
; LINDEX = radial inverse power law index of absorbed luminosity,
; (default = 2, which is for optically thin dust).
; Must be vector of length Ncomponent, giving absorbed luminosity
; index for graphite and silicates, etc.
;
; EMISSINX = matrix of size (2,Ncomp) giving emissivity indices.
; default: [ [2.,1.], $ ;graphite emissivity indices < & > Tcut.
; [2.,0.] ] ;silicate emissivity indices < & > Tcut.
;
; TCUT = vector (Ncomp) of temperatures at which emissivity index
; transition occurs, default = 70K (graphite), 150K (silicates).
;
; TMAX = matrix (Ncomp,Nphase) of dust sublimation temperatures
; to use for the power law temperature distribution.
; (default is to use values from Dust_Kappa.Tmax).
;
; OUTPUTS:
; Function returns a structure containing following tags:
; Wavelen = wavelengths for spectra, in microns.
; Tdust = single dust temperature.
; Dspec1 = emission spectrum from dust at single temperature.
; Tdmin = minimum dust temperature of power-law distributions.
; Tdmax = maximum dust temperature of power-law distributions.
; Dspec_Td = emission spectrum from dust having
; a power law distribution of temperatures.
;
; If FULL_STRUCT is negative then Dspec1 or Dspec_Td are replace by Dspec.
;
; COMMON BLOCKS:
; common dust_emit_spec, dust_spec ;to get results of function Dust_Temp
; common dust_Tmin_spec, dtd_spec ;to get results of function Dust_Tmin
; EXTERNAL CALLS:
; function Dust_Temp
; function Dust_Tmin
; function Trapez
; PROCEDURE:
; HISTORY:
; Written: Frank Varosi, HSTX @ NASA/GSFC, 1995.
; F.V. 1997, compute temperature of multi-phase medium: ICM and clumps.
; F.V. 1999, added keyword KAPPA_ICM_DUST for different dust in ICM.
; F.V. 2000, POWER_LAW_TD option now handles ICM and clumps seperately.
;-
function dust_emission, Lumins, DUST_MASS=Mdust, $
KAPPA_DUST=Dust_Kappa, $
KAPPA_ICM_DUST=Dust_ICM_Kappa, $
FULL_STRUCTURE=full_struct, PLOTIT=plotit, $
POWER_LAW_TDIST=powTdist, VERBOSE=verb, $
DINDEX=dinx, LINDEX=Linx, TMAX=Tmax, $
EMISSINX=eminx, TCUT=Tcut
common dust_emit_spec, dust_spec ;to get results of function Dust_Temp
common dust_Tmin_spec, dtd_spec ;to get results of function Dust_Tmin
sz = size( Lumins.Lum_Dust )
ndc = sz(1)
if sz(0) eq 2 then begin
nphase = sz(2)
dust_phase = Lumins.dust_phase
endif else begin
nphase=1
dust_phase = "homog"
endelse
if N_elements( powTdist ) ne nphase then powTdist = replicate(0,nphase)
if N_elements( dinx ) ne nphase then dinx = [0.,0.]
if N_elements( Linx ) ne ndc then Linx = [2.,2.]
if N_elements( Tcut ) ne ndc then Tcut = [ 70., 150. ]
sz = size( eminx )
if (sz(0) ne 2) or (sz(1) ne ndc) or (sz(2) ne 2) then $
eminx = [ [2.,1.], $ ;graphite emissivity indices < & > Tcut.
[2.,0.] ] ;silicate emissivity indices < & > Tcut.
Tdust = fltarr( ndc, nphase )
Tdmin = fltarr( ndc, nphase )
Tdav = fltarr( ndc, nphase )
Lum_spec = fltarr( ndc, nphase )
sz = size( Tmax )
if (sz(0) ne 2) or (sz(1) ne ndc) or (sz(1) ne nphase) then begin
Tdmax = Dust_Kappa.Tmax
Tdmax = [[Tdmax],[Tdmax]]
endif else Tdmax = Tmax
wave = Dust_Kappa.wavelen
wir = where( wave GT 0.5, Nwave ) ;ignore shorter wavelengths.
wave = wave(wir)
wang = wave * 1e4
Dspec = fltarr( Nwave, ndc, nphase )
Dspec_Td = fltarr( Nwave, ndc, nphase )
md = Mdust * (!cv.Msun/!cv.Lsun)
if keyword_set( plotit ) then begin
w = where( Lumins.Lum_dust GT 1e-33, nw )
if (nw gt 2) then !P.multi=[0,ndc,nphase] else !P.multi=[0,1,nw]
endif
for k=0,nphase-1 do begin
ptdist = powTdist(k)
for i=0,ndc-1 do begin
Ldust = Lumins.Lum_dust(i,k)
if (Ldust GT 1e-33) AND (md(i,k) GT 1e-33) then begin
if (k eq 0) and N_struct( Dust_ICM_Kappa ) then $
Kemit = Dust_ICM_Kappa.Kabs(wir,i) * md(i,k) $
else Kemit = Dust_Kappa.Kabs(wir,i) * md(i,k)
Tdust(i,k) = Dust_Temp( Ldust, wave, Kemit )
Dspec(*,i,k) = dust_spec
if keyword_set( ptdist ) then begin
Tdmin(i,k) = Dust_Tmin( Ldust, Tdust(i,k), Tcut(i), Tdmax(i,k),$
WAV=wave, KEMI=Kemit, $
EM=eminx(*,i), D=dinx(k), L=Linx(i) )
Dspec_Td(*,i,k) = dtd_spec
Lum_spec(i,k) = Trapez( Dspec_Td(*,i,k) > 1e-37, wang )
Tdav(i,k) = dTemp_Prob( [Tdmin(i,k),Tdmax(i,k)], TCUT=Tcut(i),$
/AV, EM=eminx(*,i),D=dinx(k),L=Linx(i) )
endif else Lum_spec(i,k) = Trapez( Dspec(*,i,k) > 1e-37, wang )
if keyword_set( plotit ) then begin
get_window,2,/SHOW
convMW = !cv.Lsun * 1d-9/!cv.c * wave^2
if keyword_set( ptdist ) then begin
Dsp = Dspec_Td(*,i,k) * convMW
ptit = Lumins.dust_comp(i) + ": " $
+ dust_phase(k) + ": Tmin=" + $
string( Tdmin(i,k), F="(I3)" ) + $
"K Tav=" + string( Tdav(i,k), F="(I3)" )+"K"
endif else begin
Dsp = Dspec(*,i,k) * convMW
ptit = Lumins.dust_comp(i) + ": " $
+ dust_phase(k) + $
": Tbb=" + string( Tdust(i,k), F="(I3)" )+"K"
endelse
ymax = ceil( alog10( max( Dsp ) ) )
plot, wave, Dsp, XTIT="Wavelength (microns)",/YTYP, /YSTY, $
YTIT="MW/m^2/Hz",YR=10^[ymax-3.0,ymax], $
XR=[1,1000],/XTYPE, TIT=ptit
oplot, wave, Dspec(*,i,k)*convMW, Line=2
empty
endif
endif
endfor
endfor
if keyword_set( plotit ) then !P.multi=0
if keyword_set( verb ) then begin
print," "
formg = "(A10,3G11.3)"
print," Mass: " + strconcat( " " + Lumins.dust_comp )
for k=0,nphase-1 do print,dust_phase(k),Mdust(*,k), F=formg
print,"total", total( Mdust, 2 ), total( Mdust ), F=formg
formi = "(A10,2(I9,' K'))"
print," Tdust:"
for k=0,nphase-1 do print,dust_phase(k),Tdust(*,k),F=formi
print," Tdmin:"
for k=0,nphase-1 do print,dust_phase(k),Tdmin(*,k),F=formi
print," Tdavg:"
for k=0,nphase-1 do print,dust_phase(k),Tdav(*,k),F=formi
print," Tdmax:"
for k=0,nphase-1 do print,dust_phase(k),Tdmax(*,k),F=formi
print," Lumin Emit:"
for k=0,nphase-1 do print,dust_phase(k),Lum_spec(*,k), F=formg
print,"total",total( Lum_spec, nphase<2 ), $
total( Lum_spec ), F=formg
endif
if keyword_set( full_struct ) then begin
if full_struct LT 0 then begin
if keyword_set( ptdist ) then begin
return, { Tdust:Tdust, Tdmin:Tdmin, Tdmax:Tdmax, Tdav:Tdav, $
Lum_spec:Lum_spec, Mdust:Mdust, $
Lum_index:Linx, Den_index:dinx, $
Emissindex:eminx, Tcut:Tcut, $
WaveLen:wave, Dspec:Dspec_Td }
endif else begin
return, { Tdust:Tdust, Tdmin:Tdmin, Tdmax:Tdmax, Tdav:Tdav, $
Lum_spec:Lum_spec, Mdust:Mdust, $
Lum_index:Linx, Den_index:dinx, $
Emissindex:eminx, Tcut:Tcut, $
WaveLen:wave, Dspec:Dspec }
endelse
endif else begin
return, { Tdust:Tdust, Tdmin:Tdmin, Tdmax:Tdmax, Tdav:Tdav, $
Lum_spec:Lum_spec, Mdust:Mdust, $
Lum_index:Linx, Den_index:dinx, $
Emissindex:eminx, Tcut:Tcut, $
WaveLen:wave, Dspec1:Dspec, Dspec_Td:Dspec_Td }
endelse
endif else begin
return, { Tdust:Tdust, Tdmin:Tdmin, Tdmax:Tdmax, Tdav:Tdav, $
Lum_spec:Lum_spec, Mdust:Mdust, $
Lum_index:Linx, Den_index:dinx, $
Emissindex:eminx, Tcut:Tcut }
endelse
end