/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*IonSilic determine ionization balance of Silicon */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#define	NDIM	15
#include "co.h"
#include "trace.h"
#include "dense.h"
#include "sil.h"
#include "ionbal.h"

void IonSilic(void)
{
	int lgPrtDebug=FALSE;
	long int i, 
	  _r;
	double save_rec;
	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static long nzUsed = -1;
	static double OldRate = 0.;
	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, atomic number 14 */
	if( !dense.lgElmtOn[13] )
	{
		sil.p1895 = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->IonSilic()\n", debug_fp );
#		endif
		return;
	}

	ion_zero(ipSILICON);

	ion_photo(ipSILICON,FALSE);
	if( nzone > 1 && OldRate > 0. )
	{
		if( nzone != nzUsed )
		{
			ionbal.PhotoRate_Shell[ipSILICON][0][4][0] = 
				(ionbal.PhotoRate_Shell[ipSILICON][0][4][0] + OldRate)/2.;
			nzUsed = nzone;
		}
		else
		{
			ionbal.PhotoRate_Shell[ipSILICON][0][4][0] = OldRate;
		}
	}
	OldRate = ionbal.PhotoRate_Shell[ipSILICON][0][4][0];

	/* find collisional ionization rates */
	ion_collis(ipSILICON);

	/* get recombination coefficients */
	/*lint -e740 unusual pointer cast */
	ion_recomb(FALSE,(double*)dicoef,(double*)dite,ditcrt,aa,bb,cc,
	  dd,ff,ipSILICON);
	/*lint +e740 unusual pointer cast */

	/* pumping to 3p of Si III */
	sil.p1895 = ionbal.PhotoRate_Shell[ipSILICON][1][2][0];

	/* >>chng 03 sep 29, synch up ion and co solvers.
	 * the co solver will have a different C/C+ balance that
	 * is strongly affected by the chemistryif we are in neutral gas
	 * in this case use co solver's C/C+ balance */
	save_rec = ionbal.RateRecomTot[ipSILICON][0];
	/*>>chng 04 apr 27, do not test for ionhigh being 1,
	 * no reason for test on upper stage of ionization,
	 * if codrive called but not evaluted then hevmol is all zero */
	/* >>chng 04 sep 10, rm check on search phase, no reason for it */
	if( dense.IonLow[ipSILICON]==0  && /*!conv.lgSearch &&*/
		co.hevmol[ipSIP] > SMALLFLOAT && 
		ionbal.RateIonizTot[ipSILICON][0]*co.hevmol[ipATSI]>0.)
	{
		/* >>chng 04 sep 11, use damped ratio, were Si+ oscillations deep in molecular cloud,
		 * where Si+ is dominant electron donor */
		ionbal.RateRecomTot[ipSILICON][0] = 
			ionbal.RateIonizTot[ipSILICON][0]*
			co.hevmol[ipATSI]/co.hevmol[ipSIP];
		/*fprintf(ioQQQ,"DEBUG sil atm\t%.2f\t%.3e\t%.3e\t%.3e\n",
			fnzone , 
			co.hevmol[ipATSI],
			co.hevmol[ipSIP],
			ionbal.RateIonizTot[ipSILICON][0]);*/
	}

	/* solve for ionization balance */
	/*if( nzone>700 ) lgPrtDebug = TRUE;*/
	ion_solver(ipSILICON,lgPrtDebug);

	/* reset the var we just hosed */
	if( save_rec > 0. )
		ionbal.RateRecomTot[ipSILICON][0] = save_rec;	

	/* get pmp rate cm-3 s-1 */
	sil.p1895 *= dense.xIonDense[ipSILICON][1]*0.85;

	if( trace.lgTrace && trace.lgHeavyBug )
	{
		fprintf( ioQQQ, "    silicon\t%.2f",fnzone );
		for( i=0; i <= 10; i++ )
		{
			fprintf( ioQQQ, "\t%.3e", dense.xIonDense[ipSILICON][i]);
		}
		fprintf( ioQQQ, "\n" );
	}
#	if 0
	fprintf( ioQQQ, "DEBUG silicon\t%.2f",fnzone );
	for( i=0; i < 2; i++ )
	{
		fprintf( ioQQQ, "\t%.3e", dense.xIonDense[ipSILICON][i]/dense.gas_phase[ipSILICON]);
	}
	fprintf( ioQQQ, "\t%.3e", co.hevmol[ipATSI]/SDIV(co.hevmol[ipSIP]));
	fprintf( ioQQQ, "\t%.3e", ionbal.RateIonizTot[ipSILICON][0]);
	fprintf( ioQQQ, "\t%.3e", dense.gas_phase[ipHYDROGEN] );
	{
#	include "grainvar.h"
	fprintf( ioQQQ, "\t%.3e", gv.GrainChTrRate[ipSILICON][1][0] );
	}
	fprintf( ioQQQ, "\n" );
#	endif


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

	return;
}

