/*HydroT2Low called to do hydrogenic level populations when temp too low for matrix */
#include "cddefines.h"
#include "taulines.h"
#include "hydrogenic.h"
#include "hcolrt.h"
#include "secondaries.h"
#include "elecden.h"
#include "hphoto.h"
#include "charexc.h"
#include "trace.h"
#include "phycon.h"
#include "thlo.h"
#include "hcaseb.h"
#include "hlteok.h"
#include "rcota.h"
#include "hionfraccom.h"

void HydroT2Low(long int ipZ, 
  double rfac)
{
	long int iup, 
	  lev, 
	  low;
	double above, 
	  below, 
	  topoff;

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

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

	/* changed to EdenHCorr
	 * edhi = eden + hi * 1.7e-4
	 * */
	if( trace.lgHBug && trace.lgTrace )
	{
		if( phycon.te <= thlo.HydTempLimit*(ipZ+1) )
		{
			fprintf( ioQQQ, " H%3ld matrix not used, te too low\n", ipZ );
		}
		if( rfac < 1e-7 )
		{
			fprintf( ioQQQ, " H%3ld matrix not used, H too neutral\n", ipZ );
		}
		if( !HlteOk.lgHlteOk[ipZ] )
		{
			fprintf( ioQQQ, " H%3ld matrix not used, Hlte not ok\n", ipZ );
		}
	}

	/* following is case b rec coef */
	hcaseb.HRecEffec[ipZ] = 
		hydro.hrec[ipZ][IP1S][ipRecRad]*
	  hydro.hrec[ipZ][IP1S][ipRecNetEsc] + hcaseb.HRecRad[ipZ] + 
	  rcota.CotaRate[ipZ];

	/* do radiative only cascades */
	topoff = 0.;
	for( lev=nhlevel; lev >= IP2S; lev-- )
	{
		/* radiative rec plus three body rec */
		above = (hydro.hrec[ipZ][lev][ipRecRad]*
			hydro.hrec[ipZ][lev][ipRecNetEsc] + 
		   hydro.hcolnc[ipZ][lev]*hcolrt.hlte[ipZ][lev]*phycon.eden)*
	  	   phycon.eden;

		topoff += above;

		for( iup=lev + 1; iup <= nhlevel; iup++ )
		{
			/* hlev(lev,low) is A * esc prob */
			above += hydro.hn[ipZ][iup]*(hydro.HRadNet[ipZ][lev][iup] + 
			  hydro.hcol[ipZ][lev][iup]*phycon.eden);
		}

		below = 0.;
		for( low=IP1S; low <= (lev - 1); low++ )
		{
			/* HRadNet is (escape+destruc)*A for upper to lower */
			below += hydro.HRadNet[ipZ][low][lev] + hydro.hcol[ipZ][low][lev]*
			  phycon.eden;
			/* hcol(up,low) is coll deec rate coef */
		}
		if( below > 1e-25 )
		{
			hydro.hn[ipZ][lev] = (float)(above/below);
		}
		else
		{
			hydro.hn[ipZ][lev] = 0.;
		}

		if( hcolrt.hlte[ipZ][lev] > 1e-25 )
		{
			hydro.hbn[ipZ][lev] = 
				hydro.hn[ipZ][lev]/(hcolrt.hlte[ipZ][lev]*
			  phycon.eden);
		}
		else
		{
			hydro.hbn[ipZ][lev] = 0.;
		}
	}

	/* do simple two level for hydrogen, then fill in HN(i) for rest
	 * above in below is tot rec */
	hydro.hn[ipZ][IP1S] = (float)((topoff + CharExc.CTHrec[MIN2(ipZ+1,2)-1][1])/
	  (HPhoto.hgamnc[ipZ][IP1S] + Secondaries.csupra + hydro.hcolnc[ipZ][IP1S]*
	  ElecDen.EdenHCorr + CharExc.CTHion[MIN2(ipZ+1,2)-1][1]));

	if( hcolrt.hlte[ipZ][IP1S] > 1e-25 )
	{
		hydro.hbn[ipZ][IP1S] = 
			hydro.hn[ipZ][IP1S]/(hcolrt.hlte[ipZ][IP1S]*
		  phycon.eden);
	}
	else
	{
		hydro.hbn[ipZ][IP1S] = 0.;
	}

	hydro.hn[ipZ][IP2S] = 
		(float)((hcaseb.HRecRad[ipZ]*0.33*phycon.eden + 
	  hydro.hn[ipZ][IP1S]*Secondaries.Hx12[MIN2(ipZ+1,2)-1][IP2S] + 
	  hydro.hn[ipZ][IP2P]*
	  hydro.hcol[ipZ][IP2S][IP2P]*phycon.eden)/
	  (HydroLines[ipZ][IP2S][IP1S].Aul + 
	  hydro.hcol[ipZ][IP1S][IP2S]*ElecDen.EdenHCorr + 
	  hydro.hcol[ipZ][IP2S][IP2P]*
	  phycon.eden*6.));

	hydro.hn[ipZ][IP2P] = (float)((hcaseb.HRecRad[ipZ]*0.67*phycon.eden + 
	  hydro.hn[ipZ][IP1S]*Secondaries.Hx12[MIN2(ipZ+1,2)-1][IP2P] + 
	  hydro.hn[ipZ][IP2S]*
	  hydro.hcol[ipZ][IP2S][IP2P]*phycon.eden*6.)/
	  (hydro.HRadNet[ipZ][IP1S][IP2P] + 
	  hydro.hcol[ipZ][IP1S][IP2P]*phycon.eden + 
	  hydro.hcol[ipZ][IP2S][IP2P]*
	  phycon.eden));

	/* were upper n levels set? */
	HIonFraccom.HIonFrac[ipZ] = 
		1./(hydro.hn[ipZ][IP1S] + 
		hydro.hn[ipZ][IP2P] + 
	  hydro.hn[ipZ][IP2S] + 
	  hydro.hn[ipZ][3] + hydro.hn[ipZ][4] + 
	  hydro.hn[ipZ][5] + hydro.hn[ipZ][6]);

	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, hydro.hn[ipZ][IP1S], hcaseb.HRecEffec[ipZ], 
		  HPhoto.hgamnc[ipZ][IP1S] );
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "       HLEVEL return, LOW TE used, HII/HI=%10.3e simple=%10.3e IonRate=%10.2e RecCo=%10.2e CharExc,d%10.2e%10.2e\n", 
		  HIonFraccom.HIonFrac[ipZ], HIonFraccom.HIonSimple[ipZ], 
		  HPhoto.hgamnc[ipZ][IP1S] + Secondaries.csupra + hydro.hcolnc[ipZ][IP1S]*
		  phycon.eden, topoff, CharExc.CTHrec[MIN2(ipZ+1,2)-1][1], CharExc.CTHion[MIN2(ipZ+1,2)-1][1] );
	}


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