/*IonCalci perform ionization balance for calcium */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#ifdef NELEM
#	undef NELEM
#endif
#define	NDIM	21
#define	NELEM	20
#include "elmton.h"
#include "timed.h"
#include "charexc.h"
#include "phycon.h"
#include "dcala.h"
#include "collidionize.h"
#include "makerecomb.h"
#include "recomrate.h"
#include "ionzer.h"
#include "theavy.h"
#include "photoionize.h"
#include "bidiag.h"
#include "ionheavy.h"

void IonCalci(void)
{
	long int i, 
	  _r;
	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static double rec[NDIM - 2]={4.5e-10,1.07e-9,2.50e-9,9.33e-9,2.45e-8,
	  3.43e-8,4.57e-8,6.53e-8,6.64e-8,2.57e-7,1.62e-7,2.05e-7,2.50e-7,
	  1.70e-7,4.12e-7,4.27e-7,9.57e-7,1.29e-6,1.35e-6};
	static double pl[NDIM - 2]={-0.900,-0.800,-0.700,-0.780,-0.840,
	  -0.820,-0.820,-0.810,-0.780,-0.900,-0.820,-0.810,-0.800,-0.730,
	  -0.800,-0.780,-0.850,-0.850,-0.830};
	static double tlow[2]={1.06e-10,-0.683};
	static double ditcrt[NDIM - 1]={6e3,2e4,4e4,5e4,7e4,8e4,8e4,3e4,
	  3e4,3e4,3e4,9e4,4e4,5e4,3e4,9e5,2e5,2e5,2e5,1e20};
	static double aa[NDIM - 1]={0.,0.,0.,0.,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.,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.,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.,0.,0.,0.,0.};
	static double ff[NDIM - 1]={0.,0.1,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.,0.,0.,0.};
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp2[] = {3.28e-4,.0584,.122,.132,.133,.126,
		  .139,.0955,.0402,.0419,.0257,.0445,.0548,.0713,.0903,.110,
		  .0205,.549,.355,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[0][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {.0907,.110,.0174,.132,.114,.162,
		  .0878,.263,.0627,.0616,2.77,2.23,2.00,1.82,.424,.243,.185,
		  .292,.275,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[1][i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {3.46e4,3.84e5,4.08e5,3.82e5,3.53e5,
		  3.19e5,3.22e5,2.47e5,2.29e5,3.73e6,9.26e5,7.96e5,6.90e5,
		  6.70e5,4.72e5,5.67e5,4.21e5,3.65e7,3.78e7,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[0][i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {1.64e4,2.45e5,4.27e5,6.92e5,8.78e5,
		  7.43e5,6.99e5,4.43e5,2.81e5,5.84e6,4.89e6,4.62e6,4.52e6,
		  3.32e6,1.37e6,4.41e6,2.27e6,7.25e6,7.68e6,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[1][i-1] = _itmp5[_r++];
			}
		}
		_aini = 0;
	}

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

	/* calcium nelem=20
	 *
	 * rates from Shull and van Steenberg, Ap.J. Sup 48, 95. */

	/* rates from Shull and van Steenberg, Ap.J. Sup 48, 95. */
	/* Pequignot and Aldrovandi Ast Ap 161, 169. */

	/* calcium, atomic number 20 */
	if( !elmton.lgElmtOn[19] )
	{
#		ifdef DEBUG_FUN
		fputs( " <->IonCalci()\n", debug_fp );
#		endif
		return;
	}

	ionzer(NELEM);

	/* last par is option to print info about photo rates */
	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);

	/* correct low temp rad rec coef */
	if( phycon.te < 1e3 )
	{
		RecomRate.RecombinRate[1][NELEM-1] += (float)((5.49e-10*(pow(phycon.te,-0.647)) - 
		  rec[1]*(pow(phycon.te,pl[1])))*phycon.eden);
	}

	/* Ly-alpha photoionization of Ca+ 
	 * this array was set to zero in ionzer call above, so this only happens
	 * one time per setup */
	CharExc.CTHeavy[0][1] += dcala.dstCala;

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

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

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