/*IonCarbo compute ionization balance for carbon */
#include "cddefines.h"
#	ifdef NDIM
#undef NDIM
#	endif
#	ifdef NELEM
#undef NELEM
#	endif
#define	NDIM	7
#define	NELEM	6
#include "heating.h"
#include "pmp2s.h"
#include "charexc.h"
#include "timed.h"
#include "trace.h"
#include "elmton.h"
#include "tepowers.h"
#include "hevmolec.h"
#include "phycon.h"
#include "coatom.h"
#include "he.h"
#include "photrate.h"
#include "ionfracs.h"
#include "collidionize.h"
#include "makerecomb.h"
#include "ionzer.h"
#include "theavy.h"
#include "photoionize.h"
#include "bidiag.h"
#include "ionheavy.h"

void IonCarbo(void)
{
	long int i, 
	  _r;
	float save;
	static double 
	  dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static double rec[NDIM - 2]={1.47e-10,8.74e-10,7.98e-9,1.34e-8,
	  1.30e-8};
	static double pl[NDIM - 2]={-0.624,-0.645,-0.803,-0.791,-0.721};
	static double tlow[2]={1.44e-10,-0.621};
	static double ditcrt[NDIM - 1]={1.2e4,1.2e4,1.1e4,4.4e5,7.0e5,1e20};
	static double aa[NDIM - 1]={.0108,1.8267,2.3196,0.,0.,0.};
	static double bb[NDIM - 1]={-0.1075,4.1012,10.7328,0.,0.,0.};
	static double cc[NDIM - 1]={.2810,4.8443,6.8830,0.,0.,0.};
	static double dd[NDIM - 1]={-0.0193,.2261,-0.1824,0.,0.,0.};
	static double ff[NDIM - 1]={0.000,0.5960,0.4101,0.1,0.1,0.};
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp2[] = {2.54e-3,6.15e-3,1.62e-3,4.78e-2,
		  3.22e-2,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[0][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {4.42e-2,5.88e-2,0.343,0.362,0.315,
		  0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[1][i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {1.57e5,1.41e5,8.19e4,3.44e6,4.06e6,
		  0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[0][i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {3.74e5,1.41e5,1.59e5,5.87e5,8.31e5,
		  0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[1][i-1] = _itmp5[_r++];
			}
		}
		_aini = 0;
	}

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

	/* carbon nelem=6
	 *
	 * real CollidRate(nelem,2)
	 * */

	/* rates from Shull and van Steenberg, Ap.J. Sup 48, 95. */

	/* rec from +3, +4 from Arnaud et al Ast Ap Sup 60 425. (1985)
	 * rec from fully ionized uses Seaton '79 in ionrat */
	/* Pequignot and Aldrovandi Ast Ap 161, 169. */


	/* in above, ff(1) has been reset to 0 since NS fit diverges at very low t
	 * better to use mean background diel
	 * */
	if( trace.lgTrace && trace.lgCarBug )
		{
		fprintf( ioQQQ, "     IonCarbo called.\n" );
		}

	if( !elmton.lgElmtOn[5] )
	{
		pmp2s.p1909 = 0.;
		pmp2s.p2326 = 0.;
		coatom.catmic = 0.;
		HeatingCom.heating[9][5] = 0.;
		
#	ifdef DEBUG_FUN
		fputs( " <->IonCarbo()\n", debug_fp );
#	endif
		return;
	}

	/* zero out ionization balance arrays */
	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);

	/* Photoproduction of 3P */
	pmp2s.p1909 = PhotRate.PhotoRate[0][1][1][5];

	/* photoproduction of C II] 2326 */
	pmp2s.p2326 = PhotRate.PhotoRate[0][1][0][5];

	CharExc.CTHeavy[0][0] += (float)(he.heii*4.17e-17*(phycon.te/TePowers.te03));

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

	/* following accounts for part of carbon that is in CO, since BiDiag ony
	 * wants atoms/ions, but [0][7] includes molecules too.  Must temperariyly
	 * convert [0][5] to atoms/molecules only */
	coatom.catmic = (float)(xIonFracs[0][5]*MAX2(1e-7,1.-hevmolec.hevmol[IPCO-1]/
	  xIonFracs[0][5]));

	/* solve for ionization balance
	 * last var is option to print stuff */
	save = xIonFracs[0][5];

	/* must trick BiDiag into conserving C */
	xIonFracs[0][5] = coatom.catmic;

	BiDiag(NELEM-1,FALSE);

	xIonFracs[0][5] = save;

	/* rate will be cm-3 s-1, into triplets only
	 * >>chng 96 nov 22, rid of carb() */
	pmp2s.p1909 *= xIonFracs[2][5]*0.62;
	pmp2s.p2326 *= xIonFracs[1][5]*0.1;

	/* charge exchange rates (s-1) for hydrogen
	 * >>chng 96 nov 22, rid of carb() */
	CharExc.CTHionAgnt[0] = (float)(xIonFracs[2][5]*1.4e-17);
	CharExc.CTHion[0][0] += CharExc.CTHionAgnt[0];
	CharExc.CTHrec[0][0] += (float)(xIonFracs[1][5]*2.18e-20*phycon.te);

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "     IonCarbo returns; CO=%10.3e", 
		  hevmolec.hevmol[IPCO-1]/xIonFracs[0][5] );
		fprintf( ioQQQ, " fracs=" );
		for( i=1; i <= 7; i++ )
		{
			fprintf( ioQQQ, " %10.3e", xIonFracs[i][5]/
			  xIonFracs[0][5] );
		}
		fprintf( ioQQQ, "\n" );
	}


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