/*IonArgon compute ionization balance of argon */
#include "cddefines.h"
#	ifdef NDIM
#undef NDIM
#	endif
#	ifdef NELEM
#undef NELEM
#	endif
#define	NDIM	19
#define	NELEM	18
#include "elmton.h"
#include "he.h"
#include "timed.h"
#include "charexc.h"
#include "collidionize.h"
#include "makerecomb.h"
#include "ionzer.h"
#include "theavy.h"
#include "photoionize.h"
#include "bidiag.h"
#include "ionheavy.h"

void IonArgon(void)
{
	long int i, 
	  _r;
	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static double rec[NDIM - 2]={3.67e-11,1.16e-9,4.1e-9,6.3e-9,8.5e-9,
	  1.15e-8,1.42e-8,6.71e-8,4.99e-8,6.5e-8,5.7e-8,5.1e-8,7.6e-8,
	  8.3e-8,5.80e-7,1.14e-6,9.45e-7};
	static double pl[NDIM - 2]={-0.497,-0.694,-0.778,-0.756,-0.724,
	  -0.706,-0.676,-0.830,-0.750,-0.755,-0.707,-0.647,-0.694,-0.682,
	  -0.819,-0.854,-0.803};
	static double tlow[2]={0.,0.};
	static double ditcrt[NDIM - 1]={2.5e4,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.,0.};
	static double bb[NDIM - 1]={0.,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.,0.};
	static double dd[NDIM - 1]={0.,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.,0.};
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp2[] = {1.00e-3,1.10e-2,3.40e-2,6.85e-2,
		  9.00e-2,6.35e-2,2.60e-2,1.70e-2,2.10e-2,3.50e-2,4.30e-2,
		  7.13e-2,9.60e-2,8.50e-2,1.70e-2,.476,.297,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[0][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {.005,.045,.057,.087,.0769,.140,
		  .120,.1,1.92,1.66,1.67,1.40,1.31,1.02,.245,.294,.277,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[1][i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {3.20e5,2.90e5,2.39e5,2.56e5,2.50e5,
		  2.10e5,1.80e5,2.70e6,8.30e5,6.95e5,6.05e5,6.68e5,6.50e5,
		  5.30e5,3.55e5,3.01e7,3.13e7,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[0][i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {3.10e5,5.50e5,6.00e5,3.81e5,3.30e5,
		  2.15e5,2.15e5,3.30e6,3.50e6,3.60e6,3.80e6,2.90e6,3.60e6,
		  2.80e6,1.10e6,6.05e6,6.54e6,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[1][i-1] = _itmp5[_r++];
			}
		}
		_aini = 0;
	}

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

	/* argon nelem=18
	 *
	 * rates from Shull and van Steenberg, Ap.J. Sup 48, 95. */

	/* rec from +15, 16, 17 from Arnauld et al 85 */
	/* Pequignot and Aldrovandi Ast Ap 161, 169. */

	/* argon, atomic number 18 */
	if( !elmton.lgElmtOn[17] )
	{
		
#	ifdef DEBUG_FUN
		fputs( " <->IonArgon()\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);

	/* CHARGE EXCHANGE, ALL RATES FROM BUTLER AND DALGARNO 1980B
	 * A+2 => A+
	 * CTHeavy(2,2) = CTHeavy(2,2) + 1E-14 * HI + HEI * 1.3E-10
	 * >>chng 96 Oct 1, move H ct into bidiag */
	CharExc.CTHeavy[1][1] += (float)(he.hei*1.3e-10);
	/* A+ => A+2 */
	CharExc.CTHeavy[0][1] += (float)(he.heii*1.1e-10);
	/* T = MIN( SQRTE , 200. )
	 * A+3 => A+2
	 * >>chng 96 Oct 1, move H ct into bidiag
	 * CTHeavy(3,2) = CTHeavy(3,2) + 4.4E-11*MIN(sqrte,200.)*hi
	 * >>chng 96 Oct 1, move H ct into bidiag
	 * CTHeavy(4,2) = CTHeavy(4,2) + 6.3E-9 * hi
	 * */
	if( timed.itime == -1 )
	{
		theavy(NELEM);
		
#	ifdef DEBUG_FUN
		fputs( " <->IonArgon()\n", debug_fp );
#	endif
		return;
	}

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


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