/*fe2ovr compute FeII overlap with Lya */
#include "cddefines.h"
#define	WLAL	1215.6845
#include "nfeii.h"
#include "doppvel.h"
#include "ionfracs.h"
#include "rfield.h"
#include "ipfe2.h"
#include "fe2pmp.h"
#include "fe2tau.h"
#include "phycon.h"
#include "hydrogenic.h"
#include "hwidth.h"
#include "h.h"
#include "zonecnt.h"
#include "pop371.h"
#include "fe2par.h"
#include "fe2ovr.h"

void fe2ovr(void)
{
	long int i;

	static long int nZoneEval;

	double Fe2Partn, 
	  displa, 
	  hopc, 
	  rate, 
	  weight;

	static double BigFeWidth, 
	  BigHWidth, 
	  oldrat;

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

	/* wavelength of Lya in Angstroms */

	/* compute efficiency of FeII emission overlapping with Ly-alpha
	 * implemented with Fred Hamann
	 *
	 * make Ly-a width monotonically increasing to avoid oscillation
	 * in deep regions of x-ray ionized clouds.
	 *
	 * do nothing if Verner FeII atom is used */
	if( FeII.lgFeIION )
	{ 

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

	if( ZoneCnt.nzone <= 1 )
	{
		BigHWidth = hwidth.HLineWidth[0];
		BigFeWidth = DoppVel.doppler[25];
		nZoneEval = ZoneCnt.nzone;
	}

	/* do not do pumping if no population,line is thin, or turned off */
	if( (xIonFracs[2][25] <= 0. || !FeII.lgLyaPumpOn) || 
		hwidth.HLineWidth[0] <= 0. )
	{
		Fe2Partn = 0.;
		oldrat = 0.;
		hydro.dstfe2lya = 0.;

		for( i=0; i < NFEII; i++ )
		{
			fe2tau.Fe2PopLte[i] = 0.;
		}
		
#		ifdef DEBUG_FUN
		fputs( " <->fe2ovr()\n", debug_fp );
#		endif
		return;
	}

	/* only evaluate this one time per zone to avoid oscillations
	 * deep in x-ray ionized clouds */
	if( nZoneEval == ZoneCnt.nzone && ZoneCnt.nzone > 1 )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->fe2ovr()\n", debug_fp );
#		endif
		return;
	}

	BigHWidth = MAX2(BigHWidth,hwidth.HLineWidth[0]);
	BigFeWidth = MAX2(BigFeWidth,DoppVel.doppler[25]);
	nZoneEval = ZoneCnt.nzone;

	/*begin sanity check */
	if( fe2pmp.fe2lam[0] <= 0. )
	{
		fprintf( ioQQQ, " Insane value for fe2lam array.\n" );
		fprintf( ioQQQ, " Was block data fe2dat linked in?\n" );
		puts( "[Stop in fe2ovr]" );
		exit(1);
	}
	/*end sanity check */

	rate = 0.;
	Fe2Partn = fe2par(phycon.te);
	for( i=0; i < NFEII; i++ )
	{
		/* this is displacement from line center in units of Lya width */
		displa = fabs(fe2pmp.fe2lam[i]-WLAL)/WLAL*3e10/BigHWidth;
		if( displa < 1.333 )
		{
			/* have variable weighting factor depending on distance away
			 * this comes form the Verner's fits to Adam's results */
			if( displa <= 0.66666 )
			{
				weight = 1.;
			}
			else
			{
				weight = MAX2(0.,1.-(displa-0.666666)/0.66666);
			}

			fe2tau.Fe2PopLte[i] = 
				(float)(fe2pmp.fe2gs[i]/Fe2Partn*rfield.ContBoltz[ipfe2Com.ipfe2[i]-1]*
				xIonFracs[2][25]);

			fe2tau.feopc[i] = (float)(fe2tau.Fe2PopLte[i]*fe2pmp.fe2osc[i]*
			  0.0150*(fe2pmp.fe2lam[i]*1e-8)/BigFeWidth);

			/* Ly-alpha line-center opacity */
			if( hydro.hn[0][IP1S]*h.hii > 0. )
			{
				hopc = 7.60e-8*hydro.hn[0][IP1S]*h.hii/DoppVel.doppler[0];
			}
			else
			{
				hopc = 7.60e-8*h.hi/DoppVel.doppler[0];
			}
			rate += (fe2tau.feopc[i]/(fe2tau.feopc[i] + hopc))*(BigFeWidth/
			  DoppVel.doppler[0])*(1. - 1./(1. + 1.6*fe2tau.Fe2TauLte[i]))*
			  weight;
		}
	}

	/* dstfe2lya is total Lya deexcitation probability due to line overlap
	 * damper is to stop feedback between here and hydrogen ionization sln */
	hydro.dstfe2lya = (float)((rate + oldrat)/2.);
	oldrat = hydro.dstfe2lya;

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