/*IonSulph compute ionization balance for sulphur */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#ifdef NELEM
#	undef NELEM
#endif
#define	NDIM	17
#define	NELEM	16
#include "elmton.h"
#include "he.h"
#include "timed.h"
#include "charexc.h"
#include "tepowers.h"
#include "collidionize.h"
#include "ionzer.h"
#include "makerecomb.h"
#include "theavy.h"
#include "photoionize.h"
#include "bidiag.h"
#include "ionheavy.h"

void IonSulph(void)
{
	long int i, 
	  _r;
	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 rec[NDIM - 2]={1.54e-10,13.4e-10,3.37e-9,7.36e-9,
	  7.64e-9,4.23e-8,2.31e-8,2.43e-8,3.84e-8,5.41e-8,6.60e-8,1.47e-7,
	  3.35e-7,6.67e-7,7.27e-7};
	static double pl[NDIM - 2]={-0.63,-0.686,-0.745,-0.755,-0.701,-0.849,
	  -0.733,-0.696,-0.711,-0.716,-0.714,-0.755,-0.806,-0.840,-0.807};
	static double tlow[2]={1.23e-10,-0.593};
	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 nelem=16
	 *
	 * rates from Shull and van Steenberg, Ap.J. Sup 48, 95. */

	/* rec from +13 14 15 from Arnauld et al 85  */
	/* Pequignot and Aldrovandi Ast Ap 161, 169. */

	/* sulphur, atomic number 16 */
	if( !elmton.lgElmtOn[15] )
	{
		
#	ifdef DEBUG_FUN
		fputs( " <->IonSulph()\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, all rates from Butler+Dalgarno 1980B
	 * S+2 => S+
	 * >>chng 96 Oct 1, moved all H ct to bidiag
	 * CTHeavy(2,2) = CTHeavy(2,2) + 1E-14 * HI
	 * S+3 => S+2
	 * CTHeavy(3,2) = CTHeavy(3,2) + 2.5E-9 * HI + HEI*SQRTE * 1.1E-11 */
	CharExc.CTHeavy[1][2] += (float)(he.hei*TePowers.sqrte*1.1e-11);
	/* S+4 => S+3
	 * CTHeavy(4,2) = CTHeavy(4,2) + 7.5E-9 * HI
	 * */
	if( timed.itime == -1 )
	{
		theavy(NELEM);
		
#	ifdef DEBUG_FUN
		fputs( " <->IonSulph()\n", debug_fp );
#	endif
		return;
	}

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

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