/*IonAlumi ionization balance for aluminum */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#ifdef NELEM
#	undef NELEM
#endif
#define	NDIM	14
#define	NELEM	13
#include "trace.h"
#include "elmton.h"
#include "timed.h"
#include "ionfracs.h"
#include "makerecomb.h"
#include "collidionize.h"
#include "ionzer.h"
#include "theavy.h"
#include "photoionize.h"
#include "bidiag.h"
#include "ionheavy.h"

void IonAlumi(void)
{
	long int i, 
	  _r;
	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static double rec[NDIM - 2]={2.14e-10,1.05e-9,6.49e-9,1.09e-8,1.31e-8,
	  2.54e-8,3.13e-8,8.85e-8,1.22e-7,7.15e-7,5.54e-7,4.12e-7};
	static double pl[NDIM - 2]={-0.659,-0.713,-0.840,-0.784,-0.742,
	  -0.752,-0.733,-0.789,-0.790,-0.927,-0.883,-0.809};
	static double tlow[2]={1.13e-10,-0.567};
	static double ditcrt[NDIM - 1]={8.0e3,7.0e3,1.2e5,7.0e4,6.5e4,5.5e4,
	  5.5e4,5.5e4,5.8e4,4.0e4,2.9e6,5.5e6,1e20};
	static double aa[NDIM - 1]={0.0219,0.7086,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.};
	static double bb[NDIM - 1]={-0.4528,-3.1083,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.};
	static double cc[NDIM - 1]={2.5427,7.0422,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.};
	static double dd[NDIM - 1]={-0.1678,0.5998,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.};
	static double ff[NDIM - 1]={0.2276,0.4194,0.1,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-3,6.0e-3,7.5e-3,5.6e-3,1.2e-2,
		  1.9e-2,1.6e-2,1.7e-2,4.4e-2,1.7e-2,3.0e-1,3.4e-0,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[0][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {9.0e4,7.5e4,8.5e5,5.0e5,4.4e5,3.8e5,
		  3.4e5,3.4e5,3.9e5,2.3e5,1.6e7,1.7e7,1e20};
		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( "<+>IonAlumi()\n", debug_fp );
#	endif

	/* aluminium nelem=13*/

	/* rates from Shull and van Steenberg, Ap.J. Sup 48, 95.
	 * use silicon since not in shull + van steenberg */

	/* rec coef from Aldrovandi and Pequignot Revista Bras de Fisica 4, 491. */
	/* Pequignot and Aldrovandi Ast Ap 161, 169. */


	/* from Nussbaumer and Storey, Mg, Al, Si: n.b., al1 uses table 2(b) */

	/* aluminium, atomic number 13 */
	if( !elmton.lgElmtOn[12] )
	{
		
#	ifdef DEBUG_FUN
		fputs( " <->IonAlumi()\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( " <->IonAlumi()\n", debug_fp );
#	endif
		return;
	}

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

	if( trace.lgTrace && trace.lgHeavyBug )
	{
		fprintf( ioQQQ, "     IonAlumi 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( " <->IonAlumi()\n", debug_fp );
#	endif
	return;
}