/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  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 "elementnames.h"
#include "radius.h"
#include "heavy.h"
#include "conv.h"
#include "rfield.h"
#include "phycon.h"
#include "mole.h"
#include "thermal.h"
#include "iso.h"
#include "dense.h"
#include "conv.h"
#include "struc.h"
#include "ionbal.h"
#include "taulines.h"

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 */
	bool lgHi_Down = false;
	bool lgHi_Up = false;
	bool lgHi_Up_enable;
	/* this will remember that lower stages trimmed up or own*/
	bool lgLo_Up = false;
	bool lgLo_Down = false;
	long int ion_lo_save = dense.IonLow[nelem],
		ion_hi_save = dense.IonHigh[nelem];
	long int ion;
	realnum trimhi , trimlo;

	/* this is debugging code that turns on print after certain number of calls */
	/*realnum xlimit = SMALLFLOAT *10.;*/
	/*static int ncall=0;
	if( nelem==5 )
		++ncall;*/

	DEBUG_ENTRY( "ion_trim()" );

	/* this routine should not be called early in the simulation, before
	 * things have settled down */
	ASSERT( conv.nTotalIoniz>2 );

	/*confirm all ionization stages 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 */
	ASSERT( dense.IonLow[nelem] <= dense.IonHigh[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 = (realnum)ionbal.trimhi * 1e-4f;
		trimlo = (realnum)ionbal.trimlo * 1e-4f;
	}
	else
	{
		trimhi = (realnum)ionbal.trimhi;
		trimlo = (realnum)ionbal.trimlo;
	}

	/* helium is special case since abundance is so high, and He+ CT with 
	 * CO is the dominant CO destruction process in molecular regions */
	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 )
		{
			if( conv.lgSearch )
			{
				/* during search phase we may be quite far from solution, in the
				 * case of a PDR sim the electron density may be vastly higher 
				 * than it will be with stable solution.  He++ can be very important
				 * for the chemistry and we want to consider it.  Make the
				 * threshold for ignoring He++ much higher than normal to prevent
				 * a premature removal of the ion */
				trimhi = MIN2( trimhi , 1e-17f ); 
			}
			else
			{
				/* similar smaller upper limit for ion*/
				trimhi = MIN2( trimhi , 1e-12f ); 
			}
		}
	}

	/* logic for PDRs, for elements included in chemistry, need stable solutions, 
	 * keep 3 ion stages in most cases - added by NA to do HII/PDR sims */
	if( !mole_global.lgNoMole )
	{
		trimlo = SMALLFLOAT;
		if(dense.IonHigh[nelem] ==2)
		{
			trimhi = MIN2(trimhi, 1e-20f);
		}
	}

	/* raise or lower highest and lowest stages of ionization to
	 * consider only significant stages
	 * IonHigh[nelem]  stage of ionization  */

	/* 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] < dense.density_low_limit || 
			dense.xIonDense[nelem][dense.IonHigh[nelem]] < dense.density_low_limit ) && 
			/* >>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]-1] +=
				dense.xIonDense[nelem][dense.IonHigh[nelem]];
			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 );
				for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
					iso_sp[ipISO][nelem].st[level].Pop() = 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] < dense.density_low_limit || 
			dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit ) && 
			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]+1] += 
				dense.xIonDense[nelem][dense.IonLow[nelem]];
			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 );
				for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
					iso_sp[ipISO][nelem].st[level].Pop() = 0.;
			}

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

	/* sanity check */
	ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );

	/* this checks on error condition where the gas is stupendously highly ionized,
	 * the low limit is already one less than high, and we need to raise the low
	 * limit further */
	if( dense.IonHigh[nelem] > 1 &&
		(dense.IonLow[nelem]==dense.IonHigh[nelem]-1) &&
		(dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit) )
	{
#		if 0
		fprintf(ioQQQ,
			"PROBLEM DISASTER the density of ion %li of element %s is too small to be computed on this cpu.\n",
			dense.IonLow[nelem]+1,
			elementnames.chElementName[nelem]);
		fprintf(ioQQQ,
			"Turn off the element with the command ELEMENT XXX OFF or do not consider "
			"gas with low density, the current hydrogen density is %.2e cm-3.\n",
			dense.gas_phase[ipHYDROGEN]);
		fprintf(ioQQQ,
			"The outer radius of the cloud is %.2e cm - should a stopping "
			"radius be set?\n",
			radius.Radius );
		fprintf(ioQQQ,
			"The abort flag is being set - the calculation is stopping.\n");
#		endif
		//lgAbort = true;
		dense.xIonDense[nelem][dense.IonLow[nelem]] = (realnum)dense.density_low_limit;
		return;
	}

	/* trim down high stages that have too small an abundance */
	if( 
		dense.IonHigh[nelem] > dense.IonLow[nelem]  && 
		( (dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] <= 
		trimhi ) ||
		(dense.xIonDense[nelem][dense.IonHigh[nelem]] <= dense.density_low_limit ) )
		) 
	{
		/* >>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.*dense.density_low_limit) )
		{
			dense.xIonDense[nelem][dense.IonHigh[nelem]-1] += 
				dense.xIonDense[nelem][dense.IonHigh[nelem]];
			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 );
				for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
					iso_sp[ipISO][nelem].st[level].Pop() = 0.;
			}
			--dense.IonHigh[nelem];
			lgHi_Down = true;
		}
	}

	/* trim up highest stages - will this be done? */
	lgHi_Up_enable = true;
	/* trimming can be turned off */
	if( !ionbal.lgTrimhiOn )
		lgHi_Up_enable = false;
	/* high stage is already fully stripped */
	if( dense.IonHigh[nelem]>=nelem+1 )
		lgHi_Up_enable = false;
	/* we have previously trimmed this stage down */
	if( lgHi_Down )
		lgHi_Up_enable = false;
	/* we are in neutral gas */
	if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]<0.9  )
		lgHi_Up_enable = false;

	if( lgHi_Up_enable )
	{
		realnum abundold = 0;
		if( nzone>2 ) 
			abundold = struc.xIonDense[nelem][dense.IonHigh[nelem]][nzone-3]/
			SDIV( struc.gas_phase[nelem][nzone-3]);
		realnum 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 */
		if( (Heavy.Valence_IP_Ryd[nelem][dense.IonHigh[nelem]] < rfield.anu[rfield.nflux-1] ||
			Heavy.Valence_IP_Ryd[nelem][dense.IonHigh[nelem]] < phycon.te_ryd*10.) &&
			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;
			/* set abundance to almost zero so that sanity check elsewhere will be ok */
			dense.xIonDense[nelem][dense.IonHigh[nelem]] = (realnum)(dense.density_low_limit*10.);
			dense.xIonDense[nelem][dense.IonHigh[nelem]-1] -= 
				dense.xIonDense[nelem][dense.IonHigh[nelem]];
			long int ipISO = nelem - dense.IonHigh[nelem];
			if (ipISO >= 0 && ipISO < NISO)
				iso_sp[ipISO][nelem].st[0].Pop() = dense.xIonDense[nelem][dense.IonHigh[nelem]];
		}
	}

	/* sanity check */
	ASSERT( dense.IonLow[nelem] <= dense.IonHigh[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]] = (realnum)dense.density_low_limit;
		dense.xIonDense[nelem][dense.IonLow[nelem]+1] -= 
			dense.xIonDense[nelem][dense.IonLow[nelem]];
		long int ipISO = nelem - dense.IonLow[nelem];
		if (ipISO >= 0 && ipISO < NISO)
			iso_sp[ipISO][nelem].st[0].Pop() = dense.xIonDense[nelem][dense.IonLow[nelem]];
	}

	/* raise lowest stage of ionization, but only if we are near illuminated face of cloud*/
	else if( nzone < 10 &&
		(dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] <= (realnum)trimlo) && 
		(dense.IonLow[nelem] < (dense.IonHigh[nelem] - 2) ) )
	{
		/* raise lowest level of ionization */
		dense.xIonDense[nelem][dense.IonLow[nelem]+1] += 
			dense.xIonDense[nelem][dense.IonLow[nelem]];
		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 );
			for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
				iso_sp[ipISO][nelem].st[level].Pop() = 0.;
		}
		++dense.IonLow[nelem];
		lgLo_Up = true;
	}
	/* test on zero */
	else if( (dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit) && 
		/*>>chng 06 may 24, from IonLow < IonHigh to <IonHigh-1 because
		 * old test would allow low to be raised too high, which then
		 * leads to insanity */
		(dense.IonLow[nelem] < dense.IonHigh[nelem]-1) )
	{
		while(dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit && 
			/* >>chng 07 feb 27 from < IonHigh to < IonHigh-1 to prevent raising
			 * low to high */
			(dense.IonLow[nelem] < dense.IonHigh[nelem]-1) )
		{
			/* raise lowest level of ionization */
			dense.xIonDense[nelem][dense.IonLow[nelem]+1] += 
				dense.xIonDense[nelem][dense.IonLow[nelem]];
			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 );
				for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
					iso_sp[ipISO][nelem].st[level].Pop() = 0.;
			}
			++dense.IonLow[nelem];
			lgLo_Up = true;
		}
	}

	/* sanity check */
	ASSERT( dense.IonLow[nelem] <= dense.IonHigh[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 dense.density_low_limit so test must be < dense.density_low_limit */
	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]] >= dense.density_low_limit ||
		dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] >= dense.density_low_limit ||
		/*>>chng 06 may 24, include this to cover case where cap is set by not allowing
		 * lower limit to come up to upper limit */
		(dense.IonLow[nelem]==dense.IonHigh[nelem]-1 ));

	{
		/* option to print out what has happened so far .... */
		enum {DEBUG_LOC=false};
		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. );

	/* option to print if any trimming occurred */
	if( lgHi_Down || lgHi_Up || lgLo_Up || lgLo_Down )
	{
		conv.lgIonStageTrimed = true;
		{
			/* option to print out what has happened so far .... */
			enum {DEBUG_LOC=false};
			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");
			}
		}
	}

	/* asserts that that densities of ion stages that are now turned 
	 * off are zero */
	for( ion=0; ion<dense.IonLow[nelem]; ++ion )
	{
		ASSERT( dense.xIonDense[nelem][ion] == 0. );
	}
	for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
	{
		ASSERT( dense.xIonDense[nelem][ion] == 0. );
	}

	for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
	{
		/* in case where ionization stage was changed the
		 * density is zero, we set it to SMALLFLOAT since a density
		 * of zero is not possible */
		dense.xIonDense[nelem][ion] = MAX2(dense.xIonDense[nelem][ion], SMALLFLOAT);
	}
	return;
}