/* 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 */
/*ion_solver solve the bi-diagonal matrix for ionization balance */
#include "cddefines.h"
#include "yield.h"
#include "prt.h"
#include "continuum.h"
#include "iso.h"
#include "dynamics.h"
#include "grainvar.h"
#include "hmi.h"
#include "mole.h"
#include "thermal.h"
#include "lapack.h"
#include "converge.h"
#include "secondaries.h"
#include "phycon.h"
#include "atmdat.h"
#include "heavy.h"
#include "elementnames.h"
#include "dense.h"
#include "radius.h"
#include "ionbal.h"

void tridiag(double *a, double *b, long int n);
void solveions(double *ion, double *rec, double *snk, double *src,
							 long int nlev, long int nmax);

void ion_solver(
	/* this is element on the c scale, H is 0 */
	long int nelem, 
	/* option to print this element when called */
	int lgPrintIt)
{
	/* use tridiag or general matrix solver? */
	int lgTriDiag=TRUE;
	long int ion, 
	  limit, 
	  IonProduced, 
	  nej, 
	  ns,
	  jmax=-1;
	int lgNegPop;
	double *amat,*achk,*src,*snk, *source,
		rateone;
	int32 nerror;
	int32 *ipiv;
	long ion_low, ion_range, i, j, ion_to , ion_from;
	static double ratio_ion_density,
		sum_dense;
	double auger[LIMELM+1] = {-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,
		-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,
		-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,
		-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,-BIGDOUBLE,
		-BIGDOUBLE};
	static double abund_old[LIMELM][LIMELM+1];

	double abund_total, 
		renorm, 
		renormnew,
		den;
	int lgHomogeneous = TRUE;
	static long int nzone_used[LIMELM];
	int lgFirstPass;
	static int lgMustInit=TRUE;
	double *xmat , *xmatsave;

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

	/*if( nelem==13 )
	{
		broken();dense.IonHigh[nelem] = 1;
	}*/
	if( lgMustInit )
	{
		/* loop is really over nelem not ion */
		for(ion=0; ion<LIMELM; ++ion )
		{
			auger[ion] = -1;
			nzone_used[ion] = -1;
		}
		for( ion=0; ion <dense.IonHigh[nelem]; ion++ )
		{
			abund_old[nelem][ion] = dense.xIonDense[nelem][ion];
		}
	}
	lgMustInit = FALSE;

	/* this is on the c scale, so H is 0 */
	ASSERT( nelem >= 0);
	ASSERT( dense.IonLow[nelem] >= 0 );
	ASSERT( dense.IonHigh[nelem] >= 0 );

	lgFirstPass = FALSE;
	if( nzone_used[nelem]!=nzone )
	{
		lgFirstPass = TRUE;
		nzone_used[nelem] = nzone;
	}

	/* H is special because its abundance spills into three routines -
	 * the ion/atom solver (this routine), the H-mole solvers (hmole), and
	 * the heavy molecule solver.  xmolecules only includes the heavy mole
	 * part for H only.  So the difference between gas_phase and xmolecules
	 * includes the H2 part of the molecular network.  This branch does
	 * this special H case, then the general case (heavy elements are
	 * not part of the H2 molecular network) */

	/* >>chng 01 dec 07, define abund_total, total atom and ion abundance here, removing molecules */
	if( nelem == ipHYDROGEN )
	{
		/* Hydrogen is a special case since hmole does not include the
		 * atom/molecules - hmole collapses its network into H = H0 + H+
		 * and then forces the solution determined here for partitioning
		 * beweeen these two */
		abund_total = dense.xIonDense[nelem][0] + dense.xIonDense[nelem][1];
	}
	else
	{
		abund_total = SDIV( dense.gas_phase[nelem] -  dense.xMolecules[nelem] );
	}

	/* protect against case where all gas phase abundances are in molecules, use the
	 * atomic and first ion density from the molecule solver 
	 * >>chng 04 aug 15, NA change from 10 0000 to 10 pre-coef on
	 * FLT_EPSILON for stability in pdr */
	if( fabs( dense.gas_phase[nelem] -  dense.xMolecules[nelem])/SDIV(dense.gas_phase[nelem] ) <
		10.*FLT_EPSILON )
	{
		/* >>chng 04 jul 31, add logic to conserve nuclei in fully molecular limit;
		 * in first calls, when searching for soln, we may be VERY far off, and sum of first ion
		 * and atom density may be far larger than difference between total gas and molecular densities,
		 * since they reflect the previous evaluation of the soln.  Do renorm to cover this case */
		/* first form sum of all atoms and ions */
		float sum = 0.;
		for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
			sum += dense.xIonDense[nelem][ion];
		/* now renorm to this sum - this should be unity, and is not if we have
		 * now conserved particles, due to molecular fraction changing */
		renorm = dense.gas_phase[nelem] / SDIV(sum + dense.xMolecules[nelem] );

		/*fprintf(ioQQQ,"DEBUG renorm %li\t%.3e\n",nelem , renorm );*/
		/*abund_total = renorm * (dense.xIonDense[nelem][0] + dense.xIonDense[nelem][1]);*/
		abund_total = renorm * sum;
	}
	if( abund_total < 0. )
	{
#		define EPS_MOLE	0.1
		/* this would mean that the molecular and ionic solvers are badly insane */
		if( fabs(abund_total) / SDIV(dense.xMolecules[nelem]) < EPS_MOLE )
		{
			/* only very slightly different, don't redo ionization, just div by 2 */
			for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
				dense.xIonDense[nelem][ion] = dense.xIonDense[nelem][ion]/2.f;
			/* print comment if not search phase */
			if(!conv.lgSearch )
				fprintf(ioQQQ,
					"PROBLEM neg net atomic abundance zero for nelem= %li, rel val= %.2e conv.nTotalIoniz=%li, fixed\n",
					nelem,
					fabs(abund_total) / SDIV(dense.xMolecules[nelem]),
					conv.nTotalIoniz );
			abund_total = -abund_total/2.;
		}
		else if( conv.lgSearch )
		{
			/* >>chng 04 sep 19, add this branch to stop abort in search phase,
			 * before molecules have come into equilibrium */
			/* this is case in search phase, this stuff can happen */
			/*fprintf(ioQQQ,
				"PROBLEM neg net atomic abundance zero for nelem= %li, rel val= %.2e conv.nTotalIoniz=%li, fixed\n",
				nelem,
				fabs(abund_total) / SDIV(dense.xMolecules[nelem]),
				conv.nTotalIoniz );*/
			abund_total = -abund_total/2.;
		}
		else
		{
			fprintf( ioQQQ, 
				"DISASTER ion_solver: negative net abundance found for nelem %li after correction for molecules.  fnzone %.2f\n",
				nelem,
				fnzone);
			fprintf( ioQQQ, "ion_solver: ratio |ion abund|/mole = %.2e\n",
				abund_total / SDIV(dense.xMolecules[nelem]) );
			fprintf( ioQQQ, "ion_solver: abund_total: %.2e, xMolecules: %.2e gas_phase: %.2e\n",
				abund_total , 
				dense.xMolecules[nelem],
				dense.gas_phase[nelem] );
			puts( "[Stop in ion_solver]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* return if IonHigh is zero, since no ionization at all */
	if( dense.IonHigh[nelem] == 0 )
	{
		/* set the atom to the total gas phase abundance */
		dense.xIonDense[nelem][0] = (float)abund_total ;
#		ifdef DEBUG_FUN
		fputs( " <->ion_solver()\n", debug_fp );
#		endif
		return;
	}

	/* >>chng 01 may 09, add option to force ionization distribution with element name ioniz */
	if( dense.lgSetIoniz[nelem] )
	{
		for( ion=0; ion<nelem+2; ++ion )
		{
			dense.xIonDense[nelem][ion] = dense.SetIoniz[nelem][ion]*(float)abund_total;
		}
#		ifdef DEBUG_FUN
		fputs( " <->ion_solver()\n", debug_fp );
#		endif
		return;
	}

	/* impossible for HIonFrac[nelem] to be zero if IonHigh(nelem)=nelem+1
	 * HIonFrac(nelem) is stripped to hydrogen */
	/* >>chng 01 oct 30, to assert */
	ASSERT( (dense.IonHigh[nelem] < nelem + 1) || iso.pop_ion_ov_neut[ipH_LIKE][nelem] > 0. );

	/* zero out the ionization and recombination rates that we will modify here,
	 * but not the iso-electronic stages which are done elsewhere,
	 * the nelem stage of ionization is he-like,
	 * the nelem+1 stage of ionization is h-like */

	/* loop over stages of ionization that we solve for here, 
	 * up through and including one less than nelem-NISO,
	 * never actually do highest NISO stages of ionization since they
	 * come from the ionization ratio from the next lower stage */
	limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);

	/* do full range of ionization - this is number of ionization stages */
	ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;

	ion_low = dense.IonLow[nelem];

	/* this will be "new" matrix, with non-adjacent coupling included */
	if( (xmat=(double*)MALLOC( (sizeof(double)*(unsigned)((ion_range)*(ion_range)) ))) == NULL )
		BadMalloc();
	if( (xmatsave=(double*)MALLOC( (sizeof(double)*(unsigned)((ion_range)*(ion_range)) ))) == NULL )
		BadMalloc();

	/* this will be used to address the 2d arrays */
#	ifdef MAT
#		undef MAT
#	endif
#	define MAT(M_,I_,J_)	(*((M_)+(I_)*(ion_range)+(J_)))

	/* zero-out loop comes before main loop since there are off-diagonal
	 * elements in the main ionization loop, due to multi-electron processes,
	 * TotIonizRate and TotRecom were already set in h-like and he-like solvers 
	 * other recombination rates were already set by routines responsible for them */
	for( ion=0; ion <= limit; ion++ )
	{
		ionbal.RateIonizTot[nelem][ion] = 0.;
		auger[ion] = 0.;
	}
	for( ion=MAX2(0,limit+1); ion <= dense.IonHigh[nelem]-1; ion++ )
	{
		auger[ion] = 0.;
	}

	/* zero out xmat */
	for( i=0; i< ion_range; ++i )
	{
		for( j=0; j< ion_range; ++j )
		{
			MAT( xmat, i, j ) = 0.;
		}
	}
#	define PRTNEW FALSE

	{
		/*@-redef@*/
		/* this sets up a fake ionization balance problem, with a trival solution,
		 * for debugging the ionization solver */
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && nelem==ipCARBON && PRTNEW )
		{
			broken();
			dense.IonLow[nelem] = 0;
			dense.IonHigh[nelem] = 3;
			abund_total = 1.;
			ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
			/* make up ionization and recombination rates */
			for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
			{
				double fac=1;
				if(ion)
					fac = 1e-10;
				ionbal.RateRecomTot[nelem][ion] = 100.;
				for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
				{
					/* direct photoionization of this shell */
					ionbal.PhotoRate_Shell[nelem][ion][ns][0] = fac;
				}
			}
		}
	}
	/* now get actual arrays of ionization and recombination processes,
	 * but only for the ions that are done as two-level systems */
	/* in two-stage system, atom + first ion, limit is zero but must
	 * include gv.GrainChTrRate[nelem][1][0] */
	/* grain charge transfer */
	if( gv.lgDustOn && ionbal.lgGrainIonRecom && gv.lgGrainPhysicsOn )
	{
		long int low;
		/* do not double count this process for atoms that are in the co network - we use
		 * a net recombination coef derived from the co solution, this includes grain ct */
		if( nelem==ipHYDROGEN ||nelem==ipCARBON ||nelem== ipOXYGEN ||nelem==ipSILICON ||nelem==ipSULPHUR 
			 ||nelem==ipNITROGEN )
		{
			 low = MAX2(1, dense.IonLow[nelem] );
		}
		else
			low = dense.IonLow[nelem];

		for( ion_to=low; ion_to <= limit; ion_to++ )
		{
			for( ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from )
			{
				/* do not do ion onto itself */
				if( ion_to != ion_from )
				{
					/* this is the rate coef for charge transfer from ion to ion_to */
					rateone = gv.GrainChTrRate[nelem][ion_from][ion_to];
					MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
					MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
				}
			}
		}
	}

	for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
	{
		/* thermal & secondary collisional ionization */
		rateone = ionbal.CollIonRate_Ground[nelem][ion][0] +
			secondaries.csupra[nelem][ion] +
			/* inner shell ionization by UTA lines */
			ionbal.xInnerShellIonize[nelem][ion];;
		ionbal.RateIonizTot[nelem][ion] += rateone;

		/* UTA ionization */
		if( ion+1-ion_low < ion_range )
		{
			/* depopulation processes enter with negative sign */
			MAT( xmat, ion-ion_low, ion-ion_low ) -= rateone;
			MAT( xmat, ion-ion_low, ion+1-ion_low ) += rateone;
		}

		/* total recombination rate */
		if( ion-1-ion_low >= 0 )
		{
			/* loss of this ion due to recom to next lower ion stage */
			MAT( xmat,ion-ion_low, ion-ion_low ) -= ionbal.RateRecomTot[nelem][ion-1];
			MAT( xmat,ion-ion_low, ion-1-ion_low ) += ionbal.RateRecomTot[nelem][ion-1];
		}

		/* >>chng 03 oct 20, add single electron ionization on grain surface */
		/* >>chng 04 apr 10, move to formal correct treatment, this term is total grain destruction
		 * rate out of current stage of ionization - evaluated in grain_ion_charge_transfer() -
		 * destroys ion, but can be to higher or lower charge */
		ionbal.RateIonizTot[nelem][ion] += ionbal.GrainDestr[nelem][ion];

		/* loop over all atomic sub-shells to include photoionization */
		for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
		{
			/* direct photoionization of this shell */
			ionbal.RateIonizTot[nelem][ion] += ionbal.PhotoRate_Shell[nelem][ion][ns][0];

			/* this is the primary ionization rate - add to diagonal element,
			 * test on ion stage is so that we don't include ionization from the very highest
			 * ionization stage to even higher - since those even higher stages are not considered
			 * this would appear as a sink - but populations of this highest level is ensured to
			 * be nearly trivial and neglecting it production of even higher ionization OK */
			if( ion+1-ion_low < ion_range )
			{
				/* this will be redistributed into charge states in following loop */ 
				MAT( xmat, ion-ion_low, ion-ion_low ) -= ionbal.PhotoRate_Shell[nelem][ion][ns][0];
			}

			/* yield.n_elec_eject[nelem][ion][ns] is total number of electrons that can
			 * possibly be freed 
			 * loop over nej, the number of electrons ejected including the primary,
			 * nej = 1 is primary, nej > 1 includes primary plus Auger 
			 * yield.frac_elec_eject is prob of nej electrons */
			for( nej=1; nej <= yield.n_elec_eject[nelem][ion][ns]; nej++ )
			{
				/* this is the ion that is produced by this ejection,
				 * limited by highest possible stage of ionization -
				 * do not want to ignore ionization that go beyond this */
				IonProduced = MIN2(ion+nej,dense.IonHigh[nelem]);
				rateone = ionbal.PhotoRate_Shell[nelem][ion][ns][0]*
					yield.frac_elec_eject[nelem][ion][ns][nej-1] ;
				/* >>chng 04 sep 06, above had included factor of nej to get rate, but
				 * actually want events into particular ion */
					/* number of electrons ejected
					*(double)nej ; */
				/* it goes into this charge state - recall upper cap due to ion stage trimming 
				 * note that compensating loss term on diagonal was done before this
				 * loop, since frac_elec_eject adds to unity */
				MAT( xmat, ion-ion_low, IonProduced-ion_low ) += rateone;
			}

			/* following loop is over number of electrons that come out of shell
			 * with multiple electron ejection,
			 * yield.n_elec_eject[nelem][ion][ns] is number of electrons, =1 for single electron photoionization
			 * this loop is only over multiple ejection */
			for( nej=2; nej <= yield.n_elec_eject[nelem][ion][ns]; nej++ )
			{
				double hold;
				/* this is the ion that is produced by this ejection,
				 * limited by highest possible stage of ionization -
				 * do not want to ignore ionization that go beyond this */
				IonProduced = MIN2(ion+nej,dense.IonHigh[nelem]);

				/* will treat this is an adjoining ionization process,
				 * from IonProduced-1 to IonProduced */
				if( IonProduced <=dense.IonHigh[nelem] &&
					dense.xIonDense[nelem][IonProduced-1] > 1e-30 )
				{
					if( nzone <3 || lgFirstPass )
					{
						/* ion points to current stage of ionization */
						ratio_ion_density = (double)(dense.xIonDense[nelem][ion])/
							(double)(dense.xIonDense[nelem][IonProduced-1]);
					}
					else
					{
						/* >>chng 04 sep 01, from 0.5 to 0.75 to stop oscillations */
						double frac_old = 0.5;
						/* ion points to current stage of ionization */
						ratio_ion_density = 
							frac_old * abund_old[nelem][ion]/abund_old[nelem][IonProduced-1] +
							(1.-frac_old) * (double)(dense.xIonDense[nelem][ion])/
							(double)(dense.xIonDense[nelem][IonProduced-1]);
					}
				}
				else
				{
					ratio_ion_density = 1.;
				}

				/* yield here is fraction removing ion electrons */
				hold = ratio_ion_density *
					ionbal.PhotoRate_Shell[nelem][ion][ns][0]*
					yield.frac_elec_eject[nelem][ion][ns][nej-1] ;
					/* >>chng 01 dec 18, number of electrons
					(double)nej ; */

				/* add this term to total Auger term */
				auger[IonProduced-1] += hold;

				/* also add to total ionization rate */
				ionbal.RateIonizTot[nelem][IonProduced-1] += hold;
			}
		}

		/* this is charge transfer ionization of this species by hydrogen and helium */
		rateone = 
			atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[ipHELIUM][1]+ 
			atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[ipHYDROGEN][1];
		ionbal.RateIonizTot[nelem][ion] += rateone;
		if( ion+1-ion_low < ion_range )
		{
			MAT( xmat, ion-ion_low, ion-ion_low ) -= rateone;
			MAT( xmat, ion-ion_low, ion+1-ion_low ) += rateone;
		}
	}
	/* after this loop, ionbal.RateIonizTot and ionbal.RateRecomTot have been defined for the
	 * stages of ionization that are done with simple */
	/* begin loop at first species that is treated with full model atom */
	j = MAX2(0,limit+1);
	/* possible that lowest stage of ionization is higher than this */
	j = MAX2( ion_low , j );
	for( ion=j; ion<=dense.IonHigh[nelem] ; ion++ )
	{
		ASSERT( ion>=0 && ion<nelem+2 );
		/* use total ionization/recombination rates */
		/* UTA ionization */
		if( ion+1-ion_low < ion_range )
		{
			/* depopulation processes enter with negative sign */
			MAT( xmat, ion-ion_low, ion-ion_low ) -= ionbal.RateIonizTot[nelem][ion];
			MAT( xmat, ion-ion_low, ion+1-ion_low ) += ionbal.RateIonizTot[nelem][ion];
		}

		if( ion-1-ion_low >= 0 )
		{
			/* loss of this ion due to recom to next lower ion stage */
			MAT( xmat,ion-ion_low, ion-ion_low ) -= ionbal.RateRecomTot[nelem][ion-1];
			MAT( xmat,ion-ion_low, ion-1-ion_low ) += ionbal.RateRecomTot[nelem][ion-1];
		}
	}

	/* >> chng 03 jan 13 rjrw, always solve using full matrix solution to allow for sources
	 * from molecular netork */
	/* last test - only do advection if we have not overrun the radius scale */
#define RJRW 1

	/* MALLOC space for the  1-d array */
	if( (amat=(double*)MALLOC( (sizeof(double)*(unsigned)((ion_range)*(ion_range)) ))) == NULL )
		BadMalloc();
	if( (achk=(double*)MALLOC( (sizeof(double)*(unsigned)((ion_range)*(ion_range)) ))) == NULL )
		BadMalloc();
	if( (src=(double*)MALLOC( (sizeof(double)*(unsigned)((ion_range)) ))) == NULL )
		BadMalloc();
	if( (snk=(double*)MALLOC( (sizeof(double)*(unsigned)((ion_range)) ))) == NULL )
		BadMalloc();
	if( (source=(double*)MALLOC( (sizeof(double)*(unsigned)((ion_range)) ))) == NULL )
		BadMalloc();
	if( (ipiv=(int32*)MALLOC( (sizeof(int32)*(unsigned)((ion_range)) ))) == NULL )
		BadMalloc();

	for( i=0; i<ion_range;i++ )
	{
		for( j=0;j<ion_range;j++ )
		{
			MAT(amat,i,j) = 0.;
		}
		src[i] = 0.;
		source[i] = 0.;
	}

	/* Set lgTriDiag = FALSE if non-tridiagonal terms are ever added to amat */
	lgTriDiag = FALSE; 

	for( i=0; i<ion_range;i++ )
	{
		/* this will be sum of source and sink terms, will be used to decide if
		 * matrix is singular */
		double totsrc = 0., totsnk = 0.;
		ion = i+ion_low;

		/* these are the external source and sink terms */
		/* source first */
		src[i] += mole.source[nelem][ion];
		/* need negative sign to get positive pops */
		source[i] -= mole.source[nelem][ion];
		totsrc += mole.source[nelem][ion];
		/* sink next */
		MAT(amat,i,i) += mole.sink[nelem][ion];
		MAT( xmat, i, i ) -= mole.sink[nelem][ion];
		totsnk += mole.sink[nelem][ion];
		/* matrix is not homogeneous if source is non-zero */
		if( totsrc != 0. /*&& totsnk != 0.*/)
			lgHomogeneous = FALSE;

		/* fprintf(ioQQQ," %li %li %.3e\n",i,i,MAT(amat,i,i)); */
		/* grain ct rec already included in RateRecomTot for these cases, for rec to atom */
		if( i==0 &&
			(nelem==ipHYDROGEN ||nelem==ipCARBON ||nelem== ipOXYGEN ||nelem==ipSILICON ||nelem==ipSULPHUR 
			 ||nelem==ipNITROGEN ) )
		{
			if( i != ion_range-1 )
			{
				MAT(amat,i,i) += ionbal.RateIonizTot[nelem][ion];
				MAT(amat,i+1,i) = -(ionbal.RateRecomTot[nelem][ion]);
				/* fprintf(ioQQQ," %li %li %.3e\n",i,i+1,MAT(amat,i,i+1)); */
			}
		}
		else
		{
			if( i != ion_range-1 )
			{
				MAT(amat,i,i) += ionbal.RateIonizTot[nelem][ion];
				MAT(amat,i+1,i) = -(ionbal.RateRecomTot[nelem][ion]+ionbal.GrainCreat[nelem][ion]);
				/* fprintf(ioQQQ," %li %li %.3e\n",i,i+1,MAT(amat,i,i+1)); */
			}
		}
		if( i==1 &&
			(nelem==ipHYDROGEN ||nelem==ipCARBON ||nelem== ipOXYGEN ||nelem==ipSILICON ||nelem==ipSULPHUR 
			 ||nelem==ipNITROGEN ) )
		{
			if( i != 0 )
			{
				MAT(amat,i,i) += (ionbal.RateRecomTot[nelem][ion-1]);
				MAT(amat,i-1,i) = -ionbal.RateIonizTot[nelem][ion-1];
				/* fprintf(ioQQQ," %li %li %.3e\n",i,i-1,MAT(amat,i,i-1)); */
			}
		}
		else
		{
			if( i != 0 )
			{
				MAT(amat,i,i) += (ionbal.RateRecomTot[nelem][ion-1]+ionbal.GrainCreat[nelem][ion-1]);
				MAT(amat,i-1,i) = -ionbal.RateIonizTot[nelem][ion-1];
				/* fprintf(ioQQQ," %li %li %.3e\n",i,i-1,MAT(amat,i,i-1)); */
			}
		}

	}

#	if 0
	/* >>chng 04 sep 03, integrate following into above with new source sink struc */
	if( nelem == ipHYDROGEN && !hmi.lgNoH2Mole)
	{
		double totsrc = 0., totsnk = 0.;
		/* for H ion_low, the lowest ionization stage, has to be zero,
		 * ion_range, the range in ionization, is 0 to 1 plus 1, so must be 2 */
		ASSERT( ion_low==0 && ion_range==2 );

		for (i=0;i<ion_range;i++) 
		{
			ion = i+ion_low;
			src[i] += mole.source[ipHYDROGEN][ion];
			source[i] += mole.source[ipHYDROGEN][ion];
			totsrc += mole.source[ipHYDROGEN][ion];
			MAT(amat,i,i) += mole.sink[ipHYDROGEN][ion];
			totsnk += mole.sink[ipHYDROGEN][ion];
		}
		if( totsrc != 0. /*&& totsnk != 0.*/)
			lgHomogeneous = FALSE;
	}
#	endif

	/* chng 03 jan 13 rjrw, add in dynamics if required here,
	 * last test - only do advection if we have not overrun the radius scale */
	if( iteration >= 2 && dynamics.lgAdvection && radius.depth < dynamics.oldFullDepth)
	{
		for( i=0; i<ion_range;i++ )
		{			 
			ion = i+ion_low;
			MAT(amat,i,i) += dynamics.Rate;
			MAT( xmat, i, i ) -= dynamics.Rate;
			src[i] += dynamics.Source[nelem][ion];
			source[i] -= dynamics.Source[nelem][ion];
			/* fprintf(ioQQQ," %li %li %.3e (%.3e %.3e)\n",i,i,MAT(amat,i,i),dynamics.Rate, dynamics.Source[nelem][ion]);*/
		}
		lgHomogeneous = FALSE;
	}
	for( i=0; i< ion_range; ++i )
	{
		for( j=0; j< ion_range; ++j )
		{
			MAT( xmatsave, i, j ) = MAT( xmat, i, j );
		}
	}


	if ( PRTNEW && nelem==13 )
	{
		fprintf(ioQQQ,"DEBUG old, matrix %li\n", nelem );
		for( i=0; i<ion_range;i++ )
		{
			for( j=0;j<ion_range;j++ )
			{
				fprintf(ioQQQ,"%e\t",MAT(amat,j,i));
				MAT(achk,j,i) = MAT(amat,j,i);
			}
			fprintf(ioQQQ,"\n");
			for( j=0;j<ion_range;j++ )
			{
				fprintf(ioQQQ,"%e\t",MAT(xmat,j,i));
			}
			fprintf(ioQQQ,"\n");
			for( j=0;j<ion_range;j++ )
			{
				fprintf(ioQQQ,"%e\t",MAT(amat,j,i)/SDIV(fabs(MAT(xmat,j,i))));
			}
			fprintf(ioQQQ,"\n\n");
		}
		fprintf(ioQQQ," source old new \n");
		for( i=0; i<ion_range;i++ )
		{
			fprintf(ioQQQ,"%e\t",src[i]);
		}
		fprintf(ioQQQ,"\n");
		for( i=0; i<ion_range;i++ )
		{
			fprintf(ioQQQ,"%e\t",source[i]);
		}
		fprintf(ioQQQ,"\n");
	}

	/* this is true if no source terms */
	if( lgHomogeneous )
	{
		double ionprod=1., recomb, scale = 0., value;
		/* Simple estimate of most abundant ion */
		jmax = 0;
		for( i=0; i<ion_range-1;i++)
		{ 
			ion = i+ion_low;
			ionprod *= ionbal.RateIonizTot[nelem][ion];
			recomb = ionbal.RateRecomTot[nelem][ion];
			/* trips if ion rate zero, so ll the gas will be more neutral than this */
			if( ionprod == 0)
				break;
			/* rec rate is zero */
			if (recomb <= 0.) 
				break;

			ionprod /= recomb;
			if (ionprod > 1.) 
			{
				/* this is peak ionization stage */
				jmax = i;
				ionprod = 1.;
			}
		}
		
		/* Matrix will be singular, so replace most abundant ion row with sum rule,
		 * scaled to maintain matrix condition. */
		scale = 0.;
		for (i=0; i<ion_range;i++)
		{
			value = MAT(amat,i,i);
			if (value > scale)
			{
				scale = value;
			}
		}
		scale *= 1e-6;
		
		for( i=0; i<ion_range;i++ )
		{
			MAT(amat,i,jmax) = 0.;
			MAT(xmat,i,jmax) = 1.;
		}
		MAT(amat,jmax,jmax) = scale;
		src[jmax] = scale;
		source[jmax] = abund_total;
	}

	if ( FALSE && nelem == ipHYDROGEN && dynamics.lgAdvection&& iteration>1 ) 
	{
		fprintf(ioQQQ,"DEBUGG Rate %.2f %.3e \n",fnzone,dynamics.Rate);
		fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateIonizTot[nelem][0], ionbal.RateIonizTot[nelem][1]);
		fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateRecomTot[nelem][0], ionbal.RateRecomTot[nelem][1]);
		fprintf(ioQQQ," %.3e %.3e %.3e\n\n", dynamics.Source[nelem][0], dynamics.Source[nelem][1], dynamics.Source[nelem][2]);
	}

	/* first branch is tridiag, following branch is just general matrix solver,
	 * this first branch is only one ever used */
	if (lgTriDiag) 
	{
		int lgLapack=FALSE , lgTriagSolver=FALSE  , lgTriOptimized=TRUE;
		/* there are three choices for the tridiag solver */
		if(lgLapack) 
		{
			/* this branch uses lapack tridiag solver, and should work 
			 * it is hardwired to not be used because not extensively tested
			 * issues - is this slower than others, and is it robust in
			 * sinrular cases? */
			int lgBad = FALSE;
			/* Use Lapack tridiagonal solver */
			double *dl, *d, *du;
			
			d = (double *) MALLOC((unsigned)ion_range*sizeof(double));
			du = (double *) MALLOC((unsigned)(ion_range-1)*sizeof(double));
			dl = (double *) MALLOC((unsigned)(ion_range-1)*sizeof(double));

			for (i=0;i<ion_range-1;i++) 
			{
				du[i] = MAT(amat,i+1,i);
				d[i] = MAT(amat,i,i);
				dl[i] = MAT(amat,i,i+1);
			}
			d[i] = MAT(amat,i,i);

			if( lgBad )
				fprintf(ioQQQ," dgtsz error.\n");
			
			free(dl);free(du);free(d);
		} 
		else if (lgTriagSolver) 
		{

			/* Use tridiagonal solver - default is not used */
			tridiag(amat,src,ion_range);

		}
		else if (lgTriOptimized)
		{
			double recrat[LIMELM+1];
			/* this is default
			 * Use tridiagonal solver re-coded to avoid rounding errors
			 * on diagonal -- uses determination of jmax for the
			 * singular case, but is otherwise independent of the array
			 * filling code above */
			
			for (i=0;i<ion_range;i++) 
			{
				src[i] = snk[i] = 0.;
			}
			if (nelem == ipHYDROGEN && !hmi.lgNoH2Mole)
			{
				for (i=0;i<ion_range;i++) 
				{
					ion = i+ion_low;
					src[i] += mole.source[ipHYDROGEN][ion];
					snk[i] += mole.sink[ipHYDROGEN][ion];
				}
			}

			/* add advection terms if we have previous solution and wind enabled */
			if ( iteration >= 2 && dynamics.lgAdvection && radius.depth < dynamics.oldFullDepth )
			{
				for (i=0;i<ion_range;i++) 
				{
					ion = i+ion_low;
					src[i] += dynamics.Source[nelem][ion];
					snk[i] += dynamics.Rate;
				}
			}
			for (i=0;i<ion_range-1;i++) 
			{
				ion = i+ion_low;
				/* need to include grain charge transfer here */
				recrat[i] = ionbal.RateRecomTot[nelem][ion] + ionbal.GrainCreat[nelem][ion];
			}
			
			solveions(ionbal.RateIonizTot[nelem]+ion_low,/*ionbal.RateRecomTot[nelem]+ion_low*/recrat ,
				snk,src,ion_range,
				/* this was init to -1 and set to >=0 if singular case encountered */
				jmax);
		}
	} 
	else 
	{
		nerror = 0;
		/* Use general matrix solver */
		getrf_wrapper(ion_range, ion_range, amat, ion_range, ipiv, &nerror);
		if( nerror != 0 )
		{
			fprintf( ioQQQ, " ion_solver: dgetrf finds singular or ill-conditioned matrix\n" );
			puts( "[Stop in ion_solver]" );
			cdEXIT(EXIT_FAILURE);
		}
		getrs_wrapper('N', ion_range, 1, amat, ion_range, ipiv, src, ion_range, &nerror);
		if( nerror != 0 )
		{
			fprintf( ioQQQ, " ion_solver: dgetrs finds singular or ill-conditioned matrix\n" );
			puts( "[Stop in ion_solver]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* get new soln */
	nerror = 0;
	/* Use general matrix solver */
	getrf_wrapper(ion_range, ion_range, xmat, ion_range, ipiv, &nerror);
	if( nerror != 0 )
	{
		fprintf( ioQQQ, 
			" ion_solver: dgetrf finds singular or ill-conditioned matrix nelem=%li ion_range=%li, limit=%li, xmat follows\n",
			nelem , ion_range,limit  );
		for( i=0; i<ion_range; ++i )
		{
			for( j=0;j<ion_range;j++ )
			{
				fprintf(ioQQQ,"%e\t",MAT(xmat,j,i));
			}
			fprintf(ioQQQ,"\n");
		}
		fprintf(ioQQQ,"source follows\n");
		for( i=0; i<ion_range;i++ )
		{
			fprintf(ioQQQ,"%e\t",source[i]);
		}
		fprintf(ioQQQ,"\n");
		puts( "[Stop in ion_solver]" );
		cdEXIT(EXIT_FAILURE);
	}
	getrs_wrapper('N', ion_range, 1, xmat, ion_range, ipiv, source, ion_range, &nerror);
	if( nerror != 0 )
	{
		fprintf( ioQQQ, " ion_solver: dgetrs finds singular or ill-conditioned matrix nelem=%li ionrange=%li\n",
			nelem , ion_range );
		puts( "[Stop in ion_solver]" );
		cdEXIT(EXIT_FAILURE);
	}

	{
		/*@-redef@*/
		/* this is to debug following failed assert */
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && (nzone >380 ) && nelem == ipHYDROGEN )
		{
			fprintf(ioQQQ,"debuggg\t%.2f\t%.4e\t%.4e\tIon\t%.3e\tRec\t%.3e\n", 
				fnzone,
				phycon.te,
				dense.eden,
				ionbal.RateIonizTot[nelem][0] , 
				ionbal.RateRecomTot[nelem][0]);
			fprintf(ioQQQ," Msrc %.3e %.3e\n", mole.source[ipHYDROGEN][0], mole.source[ipHYDROGEN][1]);
			fprintf(ioQQQ," Msnk %.3e %.3e\n", mole.sink[ipHYDROGEN][0], mole.sink[ipHYDROGEN][1]);
			fprintf(ioQQQ," Poprat %.3e nomol %.3e\n",src[1]/src[0],
				ionbal.RateIonizTot[nelem][0]/ionbal.RateRecomTot[nelem][0]);
		}
	}

	if (RJRW && 0)
	{
		/* verify that the rates are sensible */
		double test;
		for (i=0; i<ion_range; i++) {
			test = 0.;
			for (j=0; j<ion_range; j++) {
				test = test+src[j]*MAT(achk,j,i);
			}
			fprintf(ioQQQ,"%e\t",test);
		}
		fprintf(ioQQQ,"\n");

		test = 0.;
		fprintf(ioQQQ," ion %li abundance %.3e\n",nelem,abund_total);
		for( ion=dense.IonLow[nelem]; ion < dense.IonHigh[nelem]; ion++ )
		{
			if ( ionbal.RateRecomTot[nelem][ion] != 0 && src[ion-ion_low] != 0 )
				fprintf(ioQQQ," %li %.3e %.3e : %.3e\n",ion,src[ion-ion_low],
								src[ion-ion_low+1]/src[ion-ion_low],
								ionbal.RateIonizTot[nelem][ion]/ionbal.RateRecomTot[nelem][ion]);
			else
				fprintf(ioQQQ," %li %.3e [One ratio infinity]\n",ion,src[ion-ion_low]);
			test += src[ion-ion_low];
		}
		test += src[ion-ion_low];
	}

	/* 
	 * >> chng 03 jan 15 rjrw:- terms are now included for
	 * molecular sources and sinks of H and H+.
	 *
	 * When the network is not in equilibrium, this will lead to a
	 * change in the derived abundance of H and H+ when after the
	 * matrix solution -- the difference between `renorm' and 1. is a
	 * measure of the quality of the solution (it will be 1. if the
	 * rate of transfer into Ho/H+ balances the rate of transfer
	 * out, for the consistent relative abundances).
	 *
	 * We therefore renormalize to keep the total H abundance
	 * correct -- only the molecular network is allowed to change
	 * this.
	 *
	 * To do this, only the ion abundances are corrected, as the
	 * molecular abundances may depend on several different
	 * conserved species.
	 *
	 */
	if( lgHomogeneous )
	{
		dense.xIonDense[nelem][ion_low] = (float)abund_total;
		for ( i=1;i < ion_range; i++ )
		{
			dense.xIonDense[nelem][i+ion_low] = 0.;
		}
	}

	renorm = 1.;
	renormnew = 1.;
	if (iteration >= 2 && dynamics.lgAdvection && radius.depth < dynamics.oldFullDepth && 
			nelem == ipHYDROGEN && hmi.lgNoH2Mole)
	{
		/* The normalization out of the matrix solution is correct and
		 * should be retained if: dynamics is on and the total
		 * abundance of HI & H+ isn't being controlled by the
		 * molecular network */
		renorm = 1.;
		renormnew = 1.;
	}
	else
	{
		double dennew = 0.;
		sum_dense = 0.;
		den = 0.;
		
		for( i=0;i < ion_range; i++ )
		{
			ion = i+ion_low;
			sum_dense += dense.xIonDense[nelem][ion];
			den += src[i];
			dennew += source[i];
		} 
		
		if( den > 0.)
		{
			renorm = sum_dense / den;
			renormnew = sum_dense / dennew;
			/*TODO	2	renorm should == 1 when the molecules and
			 * ionization are in equilibrium.  Should monitor
			 * this figure of merit in calling routine.
			 * */
			if (0 && !lgHomogeneous)
				fprintf(ioQQQ,"Zone %li elem %ld normalization error %g\n",nzone,nelem,renorm-1.);
		}
		else
		{
			renorm = 1.;
			renormnew = 1.;
		}
	}
	if( renorm <= 0 )
	{
		fprintf(ioQQQ,"impossible value of renorm \n");
	}
	ASSERT( renorm>0 && renormnew>0 );

	/* save resulting abundances into main ionization density array, 
	 * while checking whether any negative level populations occured */
	lgNegPop = FALSE;
	{
		double diffmax; long int ionmax;
	diffmax = 0.;
	ionmax = -1;
	for( i=0; i < ion_range; i++ )
	{
		ion = i+ion_low;

		/*fprintf(ioQQQ," %li %li %.3e %.3e\n",nelem,ion,src[ion-ion_low+1],src[ion-ion_low]);
		pop_ion_ov_neut[ion] = src[ion-ion_low+1]/src[ion-ion_low];
		if( lgHomogeneous )
			renormnew = 1.;*/
		/* use old solution */
		dense.xIonDense[nelem][ion] = (float)(src[i]*renorm);
		/* use new solution */
		dense.xIonDense[nelem][ion] = (float)(source[i]*renormnew);
		if( fabs( source[i]*renormnew - src[i]*renorm ) /abund_total > diffmax )
		{
			diffmax = fabs( source[i]*renormnew - src[i]*renorm ) /abund_total;
			ionmax = ion;
		}
		/* fprintf(ioQQQ,"%ld %g [%g]\n",ion,dense.xIonDense[nelem][ion],src[i]*renorm); */
		if( dense.xIonDense[nelem][ion]< 0. )
		{
			/* early solns may have too broad an ionization distribution and so slightly
			 * negative charge states towards the edge.  set to zero during first solns */
			if( dense.xIonDense[nelem][ion]>-1e-10 )
			{
				dense.xIonDense[nelem][ion] = 0.;
			}
			else
			{
				fprintf(ioQQQ,"PROBLEM new ion neg\t%li\t%li\t%li\t%e\n",
					conv.nTotalIoniz,nelem , ion, source[i] );
				lgNegPop = TRUE;
			}
		}
	}
	/* only print this problem if we have done this enough to be close to soln,
	 * original solver used previous abundances to get new ones */
	/* >>chng 04 sep 19, require non-zero zone number to silence false alarms */
	if( nzone && !conv.lgSearch && diffmax > 0.1 )
		fprintf(ioQQQ,"PROBLEM big diff\t%li\t%.2f\t%li\t%li\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n",
			conv.nTotalIoniz,
			fnzone,nelem , ionmax, diffmax, renorm, renormnew ,sum_dense ,abund_total,
			src[0]*renorm , src[1]*renorm , 
			source[0]*renormnew , source[1]*renormnew);
	}

	if( PRTNEW && nelem==13 )
	{
		fprintf(ioQQQ,"DEBUG old, new, ratio pops,");
		if( lgHomogeneous )
			fprintf(ioQQQ," Homogeneous nelem=%li jmax=%li", nelem,jmax);
		fprintf(ioQQQ,"\n");
		for( i=0; i < ion_range; i++ )
		{
			ion = i+ion_low;
			fprintf(ioQQQ,"\t%e", src[i]*renorm );
		}
		fprintf(ioQQQ,"\n");
		for( i=0; i < ion_range; i++ )
		{
			ion = i+ion_low;
			fprintf(ioQQQ,"\t%e", source[i]*renormnew );
		}
		fprintf(ioQQQ,"\n");
		for( i=0; i < ion_range; i++ )
		{
			double fsource;
			ion = i+ion_low;
			fsource = fabs(source[i]*renorm);
			fprintf(ioQQQ,"\t%e", src[i]*renorm/SDIV(fsource) );
		}
		fprintf(ioQQQ,"\n\n");
		/*fprintf(ioQQQ,"DEBUG matrix nelem=%li \n%e\t%e\t%e\n%e\t%e\t%e\n%e\t%e\t%e\t%e\n",
			nelem,
			MAT(xmatsave,0,0),
			MAT(xmatsave,0,1),
			MAT(xmatsave,0,3),
			MAT(xmatsave,1,0),
			MAT(xmatsave,1,1),
			MAT(xmatsave,1,2),
			MAT(xmatsave,2,0),
			MAT(xmatsave,2,1),
			MAT(xmatsave,2,2),
			((source[0]*MAT(xmatsave,0,1) + source[2]*MAT(xmatsave,2,1)) /
			MAT(xmatsave,1,1) ));*/
		/*if( nelem == ipCARBON )*/
		{static int ncount = 0;++ncount;
			fprintf(ioQQQ,"ion_solver exitttt\n");
			if( ncount>10 )exit(99);
		}
	}
	
	/*fprintf(ioQQQ," Sums %.3e %.3e %.3e\n",sum1,sum2,dynamics.Rate);*/
	free(ipiv);
	free(src);
	free(source);
	free(achk);
	free(snk);
	free(amat);
	free(xmat);
	free(xmatsave);

	/* Zero levels with abundances < 1e-25 which which will suffer numerical noise */
	while( dense.IonHigh[nelem] > dense.IonLow[nelem] && 
					dense.xIonDense[nelem][dense.IonHigh[nelem]] < 1e-25*abund_total )
	{
		ASSERT( dense.xIonDense[nelem][dense.IonHigh[nelem]] >= 0. );
		/* zero out abundance and heating due to stage of ionization we are about to zero out */
		/* fprintf(ioQQQ,"Zone %ld: elem %ld losing %ld (%g %g)\n",nzone,nelem,dense.IonHigh[nelem],
		   dense.xIonDense[nelem][dense.IonHigh[nelem]],abund_total); */
		dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
		thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.;
		/* decrement counter */
		--dense.IonHigh[nelem];
	}

	/* sanity check, either offset stages of low and high ionization,
	 * or no ionization at all */
	ASSERT( (dense.IonLow[nelem] < dense.IonHigh[nelem]) ||
		(dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) );

	/* save ionization balance for this element */
	for( ion=0; ion <dense.IonHigh[nelem]; ion++ )
	{
		abund_old[nelem][ion] = dense.xIonDense[nelem][ion];
	}

	/* this can`t possibly happen */
	if( lgNegPop )
	{
		fprintf( ioQQQ, " Negative population found for abundance of ionization stage of element %4.4s, ZONE=%4ld\n", 
		  elementnames.chElementNameShort[nelem], nzone );

		fprintf( ioQQQ, " Populations were" );
		for( ion=1; ion <= dense.IonHigh[nelem]+1; ion++ )
		{
			fprintf( ioQQQ, "%9.1e", dense.xIonDense[nelem][ion-1] );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, " destroy vector =" );
		for( ion=1; ion <= dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, "%9.1e", ionbal.RateIonizTot[nelem][ion-1] );
		}
		fprintf( ioQQQ, "\n" );

		/* print some extra stuff if destroy was negative */
		if( lgNegPop )
		{
			fprintf( ioQQQ, " CTHeavy  vector =" );
			for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
			{
				fprintf( ioQQQ, "%9.1e", atmdat.HeCharExcIonOf[nelem][ion] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, " HCharExcIonOf vtr=" );
			for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
			{
				fprintf( ioQQQ, "%9.1e", atmdat.HCharExcIonOf[nelem][ion] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, " CollidRate  vtr=" );
			for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
			{
				fprintf( ioQQQ, "%9.1e", ionbal.CollIonRate_Ground[nelem][ion][0] );
			}
			fprintf( ioQQQ, "\n" );

			/* photo rates per subshell */
			fprintf( ioQQQ, " photo rates per subshell, ion\n" );
			for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
			{
				fprintf( ioQQQ, "%3ld", ion );
				for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
				{
					fprintf( ioQQQ, "%9.1e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
				}
				fprintf( ioQQQ, "\n" );
			}
		}

		/* now check out creation vector */
		fprintf( ioQQQ, " create  vector =" );
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, "%9.1e", ionbal.RateRecomTot[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		ContNegative();
		ShowMe();
		puts( "[Stop in bidiag]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* option to print ionization and recombination arrays
	 * prt flag set with "print array" command */
	if( prt.lgPrtArry || lgPrintIt )
	{
		/* say who we are, what we are doing .... */
		fprintf( ioQQQ, 
			"\n %s ion_solver DEBUG ion/rec rt [s-1] %s nz%.2f Te%.4e ne%.4e Tot abun:%.3e ion abun%.2e mole%.2e\n", 
			elementnames.chElementSym[nelem],
			elementnames.chElementName[nelem],
			fnzone,
			phycon.te , 
			dense.eden,
			dense.gas_phase[nelem],
			abund_total ,
			dense.xMolecules[nelem] );
		/* total ionization rate, all processes */
		fprintf( ioQQQ, " %s Ioniz total " ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", ionbal.RateIonizTot[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );
		if( nelem==ipHYDROGEN )
			fprintf(ioQQQ," H  mole sink:   %9.2e %9.2e\n",
			mole.sink[ipHYDROGEN][0], mole.sink[ipHYDROGEN][1]);


		/* sum of all creation processes */
		fprintf( ioQQQ, " %s Recom total " ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", ionbal.RateRecomTot[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );
		if( nelem==ipHYDROGEN )
			fprintf(ioQQQ," H  mole source: %9.2e %9.2e\n",
			mole.source[ipHYDROGEN][0] , mole.source[ipHYDROGEN][1]);

		/* collisional ionization */
		fprintf( ioQQQ, " %s Coll ioniz  " ,elementnames.chElementSym[nelem] );
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", ionbal.CollIonRate_Ground[nelem][ion][0] );
		}
		fprintf( ioQQQ, "\n" );

		/* UTA ionization */
		fprintf( ioQQQ, " %s UTA ioniz   " ,elementnames.chElementSym[nelem] );
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", ionbal.xInnerShellIonize[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* photo ionization */
		fprintf( ioQQQ, " %s Phot ioniz  " ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", 
				ionbal.PhotoRate_Shell[nelem][ion][Heavy.nsShells[nelem][ion]-1][0] );
		}
		fprintf( ioQQQ, "\n" );

		/* auger ionization */
		fprintf( ioQQQ, " %s Auger ioniz " ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", 
				auger[ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* secondary ionization */
		fprintf( ioQQQ, " %s Secon ioniz " ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", 
				secondaries.csupra[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* grain ionization */
		fprintf( ioQQQ, " %s dest on grn "  ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", ionbal.GrainDestr[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* charge exchange ionization */
		fprintf( ioQQQ, " %s chr trn ion " ,elementnames.chElementSym[nelem] );
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			double sum = atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[ipHELIUM][1]+ 
				atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[ipHYDROGEN][1];

			if( nelem==ipHELIUM && ion==0 )
			{
				sum += atmdat.HeCharExcIonTotal;
			}
			else if( nelem==ipHYDROGEN && ion==0 )
			{
				sum += atmdat.HCharExcIonTotal;
			}
			fprintf( ioQQQ, " %9.2e", sum );
		}
		fprintf( ioQQQ, "\n" );

		/* radiative recombination */
		fprintf( ioQQQ, " %s radiati rec "  ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", dense.eden*ionbal.RadRecomRateCoef[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* grain recombination */
		fprintf( ioQQQ, " %s form on grn "  ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", ionbal.GrainCreat[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );

		/* charge exchange recombination */
		fprintf( ioQQQ, " %s chr trn rec "  ,elementnames.chElementSym[nelem]);
		for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
		{
			double sum = 
				atmdat.HCharExcRecTo[nelem][ion]*
				iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]*
				dense.xIonDense[ipHYDROGEN][1] +

				atmdat.HeCharExcRecTo[nelem][ion]*
				iso.Pop2Ion[ipHE_LIKE][ipHELIUM][ipHe1s1S]*
				dense.xIonDense[ipHELIUM][1];

			if( nelem==ipHELIUM && ion==0 )
			{
				sum += atmdat.HeCharExcRecTotal;
			}
			else if( nelem==ipHYDROGEN && ion==0 )
			{
				sum += atmdat.HCharExcRecTotal;
			}
			fprintf( ioQQQ, " %9.2e", sum );
		}
		fprintf( ioQQQ, "\n" );

		/* the "new" abundances the resulted from the previous ratio */
		fprintf( ioQQQ, " %s Abun [cm-3] " ,elementnames.chElementSym[nelem] );
		for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
		{
			fprintf( ioQQQ, " %9.2e", dense.xIonDense[nelem][ion] );
		}
		fprintf( ioQQQ, "\n" );
	}

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


/* 

	 Solve an ionization level system with specified ionization and
	 recombination rates between neighbouring ions, and additional sink
	 and source terms.  The sink array is overwritten, and the results
	 appear in the source array.

	 Written in matrix form, the algorithm is equivalent to the
	 tridiagonal algorithm in Numerical Recipes applied to:

	 / i_0+a_0     -r_0          .           .    .  \ / x_0 \   / s_0 \
	 |  -i_0    i_1+a_1+r_0    -r_1          .    .  | | x_1 |   | s_1 |
	 |    .        -i_1      i_2+a_2+r_1   -r_2   .  | | x_2 |   | s_2 |
     |    .          .       (etc....)               | | ... | = | ... |
     \    .          .          .                    / \     /   \     /

	 where i, r are the ionization and recombination rates, s is the
	 source rate and a is the sink rate.

	 This matrix is diagonally dominant only when the sink terms are
	 large -- the alternative method coded here prevents rounding error
	 in the diagonal terms disturbing the solution when this is not the
	 case.

*/

/* solveions tridiagonal solver but optimized for structure of balance matrix */
void solveions(double *ion, double *rec, double *snk, double *src,
	       long int nlev, long int nmax)
{
  double kap, bet;
  long int i;

  if( nmax != -1) 
  {
    /* Singular case */
    src[nmax] = 1.;
    for (i=nmax;i<nlev-1;i++)
      src[i+1] = src[i]*ion[i]/rec[i];
    for (i=nmax-1;i>=0;i--)
      src[i] = src[i+1]*rec[i]/ion[i];
  } 
  else 
  {
	  /* this is usual, non-singular case */
    i = 0;
    kap = snk[0];    
    for (i=0;i<nlev-1;i++) 
    {
      bet = ion[i]+kap;
      if( bet == 0. )
      {
			fprintf(ioQQQ,"Ionization solver error\n");
			puts("[Stop in solveions]");
			cdEXIT(EXIT_FAILURE);
      }
      bet = 1./bet;
      src[i] *= bet;
      src[i+1] += ion[i]*src[i];
      snk[i] = bet*rec[i];
      kap = kap*snk[i]+snk[i+1];
    }
	bet = kap;
	if (bet == 0.)
	{
		fprintf(ioQQQ,"Ionization solver error\n");
		puts("[Stop in solveions]");
		cdEXIT(EXIT_FAILURE);
	}
	src[i] /= bet;
    
    for (i=nlev-2;i>=0;i--)
    {
      src[i] += snk[i]*src[i+1];
    }
  }
}

#ifdef MAT
#	undef MAT
#	endif
#define MAT(M_,I_,J_)	(*((M_)+(I_)*(n)+(J_)))
/*tridiag - general tridiag solver - */
void tridiag(double *a, double *b, long int n)
{
	long j;
	double c,*g,t;
	
	g = (double *) MALLOC((unsigned)n*sizeof(double));
	if (MAT(a,0,0) == 0.) 
	{
		fprintf(ioQQQ,"Error 1 in tridiag\n");
		exit(-1);
	}
	c = MAT(a,0,0);
	b[0] /= c;
	for (j=1; j<n; j++) 
	{
		g[j] = MAT(a,j,j-1)/c;
		t = MAT(a,j-1,j);
		c = MAT(a,j,j)-t*g[j];
		if (c == 0.) {
			fprintf(ioQQQ,"Error 2 in tridiag\n");
			exit(-1);
		}
		b[j] = (b[j]-t*b[j-1])/c;
	}
	for (j=n-1;j>0;j--)
		b[j-1] -= g[j]*b[j];
	free(g);
}

