/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*IonNitro ionization balance for nitrogen */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#define	NDIM	8
#include "dense.h"
#include "thermal.h"
#include "atoms.h"
#include "opacity.h"
#include "trace.h"
#include "conv.h"
#include "gammas.h"
#include "ionbal.h"

void IonNitro(void)
{
	long int i, 
	  _r;
	double save_rec;
	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static double ditcrt[NDIM - 1]={1.8e4,1.8e4,2.4e4,1.5e4,6.8e5,1.0e6,
	  1e20};
	static double aa[NDIM - 1]={0.0,0.0320,-0.8806,0.4134,0.,0.,0.};
	static double bb[NDIM - 1]={0.6310,-0.6624,11.2406,-4.6319,0.,0.,
	  0.};
	static double cc[NDIM - 1]={0.1990,4.3191,30.7066,25.9172,0.,0.,
	  0.};
	static double dd[NDIM - 1]={-0.0197,0.0003,-1.1721,-2.2290,0.,0.,
	  0.};
	static double ff[NDIM - 1]={0.4398,0.5946,0.6127,0.2360,0.1,0.,
	  0.};
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp2[] = {2.98e-3,7.41e-3,1.13e-2,2.62e-3,
		  7.5e-2,4.61e-2,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[0][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {0.,.0764,.164,.243,.35,.309,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[1][i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {2.2e5,2.01e5,1.72e5,1.02e5,4.75e6,
		  5.44e6,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[0][i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {0.,7.37e4,2.25e5,1.25e5,8.35e5,
		  1.14e6,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[1][i-1] = _itmp5[_r++];
			}
		}
		_aini = 0;
	}

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

	/* nitrogen, atomic number 7 */
	if( !dense.lgElmtOn[6] )
	{
		
#		ifdef DEBUG_FUN
		fputs( " <->IonNitro()\n", debug_fp );
#		endif
		return;
	}

	ion_zero(ipNITROGEN);

	ion_photo(ipNITROGEN,FALSE);

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

	/* get recombination coefficients */
	/*lint -e740 unusual pointer casts */
	/*lint -e64 type mismatch */
	ion_recomb(FALSE,(double*)dicoef,(double*)dite,ditcrt,aa,bb,cc,
	  dd,ff,ipNITROGEN);
	/*lint +e740 unusual pointer casts */
	/*lint +e64 type mismatch */

		/* >>chng 04 jul 09, synch up ion and co solvers.
	 * the co solver will have a different N/N+ balance that
	 * is strongly affected by the chemistry if we are in neutral gas
	 * in this case use co solver's N/N+ balance */
	save_rec = ionbal.RateRecomTot[ipNITROGEN][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 0
	if( dense.IonLow[ipNITROGEN]==0 && /*!conv.lgSearch &&*/
		co.hevmol[ipNP] > SMALLFLOAT && 
		ionbal.RateIonizTot[ipNITROGEN][0]*co.hevmol[ipATN]> 0. )
	{
		ionbal.RateRecomTot[ipNITROGEN][0] = 
			ionbal.RateIonizTot[ipNITROGEN][0]*
			co.hevmol[ipATN]/co.hevmol[ipNP];
	}
#	endif


	/* photoionization from 2D of NI, is atomic nitrogen present? */
	if( dense.xIonDense[ipNITROGEN][0] > 0. )
	{
		atoms.d5200r = (float)GammaK(opac.in1[0],opac.in1[1],opac.in1[2],1.);
		/*aeff = atoms.d5200r + 1.28e-5;*/
		/* >>chng 01 sep 02, use values from coolnitr */
		/*cs5200 = 1.32e-4*phycon.te/(phycon.te10*phycon.te01);*/
		/* uses 1 for total pop of atom - so is normalized to unity */
		/*atoms.p2nit = (float)(atom_pop2(cs5200,4.,10.,aeff,2.769e4,1.)/aeff);*/
		/* if very first call then set to zero */
		if( !conv.nTotalIoniz ) 
			atoms.p2nit = 0.;
		/* valence shell photoionization, [0][6] => atomic nitrogen */
		ionbal.PhotoRate_Shell[ipNITROGEN][0][2][0] = ionbal.PhotoRate_Shell[ipNITROGEN][0][2][0]*
		  (1. - atoms.p2nit) + atoms.p2nit*atoms.d5200r;
		ionbal.PhotoRate_Shell[ipNITROGEN][0][2][1] = ionbal.PhotoRate_Shell[ipNITROGEN][0][2][1]*
		  (1. - atoms.p2nit) + thermal.HeatNet*atoms.p2nit;
	}
	else
	{
		atoms.p2nit = 0.;
		atoms.d5200r = 0.;
	}

	/* solve for ionization balance */
	ion_solver(ipNITROGEN,FALSE);

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

	if( trace.lgTrace && trace.lgHeavyBug )
	{
		fprintf( ioQQQ, "     IonNitro retun; frac=" );
		for( i=1; i <= 8; i++ )
		{
			fprintf( ioQQQ, "%10.3e", dense.xIonDense[ipNITROGEN][i-1]/
			  dense.gas_phase[ipNITROGEN] );
		}
		fprintf( ioQQQ, "\n" );
	}

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

	return;
}

