/* 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 */
/*IonSulph compute ionization balance for sulphur */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#define	NDIM	17
#include "dense.h"
#include "ionbal.h"
#include "co.h"
#include "converge.h"

void IonSulph(void)
{
	int lgPrtDebug=FALSE;
	long int i, 
	  _r;
	double save_rec;
	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static double aa[NDIM - 1]={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.};
	static double cc[NDIM - 1]={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.};
	static double ff[NDIM - 1]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.};
	static double ditcrt[NDIM - 1]={2.2e4,1.2e4,1.4e4,1.5e4,1.4e4,2.9e5,
	  1.3e5,1.1e5,9.0e4,9.0e4,9.0e4,8.3e4,6.0e4,5.0e6,9.0e6,1e20};
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp2[] = {1.62e-3,1.09e-2,3.35e-2,3.14e-2,
		  1.27e-2,1.47e-2,1.34e-2,2.38e-2,3.19e-2,7.13e-2,.08,.0796,
		  .0134,.402,.241,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[0][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {0.,1.20e-2,6.59e-2,6.89e-2,.187,
		  .129,1.04,1.12,1.40,1.00,.555,1.63,.304,.298,.281,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[1][i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {1.25e5,1.92e5,1.89e5,1.68e5,1.38e5,
		  1.80e6,6.90e5,5.84e5,5.17e5,6.66e5,6.00e5,5.09e5,2.91e5,
		  2.41e7,2.54e7,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[0][i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {0.,1.80e4,1.59e5,8.04e4,1.71e5,
		  1.75e6,2.15e6,2.59e6,2.91e6,2.32e6,2.41e6,6.37e6,1.04e6,
		  4.67e6,5.30e6,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[1][i-1] = _itmp5[_r++];
			}
		}
		_aini = 0;
	}

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

	/* sulphur, atomic number 16 */
	if( !dense.lgElmtOn[ipSULPHUR] )
	{
		
#		ifdef DEBUG_FUN
		fputs( " <->IonSulph()\n", debug_fp );
#		endif
		return;
	}

	ion_zero(ipSULPHUR);

	ion_photo(ipSULPHUR,lgPrtDebug);

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

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

	/* >>chng 03 sep 29, synch up ion and co solvers.
	 * the co solver will have a different S/S+ balance that
	 * is strongly affected by the chemistry if we are in neutral gas
	 * (hence the test that the highest stage of ionization is <=2 on 
	 * physics scale)
	 * in this case use co solver's S/S+ balance */
	save_rec = ionbal.RateRecomTot[ipSULPHUR][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[ipSULPHUR]==0  && /*!conv.lgSearch &&*/
		co.hevmol[ipSP] > SMALLFLOAT && 
		ionbal.RateIonizTot[ipSULPHUR][0]*co.hevmol[ipATS]> 0. )
	{
		ionbal.RateRecomTot[ipSULPHUR][0] = 
			ionbal.RateIonizTot[ipSULPHUR][0]*
			co.hevmol[ipATS]/co.hevmol[ipSP];
	}

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

	/* reset the var we just hosed */
	/* if ion just trimed down, old rec may be zero, must not use it*/
	if( save_rec > 0. )
		ionbal.RateRecomTot[ipSULPHUR][0] = save_rec;

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

	return;
}

