;+
; NAME:
;	Pinteract
;
; PURPOSE:
;	Compute and return the probability that a photon entering a sphere
;	of absorbers and/or scatterers (e.g. dust or gas) at a random angle
;	will interact with an absorber/scatterer. In otherwords, it gives
;	the interacting fraction of a collection of randomly impacting photons.
;	If the scattering albedo is given, then the actual absorbed fraction is
;	computed approximately. When scattering is specified to be purely
;	forward then an exact formula used.
;
; CALLING:
;	Pint = Pinteract( tau )
;
; INPUTS:
;	tau = optical radius of sphere (absorption and/or scattering).
;
;	albedo = scattering albedo (optional), if specified then the eventual
;		escape of photons due to scattering is accounted for approx.,
;		in which case Pinteract returns the actual absorbed fraction.
;		If scattering is purely forward then exact formula is used.
; KEYWORDS:
;	/FORWARD : indicates purely forward scattering (when albedo > 0).
; OUTPUTS:
;	Returns probability of interaction, same dimensions as input.
; PROCEDURE:
;	Exact formula is derived by averaging transmission over all
;	impact parameters. Formula was given in Neufeld, ApJL., 1991, and
;	Hobson & Padman, MNRAS, 1993.
;	Exact formula for purely forward scattering is using tau*(1-albedo),
;	approx formula for non-forward is by using escape probability function.
; EXTERNAL CALLS:
;	function Pinteract	(recursively)
;	function escape_prob	(when scattering albedo is specified).
; HISTORY:
;	Written: Frank Varosi, HSTX @ NASA/GSFC, 1997.
;-

function Pinteract, tau, albedo, FORWARD=forw

	if N_params() ge 2 then begin
	    if keyword_set( forw ) then return, Pinteract( tau*(1-albedo) )
	    return, Pinteract( tau ) * (1 - albedo * escape_prob( tau, albedo ))
	 endif

	Pint = 4 * (tau/3.)			;optically thin case.
	w = where( tau GE 5e-5, nw )		;numerical limitation.

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

return, Pint
end