Viewing contents of file '../idllib/astron/contrib/varosi/vlibm/allpro/mg_ep_2phase.pro'
;+
; NAME:
; MG_EP_2phase
; PURPOSE:
; Approximate the radiative transfer of photons in a
; two-phase clumpy medium (dust & gas) of spherical/disk geometry, by
; calling MGEP_Rad_Trans, and then compute equilibrium dust temperatures
; and infrared emission using function Dust_Emission.
; Can optionally model an environment where the SED of stars in clumps
; is different that stars in ICM (see keywords SED_CLUMPS and SRC_CLUMPS),
; or case of the dust in ICM different than in clumps (KAPPA_ICM_DUST).
; CALLING:
; MG_EP_2phase, SED_wave, SED_source, SED_esc, Lumins, dust_emiss
; INPUTS:
; SED_wave = array of SED wavelengths, in microns.
;
; SED_source = Spectral Energy Distribution (SED) of source emission
; vs. wavelength, in solar luminosities per Angstrom (Lsun/A).
; OUTPUTS:
; SED_esc = escaping SED vs. wavelength (Lsun/A).
;
; Lumins = structure variable accounting of all Luminosities emitted,
; escaped, and absorbed in dust (tag LUM_DUST) or gas.
;
; dust_emiss = structure variable containing the temperature of dust
; in each phase, and optionally the infrared emission spectra.
;
; 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 mass)
; Dust_Kappa.Ksca (Nwave,Ncomp) (scattering per mass)
; Dust_Kappa.gcos (Nwave,Ncomp) (asymmetry params)
; Dust_Kappa.dust_comp (Ncomp) (eg. ["graphite","silicates"])
; Dust_Kappa.Tmax (Ncomp) (sublimation temperatures)
;
; (note: the following 3 structures are defined by IDL> @def_mgep)
;
; GEOMETRY = structure variable:
; geometry.R = radius of sphere/disk/cylinder (in pc).
; geometry.H = height of disk/cylinder (if H=0 then sphere).
;
; CLUMPY = structure variable:
; clumpy.fill_fact (filling factor of clumps)
; clumpy.Clump_ICM_Ratio (clump/ICM density ratio)
; clumpy.Rcl (radius of clumps, in pc)
;
; MEDIUM = structure variable, specify:
; medium.fsil (ratio of silicates to total dust mass)
; medium.Dust_Mass_Dens (gm/cm^3)
; medium.gas_dens (n(H)/cm^3)
;
; SOURCE = string (upper or lowercase) indicating type of photon source:
; "C" = central point source,
; "X" = uniformly illuminating external source,
; "U" = uniformly distributed internal sources (default case).
;
; SED_CLUMPS = optional SED of stellar emission in clumps,
; and then input array SED_source is for stars in the ICM.
; Output array SED_esc and structure Lumins is then the
; sum of effects from SED_source and SED_clumps.
;
; SRC_CLUMPS = string indicating type of source emission in clumps:
; "C" = central point source,
; "U" = uniformly distributed sources (default case).
; Must specify SED_CLUMPS and then SOURCE=type of emission in ICM.
;
; /PLOT : plot intermediate results.
; /VERBOSE : print informational messages about results.
; /FULL_STRUCT : include all results in output structure variable "Lumins"
;
; KEYWORD OUTPUTS (optional):
;
; SED_ABS_DUST = 3D matrix of flux vs. wavelength absorbed by each
; component and phase of the dust: (Nwave,Ncomp,Nphase) (Lsun/A).
; SED_ESC_CL = escaping flux from source just in clumps.
; LUM_CLUMPS = absorbed luminosities (structure) from source in clumps.
;
; EXTERNAL CALLS:
; pro MGEP_Rad_Trans
; pro add_struct
; function Dust_Emission
; HISTORY:
; Written: Frank Varosi, RSTX @ NASA/GSFC, 1997
; F.V. 1999, added keyword KAPPA_ICM_DUST for different dust in ICM.
; F.V. 2000, added option for separate SED for stars in clumps.
; F.V. 2000, added keyword output SED_ABS_DUST.
;-
pro MG_EP_2phase, SED_wave, SED_source, SED_esc, Lumins, dust_emiss, $
SED_ABS_DUST=SED_abs_dust, SOURCE=srctyp, $
SED_CLUMPS=SED_clumps, SRC_CLUMPS=srctcl, $
SED_ESC_CL=SED_cl_esc, LUM_CLUMPS=Lumcl, $
KAPPA_DUST=Dust_Kappa, KAPPA_ICM_DUST=Dust_ICM_Kappa, $
FULL_STRUCT=fulls, VERBOSE=verbose, PLOTIT=plotit,$
POWER_LAW_TDIST=ptdist, TDIST_PARM=tdistp, TMAX=Tmax, $
GEOMETRY=geometry, CLUMPY=clumpy, MEDIUM=medium
if( geometry.H LE 0 ) then begin ;case of just simple sphere:
geometry.Volume = (4*!pi/3) * geometry.R^3
geometry.Reff = geometry.R
endif else begin ;case of disk/cylinder
geometry.Volume = geometry.H * !pi * geometry.R^2
geometry.Reff = 1.5 * geometry.R * geometry.H /(geometry.R + geometry.H)
endelse
ffc = clumpy.fill_fact
clumpy.Ncl = round( -alog( 1-ffc ) * 3*geometry.Volume/(4*!PI*clumpy.Rcl^3) )
if max( tag_names( clumpy ) EQ "CLUMP_ICM_RATIO" ) then $
clumpy.dens_ratio = 1./clumpy.Clump_ICM_ratio
medium.dust_mass_dens_CL = $
medium.dust_mass_dens/( ffc + (1-ffc)*clumpy.dens_ratio )
medium.dust_mass_dens_ICM = medium.dust_mass_dens_CL * clumpy.dens_ratio
dfac = !cv.pc * !cv.pc * (!cv.pc/!cv.msun)
medium.dust_mass_CL = medium.dust_mass_dens_CL * ffc * geometry.Volume * dfac
medium.dust_mass_ICM = medium.dust_mass_dens_ICM *(1-ffc)* geometry.Volume *dfac
medium.dust_mass = geometry.Volume * medium.dust_mass_dens * dfac
medium.dust_gas_ratio = medium.dust_mass_dens/(medium.gas_dens*1.4*!cv.mh)
medium.gas_dens_CL = medium.dust_mass_dens_CL/(medium.dust_gas_ratio*1.4*!cv.mh)
medium.gas_dens_ICM = medium.gas_dens_CL * clumpy.dens_ratio
medium.fsil = ( medium.fsil < 1 ) > 0
medium.fc = (1-medium.fsil)
dust_fracs = [ medium.fc, medium.fsil ]
dmass_matrix = dust_fracs # [ medium.dust_mass_ICM, medium.dust_mass_CL ]
if max( tag_names( medium ) eq "FSIL_ICM" ) then begin
medium.fsil_ICM = ( medium.fsil_ICM < 1 ) > 0
medium.fc_ICM = (1-medium.fsil_ICM)
dust_fracs_ICM = [ medium.fc_ICM, medium.fsil_ICM ]
dmass_matrix = [[ dust_fracs_ICM * medium.dust_mass_ICM ], $
[ dust_fracs * medium.dust_mass_CL ] ]
endif
if N_elements( srctyp ) ne 1 then srctyp = "U" $
else if strlen( strtrim( srctyp, 2 ) ) LE 0 then srctyp = "U"
if keyword_set( verbose ) then begin
print," "
print_struct, clumpy, FORM=['g','11','3']
print," "
print_struct, geometry, FORM=['g','11','3']
print," "
print_struct, medium, TR=[0,6], FORM=['g','11','3']
print," "
print_struct, medium, TR=[7,20], FORM=['g','11','3']
print," "
print," SOURCE: ",srctyp
endif
MGEP_Rad_Trans, SED_wave, SED_source, SED_esc, Lumins, FDUST=frac_abs_dust, $
FILL_FACT=ffc, RATIO=clumpy.dens_ratio, $
CLUMP_RADIUS = clumpy.Rcl * !cv.pc, $
RADIUS = geometry.Reff*!cv.pc, SOURCE=srctyp, $
DUST_MASS_DEN=medium.dust_mass_dens, $
KAPPA_DUST=Dust_Kappa, KAPPA_ICM_DUST=Dust_ICM_Kappa, $
DUST_FRAC=dust_fracs, DUST_ICM_FRAC=dust_fracs_ICM, $
GAS_DEN=medium.gas_dens, GAS_TEMP=medium.gas_Temp, $
PLOTIT=plotit, VERB=verbose, FULL=fulls
if Arg_Present( SED_abs_dust ) then begin
SED_abs_dust = make_array( SIZE=size( frac_abs_dust ) )
SED_abs = SED_source - SED_esc
sz = size( SED_abs_dust )
ndc = sz(2)
nphase = sz(3)
for j=0,nphase-1 do begin
for i=0,ndc-1 do SED_abs_dust(*,i,j) = frac_abs_dust(*,i,j) * SED_abs
endfor
frac_abs_dust = 0
endif
if N_elements( SED_clumps ) eq N_elements( SED_wave ) then begin
if N_elements( srctcl ) ne 1 then srctcl = "U"
print," "
print," SOURCE in Clumps: ",srctcl
MGEP_Rad_Trans, SED_wave, SED_clumps, SED_cl_esc, Lumcl, $
FILL_FACT=ffc, RATIO=clumpy.dens_ratio, $
CLUMP_RADIUS = clumpy.Rcl * !cv.pc, $
DUST_MASS_DEN=medium.dust_mass_dens, $
RADIUS = geometry.Reff*!cv.pc, SOURCE="M"+srctcl, $
KAPPA_DUST=Dust_Kappa, KAPPA_ICM_DUST=Dust_ICM_Kappa, $
DUST_FRAC=dust_fracs, DUST_ICM_FRAC=dust_fracs_ICM, $
GAS_DEN=medium.gas_dens, GAS_TEMP=medium.gas_Temp, $
PLOT=plotit, VERB=verbose, FULL=fulls, FD=frac_abs_dust
if Arg_Present( SED_abs_dust ) then begin
SED_abs = SED_clumps - SED_cl_esc
for j=0,nphase-1 do begin
for i=0,ndc-1 do SED_abs_dust(*,i,j) = $
SED_abs_dust(*,i,j) + frac_abs_dust(*,i,j) * SED_abs
endfor
frac_abs_dust = 0
endif
SED_esc = SED_esc + SED_cl_esc
add_struct, Lumins, Lumcl, TAG_MATCH="LUM_"
add_struct, Lumins, Lumcl, TAG_MATCH="SED_E"
add_struct, Lumins, Lumcl, TAG_MATCH="SED_S"
endif
if N_struct( tdistp ) eq 1 then begin
Lindex = [ tdistp.Lindex_C, tdistp.Lindex_Sil ]
dindex = [ tdistp.dindex_ICM, tdistp.dindex_CL ]
endif
if N_elements( ptdist ) ne 2 then begin
ptdist = replicate( strupcase( srctyp ) EQ "C", 2 )
if N_elements( SED_clumps ) eq N_elements( SED_wave ) then begin
ptdist(1) = strupcase( srctcl ) EQ "C"
endif
endif
dust_emiss = Dust_Emission( Lumins, TMAX=Tmax, DUST_MASS=dmass_matrix, $
KAPPA_DUST=Dust_Kappa, $
KAPPA_ICM_DUST=Dust_ICM_Kappa, $
PLOTIT=plotit, FULL=fulls, VERB=verbose, $
LINDEX=Lindex, DINDEX=dindex, POWER=ptdist )
end