/* 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 */
/*HeLikeLevel level populations  */
#include "cddefines.h" 
#include "lapack.h"
#include "trace.h"
#include "secondaries.h"
#include "atmdat.h"
#include "phycon.h"
#include "ionbal.h"
#include "dense.h"
#include "rfield.h"
#include "he.h"
#include "iso.h"
#include "helike.h"
#include "co.h"
#include "taulines.h"
#include "physconst.h" 
#include "dynamics.h"
#include "continuum.h"

/*lint -e662 out of bounds pointer */
/*===================================================================================*/
/* solve for level populations  */
void HeLikeLevel( long nelem)
{
	int32 nerror;
	long int n, i, ipHi, ipLo, j, 
		level;
	double HighestPColOut[2] = {0.,0.}, 
		sum_popn_ov_ion,
		TooSmall;
	int lgNegPop;
	static int lgAMAT_alloc=FALSE;
	static double *amat;
	static double SaveHe23S_photorate=0.;
	static int32 *ipiv ; /* malloc out to [helike.numLevels[nelem]] */
	/* this block of variables will be obtained and freed here */
	static double
		*creation ,
		*error/*[helike.numLevels[nelem]]*/, 
		*work/*[helike.numLevels[nelem]]*/, 
		**SaveZ/*[helike.numLevels[nelem]+1][iso.numLevels[ipHE_LIKE]]*/, 
		**z/*[helike.numLevels[nelem]][iso.numLevels[ipHE_LIKE]]*/;
	static double
		*CollisionsGoingUp,	*CollisionsGoingDown,
		*OutOfNGoingUp,	*OutOfNGoingDown;
	static int lgSpaceMalloc=FALSE;
	static FILE *matrix;

	/* this flag true means to malloc and free all space each time */
#	define PARALLEL	FALSE
	
#	ifdef DEBUG_FUN
	fputs( "<+>HeLikeLevel()\n", debug_fp );
#	endif

	ASSERT( nelem < LIMELM );

	/* malloc some scratch space */
	if( PARALLEL )
	{
		if( (ipiv = (int32 *)MALLOC(sizeof(int32)*(unsigned)helike.numLevels[nelem] ) )==NULL )
			BadMalloc();

		if( (creation = (double *)MALLOC(sizeof(double)*(unsigned)helike.numLevels[nelem] ) )==NULL )
			BadMalloc();

		if( (CollisionsGoingUp = (double *)MALLOC(sizeof(double)*(unsigned)max_num_levels ) )==NULL )
			BadMalloc();

		if( (CollisionsGoingDown = (double *)MALLOC(sizeof(double)*(unsigned)max_num_levels ) )==NULL )
			BadMalloc();

		if( (OutOfNGoingUp = (double *)MALLOC(sizeof(double)*(unsigned)(max_n + 1) ) )==NULL )
			BadMalloc();

		if( (OutOfNGoingDown = (double *)MALLOC(sizeof(double)*(unsigned)(max_n + 1) ) )==NULL )
			BadMalloc();

		if( (error = (double *)MALLOC(sizeof(double)*(unsigned)helike.numLevels[nelem] ) )==NULL )
			BadMalloc();

		if( (work = (double *)MALLOC(sizeof(double)*(unsigned)helike.numLevels[nelem] ) )==NULL )
			BadMalloc();

		/* now do the 2D arrays */
		if( (SaveZ = (double **)MALLOC(sizeof(double *)*(unsigned)helike.numLevels[nelem] ) ) ==NULL )
			BadMalloc();

		if(  (z = (double **)MALLOC(sizeof(double *)*(unsigned)helike.numLevels[nelem] ))==NULL  )
			BadMalloc();
		
		/* now do the second dimension */
		for( i=0; i<(helike.numLevels[nelem]); ++i )
		{
			if( (SaveZ[i] = (double *)MALLOC(sizeof(double)*(unsigned)helike.numLevels[nelem] ) )==NULL )
				BadMalloc();

			if( (z[i] = (double *)MALLOC(sizeof(double)*(unsigned)helike.numLevels[nelem] ) )==NULL )
				BadMalloc();
		}
	}
	else if( !lgSpaceMalloc )
	{
		/* space has not been malloced yet, but will only do it one time */
		lgSpaceMalloc = TRUE;

		if( (ipiv = (int32 *)MALLOC(sizeof(int32)*(unsigned)(max_num_levels) ) )==NULL )
			BadMalloc();

		if( (creation = (double *)MALLOC(sizeof(double)*(unsigned)(max_num_levels) ) )==NULL )
			BadMalloc();

		if( (CollisionsGoingUp = (double *)MALLOC(sizeof(double)*(unsigned)
			(max_num_levels) ) )==NULL )
			BadMalloc();

		if( (CollisionsGoingDown = (double *)MALLOC(sizeof(double)*(unsigned)
			(max_num_levels) ) )==NULL )
			BadMalloc();

		if( (OutOfNGoingUp = (double *)MALLOC(sizeof(double)*(unsigned)(max_n + 1) ) )==NULL )
			BadMalloc();

		if( (OutOfNGoingDown = (double *)MALLOC(sizeof(double)*(unsigned)(max_n + 1) ) )==NULL )
			BadMalloc();

		if( (error = (double *)MALLOC(sizeof(double)*(unsigned)(max_num_levels) ) )==NULL )
			BadMalloc();

		if( (work = (double *)MALLOC(sizeof(double)*(unsigned)(max_num_levels) ) )==NULL )
			BadMalloc();

		/* now do the 2D arrays */
		if( (SaveZ = (double **)MALLOC(sizeof(double *)*(unsigned)max_num_levels ) ) ==NULL )
			BadMalloc();

		if(  (z = (double **)MALLOC(sizeof(double *)*(unsigned)max_num_levels ))==NULL  )
			BadMalloc();

		/* now do the second dimension */
		for( i=0; i<max_num_levels; ++i )
		{
			if( (SaveZ[i] = (double *)MALLOC(sizeof(double)*(unsigned)(max_num_levels+2) ) )==NULL )
				BadMalloc();

			if( (z[i] = (double *)MALLOC(sizeof(double)*(unsigned)(max_num_levels+2) ) )==NULL )
				BadMalloc();
		}
	}

	for( i=0; i < max_num_levels; ++i )
	{
		CollisionsGoingUp[i] = 0.;
		CollisionsGoingDown[i] = 0.;
	}
	sum_popn_ov_ion = 0.;

	for( i=0; i <= max_n; ++i )
	{
		OutOfNGoingUp[i] = 0.;
		OutOfNGoingDown[i] = 0.;
	}

	/* >>chng 03 may 01, move defn of these two up here, before entering matrix */
	/* >>chng 01 jun 21, add */
	/* charge transfer with helium itself, onto everything heavier */
	if( nelem == ipHELIUM )
	{
		/* >>chng 04 mar 09, add this average of old and new */
		/* take average of old and new 23S photo rate - can be dominated by Lya and a bit noisy 
		 * until ots rates have settled down */
		if( nzone>1 )
		{
			double frac=0.5;
			iso.gamnc[ipHE_LIKE][nelem][ipHe2s3S] = frac*iso.gamnc[ipHE_LIKE][nelem][ipHe2s3S] +
				(1.-frac)*SaveHe23S_photorate;
		}
		SaveHe23S_photorate = iso.gamnc[ipHE_LIKE][nelem][ipHe2s3S];
		}

	for( n=0; n < helike.numLevels[nelem]; n++ )
	{
		/* total recombination to level n rate [s-1] */
		creation[n] = iso.RateCont2Level[ipHE_LIKE][nelem][n];
	}

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

	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && nelem==1 )
		{
			fprintf(ioQQQ,"hebugggg\t%li\t%.2e\t%.2e\t%.2e\n",
				nzone,
				iso.xIonSimple[ipHE_LIKE][nelem], 
				ionbal.RateIonizTot[nelem][nelem-1] ,
				ionbal.RateRecomTot[nelem][nelem-1] );
		}
	}

	/* >>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==ipHELIUM )
		TooSmall = 1e-20;
	else
		TooSmall = 1e-10;

	/* branch for very little ionization, use simple estimate but do not do level populations
	 * never set to zero, use simple two-level system, since He+ important for chemistry */
	if( (strcmp( iso.chTypeAtomSet[ipHE_LIKE] , "LOWT" )==0) ||
		iso.xIonSimple[ipHE_LIKE][nelem] < TooSmall )
	{
		strcpy( iso.chTypeAtomUsed[ipH_LIKE], "LowT " );
		/* total ionization will just be the ground state */
		ionbal.RateIonizTot[nelem][nelem-1] = iso.RateLevel2Cont[ipHE_LIKE][nelem][0];
		/* >>chng 01 nov 23, increase upper limit by 2 */
		lgNegPop = FALSE;
		/* >>chng 02 apr 10, had set to zero,
		 * we will leave the level populations at zero but define an ion ratio,
		 * since He+ abundance cannot go to zero in a molecular cloud,
		 * - He+ charge transfer is the main CO destruction mechanism */
		iso.pop_ion_ov_neut[ipHE_LIKE][nelem] = iso.xIonSimple[ipHE_LIKE][nelem];
		iso.Pop2Ion[ipHE_LIKE][nelem][0] =  1./SDIV(iso.pop_ion_ov_neut[ipHE_LIKE][nelem]);
		for( n=1; n < helike.numLevels[nelem]; n++ )
		{
			iso.Pop2Ion[ipHE_LIKE][nelem][n] =  0.;
		}
	}
	else
	{
		strcpy( iso.chTypeAtomUsed[ipHE_LIKE], "level" );
		/* master balance equation, use when significant population */
		for( level=ipHe1s1S; level < helike.numLevels[nelem]; level++ )
		{
			/* all process depopulating level */

			/* there is no attempt to use scaled csupra here, since this is mostly 
			 * ionic and csupra will be zero.  For He itself, other routine is used.
			 * when this routine becomes fundamental for He a scale factor should
			 * use used.  TH85 Tab 10 suggests should be 0.5, but this is for
			 * their background rate - the rate derived here is self-consistent with unity */
			/*z[level][level] = iso.gamnc[ipHE_LIKE][nelem][level] + 
			  iso.ColIoniz[ipHE_LIKE][nelem][level]* dense.EdenHCorr + secondaries.csupra[nelem][nelem-1];*/
			z[level][level] = iso.RateLevel2Cont[ipHE_LIKE][nelem][level];
			
			/* this saves rate level goes to continuum, DO NOT want to include collisional
			 * and radiative bounb - bound transitions that come next 
			iso.RateLevel2Cont[ipHE_LIKE][nelem][level] = z[level][level]; */
			
			
			if( (nelem == ipHELIUM) && (level == ipHe2s3S) )
				helike.qTot2TripS = iso.ColIoniz[ipHE_LIKE][nelem][level]*dense.eden;

			/* all processes populating level from below */
			for( ipLo=ipHe1s1S; ipLo < level; ipLo++ )
			{
				double coll_down = EmisLines[ipHE_LIKE][nelem][level][ipLo].ColUL * dense.EdenHCorr;
				double pump = (double)EmisLines[ipHE_LIKE][nelem][level][ipLo].pump *
					KillIfBelowPlasma(EmisLines[ipHE_LIKE][nelem][level][ipLo].EnergyWN/WAVNRYD);
				double RadDecay = EmisLines[ipHE_LIKE][nelem][level][ipLo].Aul*
					(EmisLines[ipHE_LIKE][nelem][level][ipLo].Pesc + 
					EmisLines[ipHE_LIKE][nelem][level][ipLo].Pelec_esc + 
					EmisLines[ipHE_LIKE][nelem][level][ipLo].Pdest)*
					KillIfBelowPlasma(EmisLines[ipHE_LIKE][nelem][level][ipLo].EnergyWN/WAVNRYD);

				if( helike.lgRandErrGen )
				{
					coll_down *= helike.ErrorFactor[nelem][level][ipLo][ipCollis];
					RadDecay *= helike.ErrorFactor[nelem][level][ipLo][ipRad];
					pump *= helike.ErrorFactor[nelem][level][ipLo][ipRad];
				}

				z[ipLo][level] = 
					-coll_down* 
					(double)EmisLines[ipHE_LIKE][nelem][level][ipLo].gHi/
					(double)EmisLines[ipHE_LIKE][nelem][level][ipLo].gLo*
					iso.Boltzmann[ipHE_LIKE][nelem][level][ipLo]
					-pump;

				/* pumping out of here to lower level */
				/* >>chng 02 feb 16, add casts to double */
				z[level][level] += pump *
				  (double)EmisLines[ipHE_LIKE][nelem][level][ipLo].gLo/
				  (double)EmisLines[ipHE_LIKE][nelem][level][ipLo].gHi;

				/* collisions out of here to lower level */
				z[level][level] += coll_down;

				if( iso.quant_desig[ipHE_LIKE][nelem][level].n > 
					iso.quant_desig[ipHE_LIKE][nelem][ipLo].n )
					CollisionsGoingDown[level] += coll_down;

				/* radiative decays out of here to lower level */
				z[level][level] += RadDecay;

				if( (nelem == ipHELIUM) && iso.quant_desig[ipHE_LIKE][nelem][level].n == helike.n_HighestResolved[nelem]
					&& iso.quant_desig[ipHE_LIKE][nelem][level].l == 1 )
				{
					HighestPColOut[iso.quant_desig[ipHE_LIKE][nelem][level].s] += coll_down;
				}

				if( (nelem == ipHELIUM) && (level == ipHe2s3S) && (iso.quant_desig[ipHE_LIKE][nelem][ipLo].s==0)){
					helike.qTot2TripS += coll_down;
				}

			}

			/* all processes populating level from above */
			/* >>chng 01 nov 23, increase upper limit by two */
			for( ipHi=level + 1; ipHi < helike.numLevels[nelem]; ipHi++ )
			{
				double coll_down = EmisLines[ipHE_LIKE][nelem][ipHi][level].ColUL * dense.EdenHCorr;
				double coll_up;
				double pump = (double)EmisLines[ipHE_LIKE][nelem][ipHi][level].pump *
					KillIfBelowPlasma(EmisLines[ipHE_LIKE][nelem][ipHi][level].EnergyWN/WAVNRYD);
				double RadDecay = EmisLines[ipHE_LIKE][nelem][ipHi][level].Aul*
					(EmisLines[ipHE_LIKE][nelem][ipHi][level].Pesc + 
					EmisLines[ipHE_LIKE][nelem][ipHi][level].Pelec_esc +
					EmisLines[ipHE_LIKE][nelem][ipHi][level].Pdest)*
					KillIfBelowPlasma(EmisLines[ipHE_LIKE][nelem][ipHi][level].EnergyWN/WAVNRYD);

				if( helike.lgRandErrGen )
				{
					coll_down *= helike.ErrorFactor[nelem][ipHi][level][ipCollis];
					RadDecay *= helike.ErrorFactor[nelem][ipHi][level][ipRad];
					pump *= helike.ErrorFactor[nelem][ipHi][level][ipRad];
				}

				coll_up = coll_down * iso.Boltzmann[ipHE_LIKE][nelem][ipHi][level] *
					(double)EmisLines[ipHE_LIKE][nelem][ipHi][level].gHi / 
					(double)EmisLines[ipHE_LIKE][nelem][ipHi][level].gLo;

				z[ipHi][level] = 
					- RadDecay
					- pump*(double)EmisLines[ipHE_LIKE][nelem][ipHi][level].gLo/
					(double)EmisLines[ipHE_LIKE][nelem][ipHi][level].gHi
					- coll_down;

				/* pumping out of here to upper level */
				z[level][level] += pump;

				/* collisions out of here to upper level */
				z[level][level] += coll_up;

				if( iso.quant_desig[ipHE_LIKE][nelem][level].n <
					iso.quant_desig[ipHE_LIKE][nelem][ipHi].n )
				{
					CollisionsGoingUp[level] += coll_up;
				}
				
				if( (nelem == ipHELIUM) && iso.quant_desig[ipHE_LIKE][nelem][level].n == helike.n_HighestResolved[nelem]
					&& iso.quant_desig[ipHE_LIKE][nelem][level].l == 1 )
				{
					HighestPColOut[iso.quant_desig[ipHE_LIKE][nelem][level].s] += coll_up;
				}

				if( (nelem == ipHELIUM) && (level == ipHe2s3S) && (iso.quant_desig[ipHE_LIKE][nelem][ipHi].s==0))
				{
					helike.qTot2TripS += coll_up;
				}
			}
		}

		/* >>chng 02 jul 22, add induced 2-nu to level populations */
		/* induced two photon emission - special because upward and downward are
		 * not related by ratio of statistical weights */
		/* iso.lgInd2nu_On controls on or off, set with SET IND2 ON/OFF command */
		z[ipHe2s1S][ipHe1s1S] -= iso.TwoNu_induc_dn[ipHE_LIKE][nelem]*iso.lgInd2nu_On ;
		z[ipHe1s1S][ipHe2s1S] -= iso.TwoNu_induc_up[ipHE_LIKE][nelem]*iso.lgInd2nu_On ;

		/* rates out of 1s, and out of 2s */
		z[ipHe1s1S][ipHe1s1S] += iso.TwoNu_induc_up[ipHE_LIKE][nelem]*iso.lgInd2nu_On ;
		z[ipHe2s1S][ipHe2s1S] += iso.TwoNu_induc_dn[ipHE_LIKE][nelem]*iso.lgInd2nu_On ;

		/* >>chng 04 apr 12, this was a mistake */
		/*iso.RateLevel2Cont[ipHE_LIKE][nelem][ipHe1s1S] += iso.TwoNu_induc_up[ipHE_LIKE][nelem]*iso.lgInd2nu_On ;*/

		/* add on charge transfer recombination*/
		/* add it to 1s for now, can't really do any better */
		/* >>chng 01 apr 28, had been added to ground, energy conservation off?? */
		/* in current data set (june, 01) next four ct reactions are all zero */
		/* >>chng 01 nov 23, upper limit increased by 2 */
		/*z[helike.numLevels[nelem]][ipHe1s1S] += */
		/* HCharExcIonOf[nelem][ion]*hii  is ion => ion+1 for nelem */
		/* HCharExcRecTo[nelem][ion]*hi is ion+1 => ion of nelem */
		/* >>chng 02 may 02, several logical errors in which ct elements are added */

		/* >>chng 02 Sep 06 rjrw -- all elements have these terms */
		/*>>>chng 02 oct 01, only include if lgAdvection is set */
		if( dynamics.lgAdvection )
		{
			/* add in advection - these terms normally zero */
			/* assume for present that all advection is into ground state */
			creation[ipHe1s1S] += dynamics.Source[nelem][nelem-1]/SDIV(dense.xIonDense[nelem][nelem])*
				dynamics.lgISO[ipHE_LIKE];
			/* >>chng 02 Sep 06 rjrw -- advective term not recombination */
			/* can sink from all components (must do, for conservation) */
			for( ipLo=ipHe1s1S; ipLo<helike.numLevels[nelem]; ++ipLo )
			{
				z[ipLo][ipLo] += dynamics.Rate*dynamics.lgISO[ipHE_LIKE];
			}
			/* >>chng 04 feb 12, should not add the following since treated in ion_solver */
			/*iso.RateLevel2Cont[ipHE_LIKE][nelem][ipHe1s1S] += dynamics.Rate*dynamics.lgISO[ipHE_LIKE];*/
		}

		/* error in how Hx12 is allocated, sec dim has only LIMELM elements,
		 * will crash if following code ever used */
		/* TODO Can't be used for helium sequence because suprathermals rate has not been split up into
		 * individual levels for helium sequence.  So Secondaries structure must be extended to include
		 * helium levels before the below code can be implemented.  This is not a pressing correction
		 * because total suprathermal is already included for the sequence, but this should not be put
		 * off for too long because it should not be too difficult to do correctly.	*/
#		if 0
		if( Secondaries.Hx12[0][1] > 0. )
		{
			/* now add on supra thermal excitation */
			/* >>chng 01 nov 23, increase upper limit by 2 */
			for( level=1; level < helike.numLevels[nelem]; level++ )
			{
				double RateUp , RateDown;

				RateUp = Secondaries.Hx12[MIN2(nelem,1)][level];
				RateDown = RateUp * (double)EmisLines[ipHE_LIKE][nelem][level][ipH1s].gLo /
					(double)EmisLines[ipHE_LIKE][nelem][level][ipH1s].gHi;

				/* stuff in min after Hx12 evaluates to 0 for atom, or 1 for ion */
				/* total rate out of lower level */
				/* TODO, these are all wrong -- see comment above*/
				z[ipH1s][ipH1s] += RateUp;

				/* rate from the upper level to ground */
				z[level][ipH1s] -= RateDown ;

				/* rate from ground to upper level */
				z[ipH1s][level] -= RateUp ;

				z[level][level] += RateDown;  
			}
		}
#		endif

		if( trace.lgTrace && 
			trace.lgIsoTraceFull[ipHE_LIKE] && (nelem == trace.ipIsoTrace[ipHE_LIKE]) )
		{
			fprintf( ioQQQ, "  pop level     others => (HeLikeLevel)\n" );
			/* >>chng 01 nov 23 increase upper limit by 2 */
			for( n=0; n < helike.numLevels[nelem]; n++ )
			{
				fprintf( ioQQQ, "       HII%2ld", n );
				/* limit is <= since last has rec coef */
				/* >>chng 01 nov 23 increase upper limit by 2 */
				for( j=0; j < helike.numLevels[nelem]; j++ )
				{
					fprintf( ioQQQ,"\t%.9e", z[j][n] );
				}
				fprintf( ioQQQ, "\n" );
			}
			fprintf(ioQQQ," recomb          ");
			for( n=0; n < helike.numLevels[nelem]; n++ )
			{
				fprintf( ioQQQ,"\t%.9e", creation[n] );
			}
			fprintf( ioQQQ, "\n" );
			fprintf(ioQQQ," recomb ct %.2e co %.2e hectr %.2e hctr %.2e\n",
				atmdat.HeCharExcRecTotal,
				co.hevmol[ipCO] ,
				atmdat.HeCharExcRecTo[nelem][nelem-1]*dense.xIonDense[ipHELIUM][0] ,
				atmdat.HCharExcRecTo[nelem][nelem-1]*dense.xIonDense[ipHYDROGEN][0] );
		}

		/* save matrix */
		for( ipLo=ipHe1s1S; ipLo < helike.numLevels[nelem]; ipLo++ )
		{
			for( n=ipHe1s1S; n < helike.numLevels[nelem]; n++ )
			{
				SaveZ[n][ipLo] = z[n][ipLo];
			}
		}

		/* this is the default matrix solver */
#		ifdef AMAT
#			undef AMAT
#		endif
		/* >>chng 01 nov 23, increase upper limit by 2 */
#			define AMAT(I_,J_)	(*(amat+(I_)*(helike.numLevels[nelem])+(J_)))
		/* MALLOC space for the  1-d array */
		/* >>chng 02 jun 19, do not always malloc */
		if( PARALLEL )
		{
			if( (amat=(double*)MALLOC( (sizeof(double)*(unsigned)(helike.numLevels[nelem]*helike.numLevels[nelem]) ))) == NULL )
				BadMalloc();
		}
		else if( !lgAMAT_alloc )
		{
			/* this branch - only malloc one time */
			lgAMAT_alloc = TRUE;
			if( (amat=(double*)MALLOC( (sizeof(double)*(unsigned)((max_num_levels)*(max_num_levels)) ))) == NULL )
				BadMalloc();
		}

		/* this is the default */
		for( j=0; j < helike.numLevels[nelem]; j++ )
		{
			for( n=0; n < helike.numLevels[nelem]; n++ )
			{
				AMAT(n,j) = z[n][j];
			}
		}

		nerror = 0;

		getrf_wrapper(helike.numLevels[nelem],helike.numLevels[nelem],
			      amat,helike.numLevels[nelem],ipiv,&nerror);

		/*DGETRF(helike.numLevels[nelem],helike.numLevels[nelem],*/
		/*  amat,helike.numLevels[nelem],ipiv,&nerror1);*/

		getrs_wrapper('N',helike.numLevels[nelem],1,amat,helike.numLevels[nelem],ipiv,
			      creation,helike.numLevels[nelem],&nerror);

		if( nerror != 0 )
		{
			fprintf( ioQQQ, " HeLikeLevel: dgetrs finds singular or ill-conditioned matrix\n" );
			puts( "[Stop in helike]" );
			cdEXIT(EXIT_FAILURE);
		}

		if( PARALLEL )
			free( amat);

		lgNegPop = FALSE;
		for( level=0; level < helike.numLevels[nelem]; level++ )
		{
			iso.Pop2Ion[ipHE_LIKE][nelem][level] = creation[level];

			/* check for negative level populations */
			if( iso.Pop2Ion[ipHE_LIKE][nelem][level] < 0. )
				lgNegPop = TRUE;
		}

		/* zero populations of unused levels. Why are these different? */
		for( level=helike.numLevels[nelem]; level < iso.numLevels[ipHE_LIKE][nelem]; level++ )
		{
			iso.Pop2Ion[ipHE_LIKE][nelem][level] = 0.;
		}

		/* >>chng 02 nov 06, have all levels included in RateIonizTot */
		/* >>chng 02 nov 03, add ct to this term */
		/* this is total ionization of this species referenced to the ground state */
		ionbal.RateIonizTot[nelem][nelem-1] = 0.;
		sum_popn_ov_ion = 0.;

		/* this will eventually become the ratio of ionized to neutral hydrogenic 
		 * species, create sum of level pops per ion first */
		iso.pop_ion_ov_neut[ipHE_LIKE][nelem] = 0.;

		/* will redo recombination rate here, since limit to number of levels
		 * man have changed */
		ionbal.RateRecomTot[nelem][nelem-1] = 0.;
		/* >>chng 01 nov 23, increase upper limit by 2 */
		for( level=0; level < helike.numLevels[nelem]; level++ )
		{

			/* sum of recombination processes */
			ionbal.RateRecomTot[nelem][nelem-1] += iso.RateCont2Level[ipHE_LIKE][nelem][level];

			/* sum of all ionization processes from this atom to ion */
			ionbal.RateIonizTot[nelem][nelem-1] += 
				iso.Pop2Ion[ipHE_LIKE][nelem][level] * iso.RateLevel2Cont[ipHE_LIKE][nelem][level];

			/* create sum of populations,
			 * iso.Pop2Ion is actually populations relative to ion */
			sum_popn_ov_ion += iso.Pop2Ion[ipHE_LIKE][nelem][level];

			if( iso.PopLTE[ipHE_LIKE][nelem][level] > 0. )
			{
				/* the LTE departure coefficients */
				iso.DepartCoef[ipHE_LIKE][nelem][level] = 
					(iso.Pop2Ion[ipHE_LIKE][nelem][level]/
					(iso.PopLTE[ipHE_LIKE][nelem][level]* dense.eden) );
			}
			else
			{
				iso.DepartCoef[ipHE_LIKE][nelem][level] = 0.;
			}
		}

		/* we reset the total recomb rate above, so simple ionization may have changed -
		 * update it */
		iso.xIonSimple[ipHE_LIKE][nelem] = ionbal.RateIonizTot[nelem][nelem-1] / ionbal.RateRecomTot[nelem][nelem-1];

		{
			/*@-redef@*/
			enum {DEBUG_LOC=FALSE};
			/*@+redef@*/
			if( DEBUG_LOC && nzone > 700 && nelem==ipHELIUM )
			{
				fprintf(ioQQQ,"DEBUG hedest\t%.2f",fnzone);
				for( level=0; level < helike.numLevels[nelem]; level++ )
				{
					if( iso.Pop2Ion[ipHE_LIKE][nelem][level] * iso.RateLevel2Cont[ipHE_LIKE][nelem][level] /
						SDIV(ionbal.RateIonizTot[nelem][nelem-1] ) > 0.1 )
						fprintf(ioQQQ,"\t%li\t%.3f",level,
						iso.Pop2Ion[ipHE_LIKE][nelem][level] * iso.RateLevel2Cont[ipHE_LIKE][nelem][level] /
						SDIV(ionbal.RateIonizTot[nelem][nelem-1] ));
				}
				fprintf(ioQQQ,"\n");
			}
		}
		/* convert back to scaled from ground */
		if( ( ionbal.RateIonizTot[nelem][nelem-1]/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-1);
			cdEXIT(EXIT_FAILURE);
		}
		else
			ionbal.RateIonizTot[nelem][nelem-1] /= MAX2(SMALLDOUBLE , sum_popn_ov_ion);

		/* >>chng 02 jan 22, actually use real ratios, not just ground */
		/* this is sum of all populations in model, div by ion 
		 * create ratio of ion pops to total atom */
		if( sum_popn_ov_ion > SMALLDOUBLE )
		{
			iso.pop_ion_ov_neut[ipHE_LIKE][nelem] = 1./sum_popn_ov_ion;
		}
		else
		{
			iso.pop_ion_ov_neut[ipHE_LIKE][nelem] = 0.;
		}
	}/* end of branch that does actual level populations */
	/* end of branch that does full solution with level populations */
	ASSERT( ionbal.RateIonizTot[nelem][nelem-1] >= 0. );
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && nelem==1 )
		{
			fprintf(ioQQQ,"hebugggg2\t%li\t%.2e\t%.2e\t%.2e\n",
				nzone,
				iso.pop_ion_ov_neut[ipHE_LIKE][nelem], 
				ionbal.RateIonizTot[nelem][nelem-1] ,
				ionbal.RateRecomTot[nelem][nelem-1] );
		}
	}

	/* >> chng 02 Sep 20 rjrw: RateRecomTot should be from total not ratio, switched off for time being */
	/* >>chng 02 oct 01, back to old way when dynamics not being done.  Then results of large
	 * atom will feed back into solution in bidiag */
	/*fprintf(ioQQQ,"DEBUG assert\t%li\t%s\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n", 
		nelem, 
		iso.chTypeAtomUsed[ipH_LIKE],
		iso.pop_ion_ov_neut[ipHE_LIKE][nelem] ,
		ionbal.RateIonizTot[nelem][nelem-1] / 
		SDIV(ionbal.RateRecomTot[nelem][nelem-1]),
		iso.xIonSimple[ipHE_LIKE][nelem],
		ionbal.RateIonizTot[nelem][nelem-1],
		ionbal.RateRecomTot[nelem][nelem-1] );*/

	 ASSERT( helike.lgFSM || dynamics.lgAdvection ||
		 iso.pop_ion_ov_neut[ipHE_LIKE][nelem]<SMALLFLOAT ||
		fabs(iso.pop_ion_ov_neut[ipHE_LIKE][nelem] - 
		ionbal.RateIonizTot[nelem][nelem-1] / 
		SDIV(ionbal.RateRecomTot[nelem][nelem-1]) )/
		SDIV(iso.pop_ion_ov_neut[ipHE_LIKE][nelem]) < 0.001);

	/* check on the sum of the populations */
	if( lgNegPop || iso.pop_ion_ov_neut[ipHE_LIKE][nelem] < 0.)
	{
		fprintf( ioQQQ, 
			" HeLikeLevel finds negative He-like ion fraction for nelem=%ld using routine HeLikeLevel, IonFrac= %10.3e simple=%10.3e TE=%10.3e ZONE=%4ld\n", 
		  nelem, 
		  iso.pop_ion_ov_neut[ipHE_LIKE][nelem], 
		  iso.xIonSimple[ipHE_LIKE][nelem], 
		  phycon.te, 
		  nzone );

		fprintf( ioQQQ, " level pop are:\n" );
		for( i=0; i < helike.numLevels[nelem]; i++ )
		{
			fprintf( ioQQQ," %8.1e;", iso.Pop2Ion[ipHE_LIKE][nelem][i] );
			if( (i!=0) && !(i%10) ) fprintf( ioQQQ,"\n" );
		}
		fprintf( ioQQQ, "\n" );
		ContNegative();
		ShowMe();
		puts( "[Stop in helike]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* get level populations, two cases, 
	 * first, may be zero since all cases drop down to here, 
	 * this branch is for trivial abundance, so zero out species */
	/* >>chng 01 apr 21, add second test on high stage being he-like -
	 * there are tests where this is hydrogen like, and jump would be over
	 * two stages of ionization - do not do this */
	/* >>chng 03 apr 30, do not do this test at all, since pops were filled in, for
	 * low ionization limit, above.  important to let helium remain ionized
	 * somewhat in molecular gas */
#	if 0
	if( FALSE && iso.pop_ion_ov_neut[ipHE_LIKE][nelem] >= 0. && 
		iso.pop_ion_ov_neut[ipHE_LIKE][nelem] < TooSmall  &&
		dense.IonHigh[nelem] == nelem )
	{

		/* >>chng 02 apr 10, if actually zero then set to zero, but 
		 * do not kill a very small ion abundance since must let He+
		 * remain despite very low abundance, since important CO destruction
		 * process in molecular clouds */
		if( iso.pop_ion_ov_neut[ipHE_LIKE][nelem] == 0. && !dense.lgSetIoniz[nelem] )
		{
			iso.pop_ion_ov_neut[ipHE_LIKE][nelem] = 0.;
			/* decrement highest stage of ionization, must not
			 * go up to He-like being ionized again */
			dense.IonHigh[nelem] = nelem-1;
			/* >>chng 01 apr 21, added second drop on IonLow -
			 * there was a model where light atoms where stripped and he-like test triggered,
			 * which resulted in insanity downstream since IonHigh and IonLow were equal */
			dense.IonLow[nelem] = MIN2( dense.IonLow[nelem] , dense.IonHigh[nelem]-1 );
		}

		/* now zero out high stages we just cleared out */
		for( ipHi=ipHe2s3S; ipHi < helike.numLevels[nelem]; ipHi++ )
		{
			iso.Pop2Ion[ipHE_LIKE][nelem][ipHi] = 0.;
			for( ipLo=ipHe1s1S; ipLo < ipHi; ipLo++ )
			{
				/* population of lower level rel to ion */
				EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopLo = 0.;

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

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

				/* local ots destruction, cm^-3 s^-1 */
				EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ots =  0.;
			}
		}

		ionbal.RateIonizTot[nelem][nelem-1] = 0.;

		if( trace.lgTrace && trace.lgHeBug )
		{
			fprintf(ioQQQ,"  HeLike for Z=%li finds small population set to zerozero.\n",nelem);
		}
	}
	/* this is the main branch that applies when we have non-trivial abundance */
	else
#	endif
	/* 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 */

	/* now find emission in each line */
	for( ipHi=0; ipHi<helike.numLevels[nelem]; ++ipHi )
	{
		for( ipLo=0; ipLo<ipHi; ++ipLo )
		{
			EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopLo = iso.Pop2Ion[ipHE_LIKE][nelem][ipLo];
			EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi = iso.Pop2Ion[ipHE_LIKE][nelem][ipHi];
			EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopOpc = (iso.Pop2Ion[ipHE_LIKE][nelem][ipLo] - iso.Pop2Ion[ipHE_LIKE][nelem][ipHi]*
				EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gLo/EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gHi);
			EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopOpc = MAX2(SMALLDOUBLE,EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopOpc );
		}
	}

	if( !helike.lgSetBenjamin )
	{
		/* Along the iso-sequence, two of the 2^3Pj terms are always very nearly degenerate, these
		 * should have PopOpc equal to the sum of the two, the third should be left alone.	*/
		for( ipHi=ipHe2p1P; ipHi < iso.numLevels[ipHE_LIKE][nelem]; ipHi++ )
		{
			double sum;
			long OneOfDegenerate, OtherDegenerate;

			OneOfDegenerate = ipHe2p3P1;
				
			if( nelem == ipHELIUM )
				OtherDegenerate = ipHe2p3P2;
			else
				OtherDegenerate = ipHe2p3P0;
			
			sum = EmisLines[ipHE_LIKE][nelem][ipHi][OneOfDegenerate].PopOpc + 
					EmisLines[ipHE_LIKE][nelem][ipHi][OtherDegenerate].PopOpc;

			EmisLines[ipHE_LIKE][nelem][ipHi][OneOfDegenerate].PopOpc = sum; 
			EmisLines[ipHE_LIKE][nelem][ipHi][OtherDegenerate].PopOpc = sum;
		}
	}

	/* helike.lgCritDensLMix is a flag used to print warning if density is
	 * too low for first collapsed level to be l-mixed.  Check is if l-mixing collisions
	 * out of highest resolved singlet P are greater than sum of transition probs out.	*/
 	if( ( (HighestPColOut[0] < 1./helike.Lifetime[ipHELIUM][QuantumNumbers2Index[ipHELIUM][iso.n_HighestResolved[ipHE_LIKE][ipHELIUM]][1][0]] )
		|| (HighestPColOut[1] < 1./helike.Lifetime[ipHELIUM][QuantumNumbers2Index[ipHELIUM][iso.n_HighestResolved[ipHE_LIKE][ipHELIUM]][1][1]] ) )
		&& nelem==ipHELIUM ) 
	{
		helike.lgCritDensLMix = FALSE;
	}

	for( n=0; n < helike.numLevels[nelem]; n++ )
	{
		OutOfNGoingUp[iso.quant_desig[ipHE_LIKE][nelem][n].n] += 
			CollisionsGoingUp[n]*iso.stat[ipHE_LIKE][nelem][n]/(4.*N_(n)*N_(n));
		OutOfNGoingDown[iso.quant_desig[ipHE_LIKE][nelem][n].n] += 
			CollisionsGoingDown[n]*iso.stat[ipHE_LIKE][nelem][n]/(4.*N_(n)*N_(n));
	}

#if	0
	/* TODO at the moment can only assert these for the atom, and for Vriens gbar.	*/
	if( iso.lgColl_excite[ipHE_LIKE] && !helike.lgSetBenjamin && helike.lgCS_Vriens 
		&& !helike.lgRandErrGen && dense.eden > 100. )
	{
		for( n = 4; n < helike.n_HighestResolved[nelem]; n++  )
		{
			double DecayNPlusOne = 1./helike.Lifetime[nelem][QuantumNumbers2Index[nelem][n+1][1][0]];
			
			/* This is not monotonic for the highest two levels...the highest because you can't go up from there!
			 * and the second highest because there is only one level it can go up to, which is too few to make 
			 * the total greater than the third highest!	*/
			/* This doesn't work well for very low density or temperature plasmas, but we don't
			 * particularly care since those have very weak collisional rates.  Only check
			 * these trends if the rates are not much smaller than the corresponding decay rate.	*/

			if( (OutOfNGoingUp[n]>(1E-8*DecayNPlusOne)) && n!=helike.n_HighestResolved[nelem] )
			{
				
				ASSERT( (OutOfNGoingUp[n-1]<OutOfNGoingUp[n]) );
				ASSERT( (OutOfNGoingDown[n]<OutOfNGoingDown[n+1]) );
			}
		}
	}
#endif

	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && nelem == ipHELIUM )
		{
			/* limit output to one run.	*/
			static long RUNTHISONCE = FALSE;

			if( !RUNTHISONCE )
			{
#if	0				
				fprintf( ioQQQ,"n\tl\ts\tCollisionsGoingUp\tCollisionsGoingDown\n");
				
				for( n=0; n < helike.numLevels[nelem]; n++ )
				{
					fprintf( ioQQQ,"%2ld\t%2ld\t%2ld\t%.9e\t%.9e\n",
						iso.quant_desig[ipHE_LIKE][nelem][n].n,
						iso.quant_desig[ipHE_LIKE][nelem][n].l,
						iso.quant_desig[ipHE_LIKE][nelem][n].s,
						CollisionsGoingUp[n]*iso.Pop2Ion[ipHE_LIKE][nelem][n],
						CollisionsGoingDown[n]*iso.Pop2Ion[ipHE_LIKE][nelem][n]	);
				}
#endif

				fprintf( ioQQQ,"n\tOutOfNGoingUp\tOutOfNGoingDown\n");
				
				for( n=1; n <= helike.n_HighestResolved[nelem] + helike.nCollapsed[nelem]; n++ )
				{
					fprintf( ioQQQ,"%2ld\t%.9e\t%.9e\n",
						n, OutOfNGoingUp[n], OutOfNGoingDown[n] );
				}

				fprintf( ioQQQ, "\n" );
			}

			RUNTHISONCE = TRUE;
		}
	}

	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && nelem == ipHELIUM )
		{

			/* limit output to one run.	*/
			static long RUNONCE = FALSE, nlo, nhi;
			/*static double **zPerN;*/

			if( !RUNONCE )
			{
				if ( (matrix = fopen("RateMatrix.txt","w")) == NULL )
					return;

#if	0
				if( (zPerN = (double **)MALLOC(sizeof(double*)*(unsigned)
					(helike.n_HighestResolved[nelem] + helike.nCollapsed[nelem] + 1) ) )==NULL )
					BadMalloc();

				for( i=0; i <= (helike.n_HighestResolved[nelem] + helike.nCollapsed[nelem]); ++i )
				{
					if( (zPerN[i] = (double *)MALLOC(sizeof(double)*(unsigned)
						(helike.n_HighestResolved[nelem] + helike.nCollapsed[nelem] + 1) ) )==NULL )
						BadMalloc();
				}

				/* first zero array zPerN	*/
				for( nlo = 1; nlo <= helike.n_HighestResolved[nelem] + helike.nCollapsed[nelem]; nlo++ )
				{
					for( nhi = 1; nhi <= iso.n_HighestResolved[ipHE_LIKE][ipHELIUM] + iso.nCollapsed[ipHE_LIKE][ipHELIUM]; nhi++ )
					{		
						zPerN[nlo][nhi] = 0.;
					}
				}

				/* then stuff it with values from original z matrix	*/
				for( level=ipHe1s1S; level < helike.numLevels[nelem]; level++ )
				{
					for( ipLo=ipHe1s1S; ipLo < level; ipLo++ )
					{
						zPerN[iso.quant_desig[ipHE_LIKE][nelem][ipLo].n][iso.quant_desig[ipHE_LIKE][nelem][level].n] 
							+= z[ipLo][level]*iso.Pop2Ion[ipHE_LIKE][nelem][ipLo];
					}

					for( ipHi=level; ipHi < helike.numLevels[nelem]; ipHi++ )
					{
						zPerN[iso.quant_desig[ipHE_LIKE][nelem][ipHi].n][iso.quant_desig[ipHE_LIKE][nelem][level].n] 
							+= z[ipHi][level]*iso.Pop2Ion[ipHE_LIKE][nelem][ipHi];
					}
				}


				/* then print it!	*/
				for( nlo = 1; nlo <= helike.n_HighestResolved[nelem] + helike.nCollapsed[nelem]; nlo++ )
				{
					for( nhi = 1; nhi <= helike.n_HighestResolved[nelem] + helike.nCollapsed[nelem]; nhi++ )
					{		
						double f=fabs(zPerN[nlo][nhi]);
						fprintf( matrix, "%.2e\t", log10(SDIV(f)) );
					}
					fprintf( matrix, "\n" );
				}
#endif
				/* then print it!	*/
				for( nlo = 0; nlo < helike.numLevels[nelem]; nlo++ )
				{
					for( nhi = 0; nhi < helike.numLevels[nelem]; nhi++ )
					{	
						fprintf( matrix, "%.2e\t", z[nlo][nhi] );
#if	0
						/*lint -e777 float test equality */
						if( z[nlo][nhi] == SaveZ[nlo][nhi] )
							fprintf( matrix, "%.2e\t", 1.0 );
						else if( SaveZ[nlo][nhi] == 0. )
							fprintf( matrix, "%.2e\t", 0.0 );
						else
							fprintf( matrix, "%.2e\t", z[nlo][nhi]/SaveZ[nlo][nhi] );
						/*lint +e777 float test equality */
#endif
					}
					fprintf( matrix, "\n" );
				}
			
				fprintf( matrix, "\n" );
				fclose( matrix );
			}
			
			RUNONCE = TRUE;
		}
	}

	{
		/* following should be set true to print information */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC /*&& (iteration==2)*/ )
		{
			if( nelem==ipHELIUM )
				fprintf(ioQQQ,"\n");
			fprintf(ioQQQ,"%li\t%.4e", nelem,
				iso.pop_ion_ov_neut[ipHE_LIKE][nelem] );
			if( dense.xIonDense[nelem][nelem]>SMALLFLOAT )
			{
				float ratio = dense.xIonDense[nelem][nelem]/ dense.xIonDense[nelem][nelem-1];

				fprintf(ioQQQ,"\t%.4e\t%.4e", ratio ,
					(ratio-iso.pop_ion_ov_neut[ipHE_LIKE][nelem])/ratio );
			}
			fprintf(ioQQQ,"\n");
		}

	}
	{
		/* following should be set true to print information */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && (nelem==ipCARBON)/**/ )
		{
			{
				fprintf( ioQQQ, "     He-like Z=%2ld H2OVH1=",nelem);
				fprintf( ioQQQ,"%10.3e", iso.pop_ion_ov_neut[ipHE_LIKE][nelem]);
				fprintf( ioQQQ, " simple=");
				fprintf( ioQQQ,"%10.3e", iso.xIonSimple[ipHE_LIKE][nelem]);
				fprintf( ioQQQ," dest, creat rates %.2e %.2e", ionbal.RateIonizTot[nelem][nelem-1] , ionbal.RateRecomTot[nelem][nelem-1]);
				fprintf( ioQQQ, "\n"); 
			}
		}
	}

	/*if( nzone>700 )
		fprintf(ioQQQ,"DEBUG He1\t%.2f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
		fnzone,
		z[ipHe1s1S][ipHe1s1S] ,
		ionbal.RateRecomTot[nelem][nelem-1],
		ionbal.RateIonizTot[nelem][nelem-1],
		iso.Pop2Ion[ipHE_LIKE][nelem][0] , iso.RateLevel2Cont[ipHE_LIKE][nelem][0]);*/
	if( nzone>0 && nelem==ipHELIUM )
	{
		/* find fraction of He0 destructions due to photoionization of 2 ^S */
		double ratio = iso.Pop2Ion[ipHE_LIKE][nelem][ipHe2s3S] * iso.RateLevel2Cont[ipHE_LIKE][nelem][ipHe2s3S] /
			SDIV(iso.Pop2Ion[ipHE_LIKE][nelem][ipHe1s1S]*ionbal.RateIonizTot[nelem][nelem-1]);
		if( ratio > he.frac_he0dest_23S )
		{
			/* remember zone where this happended and fraction, and frac due to photoionization */
			he.nzone = nzone;
			he.frac_he0dest_23S = (float)ratio;
			he.frac_he0dest_23S_photo = (float)(iso.Pop2Ion[ipHE_LIKE][nelem][ipHe2s3S] *iso.gamnc[ipHE_LIKE][nelem][ipHe2s3S]/
				SDIV(iso.Pop2Ion[ipHE_LIKE][nelem][ipHe1s1S]*ionbal.RateIonizTot[nelem][nelem-1]));
		}
	}
	/*fprintf(ioQQQ,"DEBUG heoscil\t%.2f\t%.3e\t%.3e\t%.3e\n",
		fnzone ,
		iso.Pop2Ion[ipHE_LIKE][nelem][1] , 
		iso.RateLevel2Cont[ipHE_LIKE][nelem][1], 
		iso.gamnc[ipHE_LIKE][nelem][ipHe2s3S]);*/
	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "     He-like %.2f Z:%2ld H2OVH1:",fnzone,nelem);
		fprintf( ioQQQ,"%10.3e", iso.pop_ion_ov_neut[ipHE_LIKE][nelem]);
		fprintf( ioQQQ, " simple: %.3e", iso.xIonSimple[ipHE_LIKE][nelem]);
		fprintf( ioQQQ," rate dest:%.2e creat:%.2e", SaveZ[0][0] , ionbal.RateRecomTot[nelem][nelem-1]);
		fprintf(ioQQQ," photo:%.3e coll:%.3e sec:%.3e ",
		iso.gamnc[ipHE_LIKE][nelem][0],
		iso.ColIoniz[ipHE_LIKE][nelem][0]* dense.EdenHCorr,
		secondaries.csupra[nelem][nelem-1]);
		fprintf( ioQQQ, "\n"); 
	}

	/* only release this if option for parallel processing */
	if( PARALLEL )
	{
		free( ipiv );
		free( creation );
		free( OutOfNGoingUp );
		free( OutOfNGoingDown );
		free( CollisionsGoingUp );
		free( CollisionsGoingDown );
		free( error );
		free( work );

		/* now do the second dimension */
		for( i=0; i<helike.numLevels[nelem]; ++i )
		{
			free( SaveZ[i] );

			free( z[i] );
		}

		free( SaveZ );
		free( z );
	}

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

void HeLikeError( long nelem )
{

	long ipHi, ipLo, typeOfRate;
	float ErrorToPut;

	/* put recombination errors into helike.Error array */
	for( ipHi=ipHe1s1S; ipHi<iso.numLevels[ipHE_LIKE][nelem]; ipHi++ )
	{
		if( N_(ipHi)<=5 )
		{
			ErrorToPut = 0.02f;
		}
		else if( N_(ipHi)>iso.n_HighestResolved[ipHE_LIKE][nelem] )
		{
			ErrorToPut = 0.02f;
		}
		else if( L_(ipHi)>=3 )
		{
			ErrorToPut = 0.005f;
		}
		else if( L_(ipHi)==2 )
		{
			ErrorToPut = 0.01f;
		}
		else if( L_(ipHi)==1 )
		{
			ErrorToPut = 0.025f;
		}
		else
		{
			ErrorToPut = 0.05f;
		}
		putError(nelem,iso.numLevels[ipHE_LIKE][nelem],ipHi,ipRad,0.05f);
	}
	/* Error for total recombination */
	putError(nelem,iso.numLevels[ipHE_LIKE][nelem],iso.numLevels[ipHE_LIKE][nelem],ipRad,0.05f);

	/* put collisional errors into helike.Error array */
	/* TODO */

	for( ipHi=ipHe2s3S; ipHi<= iso.numLevels[ipHE_LIKE][nelem]; ipHi++ )
	{
		for( ipLo=ipHe1s1S; ipLo<= iso.numLevels[ipHE_LIKE][nelem]; ipLo++ )
		{
			for( typeOfRate=0; typeOfRate<=1; typeOfRate++ )
			{
				if( helike.Error[nelem][ipHi][ipLo][typeOfRate] >= 0. )
				{
					helike.ErrorFactor[nelem][ipHi][ipLo][typeOfRate] =  
						(float)MyGaussRand( helike.Error[nelem][ipHi][ipLo][typeOfRate] );
					ASSERT( helike.ErrorFactor[nelem][ipHi][ipLo][typeOfRate] > 0. );
				}
				else
				{
					helike.ErrorFactor[nelem][ipHi][ipLo][typeOfRate] = 1.0f;
				}
			}
		}
	}

	return;
}

double MyGaussRand( double PctUncertainty )
{
	double logStdDev;
	double result;

	/*logStdDev = log10( PctUncertainty );*/
	logStdDev = PctUncertainty;

	do
	{
		/*result = pow( 10., RandGauss( 0., logStdDev ) );*/
		result = 1.+RandGauss( 0., logStdDev );
	}
	/* This loop is to exclude results that are different from the mean by 
	 * more than two standard deviations.  This will give us a result that is 
	 * less than or equal to the percent uncertainty about 70% of the time, 
	 * and less than twice the uncertainty the rest of the time. */
	/*while( (result < 1.-2.*PctUncertainty) || (result > 1+2.*PctUncertainty) );*/
	while( (result < 0.) || (result > 2.) );
	
	ASSERT( result>0. && result<2. );

	return result;
}

	

/*lint +e662 out of bounds pointer */

