/*lgConverg check whether ionization of element nelem has converged, called by ionize */
#include "cddefines.h"
#include "ionfracs.h"
#include "zonecnt.h"
#include "lgconverg.h"

int lgConverg(
	/* atomic number on the C scale, 0 for H, 25 for Fe */
	long int nelem ,
	/* this is allowed error as a fractional change.  Most are 0.15 */
	double delta )
{
	int lgConverg_v;
	long int nstage;
	double Abund, 
	  change ,
	  one ;

	/* this is fractions [ion stage][nelem], ion stage = 1 for atom, nelem=0 for H*/
	static float OldFracs[LIMELM + 2][LIMELM];

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

	/* args are atomic number, ionization stage
	 * OldFracs[0][nelem] is abundances of nelem (cm^-3) */

	/* this function returns true if ionization of element
	 * with atomic number nelem+ has not changed by more than delta,*/

	/* check whether abundances exist yet, only do this for after first zone */
	if( ZoneCnt.nzone > 0 )
	{
		change = 0.;
		Abund = xIonFracs[0][nelem];

		/* loop is from 1 since 0 stores total abundance */
		for( nstage=1; nstage <= (nelem + 2); nstage++ )
		{
			if( OldFracs[nstage][nelem]/Abund > 1e-4 && 
				xIonFracs[nstage][nelem]/Abund > 1e-4 )
			{
				/* check change in old to new value */
				one = fabs(xIonFracs[nstage][nelem]-OldFracs[nstage][nelem])/
					OldFracs[nstage][nelem];
				change = MAX2(change, one );
			}
			/* now set new value */
			OldFracs[nstage][nelem] = xIonFracs[nstage][nelem];
		}

		if( change < delta )
		{
			lgConverg_v = TRUE;
		}
		else
		{
			lgConverg_v = FALSE;
		}
	}
	else
	{
		for( nstage=1; nstage <= (nelem + 2); nstage++ )
		{
			OldFracs[nstage][nelem] = xIonFracs[nstage][nelem];
		}

		lgConverg_v = TRUE;
	}


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