/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*ion_trim raise or lower most extreme stages of ionization considered,
 * called by ConvBase - ion limits were originally set by  */
#include "cddefines.h"
#include "heavy.h"
#include "conv.h"
#include "rfield.h"
#include "thermal.h"
#include "iso.h"
#include "dense.h"
#include "struc.h"
#include "ionbal.h"
/*lint -e661 Possible access of out-of-bounds pointer*/

void ion_trim(
	/* nelem is on the C scale, 0 for H, 5 for C */
	long int nelem )
{

	/* this will remember that higher stages trimed up or down */
	int lgHi_Down = FALSE;
	int lgHi_Up = FALSE;
	int lgHi_Up_enable;
	/* this will remember that lower stages trimed up or own*/
	int lgLo_Up = FALSE;
	int lgLo_Down = FALSE;
	long int ion_lo_save = dense.IonLow[nelem],
		ion_hi_save = dense.IonHigh[nelem];
	float trimhi , trimlo;
	float xlimit = SMALLFLOAT *10.;
	/*static int ncall=0;
	if( nelem==5 )
		++ncall;*/

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

	/*confirm all numbers are within their range of validity */
	ASSERT( nelem >= ipHELIUM && nelem < LIMELM );
	ASSERT( dense.IonLow[nelem] >= 0 );
	ASSERT( dense.IonHigh[nelem] <= nelem+1 );
	/* IonHigh can be equal to IonLow if both are zero - so no ionization,
	 * or is ionization has been set with "element ionization" command */
	ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] || 
		(dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 )  || 
		dense.lgSetIoniz[nelem] );

	/* during search phase of mostly neutral matter the electron density
	 * can be vastly too large, and the ionization suppressed as a result.
	 * during search do not trim down or up as much */
	if( conv.lgSearch )
	{
		trimhi = (float)ionbal.trimhi * 1e-4f;
		trimlo = (float)ionbal.trimlo * 1e-4f;
	}
	else
	{
		trimhi = (float)ionbal.trimhi;
		trimlo = (float)ionbal.trimlo;
	}

	/* >>chng 04 jun 04, move he trim here from con ion do */
	/* helium is special case since abundance is so high */
	if( nelem == ipHELIUM )
	{
		/* never want to trip up a lower stage of ionization */
		trimlo = SMALLFLOAT;

		/* if He+ is highest stage of ionization, probably want to keep it
		 * since important for CO chemistry in molecular clouds */
		if( dense.IonHigh[ipHELIUM] == 1 )
		{
			trimhi = MIN2( trimhi , 1e-20f );
		}
		else if( dense.IonHigh[ipHELIUM] == 2 )
		{
			/* >>chng 04 jul 27, add smaller upper limit for ion*/
			trimhi = MIN2( trimhi , 1e-12f ); 
			/* >>chng 05 jul 09, use same high limit for He++ due to large
			 * jump in chem when turned off 
			TestCode();*//* comment out 3 line above, put in following
			trimhi = MIN2( trimhi , 1e-20f ); */
		}
	}

	/* >>chng 04 aug 15, add logic for PDRs, these elements are
	 * included in chemistry, need stable solns, keep 3 ion stages
	 * in most cases - added by NA to do HII/PDR sims */
	else if (nelem == ipCARBON)
	{
		trimlo = SMALLFLOAT;
		if (dense.IonHigh[ipCARBON] ==2)
		{
			trimhi = MIN2(trimhi, 1e-20f);
		}
	}

	else if (nelem == ipOXYGEN)
	{
		trimlo = SMALLFLOAT;
		if (dense.IonHigh[ipOXYGEN] ==2)
		{
			trimhi = MIN2(trimhi, 1e-20f);
		}
	}

	else if (nelem == ipNITROGEN)
	{
		trimlo = SMALLFLOAT;
		if (dense.IonHigh[ipNITROGEN] ==2)
		{
			trimhi = MIN2(trimhi, 1e-20f);
		}
	}

	else if (nelem == ipSILICON)
	{
		trimlo = SMALLFLOAT;
		if (dense.IonHigh[ipSILICON] ==2)
		{
			trimhi = MIN2(trimhi, 1e-20f);
		}
	}

	else if (nelem == ipSULPHUR)
	{
		trimlo = SMALLFLOAT;
		if (dense.IonHigh[ipSULPHUR] ==2)
		{
			trimhi = MIN2(trimhi, 1e-20f);
		}
	}
	/* will we consider raising the highest stage of ionization?  */
	lgHi_Up_enable = FALSE;
	if( nzone > 5 )
	{
		float abundnew = dense.xIonDense[nelem][dense.IonHigh[nelem]]/SDIV( dense.gas_phase[nelem]);
		float abundold = struc.xIonDense[nelem][dense.IonHigh[nelem]][nzone-3]/SDIV( struc.gas_phase[nelem][nzone-3]);
		if( abundnew / SDIV( abundold ) > 1. )
		{
			lgHi_Up_enable = TRUE;
			/*if( nelem==1) 
			{
				fprintf(ioQQQ,"DEBUG\t%e\t%e\t%e\t%e\t%e\n",
					abundnew , abundold , abundnew / SDIV( abundold ),
					dense.xIonDense[ipHELIUM][1],dense.xIonDense[ipHELIUM][2] );
			}*/
		}
	}

	/* raise or lower highest and lowest stages of ionization to
	 * consider only significant stages
	 * IonLow[nelem] lowest stage of ionization on C scale, atom == 0
	 * IonHigh[nelem]  stage of ionization
	 *
	 * NB  NB  NB  NB  dense.xIonDense must have 1 added to above to get over tot abund */

	/* this is a special block for initialization only - it checks on absurd abundances
	 * and can trim multiple stages of ionization at one time. */
	if( conv.lgSearch )
	{
		/* special - trim down higher stages if they have essentially zero abundance */
		while( 
			(dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] < xlimit || 
			dense.xIonDense[nelem][dense.IonHigh[nelem]] < xlimit ) && 
			/* >>chng 02 may 12, rm +1 since this had effect of not allowing fully atomic */
			dense.IonHigh[nelem] > dense.IonLow[nelem] )
		{
			/* dense.xIonDense[nelem][i] is density of i+1 th ionization stage (cm^-3)
			 * the -1 is correct for the heating, -1 since heating is caused by ionization of stage below it */
			dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
			thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.;
			if( dense.IonHigh[nelem] == nelem+1-NISO )
			{
				long int ipISO = nelem - dense.IonHigh[nelem];
				ASSERT( ipISO>=0 && ipISO<NISO );
				iso.Pop2Ion[ipISO][nelem][0] = 0.;
			}

			/* decrement high stage limit */
			--dense.IonHigh[nelem];
			ASSERT( dense.IonHigh[nelem] >= 0);
			/* remember this happened */
			lgHi_Down = TRUE;
		}

		/* special - trim up lower stages trim if they have essentially zero abundance */
		while( 
			(dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] < xlimit || 
			dense.xIonDense[nelem][dense.IonLow[nelem]] < xlimit ) && 
			dense.IonLow[nelem] < dense.IonHigh[nelem] - 1 )
		{
			/* dense.xIonDense[nelem][i] is density of ith ionization stage (cm^-3)
			 * there is no-1 since we are removing the agent that heats */
			dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
			/* >>chng 01 aug 04, remove -1 which clobbers thermal.heating when IonLow == 0 */
			/*thermal.heating[nelem][dense.IonLow[nelem]-1] = 0.;*/
			thermal.heating[nelem][dense.IonLow[nelem]] = 0.;
			if( dense.IonLow[nelem] == nelem+1-NISO )
			{
				long int ipISO = nelem - dense.IonLow[nelem];
				ASSERT( ipISO>=0 && ipISO<NISO );
				iso.Pop2Ion[ipISO][nelem][0] = 0.;
			}

			/* increment low stage limit */
			++dense.IonLow[nelem];
			lgLo_Up = TRUE;
		}
	}

	/* sanity check */
	/* IonHigh can be equal to IonLow if both are zero - no ionization*/
	ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] || 
		(dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 )  || 
		dense.lgSetIoniz[nelem] );

	/* trim down high stages that have too small an abundance */
	/* >>chng 01 apr 21, last comparison had been IonHigh > 2, change to relative to IonLow */
	if( 
		/*!ionbal.lgNoDec && dense.IonHigh[nelem] > dense.IonLow[nelem]+1  && */
		/* >>chng 03 sep 19, from > IonLow+1, which was one too many, when totally neutral,
		 * IonHigh and IonLow are equal, and both 0, since loop is <=IonHigh */
		!ionbal.lgNoDec && dense.IonHigh[nelem] > dense.IonLow[nelem]  && 
		( (dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] <= 
		trimhi ) ||
		(dense.xIonDense[nelem][dense.IonHigh[nelem]] <= xlimit ) )
		) 
	{
		/* >>chng 03 sep 30, the atom and its first ion are a special case
		 * since we want to compute even trivial ions in molecular clouds */
		if( dense.IonHigh[nelem]>1 ||
			(dense.IonHigh[nelem]==1&&dense.xIonDense[nelem][1]<100.*xlimit) )
		{
			dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
			thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.;
			if( dense.IonHigh[nelem] == nelem+1-NISO )
			{
				long int ipISO = nelem - dense.IonHigh[nelem];
				ASSERT( ipISO>=0 && ipISO<NISO );
				iso.Pop2Ion[ipISO][nelem][0] = 0.;
			}
			--dense.IonHigh[nelem];
			lgHi_Down = TRUE;
		}
	}

	/* >>chng 02 marf 24, this option added - previously code did not do increasing ionization*/
	/* option to trim up highest stages */
	/* only consider this if we are fairly deep into the model, there are
	 * unconsidered ionization stages above the current limit, and the upper
	 * stage was not just trimmed down,
	 * last test is to avoid this logic in vastly deep neutral zone */
	if( lgHi_Up_enable && ionbal.lgTrimhiOn && (nzone > 10 )
		/* >>chng 02 nov 27, had been nelem+1 which allowed array overwrite.
		 * should be nelem since only want to increment bounds if one below highest */
		/*&& (dense.IonHigh[nelem]<nelem+1 )*/
		&& (dense.IonHigh[nelem]<nelem+1 )
		&& (!lgHi_Down )
		&& (dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]>0.9 ) 
		/* >>chng 03 oct 14, do not keep changing the trim after the first two calls within
		 * this zone - doing so introduced very small level of noise as some stages
		 * went up and down - O+2 - O+3 was good example, when small H+3 after He+ i-front */
		 /* >>chng 03 nov 30, limit to one increase per element per zone */
		/*&& conv.nPres2Ioniz < 10*/
		&& !conv.nPres2Ioniz
		)
	{
		float abundold = struc.xIonDense[nelem][dense.IonHigh[nelem]][nzone-3]/SDIV( struc.gas_phase[nelem][nzone-3]);
		float abundnew = dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem];
		/* only raise highest stage if ionization potential of next highest stage is within
		 * continuum array and the abundance of the highest stage is significant */
		/* in this one, IonHigh+1 since want next higher stage of ionization */
		if( Heavy.Valence_IP_Ryd[nelem][dense.IonHigh[nelem]+1] < rfield.anu[rfield.nflux-1] &&
			/* in this one, IonHigh+1 because of offset of abundance scale, with [0] being atom*/
			/* >>chng 03 nov 30, add 100x to trim limit */
			dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] > 1e-4f &&
			/* this checks that abundance of highest stage is increasing */
			abundnew > abundold*1.01 )
		{
			/*fprintf(ioQQQ,"uuppp %li %li \n", nelem, dense.IonHigh[nelem] );*/
			/* raise highest level of ionization */
			++dense.IonHigh[nelem];
			lgHi_Up = TRUE;
			/* >>chng 01 aug 02, set this to zero so that sanity check elsewhere will be ok */
			dense.xIonDense[nelem][dense.IonHigh[nelem]] = xlimit*10.f;
		}
	}

	/* sanity check */
	/* IonHigh can be equal to IonLow if both are zero - no ionization*/
	ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] || 
		(dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 )  || 
		dense.lgSetIoniz[nelem] );

	/* lower lowest stage of ionization if we have significant abundance at current lowest */
	if( dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] > 1e-3f && 
		dense.IonLow[nelem] > 0 )
	{
		/* lower lowest level of ionization */
		--dense.IonLow[nelem];
		lgLo_Down = TRUE;
		/* >>chng 01 aug 02, set this to zero so that sanity check elsewhere will be ok */
		dense.xIonDense[nelem][dense.IonLow[nelem]] = xlimit;
	}

	/* raise lowest stage of ionization, but only if we are near illuminated face of cloud*/
	/* >>chng 01 aug 08, add test on nzones */
	/* >>chng 02 feb 07, this logic was inverted!  change from nzone > 1 to nzone < 10 */
	/*else if( nzone > 1 &&*/
	else if( nzone < 10 &&
		(dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] <= (float)trimlo) && 
		(dense.IonLow[nelem] < (dense.IonHigh[nelem] - 2) ) )
	{
		/* raise lowest level of ionization */
		dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
		/* no minus one in below since this is low bound, already bounds at atom */
		thermal.heating[nelem][dense.IonLow[nelem]] = 0.;
		if( dense.IonLow[nelem] == nelem+1-NISO )
		{
			long int ipISO = nelem - dense.IonLow[nelem];
			ASSERT( ipISO>=0 && ipISO<NISO );
			iso.Pop2Ion[ipISO][nelem][0] = 0.;
		}
		++dense.IonLow[nelem];
		lgLo_Up = TRUE;
	}
	/* >>chng 03 oct 14, add this test on zero */
	else if( dense.xIonDense[nelem][dense.IonLow[nelem]] < xlimit && (dense.IonLow[nelem] < dense.IonHigh[nelem]) )
	{
		while(dense.xIonDense[nelem][dense.IonLow[nelem]] < xlimit && (dense.IonLow[nelem] < dense.IonHigh[nelem]) )
		{
			/* raise lowest level of ionization */
			dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
			/* no minus one in below since this is low bound, already bounds at atom */
			thermal.heating[nelem][dense.IonLow[nelem]] = 0.;
			if( dense.IonLow[nelem] == nelem+1-NISO )
			{
				long int ipISO = nelem - dense.IonLow[nelem];
				ASSERT( ipISO>=0 && ipISO<NISO );
				iso.Pop2Ion[ipISO][nelem][0] = 0.;
			}
			++dense.IonLow[nelem];
			lgLo_Up = TRUE;
		}
	}

	/* sanity check */
	/* IonHigh can be equal to IonLow if both are zero - no ionization*/
	ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] || 
		(dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 )  || 
		dense.lgSetIoniz[nelem] );

	/* these are standard bounds checks that appear throughout this routine
	 * dense.xIonDense[IonLow] is first one with positive abundances
	 * 
	 * in case where lower ionization stage was just lowered the abundance
	 * was set to xlimit so test must be < xlimit */
	ASSERT( dense.IonLow[nelem] <= 1 ||
		dense.xIonDense[nelem][dense.IonLow[nelem]-1] == 0. );

	ASSERT( (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) || lgLo_Up ||
		dense.xIonDense[nelem][dense.IonLow[nelem]] >= xlimit ||
		dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] >= xlimit );

	{
		/* option to print out what has happened so far .... */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/

		if ( DEBUG_LOC && nelem == ipHELIUM )
		{
			if( lgHi_Down ||lgHi_Up ||lgLo_Up ||lgLo_Down )
			{
				fprintf(ioQQQ,"DEBUG TrimZone\t%li\t",nzone );
				if(  lgHi_Down )
				{
					fprintf(ioQQQ,"high dn %li to %li",
						ion_hi_save , 
						dense.IonHigh[nelem] );
				}
				if( lgHi_Up )
				{
					fprintf(ioQQQ,"high up %li to %li",
						ion_hi_save , 
						dense.IonHigh[nelem] );
				}
				if( lgLo_Up )
				{
					fprintf(ioQQQ,"low up %li to %li",
						ion_lo_save , 
						dense.IonLow[nelem] );
				}
				if( lgLo_Down )
				{
					fprintf(ioQQQ,"low dn %li to %li",
						ion_lo_save , 
						dense.IonLow[nelem] );
				}
				fprintf(ioQQQ,"\n" );
			}
		}
	}

	/* only assert that the stage above the highest has a population of
	 * zero if that stage exists */
	if(dense.IonHigh[nelem] < nelem+1 )
		ASSERT( dense.xIonDense[nelem][dense.IonHigh[nelem]+1] == 0. );

	/* >>chng 03 nov 30, conv.lgIonStageTrimed had not been set!! */
	if( lgHi_Down || lgHi_Up || lgLo_Up || lgLo_Down )
	{
		conv.lgIonStageTrimed = TRUE;
		{
			/* option to print out what has happened so far .... */
			/*@-redef@*/
			enum {DEBUG_LOC=FALSE};
			/*@+redef@*/
			if ( DEBUG_LOC && nelem==ipHELIUM )
			{
				fprintf(ioQQQ,"DEBUG ion_trim zone\t%.2f \t%li\t", fnzone, nelem);
				if( lgHi_Down  )
					fprintf(ioQQQ,"\tHi_Down");
				if( lgHi_Up )
					fprintf(ioQQQ,"\tHi_Up");
				if( lgLo_Up )
					fprintf(ioQQQ,"\tLo_Up");
				if( lgLo_Down )
					fprintf(ioQQQ,"\tLo_Down");
				fprintf(ioQQQ,"\n");
			}
		}
	}
	/* following block only set of asserts */
#	if !defined(NDEBUG)
	{
		long int ion;

		/* check that proper abundances are either positive or zero */
		for( ion=0; ion<dense.IonLow[nelem]; ++ion )
		{
			ASSERT( dense.xIonDense[nelem][ion] == 0. );
		}
		/*if( nelem==5 ) fprintf(ioQQQ,"carbbb\t%li\n", dense.IonHigh[nelem]);*/
		for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
		{
			/* >>chng 02 feb 06, had been > o., chng to > SMALLFLOAT to
			 * trip over VERY small floats that failed on alphas, but not 386
			 * 
			 * in case where lower ionization stage was just lowered or
			 * trimmed down the abundance
			 * was set to SMALLFLOAT so test must be < SMALLFLOAT */
			/* >>chng 02 feb 19, add check for search phase.  During this search
			 * models with extreme ionization (all neutral or all ionized) can
			 * have extreme but non-zero abundances far from the ionization peak for
			 * element with lots of electrons.  These will go away once the model
			 * becomes stable */
			/* >>chng 03 dec 01, add check on whether ion trim was called 
			 * conserve.in threw assert when iontrim not called and abund grew small */
			dense.xIonDense[nelem][ion] = MAX2(dense.xIonDense[nelem][ion], SMALLFLOAT);
		}
		for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
		{
			ASSERT( dense.xIonDense[nelem][ion] == 0. );
		}
	}
#	endif

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