/*BiDiag solve the bi-diagonal matrix for ionization balance */
#include <string.h>
#include "cddefines.h"
#include "yield.h"
#include "parray.h"
#include "ionrange.h"
#include "recomrate.h"
#include "photrate.h"
#include "ionfracs.h"
#include "h.h"
#include "zonecnt.h"
#include "destcrt.h"
#include "charexc.h"
#include "nsshells.h"
#include "elementnames.h"
#include "collionrate.h"
#include "converge.h"
#include "hionfraccom.h"
#include "showme.h"
#include "negcon.h"
#include "bidiag.h"

void BiDiag(
	/* this is element on the c scale, H is 0 */
	long int nelem, 
	/* option to print this element when called */
	int lgPrintIt)
{
	int lgDone, 
	  lgNegPop;
	long int LoopBig, 
	  i, 
	  limit, 
	  MaxVal, 
	  nej, 
	  nelec, 
	  ns;
	float abundOld;
	double csav[32], 
	  ratio, 
	  sum;

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

	/* this is on the c scale */
	assert( nelem >= 0);

	/*begin sanity check
	 * impossible for HIonFrac[nelem] to be zero if IonHigh(nelem)=nelem+1
	 * HIonFrac(nelem) is stripped to hydrogen */
	if( (IonRange.IonHigh[nelem] == nelem + 2) && 
		HIonFraccom.HIonFrac[nelem] <= 0 )
	{
		fprintf( ioQQQ, " BiDiag: nelem=%3ld IonHigh=%3ld HIonFrac=%10.2e\n", 
		  nelem, IonRange.IonHigh[nelem], HIonFraccom.HIonFrac[nelem] );
		ShowMe();
	}
	/*end sanity check */


	/* this will be a sanity check */
	lgNegPop = FALSE;

	LoopBig = 0;
	lgDone = FALSE;
	while( LoopBig < 5 && (!lgDone) )
	{
		for( i=0; i < nelem; i++ )
		{
			DestCrt.destroy[i][nelem] = 0.;
			DestCrt.create[i][nelem] = 0.;
		}

		/* get total destruction rates */
		limit = MIN2(nelem,IonRange.IonHigh[nelem]-1);
		for( i=IonRange.IonLow[nelem]-1; i < limit; i++ )
		{
			assert( i>=0 );
			/* collisional ionization */
			DestCrt.destroy[i][nelem] += CollIonRate.CollidRate[0][i][nelem];

			/* photoionization loop */
			for( ns=0; ns < nsShellsCom.nsShells[i][nelem]; ns++ )
			{
				DestCrt.destroy[i][nelem] += PhotRate.PhotoRate[0][ns][i][nelem];

				/* following is most number of electrons freed */
				nelec = yield.nyield[ns][i][nelem];

				/* following loop is over number of electrons that come out of shell
				 * with multiple electron ejection */
				for( nej=2; nej <= nelec; nej++ )
				{
					/* this is highest possible stage of ionization -
					 * do not want to ignore ionization that go beyond this */

					MaxVal = MIN2(i+nej,IonRange.IonHigh[nelem]-1);
					/* if( xIonFracs(nelem,i+nej-1).gt.1e-30 ) then */
					if( xIonFracs[MaxVal][nelem] > 1e-30 )
					{
						/* i points to current stage of ionization */
						ratio = (double)(xIonFracs[i+1][nelem])/
						  (double)(xIonFracs[MaxVal][nelem]);
					}
					else
					{
						ratio = 1.;
					}
					/* yield here is fraction removing ion electrons */
					DestCrt.destroy[MaxVal-1][nelem] += PhotRate.PhotoRate[0][ns][i][nelem]*
					  yield.vyield[nej-1][ns][i][nelem]*ratio;
				}
			}
			/* this is charge transfer destruction */
			DestCrt.destroy[i][nelem] += CharExc.CTHeavy[0][i] + 
			  CharExc.HCharExcIon[nelem][i]*h.hii;

			/* recombination rates */
			DestCrt.create[i][nelem] = RecomRate.RecombinRate[i][nelem];

			/* charge transfer recombination */
			DestCrt.create[i][nelem] += CharExc.CTHeavy[1][i] + 
			  CharExc.HCharExcRec[nelem][i]*h.hi;
		}

		/* option to print ionization and recombination arrays
		 * parray flag set with "print array" command */
		if( parray.lgPrtArry || lgPrintIt )
		{
			/* total ionization rate, all processes */
			fprintf( ioQQQ, "\n %4.4s", elementnames.chElementNameShort[nelem] );
			fprintf( ioQQQ, " It" );
			for( i=0; i < (IonRange.IonHigh[nelem] - 1); i++ )
			{
				fprintf( ioQQQ, " %8.1e", DestCrt.destroy[i][nelem] );
			}
			fprintf( ioQQQ, "\n" );

			/* collisional ionization */
			fprintf( ioQQQ, " %4.4s", elementnames.chElementNameShort[nelem] );
			fprintf( ioQQQ, " C " );
			for( i=0; i < (IonRange.IonHigh[nelem] - 1); i++ )
			{
				fprintf( ioQQQ, " %8.1e", CollIonRate.CollidRate[0][i][nelem] );
			}
			fprintf( ioQQQ, "\n" );

			/* sum of all creation processes */
			fprintf( ioQQQ, " %4.4s", elementnames.chElementNameShort[nelem]);
			fprintf( ioQQQ, " R " );
			for( i=0; i < (IonRange.IonHigh[nelem] - 1); i++ )
			{
				fprintf( ioQQQ, " %8.1e", DestCrt.create[i][nelem] );
			}
			fprintf( ioQQQ, "\n" );

			/* charge exchange ionization */
			fprintf( ioQQQ, " char I " );
			for( i=0; i < (IonRange.IonHigh[nelem] - 1); i++ )
			{
				fprintf( ioQQQ, " %8.1e", CharExc.HCharExcIon[nelem][i]*
				  h.hii );
			}
			fprintf( ioQQQ, "\n" );

			/* charge exchange recombination */
			fprintf( ioQQQ, " char R " );
			for( i=0; i < (IonRange.IonHigh[nelem] - 1); i++ )
			{
				fprintf( ioQQQ, " %8.1e", CharExc.HCharExcRec[nelem][i]*
				  h.hi );
			}
			fprintf( ioQQQ, "\n" );
		}

		/* invert and solve bidiagonal matrix */
		/* >>chng 99 may 1, this had been set to 1 for all previous versions of the code,
		 * and range of ionization was limited by logic dealing with ionization parameter.
		 * that logic was removed (caused other problems) and now comphi.in crashes with
		 * overflow in loop below.  dur to csav getting larger than biggest double.
		 * change 1 to 1e-200 to get more range, could cause problems with very low
		 * ionization solns??? */
		/*csav[IonRange.IonLow[nelem]-1] = 1.e0;*/
		/*sum = 1.e0;*/
#		define CINIT 1e-200
		csav[IonRange.IonLow[nelem]-1] = CINIT;
		sum = CINIT;

		limit = MIN2(nelem+1,IonRange.IonHigh[nelem]);
		for( i=IonRange.IonLow[nelem]; i < limit; i++ )
		{
			assert( i>= 1);
			csav[i] = csav[i-1]*DestCrt.destroy[i-1][nelem]/DestCrt.create[i-1][nelem];
			sum += csav[i];
		}

		/* this is the test for hydrogenic species that appears throughout the code */
		if( IonRange.IonHigh[nelem] == nelem + 2 )
		{
			/* extends up to hydrogenic - use hydrogenic solutions,
			 * limit is the top of the loop above */
			i = limit;
			csav[i] = csav[i-1]*HIonFraccom.HIonFrac[nelem];
			sum += csav[i];
		}

		/* this is the factor needed to bring the above arry up to the total
		 * abundance of the elements, stored in xIOnFracs[0][nelem] */
		sum = (double)(xIonFracs[0][nelem])/sum;

		/* zero out abundances of ions below lowest stage we will consider */
		for( i=0; i < (IonRange.IonLow[nelem] - 1); i++ )
		{
			xIonFracs[i+1][nelem] = 0.;
		}

		/* zero out abundances of ions above highest stage we will consider */
		for( i=IonRange.IonHigh[nelem]; i < (nelem + 2); i++ )
		{
			xIonFracs[i+1][nelem] = 0.;
		}

		xIonFracs[IonRange.IonHigh[nelem]][nelem] = 
			(float)(csav[IonRange.IonHigh[nelem]-1]*sum);

		if( xIonFracs[IonRange.IonHigh[nelem]][nelem] < 0. )
		{
			lgNegPop = TRUE;
		}

		lgDone = TRUE;
		for( i=IonRange.IonLow[nelem]-1; i < (IonRange.IonHigh[nelem]-1); i++ )
		{
			if( conv.lgSearch && (xIonFracs[i+1][nelem] > 0.) && (csav[i]*sum> 0.) )
			{
				/* this branch used during search phase */
				/* >>chng 96 oct 27, sulphur and above oscillated, corrected by
				 * using log mean during search phase */
				abundOld = xIonFracs[i+1][nelem];

				/* >>chng 99 Jul 02, the log in the following is protected against 
				 * non-pos value by if branch above */
				xIonFracs[i+1][nelem] = 
					(float)(pow(10.,0.5*(log10(xIonFracs[i+1][nelem]) + 
				  log10(csav[i]*sum))));

				/* >>chng 97 jul 29, set flag to redo if changes happen */
				if( xIonFracs[i+1][nelem] > 1e-15 )
				{
					if( fabs(abundOld/xIonFracs[i+1][nelem]- 1.) > 0.2 )
					{
						lgDone = FALSE;

						/* these are printed as the old and new values that forced 
						 * declaration of no convergence */
						conv.BadConvIoniz[0] = abundOld;
						conv.BadConvIoniz[1] = xIonFracs[i+1][nelem];
					}
				}
			}
			else
			{
				/* this branch used for usual search for soln deep in model */
				abundOld = xIonFracs[i+1][nelem];
				xIonFracs[i+1][nelem] = (float)(csav[i]*sum);
				if( xIonFracs[i+1][nelem] > 1e-10 )
				{
					if( fabs(abundOld/xIonFracs[i+1][nelem]-1.) > 0.2 )
					{
						lgDone = FALSE;

						/* these are printed as the old and new values that forced 
						 * declaration of no convergence */
						conv.BadConvIoniz[0] = abundOld;
						conv.BadConvIoniz[1] = xIonFracs[i+1][nelem];
					}
				}
			}
			if( xIonFracs[i+1][nelem] < 0. )
			{
				lgNegPop = TRUE;
			}
		}

		if( parray.lgPrtArry || lgPrintIt )
		{
			/* this is the array that was inverted to produce abundances */
			fprintf( ioQQQ, " csav[] " );
			for( i=1; i <= IonRange.IonHigh[nelem]; i++ )
			{
				fprintf( ioQQQ, " %8.1e", csav[i-1] );
			}
			fprintf( ioQQQ, "\n" );
			/* these are the predicted abundances */
			fprintf( ioQQQ, " AbundX " );
			for( i=1; i <= IonRange.IonHigh[nelem]; i++ )
			{
				fprintf( ioQQQ, " %8.1e", xIonFracs[i][nelem] );
			}
			fprintf( ioQQQ, "\n" );
		}
		LoopBig += 1;
	}

	if( !lgDone )
	{
		conv.lgConvIoniz = FALSE;
		strcpy( conv.chConvIoniz,  "BiDiag:" );

		/* this is two-char short form of element name */
		strcat( conv.chConvIoniz, elementnames.chElementSym[nelem] );
	}

	/* 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], ZoneCnt.nzone );

		fprintf( ioQQQ, " Populations were" );
		for( i=1; i <= IonRange.IonHigh[nelem]; i++ )
		{
			fprintf( ioQQQ, "%9.1e", xIonFracs[i][nelem] );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, " destroy vector =" );
		for( i=1; i <= (IonRange.IonHigh[nelem] - 1); i++ )
		{
			fprintf( ioQQQ, "%9.1e", DestCrt.destroy[i-1][nelem] );
		}
		fprintf( ioQQQ, "\n" );

		/* check whether any negative level populations occured */
		lgNegPop = FALSE;
		for( i=0; i < (IonRange.IonHigh[nelem] - 1); i++ )
		{
			if( DestCrt.destroy[i][nelem] < 0. )
				lgNegPop = TRUE;
		}

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

			fprintf( ioQQQ, " HCharExcIon vtr=" );
			for( i=0; i < (IonRange.IonHigh[nelem] - 1); i++ )
			{
				fprintf( ioQQQ, "%9.1e", CharExc.HCharExcIon[nelem][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, " CollidRate  vtr=" );
			for( i=0; i < (IonRange.IonHigh[nelem] - 1); i++ )
			{
				fprintf( ioQQQ, "%9.1e", CollIonRate.CollidRate[0][i][nelem] );
			}
			fprintf( ioQQQ, "\n" );

			/* photo rates per subshell */
			fprintf( ioQQQ, " photo rates per subshell, ion\n" );
			for( i=0; i < (IonRange.IonHigh[nelem] - 1); i++ )
			{
				fprintf( ioQQQ, "%3ld", i );
				for( ns=0; ns < nsShellsCom.nsShells[i][nelem]; ns++ )
				{
					fprintf( ioQQQ, "%9.1e", PhotRate.PhotoRate[0][ns][i][nelem] );
				}
				fprintf( ioQQQ, "\n" );
			}
		}

		/* now check out creation vector */
		fprintf( ioQQQ, " create  vector =" );
		for( i=0; i < (IonRange.IonHigh[nelem] - 1); i++ )
		{
			fprintf( ioQQQ, "%9.1e", DestCrt.create[i][nelem] );
		}
		fprintf( ioQQQ, "\n" );

		negcon();
		ShowMe();
		puts( "[Stop in bidiag]" );
		exit(1);
	}

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