/*ChkRate called by nextdr to check how rates of destruction of various species changes */
#include "cddefines.h"
#include "zonecnt.h"
#include "ionfracs.h"
#include "destcrt.h"
#include "elmton.h"
#include "chkrate.h"

void ChkRate(long int nelem, 
	  double *dDestRate, 
	  long int *istage)
{
	long int i;

	double average, 
	  dDest;

	static double OldDest[LIMELM][LIMELM];

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

	/* return if this element is not turned on */
	if( !elmton.lgElmtOn[nelem-1] )
	{
		*dDestRate = 1e-3;
		
#	ifdef DEBUG_FUN
		fputs( " <->ChkRate()\n", debug_fp );
#	endif
		return;
	}

	*dDestRate = 0.;
	if( ZoneCnt.nzone <= 1 )
	{
		for( i=0; i < nelem; i++ )
		{
			OldDest[i][nelem-1] = DestCrt.destroy[i][nelem-1];
		}
	}

	else if( xIonFracs[1][nelem-1]/xIonFracs[0][nelem-1] <  0.9 )
	{
		/* do not use this method if everything is atomic */
		for( i=0; i < (nelem - 1); i++ )
		{
			/* last check below, .5 chosen so that do not key off
			 * predominantly neutral species where self-opacity
			 * could cause oscillations */
			if( ((xIonFracs[i+1][nelem-1]/xIonFracs[0][nelem-1] > 0.01 && 
				xIonFracs[i+1][nelem-1]/xIonFracs[0][nelem-1] < 0.9) && 
				xIonFracs[i+2][nelem-1]/xIonFracs[0][nelem-1] > .05) && 
				OldDest[i][nelem-1] > 0. )
			{
				/* last check on old dest in case just lowered ionization
				 * stage, so no history */
				/* check that rate is positive */
				if( DestCrt.destroy[i][nelem-1] <= 0. )
				{
					fprintf( ioQQQ, " ChkRate gets insane destruction rate for ion%4ld%4ld%10.2e\n", 
					  nelem, i, DestCrt.destroy[i][nelem-1] );
					puts( "[Stop in chkrate]" );
					exit(1);
				}

				/* do not consider unless of middling ionization, and
				 * rate is going down (to prevent dr osciallating)
				 * no absolute value in following since do not want to
				 * consider cases where ionization rate increases */
				average = (OldDest[i][nelem-1] + DestCrt.destroy[i][nelem-1])* 0.5;

				dDest = (OldDest[i][nelem-1] - (float)(DestCrt.destroy[i][nelem-1]))/ average;
				/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ term + if rate going down */

				if( *dDestRate < dDest )
				{
					*dDestRate = dDest;
					*istage = i+1;
				}

			}
			OldDest[i][nelem-1] = DestCrt.destroy[i][nelem-1];
		}
	}

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