/*humla fit Hummer and Kunasz escape probability for hydrogen atom Lya */
#include <string.h>
#include "cddefines.h"
#include "physconst.h"
#include "escdst.h"
#include "lamase.h"
#include "opac.h"
#include "humla.h"

void humla(double taume, 
  double beta, 
  float *esc, 
  float *dest,
  /* position of line in frequency array on c scale */
  long ipLine )
{
	double esc0, 
	  fac, 
	  fac1, 
	  fac2, 
	  tau, 
	  taucon, 
	  taulog;

	/* DEST0 is the smallest destruction probability to return
	 * in high metallicity models */
	const double DEST0=1e-8 ;

#	ifdef DEBUG_FUN
	fputs( "<+>humla()\n", debug_fp );
#	endif

	/* fits to numerical results of Hummer and Kunasz Ap.J. 80 */
	tau = taume*SQRTPI;

	if( strcmp(escdst.chEscFun,"SIMP") == 0 )
	{
		/* this set with "escape simple" command, used for debugging */
		esc0 = 1./(1. + taume);
	}
	else
	{
		/* this is the real escape probability */
		esc0 = 1./((0.6451 + tau)*(0.47 + 1.08/(1. + 7.3e-6*tau)));
	}

	esc0 = MAX2(0.,esc0);
	esc0 = MIN2(1.,esc0);

	if( tau > 0. )
	{
		taulog = log10(MIN2(1e8,tau));
	}
	else
	{
		/* the line mased */
		LaMase.xLaMase = (float)MIN2(LaMase.xLaMase,tau);
		taulog = 0.;
		*dest = 0.;
		*esc = (float)esc0;
	}

	if( beta > 0. )
	{
		taucon = MIN2(2.,beta*tau);

		if( taucon > 1e-3 )
		{
			fac1 = -1.25 + 0.475*taulog;
			fac2 = -0.485 + 0.1615*taulog;
			fac = -fac1*taucon + fac2*POW2(taucon);
			fac = pow(10.,fac);
			fac = MIN2(1.,fac);
		}
		else
		{
			fac = 1.;
		}

		*esc = (float)(esc0*fac);
		/* MIN puts cat at 50 */
		*dest = (float)(beta/(0.30972 - MIN2(.28972,0.03541667*taulog)));
	}

	else
	{
		*dest = 0.;
		*esc = (float)esc0;
	}

	*dest = MIN2(*dest,1.f-*esc);
	*dest = MAX2(0.f,*dest);

	/* chng>> 99 apr 12, limit destruction prob in case where gas dominated by scattering.
	 * in this case scattering is much more likely than absorption on this event */
	*dest = (float)( (1. - opac.albedo[ipLine]) * *dest + 
		opac.albedo[ipLine]*DEST0);
	/* this is for debugging H Lya */
	{
		enum {BUG=FALSE};
		if( BUG )
		{
			if( ipLine < 600 )
			{
				fprintf(ioQQQ," %.2e %.2e %.2e %.2e %.2e \n",
				taume,
				beta, 
				(1. - opac.albedo[ipLine]), 
				opac.albedo[ipLine] ,
				*dest 
				);
			}
		}
	}

#	ifdef DEBUG_FUN
	fputs( " <->humla()\n", debug_fp );
#	endif
	return;
}