/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*ion_collis fill in collisional ionization rates, and resulting cooling */
#include "cddefines.h"
#include "physconst.h"
#include "phycon.h"
#include "rfield.h"
#include "heavy.h"
#include "dense.h"
#include "atmdat.h"
#include "ionbal.h"

void ion_collis(
	/* element number on c scale, H is 0 */
	long int nelem )
{
	long int ion,
		limit;
	float DimaRate;
	double crate;

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

	/* compute collisional ionization rate */

	/* CollidRate[nelem][ion][0] is collisional ionization rate, s-1
	 * CollidRate[nelem][ion][1] is collisional ionization cooling, erg/s
	 */

	/* zero out rates below lowest ionization stage we will consider */
	for( ion=0; ion < (dense.IonLow[nelem] - 1); ion++ )
	{
		ionbal.CollIonRate_Ground[nelem][ion][0] = 0.;
		ionbal.CollIonRate_Ground[nelem][ion][1] = 0.;
	}

	/* chng logic changed to be precisely same as bidiag */
	/* >>chng 02 nov 08, change 2 to NISO */
	/*limit = MIN2(nelem-2,dense.IonHigh[nelem]-1);*/
	limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);

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

		crate = DimaRate*dense.EdenHCorr;

		/* total collisional ionization rate 
		 * with only thermal suprathermal electrons */
		ionbal.CollIonRate_Ground[nelem][ion][0] = crate ;

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

	for( ion=dense.IonHigh[nelem]; ion <= nelem; ion++ )
	{
		ionbal.CollIonRate_Ground[nelem][ion][0] = 0.;
		ionbal.CollIonRate_Ground[nelem][ion][1] = 0.;
	}

	/* check not rates are negative - in release mode this loop will optimize out */
	for( ion=0; ion <= nelem; ion++ )
	{
		/* there can be no negative rates */
		ASSERT( ionbal.CollIonRate_Ground[nelem][ion][0] >= 0. );
	}

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