/*IonNeon ionization balance for neon */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#ifdef NELEM
#	undef NELEM
#endif
#define	NDIM	11
#define	NELEM	10
#include "elmton.h"
#include "trace.h"
#include "timed.h"
#include "charexc.h"
#include "ionfracs.h"
#include "collidionize.h"
#include "ionzer.h"
#include "makerecomb.h"
#include "theavy.h"
#include "photoionize.h"
#include "bidiag.h"
#include "ionheavy.h"

void IonNeon(void)
{
	long int i, 
	  _r;
	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static double rec[NDIM - 2]={3.07e-10,8.87e-10,2.21e-9,4.28e-9,
	  8.17e-9,1.51e-8,3.21e-8,9.29e-8,1.02e-7};
	static double pl[NDIM - 2]={-0.759,-0.693,-0.675,-0.668,-0.684,
	  -0.704,-0.742,-0.803,-0.769};
	static double tlow[2]={0.,0.};
	static double aa[NDIM - 1]={0.,0.0129,3.6781,-0.0254,-0.0141,19.9280,
	  5.4751,0.,0.,0.};
	static double bb[NDIM - 1]={0.,-0.1779,14.1481,5.5365,33.8479,235.0536,
	  203.9751,0.,0.,0.};
	static double cc[NDIM - 1]={0.,0.9353,17.1175,17.0727,43.1608,152.5096,
	  86.9016,0.,0.,0.};
	static double dd[NDIM - 1]={0.,-0.0682,-0.5017,-0.7225,-1.6072,
	  9.1413,-7.4568,0.,0.,0.};
	static double ff[NDIM - 1]={0.1,0.4516,0.2313,0.1702,0.1942,0.1282,
	  2.5145,0.,0.,0.};
	static double ditcrt[NDIM - 1]={3.0e4,3.3e4,3.3e4,3.5e4,3.6e4,3.6e4,
	  2.9e4,1.5e6,3.8e6,1e20};
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp2[] = {9.77e-4,2.65e-3,3.69e-3,1.12e-2,
		  2.44e-2,3.02e-2,6.10e-3,2.52e-1,1.44e-1,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[0][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {.073,.242,1.01,.391,2.52,.445,.254,
		  .304,.296,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[1][i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {3.11e5,2.84e5,2.24e5,2.7e5,3.09e5,
		  2.83e5,1.68e5,1.4e7,1.5e7,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[0][i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {2.06e5,3.07e5,2.94e5,5.50e5,9.91e5,
		  1.73e6,6.13e5,1.80e6,2.24e6,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[1][i-1] = _itmp5[_r++];
			}
		}
		_aini = 0;
	}

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

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

	/* rec from +7, +8 fro Arnaud et al 85 */
	/* Pequignot and Aldrovandi Ast Ap 161, 169. */

	/* neon, atomic number 10 */
	if( !elmton.lgElmtOn[9] )
	{
		
#	ifdef DEBUG_FUN
		fputs( " <->IonNeon()\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, first 2 from Butler, Heil, Dalgarno
	 * Ne+2 => Ne+
	 * >>chng 96 Oct 1, all H ct moved to bidiag
	 * CTHeavy(2,2) = CTHeavy(2,2) + 1E-14 * HI
	 * Ne+3 => Ne+2
	 * CTHeavy(3,2) = CTHeavy(3,2) + 5.7E-11*SQRTE*HI
	 * Ne+4 => Ne+3, Butler and Dalgarno 1980B
	 * CTHeavy(4,2) = CTHeavy(4,2) + 6.6E-9*HI
	 * */
	if( timed.itime == -1 )
	{
		theavy(NELEM);
		
#	ifdef DEBUG_FUN
		fputs( " <->IonNeon()\n", debug_fp );
#	endif
		return;
	}

	if( trace.lgTrace && trace.lgNeonBug )
	{
		fprintf( ioQQQ, " IonNeon;" );
		for( i=1; i <= 10; i++ )
		{
			fprintf( ioQQQ, "%10.2e", CharExc.CTHeavy[0][i-1] );
		}
		for( i=1; i <= 10; i++ )
		{
			fprintf( ioQQQ, "%10.2e", CharExc.CTHeavy[1][i-1] );
		}
		fprintf( ioQQQ, "\n" );
	}

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

	if( trace.lgTrace && trace.lgHeavyBug )
	{
		fprintf( ioQQQ, "     IonNeon   returns; frac=" );
		for( i=1; i <= 10; i++ )
		{
			fprintf( ioQQQ, "%10.3e", xIonFracs[i][NELEM-1]/
			  xIonFracs[0][NELEM-1] );
		}
		fprintf( ioQQQ, "\n" );
	}


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