/*HydroPesc evaluate escape and destruction probabilities for hydrogen lines,
 * called by RTMake */
#include <string.h>
#include "cddefines.h"
#include "physconst.h"
#include "taulines.h"
#include "hydrogenic.h"
#include "strk.h"
#include "phycon.h"
#include "converge.h"
#include "touton.h"
#include "trace.h"
#include "hlife.h"
#include "escdst.h"
#include "twoway.h"
#include "ionfracs.h"
#include "wind.h"
#include "bowoi.h"
#include "rin2ph.h"
#include "p8446.h"
#include "rt.h"
#include "twopht.h"
#include "rtescprob.h"
#include "fe2ovr.h"

void HydroPesc(
	/* ipZ is on C scale, 0 for H, 1 for He, etc */
	long int ipZ ,
	/* flag saying whether we should do escape prob  and reevaluate A and opac (TRUE) or
	 * only dest probs (FALSE).  This is called with TRUE one time per zone,
	 * when optical depths are updated.  It is called FALSE many times, within ionzie,
	 * each time the opacity changes */
	int lgDoEsc )
{
	static int lgTOnSave;
	long int ipHi, 
	  ipLo, 
	  limit;
	double abund=0., 
	  coloi, 
	  es, 
	  factor, 
	  tout,
	  z4;/* physical scale z to the 4th power, used often below */
	float dest=0.f, 
	  esin;
	float FracNew ;
	/* this will be used to save old vals of destruction prob in case
	 * we need to use mean of old and new, as is necessary when ots oscillations
	 * take place */
	float DespSave[NHYDRO_MAX_LEVEL];

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

	/* check that we were called with valid charge */
	assert( ipZ >= 0);
	assert( ipZ < LIMELM );

	/* will need this for some scaling laws - physical Z to the 4th power*/
	z4 = POW2(ipZ+1.);
	z4 *= z4;

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, 
			"          HydroPesc ipZ=%ld called\n", 
		  ipZ );
	}

	/*escape dest function must be one of these three, which will
	 * evaluate to zero, so test product is zero */
	assert(	strcmp( escdst.chEscFun,"  K2") *
				strcmp( escdst.chEscFun,"INCO") *
				strcmp( escdst.chEscFun,"SIMP") == 0 ) ;

	if( lgDoEsc )
	{
		/* hydrogen lines opacity is function of A's
		 * above 15 just use actual A's for high density limit
		 * only change opacities for levels between 4 and 15 */
		limit = MIN2(15,nhlevel);
		/* do Paschen and higher lines first, 
		 * do balmer below since must separate 2s and 2p */
		for( ipHi=4; ipHi <= limit; ipHi++ )
		{
			for( ipLo=3; ipLo < ipHi; ipLo++ )
			{
				HydroLines[ipZ][ipHi][ipLo].Aul = 
					(float)(hlife.HyLife[ipHi]*
					HydroBranch(ipHi,ipLo,ipZ+1)*z4);
				assert(HydroLines[ipZ][ipHi][ipLo].Aul > 0.);

				/* make self-consistent opacity, convert new As back into opacities */
				HydroLines[ipZ][ipHi][ipLo].opacity = 
					(float)(HydroLines[ipZ][ipHi][ipLo].Aul*
				  2.2448e-26*HydroLines[ipZ][ipHi][ipLo].gHi/
				  HydroLines[ipZ][ipHi][ipLo].gLo*
				  POW3(RYDLAM/HydroLines[ipZ][ipHi][ipLo].EnergyRyd));

				/* check that results are ok */
				assert(HydroLines[ipZ][ipHi][ipLo].opacity > 0.);
			}
		}

		/* the actual branching ratio from high levels down to
		 * 2s and 2p depend on the density.  the code goes to
		 * the correct low and high density limits - I know of
		 * nothing better than can be done. */
		factor = MAX2( 0.25 , 0.32 - 0.07*phycon.eden/(phycon.eden + 1e7) );

		/* treat 2s = to 2p for HydroBranch, which returns total to 2s+2p */
		for( ipHi=4; ipHi <= limit; ipHi++ )
		{
			/* get new effective A for this density and temperature */
			HydroLines[ipZ][ipHi][IP2S].Aul = 
				(float)(hlife.HyLife[ipHi]*
				factor*HydroBranch(ipHi,2,ipZ+1)*z4);

				/* check that results are ok */
			assert(HydroLines[ipZ][ipHi][IP2S].Aul > 0.);

			/* do 2p by scaling relative to 2s */
			HydroLines[ipZ][ipHi][IP2P].Aul = (float)(
				HydroLines[ipZ][ipHi][IP2S].Aul / factor *( 1. - factor )); 

				/* check that results are ok */
			assert(HydroLines[ipZ][ipHi][IP2P].Aul > 0.);

			/* make self-consistent opaciity for 2s, from A */
			HydroLines[ipZ][ipHi][IP2S].opacity = 
				(float)(HydroLines[ipZ][ipHi][IP2S].Aul*
			  2.2448e-26*HydroLines[ipZ][ipHi][IP2S].gHi/
			  HydroLines[ipZ][ipHi][IP2S].gLo*
			  POW3(RYDLAM/HydroLines[ipZ][ipHi][IP2S].EnergyRyd));

				/* check that results are ok */
			assert(HydroLines[ipZ][ipHi][IP2S].opacity > 0.);

			/* make self-consistent opaciity for 2p, from A */
			HydroLines[ipZ][ipHi][IP2P].opacity = 
				(float)(HydroLines[ipZ][ipHi][IP2P].Aul*
			  2.2448e-26*HydroLines[ipZ][ipHi][IP2P].gHi/
			  HydroLines[ipZ][ipHi][IP2P].gLo*
			  POW3(RYDLAM/ HydroLines[ipZ][ipHi][IP2P].EnergyRyd));

				/* check that results are ok */
			assert(HydroLines[ipZ][ipHi][IP2P].opacity > 0.);
		}
	}
	/* end test branch lgDoEsc */

	/* now update escape and destruction prob */

	if( xIonFracs[ipZ+2][ipZ] > 1e-30 )
	{
		factor = xIonFracs[ipZ+2][ipZ];
	}
	else
	{
		/* case where almost no parent ion - this will make
		 * very large line opacity, so background dest small */
		factor = 1.;
	}

	/* save destruction probs for Lyman lines in case ots rates oscillate */
	for( ipHi=IP2P; ipHi<= nhlevel; ++ipHi )
	{
		DespSave[ipHi] = HydroLines[ipZ][ipHi][IP1S].Pdest;
	}

	/* first will be static solution */
	if( wind.windv == 0. )
	{
		/* hydrogenic lyman lines special since outward optical depths always set,
		 * trick routines for Lyman lines only */
		lgTOnSave = touton.lgTauOutOn;
		touton.lgTauOutOn = TRUE;

		/* first do Lya, but only if optical depths not overrun */
		tout = HydroLines[ipZ][IP2P][IP1S].TauTot*0.9 - 
			HydroLines[ipZ][IP2P][IP1S].TauIn;

		/* must temporarily make ipLnPopOpc physical */
		HydroLines[ipZ][IP2P][IP1S].PopOpc *= 
			(float)factor;

		/* generate escape prob, pumping rate, destruction prob, 
		 * inward outward fracs  */
		RTMakeStat(&HydroLines[ipZ][IP2P][IP1S] , lgDoEsc );

		/* go back to original units so that final correction ok */
		HydroLines[ipZ][IP2P][IP1S].PopOpc /= 
			(float)factor;

		/* >>>chng 99 dec 18, repair dest prob that got clobbered in call to
		 * RTMakeStat, since will not be evaluated when tout not positive */
		HydroLines[ipZ][IP2P][IP1S].Pdest = DespSave[IP2P];

		if( tout > 0. )
		{
			tout = HydroLines[ipZ][IP2P][IP1S].TauTot - 
			  HydroLines[ipZ][IP2P][IP1S].TauIn;

			abund = HydroLines[ipZ][IP2P][IP1S].PopOpc*xIonFracs[ipZ+2][ipZ];

			/* the descruction prob comes back as dest */
			HydroLines[ipZ][IP2P][IP1S].Pesc = 
				(float)(escla(&esin, &dest,abund,ipZ));

			/* this is current destruction rate */
			HydroLines[ipZ][IP2P][IP1S].Pdest = dest;

			/* chng>>> 99 apr 16, average of old and new rates to damp opacity - Lya ots
			 * rates when opacity is large, as in HeI trip opac in Balmer continuum.
			 * the pphheonly.in test case exposes this problem */
			{
#				include "zonecnt.h"
#				include "opac.h"
				enum {DEBUG=FALSE};

				if( DEBUG && ipZ==0 )
				{
					fprintf(ioQQQ,
						"z%3li Lya eval Pdest popl\t%g\tconopac\t%g\tPdest\t%g\n",
						ZoneCnt.nzone ,
						abund , 
						opac.opac[HydroLines[ipZ][IP2P][IP1S].ipCont-1],
						HydroLines[ipZ][IP2P][IP1S].Pdest);
				}
			}

			/*  this is fraction of line which is inward escaping */
			HydroLines[ipZ][IP2P][IP1S].FracInwd = twoway.fracin;
			/*if(ipZ==1)
				fprintf(ioQQQ," desp=%g\n", dest );*/
		}

		/* now do remainder of Lyman lines, skipping 2s */
		ipLo=IP1S;
		for( ipHi=3; ipHi <= nhlevel; ipHi++ )
		{
			/* must temporarily make ipLnPopOpc physical */
			HydroLines[ipZ][ipHi][ipLo].PopOpc *= 
				(float)factor;

			/* generate escape prob, pumping rate, destruction prob, 
			 * inward outward fracs  */
			RTMakeStat(&HydroLines[ipZ][ipHi][ipLo] , lgDoEsc );

			/* go back to original units so that final correction ok */
			HydroLines[ipZ][ipHi][ipLo].PopOpc /= 
				(float)factor;
		}
		/* reset the flag, so only Lyman lines forced to include outward */
		touton.lgTauOutOn = lgTOnSave;

		/* this loop for Balmer lines which are special,
		 * because must bring 2s and 2p together */
		/* now do 2s and 2p, must bring optical depths together */
		for( ipHi=3; ipHi <= nhlevel; ipHi++ )
		{
			/* these are used for saving current state of two lines */
			float tauin2s , tauin2p , tauout2s , tauout2p , opac2s , opac2p;

			/* 2s inward and total optical depths */
			tauin2s =  HydroLines[ipZ][ipHi][IP2S].TauIn;
			tauout2s = HydroLines[ipZ][ipHi][IP2S].TauTot;

			/* 2p inward and total optical depths */
			tauin2p =  HydroLines[ipZ][ipHi][IP2P].TauIn;
			tauout2p = HydroLines[ipZ][ipHi][IP2P].TauTot;

			/* add inward optical depths together */
			HydroLines[ipZ][ipHi][IP2S].TauIn += 
			HydroLines[ipZ][ipHi][IP2P].TauIn;

			HydroLines[ipZ][ipHi][IP2P].TauIn = 
			HydroLines[ipZ][ipHi][IP2S].TauIn;

			/* add both total optical depths together */
			HydroLines[ipZ][ipHi][IP2S].TauTot += 
			HydroLines[ipZ][ipHi][IP2P].TauTot;

			HydroLines[ipZ][ipHi][IP2P].TauTot = 
			HydroLines[ipZ][ipHi][IP2S].TauTot;

			/* current 2s and 2p opacities */
			opac2s = HydroLines[ipZ][ipHi][IP2S].PopOpc;
			opac2p = HydroLines[ipZ][ipHi][IP2P].PopOpc;

			/* add opacities together */
			HydroLines[ipZ][ipHi][IP2S].PopOpc += 
			HydroLines[ipZ][ipHi][IP2P].PopOpc;

			HydroLines[ipZ][ipHi][IP2P].PopOpc = 
			HydroLines[ipZ][ipHi][IP2S].PopOpc;

			/* must temporarily make ipLnPopOpc physical */
			HydroLines[ipZ][ipHi][IP2S].PopOpc *= 
				(float)factor;
			HydroLines[ipZ][ipHi][IP2P].PopOpc *= 
				(float)factor;

			/* generate escape prob, pumping rate, destruction prob, 
			 * inward outward fracs  */
			RTMakeStat(&HydroLines[ipZ][ipHi][IP2S] , lgDoEsc );
			RTMakeStat(&HydroLines[ipZ][ipHi][IP2P] , lgDoEsc );

			HydroLines[ipZ][ipHi][IP2S].PopOpc =   opac2s;
			HydroLines[ipZ][ipHi][IP2S].TauIn =   tauin2s;
			HydroLines[ipZ][ipHi][IP2S].TauTot = tauout2s;

			/* 2p inward and total optical depths */
			HydroLines[ipZ][ipHi][IP2P].PopOpc =   opac2p;
			HydroLines[ipZ][ipHi][IP2P].TauIn =   tauin2p;
			HydroLines[ipZ][ipHi][IP2P].TauTot = tauout2p;
		}

		/* now do Paschen and higher lines */
		for( ipLo=3; ipLo <= (nhlevel - 1); ipLo++ )
		{
			for( ipHi=ipLo+1; ipHi <= nhlevel; ipHi++ )
			{
				/* must temporarily make ipLnPopOpc physical */
				HydroLines[ipZ][ipHi][ipLo].PopOpc *= 
					(float)factor;

				/* generate escape prob, pumping rate, destruction prob, 
				 * inward outward fracs  */
				RTMakeStat(&HydroLines[ipZ][ipHi][ipLo] , lgDoEsc );

				/* go back to original units so that final correction ok */
				HydroLines[ipZ][ipHi][ipLo].PopOpc /= 
					(float)factor;
			}
		}

		/* 666    C90 had the following code and warning about oscillations
		 *  highest level is a mess, set escape prob to one, actually a blend
		 *  if this is removed, then oscillations in totally optically thick
		 *  blr occur, for cases where all lines are optically thick
		 *  see model badbugs/bug0.in */
		HydroLines[ipZ][nhlevel][IP1S].Pesc = (float)(
			MAX2( HydroLines[ipZ][nhlevel][IP1S].Pesc , 0.1));

		for( ipLo=IP2S; ipLo <= (nhlevel - 1); ipLo++ )
		{
			HydroLines[ipZ][nhlevel][ipLo].Pesc = 1.;
		}
	}
	else
	{
		/* hydrogenic lyman lines special since outward optical depths always set,
		 * trick routines for Lyman lines only */
		lgTOnSave = touton.lgTauOutOn;
		touton.lgTauOutOn = TRUE;

		/* windy model */
		for( ipLo=IP1S; ipLo <= (nhlevel - 1); ipLo++ )
		{
			for( ipHi=MAX2(IP2P,ipLo+1); ipHi <= nhlevel; ipHi++ )
			{
				/* must temporarily make ipLnPopOpc physical */
				HydroLines[ipZ][ipHi][ipLo].PopOpc *= 
					(float)factor;

				RTMakeWind(&HydroLines[ipZ][ipHi][ipLo] , lgDoEsc );

				/* go back to original units so that final correction ok */
				HydroLines[ipZ][ipHi][ipLo].PopOpc /= 
					(float)factor;
			}
			/* reset the flag, so only Lyman lines forced to include outward */
			touton.lgTauOutOn = lgTOnSave;
		}
	}

	/* this is option to damp out Lyman line dest probs if ots rates oscillating */
	if( conv.lgOscilOTS )
	{
		/* this is damper used to stop oscillations when now present*/
		FracNew = 0.2f;
	}
	else
	{
		/* this is damper used to stop oscillations even when not detected */
		FracNew = 0.5f;
	}
	{
		/* following should be set true to print ots contributors */
		enum {DEBUG = FALSE};
		if( DEBUG && ipZ==0 )
		{
			fprintf(ioQQQ,"Lya aver DespSave\t%g\tPdest\t%g\n",
				DespSave[IP2P],
				HydroLines[ipZ][IP2P][IP1S].Pdest);
		}
	}
	for( ipHi=IP2P; ipHi<= nhlevel; ++ipHi )
	{
		/*lint -e771 DespSave possibly not initialized */
		HydroLines[ipZ][ipHi][IP1S].Pdest = (1.f-FracNew)*DespSave[ipHi] +
			FracNew * HydroLines[ipZ][ipHi][IP1S].Pdest;
		/*lint +e771 DespSave possibly not initialized */
	}

	{
		/* following should be set true to print ots contributors */
		enum {DEBUG = FALSE};
		if( DEBUG && ipZ==0 )
		{
			fprintf(ioQQQ,"Lya ave desp\t%g\t popl\t%g\tconopac\t%g\tbeta\t%g\n",
				dest , 
				abund , 
				DespSave[IP2P],
				HydroLines[ipZ][IP2P][IP1S].Pdest);
		}
	}

	/* reset the flag */
	touton.lgTauOutOn = lgTOnSave;

	if( ipZ == 0 )
	{
		/* find Stark escape probabilities */
		strk();
		/* electron scattering escape */
		es = phycon.eden*6.65e-25;
		for( ipLo=IP1S; ipLo <= (nhlevel - 1); ipLo++ )
		{
			for( ipHi=MAX2(IP2P,ipLo+1); ipHi <= nhlevel; ipHi++ )
			{
				hydro.esesc[ipZ][ipLo][ipHi] = (float)(MIN2(1.,5.9*es/(es+
				  HydroLines[ipZ][ipHi][ipLo].PopOpc*
				  HydroLines[ipZ][ipHi][ipLo].opacity))*
				  xIonFracs[ipZ+2][ipZ]*
				  MAX2(0.,1.-HydroLines[ipZ][ipHi][ipLo].Pesc-
				  HydroLines[ipZ][ipHi][ipLo].Pdest));
				hydro.esesc[ipZ][ipLo][ipHi] = 0.;
			}
		}
		/*fixit(); in above, electrons cattering escape set to 0 since got
		           * very large in case b test */

		/* do the 8446 problem */
		p8446(&coloi);
		HydroLines[0][3][IP1S].Pesc = bowoi.pmph31/
		  HydroLines[0][3][IP1S].Aul;

		if( trace.lgTrace && trace.lgHBugFull )
		{
			fprintf( ioQQQ, "       HydroPesc calls P8446 who found pmph31=%10.2e\n", 
			  bowoi.pmph31 );
		}

		/* add on destruction of hydrogen Lya by FeII
		 * now add in FeII deexcitation via overlap,
		 * but only as loss, not photoionization, source
		 * dstfe2lya is Ly-alpha deexcit by FeII overlap - conversion into FeII em */
		/* find FeII overlap destruction rate, added to net rate in HydroSumTrans */
		fe2ovr();
#if 0
		/* this was introduced in the fort - c conversion, and is a bug since
		 * dest feii added to ots dest */
		HydroLines[0][IP2P][IP1S].Pdest = (float)(MIN2(1.,
		  HydroLines[0][IP2P][IP1S].Pdest+hydro.dstfe2lya));
#endif

		/* find the rate of induced two photon */
		twopht();
		HydroLines[0][IP2S][IP1S].pump = rin2ph.ri1s2s;
	}
	else
	{
		/* not included for non-H species */
		HydroLines[ipZ][IP2S][IP1S].pump = 0.;
	}

	/*fixit(); this was cut out of pesc and pasted here.  this needs to
				 * be added back into the total rates, even though it is not important */
#	if 0
	/* increase escape probability by scattering and stark */
	HydroLines[ipZ][ipHi][ipLo].Pesc = (float)(fmin(1.,
	  HydroLines[ipZ][ipHi][ipLo].Pesc+
	  hydro.pestrk[ipLo][ipHi]+hydro.esesc[ipZ][ipLo][ipHi]));
#	endif

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