/* 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 */
/*HydroLevel calls HydroLevelPop or HydroT2Low to solve for ionization balance 
 * and level populations of model hydrogen atom */
/*PrtHydroTrace1 print trace info for hydrogen-like species */
#include "cddefines.h"
#include "taulines.h"
#include "iso.h"
#include "dense.h"
#include "secondaries.h"
#include "trace.h"
#include "elementnames.h"
#include "phycon.h"
#include "dynamics.h"
#include "opacity.h"
#include "ionbal.h"
#include "continuum.h"
#include "hydrogenic.h"
/*lint -e661 Possible access of out-of-bounds pointer*/
/*lint -e662 Possible access of out-of-bounds pointer*/

/*PrtHydroTrace1a print trace info for hydrogen-like species */
static void PrtHydroTrace1a(long nelem )
{
	double colfrc, 
	  phtfrc, 
	  secfrc,
	  collider;

	/* >>chng 05 aug 17, must use real electron density for collisional ionization of H
	 * since in Leiden v4 pdr with its artificial high temperature coll ion can be important
	 * H on H is homonuclear and scaling laws for other elements does not apply 
	 * next two were multiplied by dense.EdenHCorr and changed to dense.eden for H,
	 * EdenHCorr for rest */
	if( nelem==ipHYDROGEN )
	{
		/* special version for H0 onto H0 */
		collider = dense.EdenHontoHCorr;
	}
	else
	{
		collider = dense.EdenHCorr;
	}

	/* identify how atom is ionized for full trace */
	if( iso.xIonSimple[ipH_LIKE][nelem] > 0. )
	{
		/* fraction of ionization due to photoionization */
		phtfrc = iso.gamnc[ipH_LIKE][nelem][ipH1s]/((dense.eden*(iso.RadRec_effec[ipH_LIKE][nelem] + 
			ionbal.CotaRate[nelem]) )*
			iso.xIonSimple[ipH_LIKE][nelem]);

		/* fraction of ionization due to collisional ionization 
		 * >>chng 05 aug 17 from EdenHCorr to collider */
		colfrc = (iso.ColIoniz[ipH_LIKE][nelem][ipH1s]*collider )/
			((dense.eden*(iso.RadRec_effec[ipH_LIKE][nelem] + 
			ionbal.CotaRate[0]) )*
			iso.xIonSimple[ipH_LIKE][nelem]);

		/* fraction of ionization due to secondary ionization */
		secfrc = secondaries.csupra[nelem][nelem]/((dense.eden*(iso.RadRec_effec[ipH_LIKE][nelem] + 
			ionbal.CotaRate[0]) )*
			iso.xIonSimple[ipH_LIKE][nelem]);
	}
	else
	{
		phtfrc = 0.;
		colfrc = 0.;
		secfrc = 0.;
	}

	fprintf( ioQQQ, "     HydroLevel Z=%2ld called, simple II/I=",nelem);
	PrintE93( ioQQQ, iso.xIonSimple[ipH_LIKE][nelem]);
	fprintf( ioQQQ," PhotFrc:");
	PrintE93( ioQQQ,phtfrc);
	fprintf(ioQQQ," ColFrc:");
	PrintE93( ioQQQ,colfrc);
	fprintf( ioQQQ," SecFrc");
	PrintE93( ioQQQ, secfrc);
	fprintf( ioQQQ,"  Te:");
	PrintE93( ioQQQ,phycon.te);
	fprintf( ioQQQ," eden:");
	PrintE93( ioQQQ,dense.eden);
	fprintf( ioQQQ,"\n"); 
	return;
}

/*PrtHydroTrace1 print trace info for hydrogen-like species */
static void PrtHydroTrace1(long nelem )
{
	long int ipHi , ipLo , i;

	fprintf( ioQQQ, 
		"       HydroLevel%3ld finds arrays, with optical depths defined? %1c induced 2ph=%12.3e\n", 
		nelem, TorF(opac.lgTauOutOn), EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s].pump );
	for( ipHi=ipH2s; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
	{
		fprintf( ioQQQ, "up:%2ld", ipHi );
		fprintf( ioQQQ, "lo" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ, "%9ld", ipLo );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " A*esc" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e",  EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul*
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pesc ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " A*ees" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e",  EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul*
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pelec_esc ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " tauin" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e",  EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauIn ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " t tot" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauTot ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " Esc  " );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pesc ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " Eesc " );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pelec_esc ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " Dest " );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e",  EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pdest) );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " A*dst" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e",  EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul*
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pdest ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " StrkE" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e",  hydro.pestrk[ipLo][ipHi] ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " B(ul)" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].pump*
				iso.stat[ipH_LIKE][nelem][ipLo]/iso.stat[ipH_LIKE][nelem][ipHi] ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " tcont" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e",  EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauCon ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "%3ld", ipHi );
		fprintf( ioQQQ, " C(ul)" );
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			fprintf( ioQQQ,PrintEfmt("%9.2e", 
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].ColUL*dense.eden ));
		}
		fprintf( ioQQQ, "\n" );

		if( ipHi == 2 )
		{
			fprintf( ioQQQ, "    FeIIo");
			fprintf( ioQQQ,PrintEfmt("%9.2e", 
				hydro.dstfe2lya* EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Aul ));
			fprintf( ioQQQ, "\n");
		}
	}

	fprintf( ioQQQ, "         " );
	for( i=1; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
	{
		fprintf( ioQQQ, "%9ld", i );
	}
	fprintf( ioQQQ, "\n" );
	return;
}


/*HydroLevel calls HydroLevelPop or HydroT2Low to solve for ionization balance 
 * and level populations of model hydrogen atom */
void HydroLevel(long int nelem)
{
	long int i, 
	  ipHi, 
	  ipLo, 
	  level;
	double 
	  sum_popn_ov_ion,
	  collider;

	double TooSmall;

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

	/* check that we were called with valid charge */
	ASSERT( nelem >= 0);
	ASSERT( nelem < LIMELM );
	/* >>chng 05 aug 17, must use real electron density for collisional ionization of H
	 * since in Leiden v4 pdr with its artificial high temperature coll ion can be important
	 * H on H is homonuclear and scaling laws for other elements does not apply 
	 * next two were multiplied by dense.EdenHCorr and changed to dense.eden for H,
	 * EdenHCorr for rest */
	if( nelem==ipHYDROGEN )
	{
		/* special version for H0 onto H0 */
		collider = dense.EdenHontoHCorr;
	}
	else
	{
		collider = dense.EdenHCorr;
	}

	/* these two collision rates must be the same or we are in big trouble,
	 * since used interchangably */
	ASSERT( ionbal.CollIonRate_Ground[nelem][nelem][0]< SMALLFLOAT ||
		fabs( iso.ColIoniz[ipH_LIKE][nelem][0]* collider /
		SDIV(ionbal.CollIonRate_Ground[nelem][nelem][0] ) - 1.) < 0.001 );

	/* option to print some rates */
	if( (trace.lgTrace && trace.lgIsoTraceFull[ipH_LIKE]) && (nelem == trace.ipIsoTrace[ipH_LIKE]) )
		PrtHydroTrace1(nelem);

	if( trace.lgHBug && trace.lgTrace )
		PrtHydroTrace1a(nelem);

	/* >>chng 03 apr 30, different limit for He itself since so important in molecular regions */
	/* >>chng 04 mar 09, above comment was in place, but there was no code, def had been just 1e-10 */
	if( nelem==ipHYDROGEN )
		TooSmall = 1e-28;
	else if( nelem==ipHELIUM )
		TooSmall = 1e-18;
	else
		TooSmall = 2e-18;

	/* for dense models close to lte all ionizations can be from
	 * excited states, and simple and actual pops are very different.
	 * code used simple all the time, caused catastrophe for Kingdon nova model.
	 * now use actual pops if we are into model */

	/* which case atom to solve??? */
	/* >>chng 03 apr 11, from ae-30 bail, to ae-35, will catch it in hydro 2 low */
	/* iso.xIonSimple[ipH_LIKE][nelem] was set in iso_ionize_recombine */
	if( iso.xIonSimple[ipH_LIKE][nelem] < 1e-35 )
	{
		/* don't bother if no ionizing radiation */
		strcpy( iso.chTypeAtomUsed[ipH_LIKE][nelem], "zero " );
		if( trace.lgHBug && trace.lgTrace )
		{
			fprintf( ioQQQ, "     HydroLevel Z=%2ld simple II/I=%10.2e so not doing equilibrium.\n", 
			  nelem, iso.xIonSimple[ipH_LIKE][nelem] );
		}

		for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
		{
			iso.DepartCoef[ipH_LIKE][nelem][i] = 1.;
			iso.Pop2Ion[ipH_LIKE][nelem][i] = 0.;
		}

		ionbal.RateRecomTot[nelem][nelem] = 0.;
		iso.xIonSimple[ipH_LIKE][nelem] = 0.;
		iso.pop_ion_ov_neut[ipH_LIKE][nelem] = 0.;

	}
	/* >>chng 99 nov 23, very dense model close to lte had very low simple
	 * ionization fraction but moderate ionization since all pops in
	 * excited states - added test for density 
	 * second change - for all cases use this routine if ionization is
	 * very low indeed */
	/* >>chng 00 aug 26, added matrix force option,
	 * logic is to override default options, "none was set in zero */

	/* NB - this test is mostly bogus since chTypeAtom is "POPU" in zero */

	/* >>chng 02 mar 13, iso.chTypeAtomSet had been set to POPU in zero, look for
	 * comment with this date.  This had the effect of killing this following test.
	 * This also meant that the code had been happily inverting these impossible
	 * matrices for quite some time.  This test changed to stringest one, in
	 * light of this */
	/* >>chng 02 jun 29 add element specific test on ionization limit - had just
	 * been any element < 1e-28, but he ii failed at 1e-21 */
	/* forced with atom h-like matrix lowt */
	/*>>chng 04 apr 17, introduce TooSmall to simplify logic, also intro
	 * chTypaAtomSet option with atom he-like matrix */
	else if( (strcmp( iso.chTypeAtomSet[ipH_LIKE] , "LOWT" )==0) ||
		(iso.xIonSimple[ipH_LIKE][nelem] < TooSmall) )
#		if 0
		( (nelem==ipHYDROGEN) && (iso.xIonSimple[ipH_LIKE][nelem] < 1e-28) ) ||
		( (nelem==ipHELIUM) && (iso.xIonSimple[ipH_LIKE][nelem] < 1e-18) ) ||
		/*>>chng 04 mar 22,break heaver than he out, change from 1 to 2 e-18,
		 * fl had neg pop at 1.05e-18 */
		( (nelem>ipHELIUM) && (iso.xIonSimple[ipH_LIKE][nelem] < 2e-18) ) )
#		endif
	{
		strcpy( iso.chTypeAtomUsed[ipH_LIKE][nelem], "Low T" );
		/* do this, since we will force answer to this */
		/*iso.xIonSimple[ipH_LIKE][nelem] = iso.xIonSimple[ipH_LIKE][nelem];*/
		/* this avoids departure coefficients, 
		 * by doing simple cascade, should always work, 
		 * but not very well at high densities, the ionization ratio WILL be equal to
		 * iso.,xIonSimple */
		HydroT2Low( nelem );
	}

	 /* >>chng 00 aug 26, added matrix force option,
	 * logic is to override default options, "none was set in zero */
	/* >>chng 01 may 09, remove all logic - always dumps to HydroLevelPop if gotten this far */
	else
	{
		/* this is the branch that uses the matrix to get level populations */
		strcpy( iso.chTypeAtomUsed[ipH_LIKE][nelem], "Popul" );
		HydroLevelPop(nelem );
	}
	/* >>chng 01 may 09, totally remove HydroLevelDep */
	/* end branch for main solve */

	/* all three cases end up down here 
	 * numLevels is set in cddefines.c, and is 26 for H and he 
	 * nTopOff is one less, or 25 for these two,
	 * numbers are ten less for heavier h-like species */
	for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
		{
			/* save old PopOpc */
			double PopOpcSave = EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc;

			/* population of lower level rel to ion */
			EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopLo = 
				iso.Pop2Ion[ipH_LIKE][nelem][ipLo];

			/* population of upper level rel to ion */
			EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi = 
				iso.Pop2Ion[ipH_LIKE][nelem][ipHi];

			/* population of lower level rel to ion, corrected for stim em */
			EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc = 
				iso.Pop2Ion[ipH_LIKE][nelem][ipLo] - EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi*
				iso.stat[ipH_LIKE][nelem][ipLo]/iso.stat[ipH_LIKE][nelem][ipHi];

			/* >>chng 03 may 07, add this case for optically thick higher
			 * transitions where stim emiss correction almost cancels out lower pop */
			/* this is special case where lower and upper pops nearly cancel in
			 * the PopOpc difference - then use mean of old and new,
			 * do not include 2s 2p since PopOpc is actually sum of
			 * both PopLo */
			/* >>chng 04 dec 07, add test on tau in > 0.01 to avoid this step
			 * when optically thin */
			if( EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauIn > 0.01 &&
				ipLo>2 && EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc /
				SDIV( EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopLo ) <0.9 )
			{
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc = (0.75*PopOpcSave+
					0.25*EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc);
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc = MIN2(
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc ,
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopLo );
				/*fprintf(ioQQQ,"DEBUG hydrolevel %e\n", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc );*/
			}
		}
	}
	/* >>chng 00 dec 18, do not let lines going to topoff level mase
	 * higest level can have artifical maser since added recombination coefficient
	 * goes into such high levels 
	 * >>chng 01 feb 24, this had been in the wrong place in hydrogenic.c, having not effect.
	 * moved here */
	/* >>chng 03 jun 09, do not include stimulated emission when upper level affected
	 * by recombination topoff */
	/*if( nelem==0 )
		fprintf(ioQQQ,"hlevelll\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e \n",
		EmisLines[ipH_LIKE][nelem][25][24].PopOpc,
		EmisLines[ipH_LIKE][nelem][25][24].PopLo,
		EmisLines[ipH_LIKE][nelem][24][23].PopOpc,
		EmisLines[ipH_LIKE][nelem][24][23].PopLo,
		EmisLines[ipH_LIKE][nelem][23][22].PopOpc,
		EmisLines[ipH_LIKE][nelem][23][22].PopLo,
		EmisLines[ipH_LIKE][nelem][24][21].PopOpc,
		EmisLines[ipH_LIKE][nelem][24][21].PopLo);*/
	for( ipHi=iso.nTopOff[ipH_LIKE][nelem]-1; ipHi<iso.numLevels_max[ipH_LIKE][nelem]; ++ipHi )
	{
		for( ipLo=ipH1s; ipLo < ipHi; ++ipLo )
		{
			/* population of lower level rel to ion, corrected for stim em */
			EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc =  
				EmisLines[ipH_LIKE][nelem][iso.numLevels_max[ipH_LIKE][nelem]-1][ipLo].PopLo;
		}
	}
	/* >>chng 01 feb 24, 2s and 2p are degenerate, so make PopOpc reflect total population. 
	 * this had major effect on radiation pressure due to n=3 to 2p transition, since 
	 * pop of 2p alone is small due to rapid escape route, so radiation pressure was large */
	for( ipHi=ipH2p+1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
	{
		/* sum of populations in 2s and 2p */
		double sum = EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].PopOpc + 
			EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].PopOpc;

		/* population of lower level rel to ion, corrected for stim em */
		/* note that this logic assumes that upper level is the same,
		 * which is it for H since only 2s and 2p are l-resolved */
		EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].PopOpc =  sum;
		EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].PopOpc =  sum;
	}

	/* both solvers end up here */
	/* sum_popn_ov_ion will become the ratio of ionized to neutral hydrogenic 
	 * species, create sum of level pops per ion first */
	sum_popn_ov_ion = 0.;

	/* this is charge transfer ionization of this specieds by hydrogen and helium */
	ionbal.RateIonizTot[nelem][nelem] = 0.;

	for( level=ipH1s; level < iso.numLevels_max[ipH_LIKE][nelem]; level++ )
	{
		/* have solution, now sum all ionization processes from each level */
		ionbal.RateIonizTot[nelem][nelem] += 
			iso.RateLevel2Cont[ipH_LIKE][nelem][level]*iso.Pop2Ion[ipH_LIKE][nelem][level];

		sum_popn_ov_ion += iso.Pop2Ion[ipH_LIKE][nelem][level];
	}

	/* convert back to scaled from ground */
	if( ( ionbal.RateIonizTot[nelem][nelem]/MAX2(SMALLDOUBLE , sum_popn_ov_ion) ) > BIGDOUBLE )
	{
		fprintf( ioQQQ, "RateIonizTot for Z=%li, ion %li is larger than BIGDOUBLE.  This is a problem.",
			nelem+1, nelem);
		cdEXIT(EXIT_FAILURE);
	}
	else
		ionbal.RateIonizTot[nelem][nelem] /= MAX2(SMALLDOUBLE , sum_popn_ov_ion);

	/* all solvers come here - this branch to catch insanity */
	if( sum_popn_ov_ion < 0. )
	{
		fprintf( ioQQQ, 
			" HydroLevel finds negative H-like ion fraction for nelem=%2ld %s using routine %s, val= %.3e simple=%.3e TE=%.3e ZONE=%4ld\n", 
		  nelem, 
		  elementnames.chElementSym[nelem],
		  iso.chTypeAtomUsed[ipH_LIKE][nelem] , 
		  sum_popn_ov_ion, 
		  iso.xIonSimple[ipH_LIKE][nelem], 
		  phycon.te, 
		  nzone );
		fprintf( ioQQQ, " level pop are:" );
		for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
		{
			fprintf( ioQQQ,PrintEfmt("%8.1e", iso.Pop2Ion[ipH_LIKE][nelem][i] ));
		}
		fprintf( ioQQQ, "\n" );
		ContNegative();
		ShowMe();
		puts( "[Stop in HydroLevel]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* all solvers come here, two cases, too small, and normal pops
	 * this branch is for trivial abundance, will zero out species */
	else if( sum_popn_ov_ion >= 0. && sum_popn_ov_ion < 1e-30 )
	{
		iso.pop_ion_ov_neut[ipH_LIKE][nelem] = 0.;

		/* >>chng 03 aug 16, zero out parent density - had not been done */
		dense.xIonDense[nelem][nelem+1] = 0.;

		/* reset pointer to one lower stage of ionization so this not
		 * considered again, hydrogenic considered if IonHigh is nelem+2 */
		dense.IonHigh[nelem] = nelem;
		ionbal.RateIonizTot[nelem][nelem] = 0.;

		/* now zero this out */
		for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
			{
				/* population of lower level rel to ion */
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopLo = 0.;

				/* population of upper level rel to ion */
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi =  0.;

				/* population of lower level rel to ion, corrected for stim em */
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc =  0.;

				/* local ots destruction, cm^-3 s^-1 */
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].ots =  0.;
			}
		}
	}
	/* following is the main branch that applies when we have non-trivial abundance */
	else
	{
		/* level populations 
		 * this brance we have significant population as sum of levels per ion,
		 * store inverse, ion per atom, here.  For non-H species this is ratio
		 * of highest to next highest stage of ionization */
		iso.pop_ion_ov_neut[ipH_LIKE][nelem] = 1./sum_popn_ov_ion;

		/* >>chng 02ionrec.RateIonizTot[nelem][nelem] Sep 20 rjrw: RateRecomTot from total not ratio */
		/* >>chng 02 oct 01, use recomtot when dynamics enabled, but the ratio of
		 * ionization rate to neutral frac for other cases.  This needs to be sorted out. */
#		ifndef NDEBUG
		{
			double pop_ion_ov_neut = ionbal.RateIonizTot[nelem][nelem] / 
				SDIV(ionbal.RateRecomTot[nelem][nelem]);

			double ratio_error = fabs(iso.pop_ion_ov_neut[ipH_LIKE][nelem] - pop_ion_ov_neut ) /
				/* >>chng 03 mar 20,
				 * first term below changed from SMALLFLOAT to 1e-8, to prevent
				 * checks on very minor species, as per Robin Williams email.
				 * with old simple bidiagonal solver, each ionization pair had
				 * very high accuracy, not related to the result itself.  Now
				 * ionization is result of soln with multiple ionization stages
				 * being directly coupled, so this level of error exists 
				 * across solution */
				MAX2(1e-7,iso.pop_ion_ov_neut[ipH_LIKE][nelem]);

			ASSERT( dynamics.lgAdvection ||
				/* following needed to pass assert on alpha */
				/* [nelem][nelem] not a bug - that arrray is on elem, ion scale */
				ionbal.RateIonizTot[nelem][nelem] < SMALLFLOAT || 
				 /* >>chng 04 apr 15, add 1e-5, to prevent crash on very trivial H ion */
				iso.pop_ion_ov_neut[ipH_LIKE][nelem] < 1e-5 ||
				ratio_error < 0.001 );
#				if 0
				||
				/* when continuum pumping is enhanced we can overpopulate highly excited levels
				 * which are then collsionally ionized */ 
				 hydro.xLymanPumpingScaleFactor > 1. );
#				endif
		}
#		endif
	}

	/* this is main trace h-like printout */
	if( (trace.lgIsoTraceFull[ipH_LIKE] && trace.lgTrace) && (nelem == trace.ipIsoTrace[ipH_LIKE]) )
	{
		fprintf( ioQQQ, "       HLEV HGAMNC" );
		PrintE93( ioQQQ, iso.gamnc[ipH_LIKE][nelem][ipH1s] );
		for( i=ipH2s; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso.gamnc[ipH_LIKE][nelem][i] ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "       HLEV TOTCAP" );
		for( i=1; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso.RateCont2Level[ipH_LIKE][nelem][i]/dense.eden ));
		}
		fprintf( ioQQQ," tot");
		fprintf( ioQQQ,PrintEfmt("%10.2e", ionbal.RateRecomTot[ipH_LIKE][nelem]/dense.eden ) );
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "       HLEV IND Rc" );
		for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso.RecomInducRate[ipH_LIKE][nelem][i]*iso.PopLTE[ipH_LIKE][nelem][i] ));
		}
		fprintf( ioQQQ, "\n" );

		/* incuded recombination rate coefficients */
		fprintf( ioQQQ, "       IND Rc LTE " );
		for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e",
				iso.gamnc[ipH_LIKE][nelem][i]*iso.PopLTE[ipH_LIKE][nelem][i] ));
		}
		fprintf( ioQQQ, "\n" );

		/* LTE level populations */
		fprintf( ioQQQ, "       HLEV   HLTE" );
		for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso.PopLTE[ipH_LIKE][nelem][i] ));
		}
		fprintf( ioQQQ, "\n" );

		/* fraction of total ionization due to collisions from given level */
		fprintf( ioQQQ, "       HLEVfr cion" );
		for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", 
				iso.ColIoniz[ipH_LIKE][nelem][i]*
			  iso.Pop2Ion[ipH_LIKE][nelem][i]*collider/MAX2(1e-30,iso.RateLevel2Cont[ipH_LIKE][nelem][i]) ) );
		}
		fprintf( ioQQQ, "\n" );

		/* fraction of total ionization due to photoionization from given level */
		if( ionbal.RateRecomTot[nelem][nelem]> 0. )
		{
			fprintf( ioQQQ, "       HLEVfrPhIon" );
			for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", 
					iso.gamnc[ipH_LIKE][nelem][i]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][i]/MAX2(1e-30,iso.RateLevel2Cont[ipH_LIKE][nelem][i]) ) );
			}
			fprintf( ioQQQ, "\n" );
		}

		fprintf( ioQQQ, "       HLEV     HN" );
		for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipH_LIKE][nelem][i] ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "       HLEV   b(n)" );
		for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
		{
			fprintf(ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipH_LIKE][nelem][i] ));
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "       HLEV   X12tot");
		fprintf(ioQQQ,PrintEfmt("%9.2e", secondaries.x12tot ));
		fprintf( ioQQQ," Grn dest:");
		fprintf(ioQQQ,PrintEfmt("%9.2e",
		  ionbal.RateIonizTot[nelem][nelem] ));
		fprintf(ioQQQ, "\n"); 
	}

	/* find totaql ionization rate while checking for non-positive level populations */
	/* >>chng 02 apr 17, only check on negative pops, alpha turned up some zero pops
	 * in extreme low ionization conditions */
	if( iso.pop_ion_ov_neut[ipH_LIKE][nelem] > 0. )
	{
		for( level=ipH1s; level < iso.numLevels_max[ipH_LIKE][nelem]; level++ )
		{

			/* >>chng 02 apr 17, only check on negative pops, alpha turned up some zero pops
			* in extreme low ionization conditions */
			/*if( iso.Pop2Ion[ipH_LIKE][nelem][level] <= 0. )*/
			if( iso.Pop2Ion[ipH_LIKE][nelem][level] < 0. )
			{
				fprintf( ioQQQ, 
					" HydroLevel finds negative hydrogen level population for %s nelem=%ld level %ld value=%.3e simple=%.3e TE=%.3e ZONE=%4ld\n", 
				  elementnames.chElementName[nelem], 
				  nelem, level, 
				  iso.Pop2Ion[ipH_LIKE][nelem][level], 
				  iso.xIonSimple[ipH_LIKE][nelem], 
				  phycon.te, nzone );
				fprintf( ioQQQ, " level pop are:" );
				for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
				{
					fprintf( ioQQQ,PrintEfmt("%8.1e", iso.Pop2Ion[ipH_LIKE][nelem][i] ));
				}
				fprintf( ioQQQ, "\n" );
				ContNegative();
				ShowMe();
				puts( "[Stop in HydroLevel]" );
				cdEXIT(EXIT_FAILURE);
			}
		}
	}
	else if( iso.pop_ion_ov_neut[ipH_LIKE][nelem] < 0.)
	{
		ionbal.RateIonizTot[nelem][nelem] = 0.;
		/* this is case where we expected ionization, but found none, or negative */
		fprintf( ioQQQ, 
			" HydroLevel finds negative hydrogen ion frac for nelem=%ld (%s) value=%.3e simple= %.3e TE=%.3e ZONE=%ld\n", 
		  nelem, 
		  elementnames.chElementSym[nelem],
		  iso.Pop2Ion[ipH_LIKE][nelem][level], 
		  iso.xIonSimple[ipH_LIKE][nelem], 
		  phycon.te, nzone );
		fprintf( ioQQQ, " level pop are:" );
		for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
		{
			fprintf( ioQQQ,PrintEfmt("%8.1e", iso.Pop2Ion[ipH_LIKE][nelem][i] ));
		}
		fprintf( ioQQQ, "\n" );
		ContNegative();
		ShowMe();
		puts( "[Stop in HydroLevel]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* the case where no ionization is present was covered above, there is no
	 * final else for these tests */

	/* >>chng 05 mar 24,
	 * renormalize the populations and emission of H atom to agree with chemistry 
	 * these were not being kept parallel with chemistry, and caused large changes in O+
	 * abundance when finally done */
	HydroRenorm();

	if( trace.lgTrace )
	{
		/* iso.RecomTotal[nelem] is gross rec coef, computed here while filling in matrix
		 * elements, all physical processes included. 
		 * RadRec_effec is total effective radiative only */
		fprintf( ioQQQ, "       HydroLevel%3ld retrn %s te=",
			nelem, 
			iso.chTypeAtomUsed[ipH_LIKE][nelem] );
		PrintE93( ioQQQ,phycon.te);
		fprintf( ioQQQ," ion/atom=%.4e",iso.pop_ion_ov_neut[ipH_LIKE][nelem]);

		fprintf( ioQQQ," simple=%.4e",iso.xIonSimple[ipH_LIKE][nelem]);

		fprintf( ioQQQ," b1=%.2e",iso.DepartCoef[ipH_LIKE][nelem][ipH1s]);

		fprintf( ioQQQ," ion rate=%.4e",ionbal.RateIonizTot[nelem][nelem]);

		fprintf( ioQQQ," TotRec=%.4e",ionbal.RateRecomTot[nelem][nelem]);

		fprintf( ioQQQ," RadRec=%.4e",iso.RadRec_effec[ipH_LIKE][nelem]);
		fprintf( ioQQQ, "\n");
	}
#if 0
	{
		static int ncall=0;
		if( nelem==0)
			++ncall;
		if( ncall > 11 )
		{
			fprintf(ioQQQ,"debugggg exit\n");
			exit(1);
		}
	}
#endif

	ASSERT( fabs( EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].PopLo - iso.Pop2Ion[ipH_LIKE][nelem][0] ) /
		SDIV( iso.Pop2Ion[ipH_LIKE][nelem][0] ) < 1e-4 || EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].PopLo < 1e-30 );
#	if !defined(NDEBUG)
	for( ipHi=ipH2s; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
	{
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			/* make sure that populations are valid, first check Pop2Ion 
			 * test on 2s and 2p due to brining PopOpc together for optical depth scale being degenerate */
			ASSERT(  (ipLo==ipH2s || ipLo==ipH2p) ||
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc <=
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopLo );
			ASSERT(  EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc <BIGFLOAT );
		}
	}
#	endif

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

