/*IonPhosi derive ionization balance for phosphorus */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#ifdef NELEM
#	undef NELEM
#endif
#define	NDIM	16
#define	NELEM	15
#include "elmton.h"
#include "timed.h"
#include "collidionize.h"
#include "ionzer.h"
#include "makerecomb.h"
#include "theavy.h"
#include "photoionize.h"
#include "bidiag.h"
#include "ionheavy.h"

void IonPhosi(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.};
	static double bb[NDIM - 1]={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.};
	static double dd[NDIM - 1]={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.};
	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};
	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};
	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,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,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,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,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,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[1][i-1] = _itmp5[_r++];
			}
		}
		_aini = 0;
	}

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

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

	/* phosphorus, atomic number 15 */
	if( !elmton.lgElmtOn[NELEM-1] )
	{
		
#	ifdef DEBUG_FUN
		fputs( " <->IonPhosi()\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( " <->IonPhosi()\n", debug_fp );
#	endif
		return;
	}

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

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