/* 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 */
/*IonCarbo compute ionization balance for carbon */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#define	NDIM	7
#include "thermal.h"
#include "carb.h"
#include "converge.h"
#include "co.h"
#include "trace.h"
#include "dense.h"
#include "phycon.h"
#include "punch.h"
#include "tfidle.h"
#include "ionbal.h"

void IonCarbo(void)
{
	long int i, 
	  _r;
	double save_rec;
	static double 
	  dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	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

	if( trace.lgTrace && trace.lgCarBug )
	{
		fprintf( ioQQQ, "     IonCarbo called.\n" );
	}

	if( !dense.lgElmtOn[ipCARBON] )
	{
		carb.p1909 = 0.;
		carb.p2326 = 0.;
		thermal.heating[ipCARBON][9] = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->IonCarbo()\n", debug_fp );
#		endif
		return;
	}

	/* zero out ionization balance arrays */
	ion_zero(ipCARBON);

	ion_photo(ipCARBON,FALSE);

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

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

	/* Photoproduction of 3P */
	carb.p1909 = ionbal.PhotoRate_Shell[ipCARBON][1][1][0];

	/* photoproduction of C II] 2326 */
	carb.p2326 = ionbal.PhotoRate_Shell[ipCARBON][0][1][0];

	/* >>chng 03 sep 29, synch up ion and co solvers.
	 * the co solver will have a different C/C+ balance that
	 * is strongly affected by the chemistry if we are in neutral gas
	 * in this case use co solver's C/C+ balance */
	save_rec = ionbal.RateRecomTot[ipCARBON][0];
	/*>>chng 04 apr 24, 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[ipCARBON]==0 && /*!conv.lgSearch &&*/
		co.hevmol[ipCP] > SMALLFLOAT && ionbal.RateIonizTot[ipCARBON][0]*co.hevmol[ipATC]>0. )
	{
		ionbal.RateRecomTot[ipCARBON][0] = 
			ionbal.RateIonizTot[ipCARBON][0]*
			co.hevmol[ipATC]/co.hevmol[ipCP];
	}

	ion_solver(ipCARBON,FALSE);
	/*fprintf(ioQQQ,"DEBUG IonCarbo\t%.2f\t%li\t%.2e\t%.2e\t%.2e\t%.2e\n",
		fnzone,
		dense.IonHigh[ipCARBON],
		ionbal.RateIonizTot[ipCARBON][0],
		ionbal.RateRecomTot[ipCARBON][0],
		dense.xIonDense[ipCARBON][1]/SDIV(dense.xIonDense[ipCARBON][0]),
		co.hevmol[ipCP]/SDIV(co.hevmol[ipATC]) );*/

	/* reset the var we just hosed */
	if( save_rec > 0. )
		ionbal.RateRecomTot[ipCARBON][0] = save_rec;

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

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "     IonCarbo returns;  fracs=" );
		for( i=1; i <= 7; i++ )
		{
			fprintf( ioQQQ, " %10.3e", dense.xIonDense[ipCARBON][i-1]/
			  dense.gas_phase[ipCARBON] );
		}
		fprintf( ioQQQ, "\n" );
	}

	{
		/* We have already checked the validity of CalcThermAveGaunt in sanitycheck.c.
			* Now we check to see if the InterpolateGff routine also works correctly.	*/
		/*@-redef@*/
		enum {AGN=FALSE};
		/* if set true there will be two problems at very low temps */
		/*@+redef@*/
		if( AGN )
		{
			phycon.te=10.;
			/* this tells ion_recomb to punch data */
			punch.lgioRecom = TRUE;
			/* this is where it will be punched */
			punch.ioRecom = ioQQQ;
			while( phycon.te<1e7 )
			{

				/* get recombination coefficients */
				ion_recomb(FALSE,(double*)dicoef,(double*)dite,ditcrt,aa,bb,cc,
					dd,ff,ipCARBON);
				phycon.te = phycon.te *2.f;
				tfidle(TRUE);
			}
			/* bail out */
			cdEXIT(EXIT_SUCCESS);
		}
	}

#	ifdef DEBUG_FUN
	fputs( " <->IonCarbo()\n", debug_fp );
#	endif
	return;
	/*lint +e740 unusual pointer cast */
}
#undef NDIM
