/*IonSilic determine ionization balance of Silicon */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#ifdef NELEM
#	undef NELEM
#endif
#define	NDIM	15
#define	NELEM	14
#include "trace.h"
#include "elmton.h"
#include "he.h"
#include "pmp2s.h"
#include "timed.h"
#include "tepowers.h"
#include "charexc.h"
#include "photrate.h"
#include "zonecnt.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 IonSilic(void)
{
	long int i, 
	  _r;
	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static long nzUsed = -1;
	static double OldRate = 0.;
	static double rec[NDIM - 2]={1.64e-10,2.00e-9,2.19e-9,1.41e-8,1.05e-8,
	  1.54e-8,1.93e-8,2.43e-8,3.76e-8,5.53e-8,2.80e-7,5.81e-8,5.99e-7};
	static double pl[NDIM - 2]={-0.601,-0.786,-0.693,-0.821,-0.735,
	  -0.716,-0.703,-0.688,-0.703,-0.714,-0.823,-0.858,-0.818};
	static double tlow[2]={1.64e-10,-0.617};
	static double ditcrt[NDIM - 1]={1.1e4,1.1e4,1.1e4,1.7e5,9.5e4,8.0e4,
	  7.4e4,6.8e4,6.6e4,6.5e4,4.5e4,3.7e6,6.3e6,1e20};
	static double aa[NDIM - 1]={-0.0219,3.2163,0.1203,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.,0.};
	static double bb[NDIM - 1]={0.4364,-12.0571,-2.6900,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.,0.,0.};
	static double cc[NDIM - 1]={0.0684,16.2118,19.1943,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.,0.,0.};
	static double dd[NDIM - 1]={-0.0032,-0.5886,-0.1479,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.,0.,0.};
	static double ff[NDIM - 1]={0.1342,0.5613,0.1118,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[] = {1.10e-3,5.87e-3,5.03e-3,5.43e-3,
		  8.86e-3,1.68e-2,2.49e-2,3.13e-2,4.25e-2,6.18e-2,1.38e-2,
		  .327,.189,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[0][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {0.,.753,.188,.450,0.,1.80,1.88,
		  2.01,1.22,.303,1.42,.306,.286,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[1][i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {7.70e4,9.63e4,8.75e4,1.05e6,1.14e6,
		  4.85e5,4.15e5,3.66e5,3.63e5,3.88e5,2.51e5,1.88e7,1.99e7,
		  0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[0][i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {0.,6.46e4,4.71e4,7.98e5,0.,1.03e6,
		  1.91e6,2.11e6,2.14e6,1.12e6,3.93e6,3.60e6,4.14e6,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[1][i-1] = _itmp5[_r++];
			}
		}
		_aini = 0;
	}

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

	/* silicon nelem=14
	 *
	 * rec fro +11, 12 13 from Arnauld et al 85 */
	/* Pequignot and Aldrovandi Ast Ap 161, 169., scaled to smoothly fit
	 * onto high temp coef */


	/* from nuss+storey, who say that =>si iv is very small */

	/* silicon, atomic number 14 */
	if( !elmton.lgElmtOn[13] )
	{
		pmp2s.p1895 = 0.;
		
#	ifdef DEBUG_FUN
		fputs( " <->IonSilic()\n", debug_fp );
#	endif
		return;
	}

	ionzer(NELEM);

	PhotoIonize(NELEM-1,FALSE);
	if( ZoneCnt.nzone > 1 && OldRate > 0. )
	{
		if( ZoneCnt.nzone != nzUsed )
		{
			PhotRate.PhotoRate[0][4][0][NELEM-1] = 
				(PhotRate.PhotoRate[0][4][0][NELEM-1] + OldRate)/2.;
			nzUsed = ZoneCnt.nzone;
		}
		else
		{
			PhotRate.PhotoRate[0][4][0][NELEM-1] = OldRate;
		}
	}
	OldRate = PhotRate.PhotoRate[0][4][0][NELEM-1];

	/* 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);

	/* pumping to 3p of Si III */
	pmp2s.p1895 = PhotRate.PhotoRate[0][2][1][NELEM-1];

	/* SI ++ - SI+ CH EXCH FROM GARGAUD ET AL., 82, AST AP 106, 197
	 * RATE = 2.53E-10 * TE/TE70*TE03
	 * >>chng 96 Oct 1, moved all H ct to bidiag
	 * CTHeavy(2,2) = CTHeavy(2,2) + HI * RATE
	 * CTHeavy(2,1) = CTHeavy(2,1) + HII*RATE*0.262 * SEXP(2.943E4/TE)
	 * SI+3 => SI+2 BUTLER AND DALGARNO 1980B
	 * CTHeavy(3,2) = CTHeavy(3,2) + 4E-10 * HI + HEI*SQRTE*TE10*TE10*
	 * >>refer Butler, S.E., &Dalgarno, A. 1980, ApJ 241, 838
	 * scale by 1.3 to bring into agreement with
	 * >>refer Fang, Z., & Kwong, H.S. 1997 ApJ 483, 527 */
	CharExc.CTHeavy[1][2] += (float)(he.hei*TePowers.sqrte*TePowers.te10*TePowers.te10*
	  1.3*1.5e-12);
	/* SI+4 => SI+3 BUTLER AND DALGARNO 1980B
	 * CTHeavy(4,2) = CTHeavy(4,2) + 2.3E-9 * HI
	 * Si+4 + Heo => Si+3 (Opradolce et al. 1985, AstAp. 148, 229.) */
	CharExc.CTHeavy[1][3] += (float)(2.54e-11*TePowers.sqrte/TePowers.te03/
	  TePowers.te01/TePowers.te01*he.hei);

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

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

	/* get pmp rate cm-3 s-1 */
	pmp2s.p1895 *= xIonFracs[2][NELEM-1]*0.85;

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