/*IonChlor ionization balance for chlorine */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#ifdef NELEM
#	undef NELEM
#endif
#define	NDIM	18
#define	NELEM	17
#include "elmton.h"
#include "timed.h"
#include "collidionize.h"
#include "ionzer.h"
#include "makerecomb.h"
#include "theavy.h"
#include "photoionize.h"
#include "bidiag.h"
#include "ionheavy.h"

void IonChlor(void)
{
	long int i, 
	  _r;
	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static double rec[NDIM - 2]={4.81e-10,1.08e-09,2.39e-09,4.33e-09,
	  6.59e-09,1.04e-08,1.63e-08,1.82e-08,2.27e-08,2.67e-08,3.17e-08,
	  3.73e-08,4.42e-08,5.24e-08,5.95e-08,6.43e-08};
	static double pl[NDIM - 2]={-0.702,-0.697,-0.678,-0.669,-0.659,
	  -0.681,-0.711,-0.670,-0.668,-0.646,-0.640,-0.629,-0.631,-0.647,
	  -0.641,-0.616};
	static double tlow[2]={0.,0.};
	static double ditcrt[NDIM - 1]={3.0e4,2.5e4,2.5e4,1.8e4,1.8e4,2.2e4,
	  5.0e5,1.6e5,1.5e5,1.5e5,1.3e5,1.3e5,1.1e5,7.6e4,6.5e6,1.4e7,
	  1e20};
	static double aa[NDIM - 1]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.};
	static double bb[NDIM - 1]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.};
	static double cc[NDIM - 1]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.};
	static double dd[NDIM - 1]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.};
	static double ff[NDIM - 1]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.};
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp2[] = {5.5e-5,1.0e-2,1.1e-2,1.0e-2,5.0e-2,
		  3.2e-2,3.4e-2,1.6e-2,2.4e-2,4.0e-2,4.0e-2,3.8e-2,6.8e-2,
		  2.6e-2,4.6e-1,11.,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[0][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {1.3e5,1.4e5,1.4e5,1.4e5,2.0e5,1.8e5,
		  2.3e5,7.2e5,6.4e5,6.0e5,4.9e5,4.6e5,5.3e5,3.2e5,2.6e7,2.8e7,
		  0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[0][i-1] = _itmp3[_r++];
			}
		}
		for( i=1; i <= (NDIM - 1); i++ )
		{
			dicoef[1][i-1] = 0.;
			}
		for( i=1; i <= (NDIM - 1); i++ )
		{
			dite[1][i-1] = 0.;
			}
		_aini = 0;
	}

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

	/* chlorine nelem=17
	 *
	 * rates from Shull and van Steenberg, Ap.J. Sup 48, 95.
	 * these are for argon, but no better can be done today */

	/* fits to Aldrovandi and Pequignot Rev Bras Fisica 4, 491 */

	/* Pequignot and Aldrovandi Ast Ap 161, 169. */

	/* clorine, atomic number 17 */
	if( !elmton.lgElmtOn[16] )
	{
		
#	ifdef DEBUG_FUN
		fputs( " <->IonChlor()\n", debug_fp );
#	endif
		return;
	}

	ionzer(NELEM);

	PhotoIonize(NELEM-1,FALSE);

	/* find collisional ionization rates */
	CollidIonize(NELEM-1);

	/* get recombination coefficients */
	MakeRecomb(rec,pl,(double*)dicoef,(double*)dite,ditcrt,aa,bb,cc,
	  dd,ff,NELEM-1,tlow);

	if( timed.itime == -1 )
	{
		theavy(NELEM);
		
#	ifdef DEBUG_FUN
		fputs( " <->IonChlor()\n", debug_fp );
#	endif
		return;
	}

	/* solve for ionization balance */
	BiDiag(NELEM-1,FALSE);

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