/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*HydroT2Low called to do hydrogenic level populations when temp too low for matrix,
 * forces total of level populations to add up to iso.xIonSimple[ipH_LIKE][nelem],
 * so results of this routine will always agree with that value */
#include "cddefines.h"
#include "taulines.h"
#include "iso.h"
#include "secondaries.h"
#include "ionbal.h"
#include "dense.h"
#include "trace.h"
#include "phycon.h"
#include "hydrogenic.h"
/*lint -e662 out of bounds pointer */

void HydroT2Low( long int nelem )
{
	long int iup, 
	  level, 
	  low;
	double above, 
	  below, 
	  sum, renorm;

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

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

	/* do radiative only cascades */
	for( level=iso.numLevels[ipH_LIKE][nelem]-1; level >= ipH2s; level-- )
	{
		/* radiative rec plus three body rec */
		/* >>chng 04 apr 12, had used these explicit terms, change over to
		 * vectors defined in iso_ionize_recombine */
		/*above = (iso.RadRecomb[ipH_LIKE][nelem][level][ipRecRad]*
			iso.RadRecomb[ipH_LIKE][nelem][level][ipRecNetEsc] + 
		   iso.ColIoniz[ipH_LIKE][nelem][level]*iso.PopLTE[ipH_LIKE][nelem][level]*dense.EdenHCorr)*
	  	   dense.eden;*/
		above = iso.RateCont2Level[ipH_LIKE][nelem][level];

		for( iup=level + 1; iup < iso.numLevels[ipH_LIKE][nelem]; iup++ )
		{
			double RadDecay;
			RadDecay =
				EmisLines[ipH_LIKE][nelem][iup][level].Aul*
			  (EmisLines[ipH_LIKE][nelem][iup][level].Pesc + 
			  EmisLines[ipH_LIKE][nelem][iup][level].Pelec_esc + 
			  EmisLines[ipH_LIKE][nelem][iup][level].Pdest);

			/* hlev(level,low) is A * esc prob */
			above += iso.Pop2Ion[ipH_LIKE][nelem][iup]*(RadDecay + 
			  EmisLines[ipH_LIKE][nelem][iup][level].ColUL*dense.EdenHCorr);
		}

		below = 0.;
		for( low=ipH1s; low <= (level - 1); low++ )
		{
			double RadDecay;
			RadDecay =
				EmisLines[ipH_LIKE][nelem][level][low].Aul*
			  (EmisLines[ipH_LIKE][nelem][level][low].Pesc + 
			  EmisLines[ipH_LIKE][nelem][level][low].Pelec_esc + 
			  EmisLines[ipH_LIKE][nelem][level][low].Pdest);

			below += RadDecay + EmisLines[ipH_LIKE][nelem][level][low].ColUL*
			  dense.EdenHCorr;
		}
		if( below > 1e-25 )
		{
			iso.Pop2Ion[ipH_LIKE][nelem][level] = above/below;
		}
		else
		{
			iso.Pop2Ion[ipH_LIKE][nelem][level] = 0.;
		}

		if( iso.PopLTE[ipH_LIKE][nelem][level] > 1e-25 )
		{
			iso.DepartCoef[ipH_LIKE][nelem][level] = 
				iso.Pop2Ion[ipH_LIKE][nelem][level]/(iso.PopLTE[ipH_LIKE][nelem][level]*
			  dense.eden);
		}
		else
		{
			iso.DepartCoef[ipH_LIKE][nelem][level] = 0.;
		}
	}

	/* following is total rec coef */
	/*ionbal.RateRecomTot[nelem][nelem] = dense.eden * (
		iso.RadRecomb[ipH_LIKE][nelem][ipH1s][ipRecRad]*
	  iso.RadRecomb[ipH_LIKE][nelem][ipH1s][ipRecNetEsc] + iso.RadRec_caseB[ipH_LIKE][nelem] + 
	  ionbal.CotaRate[nelem]);*/

	/* =================================================================== 
	 *
	 * at this point all destruction processes have been established 
	 *
	 * ==================================================================== */

	/* do simple sums for atom to ion, then fill in for rest */
	iso.Pop2Ion[ipH_LIKE][nelem][ipH1s] = ionbal.RateRecomTot[nelem][nelem]/iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];

	if( iso.PopLTE[ipH_LIKE][nelem][ipH1s] > 1e-25 )
	{
		iso.DepartCoef[ipH_LIKE][nelem][ipH1s] = 
			iso.Pop2Ion[ipH_LIKE][nelem][ipH1s]/(iso.PopLTE[ipH_LIKE][nelem][ipH1s]*
			dense.eden);
	}
	else
	{
		iso.DepartCoef[ipH_LIKE][nelem][ipH1s] = 0.;
	}

	iso.Pop2Ion[ipH_LIKE][nelem][ipH2s] = 
		(iso.RadRec_caseB[ipH_LIKE][nelem]*0.33*dense.eden + 
		iso.Pop2Ion[ipH_LIKE][nelem][ipH1s]*secondaries.Hx12[MIN2(nelem+1,2)-1][ipH2s] + 
		iso.Pop2Ion[ipH_LIKE][nelem][ipH2p]*
		EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].ColUL*dense.EdenHCorr)/
	  (EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s].Aul + 
		EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s].ColUL*dense.EdenHCorr + 
		EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].ColUL*dense.EdenHCorr*6. );

	iso.Pop2Ion[ipH_LIKE][nelem][ipH2p] = 
		(iso.RadRec_caseB[ipH_LIKE][nelem]*0.67*dense.eden + 
		iso.Pop2Ion[ipH_LIKE][nelem][ipH1s]*secondaries.Hx12[MIN2(nelem+1,2)-1][ipH2p] + 
		iso.Pop2Ion[ipH_LIKE][nelem][ipH2s]*
		EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].ColUL*dense.EdenHCorr*6.)/
		(EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Aul*
		(EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pesc + 
		EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pelec_esc +
		EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pdest) +
		EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].ColUL * dense.EdenHCorr + 
		EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].ColUL * dense.EdenHCorr);

	/* now renormalize to simple ionization ratio calculated in hydrolevel */
	sum = 0.;
	for( level=0; level < iso.numLevels[ipH_LIKE][nelem]; level++ )
	{
		sum += iso.Pop2Ion[ipH_LIKE][nelem][level];
	}
	/* rescale so that inverse of sum is iso.xIonSimple[ipH_LIKE][nelem] */
	renorm = 1. / SDIV( iso.xIonSimple[ipH_LIKE][nelem]) / SDIV( sum);
	for( level=0; level < iso.numLevels[ipH_LIKE][nelem]; level++ )
	{
		/* >>chng 04 apr 11, last index below was ipH2s rather than level - 
		 * a major bug */
		iso.Pop2Ion[ipH_LIKE][nelem][level] *= renorm;
	}

	if( trace.lgHBug && trace.lgTrace )
	{
		fprintf( ioQQQ, "       LOW TE,=%10.3e HN(1)=%10.3e rec=%10.3e Hgamnc(1s)=%10.3e\n", 
		  phycon.te, iso.Pop2Ion[ipH_LIKE][nelem][ipH1s], ionbal.RateRecomTot[nelem][nelem], 
		  iso.gamnc[ipH_LIKE][nelem][ipH1s] );
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, 
			"       HydroT2Low return, LOW TE used, HII/HI=%.3e simple=%.3e IonRate=%.3e RecCo=%.3e\n", 
		  iso.pop_ion_ov_neut[ipH_LIKE][nelem], 
		  iso.xIonSimple[ipH_LIKE][nelem], 
		  iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s], 
		  ionbal.RateRecomTot[nelem][nelem] );
	}

#	ifdef DEBUG_FUN
	fputs( " <->HydroT2Low()\n", debug_fp );
#	endif
	return;
}
/*lint +e662 out of bounds pointer */

