;+
; NAME:
;	Escape_Prob
;
; PURPOSE:
;	Compute the probability of escape for photons emitted uniformly
;	in a homogenous sphere of absorbers, possibly with scattering.
;	For no scattering (albedo = 0) the equation is exact (Osterbrock 1989).
;	With scattering the accuracy of generalized formula (Lucy 1989)
;	has been tested against Monte-Carlo simulations (Varosi 1995)
;	and was found to depend also on the angular distribution of the
;	scattering. For small optical depths ( < 1 ) the formual agrees with
;	Monte-Carlo to better than 5% for scattering ranging from isotropic
;	to slightly forward. As optical depth increases the formula
;	over-estimates the escape probability for isotropic scattering, but
;	agrees better with cases of progressively more forward scattering.
;	Formula can also be applied to a homogenous disk/cylinder/ellipsoid
;	by using effective tau = 3*tauZ * tauR / ( 2*tauZ + tauR ) where
;	tauZ and tauR are the vertical and radial optical depths to the center
;	of disk/cylinder/ellipsoid (Varosi 1995). Results are within 5% of
;	Monte-Carlo simulations for cases of moderate forward scattering.
;
; CALLING:
;	Pescape = Escape_Prob( tau, albedo )
;
; INPUTS:
;	tau	= array of optical depths (radii, from center of sphere).
;	albedo	= array of scattering albedos (ratios of scattering to total
;		cross-sections), default = 0.
;
; KEYWORDS:
;	/MATRIX : compute matrix of probabilities for all combinations
;		of tau and albedo (also the case if # taus NE # albedos).
;
; OUTPUTS:
;	Function returns vector of photon escape probabilities corresponding
;	to the arrays of optical radii and albedos if equal in number.
;	If tau and albedo arrays are different sizes or /MATRIX is set
;	then a matrix of probabilities is returned, one element (i,j) for each
;	combination of tau(i) and albedo(j).
;
; PROCEDURE:
;	The zero albedo case (exact derivation for no scattering) is from
;	Osterbrock, "Astrophysics of Gaseous Nebulae...", 1989, Appendix 2.
;	The non-zero albedo case is from Lucy et al., in "Supernovae", 1989,
;	edited by Woosley. Formula is stated with no derivation and is
;	independent of the Osterbrock formula. However, with the assumption
;	that scattering redistributes the photons uniformly,
;	any zero albedo escape probability can be applied recursively,
;	resulting in a geometric series which has a sum yielding the formula.
;
; HISTORY:
;	Written: Frank Varosi, HSTX @ NASA/GSFC, 1995.
;-


function Escape_Prob, tau, albedo, MATRIX=matrix

	ntau = N_elements( tau )
	if (ntau LE 0) then print,"SYNTAX:  pe = escape_prob( tau, albedo )"
	p0 = replicate( 1.0, ntau )
	w = where( tau GE 1e-4, nw )

	if (nw GT 0) then begin
		tauw = double( tau(w) )
		T1 = 1./tauw
		Texp = exp( -2*tauw )
		p0(w) = 0.75 * T1 * ( 1 + T1*Texp + (Texp - 1)/(2.*tauw^2) )
	   endif

	if (ntau EQ 1) then  p0 = p0(0)
	nalb = N_elements( albedo )
	if (nalb LE 0) then return,p0

	if (nalb NE ntau) OR keyword_set( matrix ) then begin

		pe = replicate( 1.0, ntau, nalb )
		for i=0,nalb-1 do pe(*,i) = p0/( 1-albedo(i)*(1-p0) )

	 endif else pe = p0/( 1 - albedo * (1-p0) )

if N_elements( pe ) eq 1 then  return, pe(0)  else  return, pe
end