;+
; NAME:
;	MGEP_Rad_Trans
; PURPOSE:
;	Approximate the radiative transfer of photons emitted in or impacting
;	a two-phase clumpy medium (dust & gas) of spherical/disk geometry,
;	using the Mega-Grains approximation of Hobson & Padman,
;	combined with the escape/interaction probability formulae:
;	Osterbrock-Lucy E.P. for uniformly distributed internal source,
;	Varosi's interaction prob. for uniformly illuminating external source,
;	or Varosi's E.P. approximation for central source.
;	See Varosi & Dwek 1999, ApJ 523, 265 for complete discussion.
;	An approximation for the ionization of Hydrogen, and resultant
;	heating of dust by Lyman-alpha photons is also attempted.
;	See keyword SOURCE for the 5 types of photon sources modeled,
;	(3 standard types and 2 special for sources only in clumps).
; CALLING:
;	MGEP_Rad_Trans, SED_wave, SED_source, SED_esc, Lumins
;
; INPUTS:
;	SED_wave = array of SED wavelengths, in microns.
;		The dust abs. & scat. coeffs. (input from KAPPA_DUST keyword)
;		are interpolated onto this SED wavelength grid.
;
;	SED_source = Spectral Energy Distribution (SED) of source emission
;		vs. wavelength, in solar luminosities per Angstrom (Lsun/A).
;
; OUTPUTS:
;	SED_esc = escaping spectral energy distribution vs. wavelength (Lsun/A).
;
;	Lumins = structure variable accounting of all Luminosities emitted,
;		escaped, and absorbed in dust/gas.
;
; KEYWORD INPUTS:	(ICM = interclump medium)
;
;	KAPPA_DUST = structure variable with dust absorption/scattering coeffs.
;	KAPPA_ICM_DUST = optional, structure with ICM dust abs. & scat. coeffs.
;
;	DUST_MASS_DEN = equivalent homogeneous mass density of dust (gm/cm^3).
;	DUST_FRAC = mass fraction of each component of the dust.
;	DUST_ICM_FRAC = optional, dust component fractions in the ICM.
;
;	GAS_DEN = homogeneous density (n(H)/cm^3) of atomic hydrogen.
;	GAS_TEMP = temperature (K) of atomic hydrogen.
;	RADIUS = effective radius of sphere/disk containing clumpy medium.
;
;	FILL_FACT = volume filling factore of spherical clumps,
;	CLUMP_RADIUS = radius of a spherical clump, same units as RADIUS.
;	RATIO_DICM_DCL = ratio of ICM density to clump density.
;
;	SOURCE = string (upper or lowercase) indicating source of photons:
;		"C" : case of central point source,
;		"X" : case of uniformly illuminating external source,
;		"U" = uniformly distributed internal sources (default case).
;		"MC" : central source only in clumps (no emission in ICM),
;		"MU" : uniform source only in clumps (no emission in ICM).
;
;	/PLOT : plot intermediate results.
;	/VERBOSE : print informational messages about results.
;	/FULL_STRUCT : include all results in output structure variable "Lumins"
;
; OPTIONAL KEYWORD OUTPUTS:
;	EPWAVE = computed photon escape probability vs. wavelength.
;	EPLYA = escape probability of Lyman alpha accounting for resonant
;		scattering in a dusty medium.
;	TAU_EFF = the effective optical depth of medium at each wavelength.
;	ALBEDO_EFF = the effective albedo of clumpy medium at each wavelength.
;	FDUST = fractions of photons absorbed by each dust component and phase.
;	FGAS = fractions of photons absorbed by each gas phase at each wavelen.
; EXTERNAL CALLS:
;	pro Mega_Grains
;	function Mega_Grain_Abs
;	function Finterpol
;	function ismeuv
;	function Tau_Eff_Scat
;	function Escape_Prob
;	function Pinteract
;	function Trapez
;	function ep_Lya
; HISTORY:
;	Written: Frank Varosi, HSTX @ NASA/GSFC, 1995
;	F.V. 1997, account for absorbtion by both ICM and clumps.
;	F.V. 1999, added keyword options for different dust in ICM.
;	F.V. 2000, added source types "MC" and "MU" (source just in clumps).
;-

pro MGEP_Rad_Trans, SED_wave, SED_source, SED_esc, Lumins, RADIUS=Rs, $
			KAPPA_DUST=Dust_Kappa, KAPPA_ICM_DUST=Dust_ICM_Kappa, $
			DUST_MASS_DEN=dust_dens, DUST_FRACTIONS=dust_fracs, $
			DUST_ICM_FRACTIONS=dust_ICM_fracs, $
			FILL_FACT=ffcl, CLUMP_RADIUS=Rcl, RATIO_DICM_DCL=drat, $
			GAS_DENSITY=gas_dens, GAS_TEMP=gas_T, PLOTIT=plotit,$
			EPWAVE=esc_prob, EPLYA=epLya, TAU_EFF=Tau_Eff, $
			ALBEDO_EFF=albedo_eff, FDUST=frac_dust, FGAS=frac_gas, $
			FULL_STRUCT=full_struct, VERBOSE=verbose, SOURCE=srctype

	if N_struct( Dust_Kappa ) ne 1 then begin
		message,"keyword KAPPA_DUST must be structure variable",/INFO
		stop
	   endif

	sz = size( Dust_Kappa.Kabs )
	ndc = sz(2)			;usually 0=graphite, 1=silicates.
	nphase = 2			;0=ICM, 1=clumps.

	if N_elements( dust_fracs ) NE ndc then begin
		message,"keyword DUST_FRACTIONS must be a " + $
			strtrim( ndc, 2 ) + " element array",/INFO
		stop
	   endif

;interpolate dust abs & scat coeffs onto wavelength grid of SED:

	Nwave = N_elements( SED_wave )
	Kabs = fltarr( Nwave, ndc )
	Kscat = fltarr( Nwave, ndc )
	gcos = fltarr( Nwave, ndc )
	r = Finterpol( dummy, alog( Dust_Kappa.WaveLen ), $
				alog( SED_wave ), /INIT,/QUIET )
	for i=0,ndc-1 do begin
		Kabs(*,i) = exp( Finterpol( alog( Dust_Kappa.Kabs(*,i) ),/QU))
		Kscat(*,i) = exp( Finterpol( alog( Dust_Kappa.Ksca(*,i) ),/QU))
		gcos(*,i) = Finterpol( Dust_Kappa.gcos(*,i),/QUIET ) < 1
	  endfor

	Kabstot = Kabs # dust_fracs
	frac_dust = fltarr( Nwave, ndc, nphase )
	for i=0,ndc-1 do frac_dust(*,i,1) = Kabs(*,i)*dust_fracs(i)/Kabstot

	Kscatot = Kscat # dust_fracs
	Ktot = Kscatot + Kabstot
	Kalb = Kscatot/Ktot
	gasp = gcos # dust_fracs

;if given, interpolate ICM dust coeffs onto SED wavelength grid:

	if N_struct( Dust_ICM_Kappa ) eq 1 then begin
	  r = Finterpol( dummy, alog( Dust_ICM_Kappa.WaveLen ), $
					alog( SED_wave ), /INIT,/QUIET )
	  for i=0,ndc-1 do begin
	     Kabs(*,i) = exp( Finterpol( alog( Dust_ICM_Kappa.Kabs(*,i) ),/QU))
	     Kscat(*,i) = exp( Finterpol( alog( Dust_ICM_Kappa.Ksca(*,i)),/QU))
	     gcos(*,i) = Finterpol( Dust_ICM_Kappa.gcos(*,i),/QUIET ) < 1
	   endfor
	  if N_elements( dust_ICM_fracs ) eq ndc then begin
		Kabstot = Kabs # dust_ICM_fracs
		for i=0,ndc-1 do $
			frac_dust(*,i,0) = Kabs(*,i)*dust_ICM_fracs(i)/Kabstot
		Kscatot = Kscat # dust_ICM_fracs
		gasp_ICM = gcos # dust_ICM_fracs
	   endif else begin
		Kabstot = Kabs # dust_fracs
		for i=0,ndc-1 do $
			frac_dust(*,i,0) = Kabs(*,i)*dust_fracs(i)/Kabstot
		Kscatot = Kscat # dust_fracs
		gasp_ICM = gcos # dust_fracs
	   endelse
	  Ktot_ICM = Kscatot + Kabstot
	  Kalb_ICM = Kscatot/Ktot
	 endif else frac_dust(*,*,0) = frac_dust(*,*,1)

	if N_elements( srctype ) ne 1 then srctype = "U"
	srctype = strupcase( srctype )
;DUST:
	Mega_Grains, kd_mg, kd_icm, albedo_eff, gp_eff, TAU_CL=tau_cd, $
				XSEC=Ktot, GP=gasp, ALB=Kalb, DEN=dust_dens, $
				ICM_X=Ktot_ICM, ICM_G=gasp_ICM, ICM_A=Kalb_ICM,$
				RADIUS_CL=Rcl, FILL=ffcl, RATIO=drat

	fd_abs_ICM = Mega_Grain_Abs( XSEC=Ktot, GP=gasp, ALB=Kalb, $
				ICM_X=Ktot_ICM, ICM_G=gasp_ICM, ICM_A=Kalb_ICM,$
				RADIUS_CL=Rcl, FILL=ffcl, RATIO=drat, $
				SOURCE=srctype, RADIUS_EFF=Rs, DEN=dust_dens )

	for i=0,ndc-1 do frac_dust(*,i,0) = frac_dust(*,i,0) * fd_abs_ICM
	for i=0,ndc-1 do frac_dust(*,i,1) = frac_dust(*,i,1) * (1 - fd_abs_ICM)
;GAS:
	wang = 1e4 * SED_wave		;convert microns to Angstroms.
	wlyc = where( wang LT 912, nlyc )

	Mega_Grains, kg_mg, kg_icm, ALB=0, RADIUS_CL=Rcl, FILL=ffcl, RAT=drat,$
			XSEC=ismeuv( wang, 1 ), DEN=gas_dens, TAU_CL=tau_cg
;Combine:
	kd_eff = kd_mg + kd_icm
	kg_eff = kg_mg + kg_icm
	frac_gas = kg_eff/(kd_eff + kg_eff)
	frac_gas = [[frac_gas],[frac_gas]]	;for two phases.

;GAS phase fractions:

	if nlyc GT 1 then begin

		fg_abs_ICM = Mega_Grain_Abs( ALB=0, GP=0, RADIUS_CL=Rcl, $
						SOURC=srctype, RATIO=drat, $
					RADIUS_EFF=Rs, FILL=ffcl, DEN=gas_dens,$
					XSEC=ismeuv( wang(wlyc) < 911.5, 1 ) )

		frac_gas(wlyc,0) = frac_gas(wlyc,0) * fg_abs_ICM
		frac_gas(wlyc,1) = frac_gas(wlyc,1) * (1 - fg_abs_ICM)
	  endif

	Tau_Eff = Rs * ( kd_eff + kg_eff )
	Tau_Clu = tau_cd + tau_cg

	CASE srctype OF

	  "C":	esc_prob = exp( -Tau_Eff_Scat( gp_eff, albedo_eff, Tau_Eff ) )

	  "X":	esc_prob = 1 - Pinteract( Tau_Eff, albedo_eff )

	  "U":	esc_prob = Escape_Prob( Tau_Eff, albedo_eff )

	  "MU":	esc_prob = Escape_Prob( Tau_Eff, albedo_eff ) $
	  					* Escape_Prob( Tau_Clu, Kalb )

	  "MC":	esc_prob = Escape_Prob( Tau_Eff, albedo_eff ) $
	  			* exp( -Tau_Eff_Scat( gasp, Kalb, Tau_Clu ) )

	  else:	esc_prob = Escape_Prob( Tau_Eff, albedo_eff )
	 ENDCASE

	if keyword_set( plotit ) then begin
		!P.multi=[0,1,2]
		!X.title="Wavelength (microns)"
		get_window,0,/SHOW
		plot, SED_wave, Tau_Eff,/XTYP,/YTYP, XTIT="", YMARG=[0,2], $
				YTIT="Effective Optical Depth", TIT=srctype
		oplot, SED_wave, Rs * (kd_mg + kg_mg), LINE=2
		oplot, SED_wave, Rs * (kd_icm + kg_icm), LINE=1
		plot, SED_wave, albedo_eff,/XTYP,YR=[0,1], YMARG=[4,2], $
					YTIT="Effective Albedo or g=<cos>"
		oplot, SED_wave, gp_eff, LINE=1
		get_window,1,/SHOW
		plot, SED_wave, frac_gas,/XLOG,YR=[0,1], XTIT="", YMARG=[2,2], $
				YTIT="Absorption Fractions", TIT=srctype,/NODATA
		for j=0,nphase-1 do begin
			oplot, SED_wave, frac_gas(*,j)
			for i=0,ndc-1 do oplot,SED_wave,frac_dust(*,i,j),LIN=2-i
		  endfor
		plot, SED_wave, esc_prob,/XTYP,YR=[0,1], $
				YTIT="Escape Probability & Frac_Abs_ICM"
		oplot, SED_wave, fd_abs_ICM, LINE=2
		empty
		!X.title=""
		!P.multi=0
	   endif

	SED_esc = SED_source * esc_prob
	SED_abs = SED_source - SED_esc

	Lum_tot = Trapez( SED_source, wang )
	Lum_esc = Trapez( SED_esc, wang )

	hc = !cv.h * !cv.c * 1e4	;units = erg-microns
	hfinv = SED_wave/hc		; = 1/(h*freq) = ergs^(-1)

	if nlyc GT 1 then begin
		Lum_Lyc = Trapez( SED_source(wlyc), wang(wlyc) )
		Nphot_Lyc = Trapez( SED_source(wlyc) * hfinv, wang(wlyc) )
	 endif else begin
		Lum_Lyc = 0.0
		Nphot_Lyc = 0.0
	  endelse

	Nphot_Lycabs = fltarr( nphase )	;absorbed by: 0=ICM, 1=clumps.
	Lum_gas = fltarr( nphase )
	Lum_dust = fltarr( ndc, nphase ) ;absorbed by: 0=graphite, 1=silicates.

	for j=0,nphase-1 do begin

		if nlyc GT 1 then begin
			abs = SED_abs(wlyc) * frac_gas(*,j)
			Lum_gas(j) = Trapez( abs, wang(wlyc) )
			Nphot_Lycabs(j) = Trapez( abs * hfinv, wang(wlyc) )
		   endif

		for i=0,ndc-1 do begin
			frac_dust(*,i,j) = frac_dust(*,i,j) * (1-frac_gas(*,j))
			Lum_dust(i,j) = Trapez( frac_dust(*,i,j)*SED_abs, wang )
		  endfor
	  endfor

; Lyman-alpha escape probability:

	if N_elements( Lyafr ) NE 1 then Lyafr = 0.68
	Lum_Lya = Lyafr * Nphot_Lycabs * hc/0.1215
	wLya = scalar( where( SED_wave GE 0.1215 ) )
	gdc  = gas_dens/( ffcl + (1-ffcl)*drat )

	if strpos( srctype, "M" ) eq 0 then begin
		epLya = ep_Lya( TG=gas_T, G=gdc, R=Rs, ALPHA=alpha, $
				TAU=Tau_Clu(wLya), ALB=Kalb(wLya) )
	 endif else begin
		epLya = ep_Lya( TG=gas_T, G=gdc * drat, R=Rs, ALPHA=alpha, $
				TAU=Tau_Eff(wLya), ALB=albedo_eff(wLya) )
	  endelse

	Lum_Lya_abs = Lum_Lya * (1-epLya)
	Lum_Lya_dust = reform( frac_dust(wLya,*,*) ) $
			* transpose( [ [Lum_Lya_abs], [Lum_Lya_abs] ] )
	Lum_dust = Lum_dust + Lum_Lya_dust
	LogLsun = alog10( !cv.Lsun )

	wV = scalar( where( SED_wave GE 0.55 ) )
	wB = scalar( where( SED_wave GE 0.44 ) )
	wU = scalar( where( SED_wave GE 0.36 ) )
	wUBV = [wU,wB,wV]

	Lumins = { 	Lum_tot: Lum_tot, $
			Lum_esc: Lum_esc, $
			Lum_Lyc: Lum_Lyc, $
			Lum_gas: Lum_gas, $
			Log_Nphot_Lyc:    alog10( Nphot_Lyc ) + LogLsun, $
			Log_Nphot_Lycabs: alog10( Nphot_Lycabs ) + LogLsun, $
			Lum_Lya:	Lum_Lya, $
			Lum_Lya_abs:	Lum_Lya_abs, $
			Lum_Lya_dust:	Lum_Lya_dust, $
			Lum_dust:	Lum_dust, $
			Lum_dust_tot:	total( Lum_dust ), $
			SrcType:	srctype, $
			dust_comp:	Dust_Kappa.dust_comp, $
			dust_phase:	["ICM","clumps"], $
			Radius_eff:	Rs, $
			Rclump:		Rcl, $
			FFCL:		ffcl, $
			RATIO_DICM_DCL:	drat, $
			DUST_MASS_DEN:	dust_dens, $
			GAS_DENSITY:	gas_dens, $
			GAS_TEMPERATURE: gas_T,  $
			Tau_UBV:	Rs * dust_dens * Ktot(wUBV), $
			Alb_UBV:	Kalb(wUBV), $
			Tau_Eff_UBV:	Tau_Eff(wUBV), $
			Alb_Eff_UBV:	Albedo_Eff(wUBV), $
			Tau_Eff_Lya:	Tau_Eff(wLya), $
			Alb_Eff_Lya:	albedo_eff(wLya), $
			ep_Lya:		epLya, $
			alpha:		alpha	}

	if keyword_set( verbose ) then begin
		print," "
		print_struct, Lumins, TR=[0,5], FORM=['g','11','3']
		print," "
		print_struct, Lumins, TR=[6,8], FORM=['g','11','3']
		print," "
		print_struct, Lumins, TR=[9,10], FORM=['g','11','3']
	   endif

	if keyword_set( full_struct ) then begin
		Lumins = create_struct( Lumins, "SED_wave", SED_wave, $
						"SED_source", SED_source, $
						"SED_esc", SED_esc, $
						"kd_mg", kd_mg, $
						"kd_icm", kd_icm, $
						"tau_cd", tau_cd, $
						"gp_eff", gp_eff, $
						"albedo_eff", albedo_eff )
		if full_struct gt 1 then $
			Lumins = create_struct( Lumins, "kg_mg", kg_mg, $
							"kg_icm", kg_icm, $
							"tau_cg", tau_cg )
	   endif

if !DEBUG then stop
end