/* This file is part of Cloudy and is copyright (C)1978-2006 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 entries, 
	  exits, 
	  sum, 
	  renorm,
	  collider;

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

	/* check that we were called with valid charge */
	ASSERT( nelem >= 0);
	ASSERT( nelem < LIMELM );
	/* >>chng 05 aug 17, use eden as collider for H and corrected eden for heavier
	 * nuclei - in hot PDR H is collisionally ionized, but not by other H0 since
	 * collision is homonuclear */
	if( nelem==ipHYDROGEN )
	{
		/* special version for H0 onto H0 */
		collider = dense.EdenHontoHCorr;
	}
	else
	{
		collider = dense.EdenHCorr;
	}

	/* do radiative and downward collisional cascades */
	for( level=iso.numLevels_max[ipH_LIKE][nelem]-1; level >= ipH1s; level-- )
	{
		/* captures into this level from the continuum,
		 * radiative rec plus three body rec, units s-1 */
		entries = iso.RateCont2Level[ipH_LIKE][nelem][level];

		/* sum over all levels between current level and top of atom */
		for( iup=level + 1; iup < iso.numLevels_max[ipH_LIKE][nelem]; iup++ )
		{
			/* radiative decays from upper level to here, s-1 */
			double 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);

			/* rate upper level (rel to continuum) decays to here, units s-1 since Pop2Ion is dimless */
			entries += iso.Pop2Ion[ipH_LIKE][nelem][iup]*(RadDecay + 
			  EmisLines[ipH_LIKE][nelem][iup][level].ColUL*collider);
		}

		/* now find rate this levels decays to all lower levels and continuum 
		 * continuum ionization rate, s-1, first */
		exits = iso.gamnc[ipH_LIKE][nelem][level];
		for( low=ipH1s; low <= (level - 1); low++ )
		{
			/* radiative decays to all lower levels - s-1 */
			double 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);

			/* add rad and collisional decays */
			exits += RadDecay + 
				EmisLines[ipH_LIKE][nelem][level][low].ColUL*collider;
		}
		if( exits > 1e-25 )
		{
			iso.Pop2Ion[ipH_LIKE][nelem][level] = entries/exits;
		}
		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.;
		}
	}

	/* =================================================================== 
	 *
	 * 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*collider)/
	  (EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s].Aul + 
		EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s].ColUL*collider + 
		EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].ColUL*collider*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*collider*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 * collider + 
		EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].ColUL * collider);

	/* now renormalize to simple ionization ratio calculated in hydrolevel */
	sum = 0.;
	for( level=0; level < iso.numLevels_max[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_max[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;
		ASSERT( iso.Pop2Ion[ipH_LIKE][nelem][level]<BIGFLOAT );
	}

	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 */

