/*eovrlp returns line destruction probability due to continuum opacity */
#include <string.h>
#include "cddefines.h"
#include "physconst.h"
#include "opac.h"
#include "converge.h"
#include "eovrlpvar.h"
#include "fhummr.h"
#include "eovrlp.h"

double eovrlp(
	  /* abundance of species */
	  double abund, 
	  /* its line absorption cross section */
	  double crsec, 
	  /* pointer to energy within continuum array, to get background opacity */
	  long int ipanu, 
	  /* line width */
	  double widl, 
	  /* escape probability */
	  double escp, 
	  /* type of redistribution function */
	  char *chCore)
{
	/* this will be the value we shall return */
	double eovrlp_v;

	/* DEST0 is the smallest destruction probability to return
	 * in high metallicity models 
	 * this was set to 1e-8 until 99nov18,
	 * in cooling flow model the actual Lya ots dest prob was 1e-16,
	 * and this lower limit of 1e-8 caused energy balance problems,
	 * since dest prob was 8 orders of magnitude too great.  
	 * >>chng 99 nov 18, to 1e-20, but beware that comments indicate that
	 * this will cause problems with high metallicity clouds(??) */
	const double DEST0=1e-20 ;

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

	/* computes "escape probability" due to continuum destruction of
	 *
	 * if esc prob gt 1 then line is masing - return small number for dest prob */
	/* >>>chng 99 apr 10, return min dest when scattering greater than abs */
	if( escp >= 1.0)
	{
		eovrlp_v = DEST0;
		
#		ifdef DEBUG_FUN
		fputs( " <->eovrlp()\n", debug_fp );
#		endif
		return( eovrlp_v );
	}

	/* find continuum opacity */
	eovrlpVar.conopc = (float)opac.opac[ipanu-1];

	/* no idea of opacity whatsoever, on very first soln for this model */
	if( !conv.nTotalIoniz )
	{
		/* opacity not yet defined on very first try, return small but finite number */
		eovrlp_v = DEST0;
		
#		ifdef DEBUG_FUN
		fputs( " <->eovrlp()\n", debug_fp );
#		endif
		return( eovrlp_v );
	}

	/* may be no population, cannot use above test since return 0 not DEST0 */
	if( abund <= 0. || eovrlpVar.conopc <= 0. )
	{
		/* do not set this to DEST0 since energy not then conserved */
		eovrlp_v = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->eovrlp()\n", debug_fp );
#		endif
		return( eovrlp_v );
	}

	/* fac of 1.7 convert to Hummer convention for line opacity */
	eovrlpVar.beta = (float)(eovrlpVar.conopc/(abund*SQRTPI*crsec/widl + eovrlpVar.conopc));
	eovrlpVar.beta = (float)MIN2(eovrlpVar.beta,1.-escp);

	if( strcmp(chCore,"INCO") == 0 )
	{
		/*  fits to 
		 *  >>>refer Hummer and Kunasz 1980 Ap.J. 236,609.
		 *  the max value of 1e-3 is so that we do not go too far
		 *  beyond what Hummer and Kunasz did */
		eovrlp_v = MIN2(1e-3,8.5*eovrlpVar.beta);
	}
	else if( strcmp(chCore,"  K2") == 0 )
	{
		/*  Doppler core only; a=0., Hummer 68 */
		eovrlp_v = fhummr(eovrlpVar.beta);
	}
	else if( strcmp(chCore,"SIMP") == 0 )
	{
		/*  this for debugging only */
		eovrlp_v = 8.5*eovrlpVar.beta;
	}
	else
	{
		fprintf( ioQQQ, " chCore of %4.4s not understood by EOVRLP.\n", 
		  chCore );
		puts( "[Stop in eovrlp]" );
		exit(1);
	}

	eovrlp_v /= 1. + eovrlp_v;
	eovrlp_v *= 1. - escp;

	/*check results in bounds */
	assert( eovrlp_v >= 0.  );
	assert( eovrlp_v <= 1.  );

	{
		/* debugging code for Lya problems */
		enum {DEBUG = FALSE};
		if( DEBUG )
		{
#			include "rfield.h"
#			include "taulines.h"
			if( rfield.anu[ipanu-1]>0.73 && rfield.anu[ipanu-1]<0.76 &&
				(float)(abund)==HydroLines[0][IP2P][IP1S].PopOpc ) 
			{
#				include "zonecnt.h"
				fprintf(ioQQQ,"%li eovrlp\t%g\n",
					ZoneCnt.nzone, eovrlp_v  );
			}
		}
	}

			/* 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 */
	eovrlp_v = (1. - opac.albedo[ipanu-1]) * eovrlp_v + 
		opac.albedo[ipanu-1]*DEST0;

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