;+
; NAME:
;	EP_Rad_Trans
; PURPOSE:
;	Approximate the radiative transfer of photons emitted in or impacting
;	a homogenous medium (dust & gas) of spherical/disk geometry,
;	using 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.
;	An approximation of absorption, ionization of gas, and resultant
;	heating of dust by Lyman-alpha from gas is also attempted,
;	but this needs more work.
; CALLING:
;	EP_Rad_Trans, SED_wave, SED_emit, 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_emit = emitted Spectral Energy Distribution (SED) vs. wavelength,
;		in solar luminosities per second per Angstrom.
;
; OUTPUTS:
;	SED_esc = escaping spectral energy distribution vs. wavelength.
;
;	Lumins = structure variable accounting of all Luminosities emitted,
;		escaped, and absorbed in dust/gas.
;
; KEYWORD INPUTS:
;
;	KAPPA_DUST = structure variable with dust absorption/scattering coeffs.
;	DUST_MASS_DEN = array of dust mass densities for each component.
;	GAS_DEN = density of atomic hydrogen, same units as RADIUS.
;	RADIUS = effective radius of sphere/disk containing clumpy medium.
;
;	SOURCE = string (upper or lowercase) indicating source of photons:
;		"C" : do case of central source,
;		"X" : do case of uniformly illuminating external source,
;			default is "U": uniformly distributed internal source.
;
;	/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.
;	FD = of the photons absorbed, the fractions going to each dust comp.
;	FG = fraction of photons absorbed by the gas at each wavelength.
; EXTERNAL CALLS:
;	function Finterpol
;	function ismeuv
;	function Tau_Eff_Scat
;	function Escape_Prob
;	function Pinteract
;	function Trapez
;	function ep_Lya
; PROCEDURE:
; HISTORY:
;	Written: Frank Varosi, HSTX @ NASA/GSFC, 1995
;-

pro EP_Rad_Trans, SED_wave, SED_emit, SED_esc, Lumins, PLOT=plotit, $
			KAPPA_DUST=Dust_Kappa, DUST_MASS_DEN=den_dustm, $
			RADIUS=Rs, GAS_DENSITY=den_gas, GAS_TEMPERATURE=Tgas, $
			EPWAVE=esc_prob, EPLYA=epLya, TAU_TOT=Tau_tot, $
			ALBEDO=albedo, FD=frac_dust, FG=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.

	if N_elements( den_dustm ) NE ndc then begin
		message,"keyword DUST_MASS_DEN 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 )
	junk = 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) )))
		Kscat(*,i) = exp( Finterpol( alog( Dust_Kappa.Ksca(*,i) )))
		gcos(*,i) = Finterpol( Dust_Kappa.gcos(*,i) ) < 1
	  endfor

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

	Kscatot = Kscat # den_dustm
	Ktot = Kscatot + Kabstot
	gcos = gcos # ( den_dustm/total( den_dustm ) )
	albedo = Kscatot/Ktot

	wang = 1e4 * SED_wave		;convert microns to Angstroms.
	wlyc = where( wang LT 912, nlyc )

	Tau_tot = Rs * ( Ktot + ismeuv( wang, den_gas ) )
	frac_gas = 1 - Rs*Ktot/Tau_tot

	if N_elements( srctype ) ne 1 then srctype = "U"
	srctype = strupcase( srctype )

	CASE srctype OF

	  "C":	esc_prob = exp( -Tau_Eff_Scat( gcos, albedo, Tau_tot ) )

	  "X":	esc_prob = 1 - Pinteract( Tau_tot, albedo )

	  else:	esc_prob = Escape_Prob( Tau_tot, albedo )
	 ENDCASE

	if keyword_set( plotit ) then begin
		!P.multi=[0,1,2]
		!X.title="Wavelength (microns)"
		get_window,0,/SHOW
		plot, SED_wave, Tau_tot,/XTYP,/YTYP, XTIT="", YMARG=[0,2], $
					YTIT="Optical Depth"
		plot, SED_wave, albedo,/XTYP,YR=[0,1], YMARG=[4,2], $
					YTIT="Albedo"
		if srctype eq "C" then oplot, SED_wave, gcos, LINE=1
		get_window,1,/SHOW
		plot, SED_wave, frac_gas,/XTYP,YR=[0,1], XTIT="", YMARG=[2,2], $
			YTIT="Absorption Fractions (Gas, graphite, silicates)"
		for i=0,1 do oplot, SED_wave, frac_dust(*,i),LINE=2-i
		plot, SED_wave, esc_prob,/XTYP,YR=[0,1], $
				YTIT="Escape Probability"
		empty
		!X.title=""
		!P.multi=0
	   endif

	SED_esc = SED_emit * esc_prob
	SED_abs = SED_emit - SED_esc

	Lum_tot = Trapez( SED_emit, 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_emit(wlyc), wang(wlyc) )
		Nphot_Lyc = Trapez( SED_emit(wlyc) * hfinv, wang(wlyc) )
	 endif else begin
		Lum_Lyc = 0.0
		Nphot_Lyc = 0.0
	  endelse

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

	Lum_dust = fltarr( ndc )	;absorbed by: 0=graphite, 1=silicates.
	abs = SED_abs * (1-frac_gas)
	for i=0,ndc-1 do Lum_dust(i) = Trapez( frac_dust(*,i)*abs, wang )

;obtain optical depth and albedo of dust at Lyman-alpha wavelength,
; and get 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 ) )

	epLya = ep_Lya( TG=Tgas, G=den_gas, R=Rs, ALPHA=alpha, $
			TAU=Tau_tot(wLya), ALB=albedo(wLya) )

	Lum_Lya_abs = Lum_Lya * (1-epLya)
	Lum_Lya_dust = reform( frac_dust(wLya,*) ) * 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, $
			Radius_eff:	Rs, $
			DUST_MASS_DEN:	den_dustm, $
			GAS_DENSITY:	den_gas, $
			GAS_TEMPERATURE: Tgas,  $
			Tau_UBV:	Rs * Ktot(wUBV), $
			Alb_UBV:	albedo(wUBV), $
			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_emit", SED_emit, $
						"SED_esc", SED_esc, $
						"Tau_tot", Tau_tot, $
						"gcos", gcos, $
						"albedo", albedo )
	   endif

if !DEBUG then stop
end