/*TrimStages raise or lower most extreme stages of ionization considered,
 * called by ionize */
#include "cddefines.h"
#include "ionfracs.h"
#include "trimstg.h"
#include "ionrange.h"
#include "heating.h"
#include "nodec.h"
#include "converge.h"
#include "trimstages.h"

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

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

	/*confirm all numbers are within their range of validity */
	assert( nelem > 0 && nelem < LIMELM );
	assert( IonRange.IonLow[nelem] > 0 );
	assert( IonRange.IonHigh[nelem] <= nelem+2 );

	/* raise or lower highest and lowest stages of ionization to
	 * consider only significant stages
	 * IonLow(nelem) lowest stage of ionization
	 * IonHigh(nelem)  stage of ionization*/

	/* special trim if upper stages have zero abundance */
	while( 
		xIonFracs[IonRange.IonHigh[nelem]][nelem]/xIonFracs[0][nelem] < SMALLFLOAT && 
		IonRange.IonHigh[nelem] > IonRange.IonLow[nelem] + 1 )
	{
		/* xIonFracs(nelem,i) is density of ith ionization stage (cm^-3) */
		xIonFracs[IonRange.IonHigh[nelem]][nelem] = 0.;
		HeatingCom.heating[IonRange.IonHigh[nelem]-2][nelem] = 0.;
		--IonRange.IonHigh[nelem];
	}

	/* only have atom and first ion, do not lower any more */
	if( IonRange.IonHigh[nelem] == 2 )
	{ 
#	ifdef DEBUG_FUN
		fputs( " <->TrimStages()\n", debug_fp );
#	endif
		return;
	}

	/* trim down high stages */
	if( 
		(!nodec.lgNoDec && 
		xIonFracs[IonRange.IonHigh[nelem]][nelem]/xIonFracs[0][nelem] <= trimstg.trimhi) && 
		IonRange.IonHigh[nelem] > 3 )
	{
		xIonFracs[IonRange.IonHigh[nelem]][nelem] = 0.;
		HeatingCom.heating[IonRange.IonHigh[nelem]-2][nelem] = 0.;
		--IonRange.IonHigh[nelem];
	}

	/* lower lowest stage of ionization */
	if( 
		xIonFracs[IonRange.IonLow[nelem]][nelem]/xIonFracs[0][nelem] > 1e-3 && 
		IonRange.IonLow[nelem] > 1 )
	{
		/* lower lowest level of ionization */
		--IonRange.IonLow[nelem];

		/* tell code to do all photo rates only in case where lowest is lowered still */
		conv.lgIonStageTrimed = TRUE;
	}

	/* raise lowest stage of ionization */
	else if( 
		xIonFracs[IonRange.IonLow[nelem]][nelem]/xIonFracs[0][nelem] <= trimstg.trimlo && 
		IonRange.IonLow[nelem] < (IonRange.IonHigh[nelem] - 2) )
	{
		/* raise lowest level of ionization */
		xIonFracs[IonRange.IonLow[nelem]][nelem] = 0.;
		HeatingCom.heating[IonRange.IonLow[nelem]-1][nelem] = 0.;
		++IonRange.IonLow[nelem];
	}


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