/*CollidIonize fill in collisional ionization rates, and resulting cooling */
#include "cddefines.h"
#include "physconst.h"
#include "secondaries.h"
#include "phycon.h"
#include "ionrange.h"
#include "rfield.h"
#include "heavy.h"
#include "collionrate.h"
#include "showme.h"
#include "cfit.h"
#include "collidionize.h"

void CollidIonize(
	/* element number on c scale, H is 0 */
	long int nelem )
{
#if !defined(NDEBUG)
	int lgInsane;
#endif
	long int nion;
	float DimaRate;
	double crate;

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

	/* compute collisional ionization rate */

	/* CollidRate(nelem,ion , 1) is collisional ionization rate, s-1
	 * CollidRate(nelem,ion , 2) is collisional ionization cooling, erg/s
	 */

	for( nion=0; nion < (IonRange.IonLow[nelem] - 2); nion++ )
	{
		CollIonRate.CollidRate[0][nion][nelem] = 0.;
		CollIonRate.CollidRate[1][nion][nelem] = 0.;
	}

	for( nion=IonRange.IonLow[nelem]-1; nion < (IonRange.IonHigh[nelem] - 1); nion++ )
	{
		/* 
		 * collisional ionization by thermal electrons
		 * >>chng 97 mar 19, to Dima's new routine using
		 * >>refer Voronov G.S., 1997, At. Data Nucl. Data Tables 65, 1
		 */
		cfit( nelem+1, nelem+1-nion , phycon.te , &DimaRate );

		crate = DimaRate*phycon.eden;

		/* total collisional ionization rate 
		 * with both thermal and suprathermal electrons */
		CollIonRate.CollidRate[0][nion][nelem] = (float)crate + Secondaries.csupra;

		/* cooling due to collisional ionization, which only includes thermal */
		CollIonRate.CollidRate[1][nion][nelem] = (float)(crate*
			rfield.anu[Heavy.ipHeavy[nion][nelem]-1]* EN1RYD);
	}

	for( nion=IonRange.IonHigh[nelem]-1; nion <= nelem; nion++ )
	{
		CollIonRate.CollidRate[0][nion][nelem] = 0.;
		CollIonRate.CollidRate[1][nion][nelem] = 0.;
	}

	/* only do this test in debug mode */
#if !defined(NDEBUG)
	/* check not rates are negative */
	lgInsane = FALSE;
	for( nion=0; nion <= nelem; nion++ )
	{
		/* there can be no negative rates */
		if( CollIonRate.CollidRate[0][nion][nelem] < 0. )
		{
			fprintf( ioQQQ, " CollidIonize finds insane coll ionize rate%10.2e nelem, ion=%3ld%3ld csupra=%10.2e te=%12.4e\n", 
			  CollIonRate.CollidRate[0][nion][nelem], nelem, nion, Secondaries.csupra, 
			  phycon.te );
			lgInsane = TRUE;
		}
	}

	if( lgInsane )
	{
		fprintf( ioQQQ, " Insanity has occurred.\n" );
		ShowMe();
		puts( "[Stop in collidionize]" );
		exit(1);
	}
#	endif

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