;+
; NAME:
;	Mega_Grains
; PURPOSE:
;	Approximate the effective absorption and scattering properties
;	of a two-phase clumpy medium of dust by considering the randomly
;	located spherical clumps as large grains:  "mega-grains".
;	Assumes that clumps are uniformly distributed, all of same radius
;	and same density relative to the inter-clump medium (ICM).
;	Based on the theory of Hobson & Padman 1993, MNRAS, 264, 161-164,
;	and extended to the case when the filling factor of the clumps > 0.3.
;	Effective albedo of the clumps (albedo_CL) is estimated using the
;	Osterbrock-Lucy escape probability instead of Hobson & Padman eq.,
;	thereby improving agreement with Monte Carlo simulations.
;	Includes new approximation for effective scattering asymmetry parameter
;	(see Varosi and Dwek, 1999 ApJ 523, 265 for all equations).
;	Optionally, the optical properties of dust in the ICM can be specified
;	with ICM_* keywords, if different than dust in clumps.
; CALLING:
;	Mega_Grains, kabscat_MG, kabscat_ICM, albedo_EFF, g_EFF
;
; KEYWORD INPUTS:
;
;	XSEC_DUST = total (absorption + scattering) cross-section of dust.
;	ALBEDO_DUST = dust scattering albedo (default = 0).
;	GP_DUST = scattering asymmetry parameter of dust.
;
;	ICM_XSEC = optional, Xsec for dust in the ICM if different than clumps.
;	ICM_ALBEDO = scattering albedo for dust in the ICM (optional).
;	ICM_GP = scatt. asymmetry param. for dust in the ICM (optional).
;
;	DENS_HOMOG = desired equivalent homogenous (average) density (default=1)
;	FILL_FACT = volume filling factor of clumps.
;	RADIUS_CLUMP = radius of a spherical clump (same units as XSEC & DENS).
;	RATIO_ICM_CL = ratio of ICM to clump densities.
;	DRATIO_CL_ICM = ratio of clump to ICM densities (overrides RATIO_ICM_CL)
;
;	GAMMA_F = exponent used in extended MG approx: (1-ff)^gamma_f
;		which is the clump radii renormalization factor.
;		Default = 1, so default renorm factor is just (1-ff).
;
;	Note: any one input keyword can be a vector (so outputs are vectors),
;		then other inputs should be scalar or vector of same size.
;
;	/NORMAL_CLUMPS : flag to use full clump optical depth instead of
;		the default due to just the overdensity.
;		This option is invoked by function mega_grain_abs.
;
;	/HP_MG : set this to use original Hobson & Padman mega-grains approx.
;	/HP_ALBEDO : set to use original Hobson & Padman approx. clump albedo,
;		default is to use Osterbrock-Lucy escape probability from clump.
; OUTPUTS:
;	kabscat_MG = effective absorption + scattering coefficient per unit
;			length in medium of the clumps only.
;
;	kabscat_ICM = absorption + scattering coefficient of inter-clump medium,
;			again per unit length.
;
;	albedo_EFF = effective scattering albedo of the medium: clumps and ICM.
;
;	g_EFF = effective scattering asymmetry parameter of the clumpy medium.
;
; KEYWORD OUTPUTS (optional):
;	ACLUMP = equivalent scattering albedo of each clump.
;	GCLUMP = equivalent scattering asymmetry parameter of each clump.
;	TAU_CL = optical depth of each clump.
;
; EXAMPLE OF USING OUTPUT:
;	Effective optical depth of sphere:
;		tau_eff = radius * ( kabscat_MG + kabscat_ICM )
;	With scattering:
;		tau_effsc = tau_eff_scat( g_EFF, albedo_EFF, tau_eff )
;	Escape probabilities:
;		central source:		ep = exp( -tau_effsc )
;		uniform source:		ep = escape_prob( tau_eff, albedo_EFF )
; EXTERNAL CALLS:
;	function Pinteract
;	function escape_prob
; PROCEDURE:
;	Based on theory of Hobson & Padman 1993, MNRAS, 264, 161-164,
;	"Radiative transfer in a clumpy medium - II:
;	 The mega-grains approximation for two-phase models".
;	Extended to filling factors > 0.3 by decreasing the effective
;	clump radii by (1-f) yielding more and smaller clumps.
;	Extended approximation to ICM/Clump density ratios > 1 (cavities) by
;	simply reversing the roles of clumps and ICM.
;	Improved approximation for effective scattering albedo.
;	Introduced new approximation for effective scattering asymmetry param.
; HISTORY:
;	Written: Frank Varosi, HSTX @ NASA/GSFC, 1996.
;	FV 1998, added code for g_EFF and keywords GAMMA_F= , /NORMAL_CLUMPS.
;	FV 1999, added keywords ICM_XSEC, ICM_ALBEDO, ICM_GP and related code.
;-

pro Mega_Grains, kabscat_MG, kabscat_ICM, albedo_EFF, g_EFF, $
					XSEC_DUST=xsec, ICM_XSEC=xsicm, $
				ALBEDO_DUST=albedo, ICM_ALBEDO=albedo_ICM, $
					GP_DUST=gp, ICM_GP=gp_ICM, $
			RADIUS_CLUMP=radius_CL, FILL_FACT=fill_fact, $
			RATIO_ICM_CL=den_rat, DRATIO_CL_ICM=dratio_CL_ICM, $
			TAU_CLUMP=tau_CL, ACLUMP=albedo_CL, GCLUMP=gclump,$
				NORMAL_CLUMPS=normc, DENS_HOMOG=den_hom, $
				HP_ALBEDO=hp_albedo, HP_MG=hp_mg, GAMMA_F=fgam
   common MG_test, fpow

	if N_elements( albedo ) LE 0 then albedo = 0
	if N_elements( den_hom ) LE 0 then den_hom = 1

	if N_elements( dratio_CL_ICM ) gt 0 then den_rat = 1.0/dratio_CL_ICM
	ndr = N_elements( den_rat )

	if keyword_set( normc ) then begin

		ffc = fill_fact
		drat = den_rat
		if (ndr GT 1) and (ndr NE N_elements( ffc )) then $
					ffc = replicate( ffc(0), ndr )
	 endif else begin

		if( ndr GT 1 ) then begin     ;case of array of density ratios:
			drat = den_rat
			ffc = fill_fact
			if ndr NE N_elements( ffc ) then $
					ffc = replicate( ffc(0), ndr )
			w = where( drat GT 1, nw )
			if (nw GT 0) then begin		;dens.clumps < dens.ICM
				drat(w) = 1.0/drat(w)
				ffc(w) = 1 - ffc(w)
			   endif
		 endif else if den_rat GT 1 then begin	;dens.clumps < dens.ICM
			ffc = 1 - fill_fact
			drat = 1.0/den_rat
		  endif else begin	 ;standard case: dens.clumps > dens.ICM
			ffc = fill_fact
			drat = den_rat
		   endelse
	  endelse

	if keyword_set( hp_mg ) then begin
		renorm = 1
		xmg = -3 * aLog( 1 - ffc )/( 4 * radius_CL )
	 endif else begin			;extension to larger fill_fact.
		if N_elements( fpow ) eq 1 then renorm = (1-ffc)^fpow $
					   else	renorm = 1 - ffc
		if N_elements( fgam ) eq 1 then renorm = (1-ffc)^fgam
		xmg = ffc * 3/( 4 * radius_CL * renorm )
	  endelse

	den_CL = den_hom/( ffc + drat*(1-ffc) )
	if keyword_set( normc ) then tau_CL = den_CL * xsec * radius_CL $
				else tau_CL = den_CL * (1-drat) * xsec*radius_CL

	kabscat_MG = xmg * Pinteract( tau_CL * renorm )

	if N_elements( xsicm ) eq N_elements( xsec ) then $
			kabscat_ICM = drat * den_CL * xsicm $
		  else 	kabscat_ICM = drat * den_CL * xsec

	if keyword_set( hp_albedo ) then $
		albedo_CL = albedo/( 1 + (1-albedo)*4*tau_CL/3 ) $	;H. & P.
	else $
	   albedo_CL = albedo * escape_prob( tau_CL * renorm, albedo )	;Varosi.

	if N_elements( albedo_ICM ) ne N_elements( albedo ) then $
					albedo_ICM = albedo

	albedo_EFF = ( albedo_CL * kabscat_MG + albedo_ICM * kabscat_ICM ) / $
				( ( kabscat_MG + kabscat_ICM ) > 1e-37 )

	if (N_elements( gp ) gt 0) and (max( albedo ) gt 0) then begin
		gg = gp > 0
		gfit = [ 1.5, 2., 3. ]
		a = gfit(0) + 4*gg^3 + 2*sqrt( gg ) * exp( -5*gg ) * albedo
		b = gfit(1) - gg*(1-gg) - 2*gg*albedo
		c = 1/( gfit(2) - sqrt(2*gg) - 2*gg*(1-gg) * albedo )
		gclump = gp - c*( 1 - ( 1+exp(-b/a) )/( 1+exp((tau_cl-b)/a) ) )
		if N_elements( gp_ICM ) ne N_elements( gp ) then gp_ICM = gp
		g_EFF = ( gclump * kabscat_MG + gp_ICM * kabscat_ICM ) / $
				( ( kabscat_MG + kabscat_ICM ) > 1e-37 )
	 endif else begin
		gclump = 0.0
		g_EFF = 0.0
	  endelse
end