/*IonNitro ionization balance for nitrogen */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#ifdef NELEM
#	undef NELEM
#endif
#define	NDIM	8
#define	NELEM	7
#include "opacpoint.h"
#include "elmton.h"
#include "he.h"
#include "heat.h"
#include "nitexc.h"
#include "timed.h"
#include "tepowers.h"
#include "ionfracs.h"
#include "charexc.h"
#include "phycon.h"
#include "trace.h"
#include "photrate.h"
#include "collidionize.h"
#include "popexc.h"
#include "ionzer.h"
#include "makerecomb.h"
#include "theavy.h"
#include "photoionize.h"
#include "gammas.h"
#include "bidiag.h"
#include "ionheavy.h"

void IonNitro(void)
{
	long int i, 
	  _r;
	double aeff, 
	  cs5200, 
	  d5200r;
	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static double rec[NDIM - 2]={1.11e-10,7.91e-10,2.53e-9,1.08e-8,
	  2.07e-8,2.90e-8};
	static double pl[NDIM - 2]={-0.608,-0.639,-0.676,-0.765,-0.780,
	  -0.750};
	static double tlow[2]={0.,0.};
	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 nelem=7
	 *
	 * rates from Shull and van Steenberg, Ap.J. Sup 48, 95. */

	/* rec from +5, +4 from Arnaud et al 85 */
	/* Pequignot and Aldrovandi Ast Ap 161, 169. */
	/* nitrogen, atomic number 7 */
	if( !elmton.lgElmtOn[6] )
	{
		
#	ifdef DEBUG_FUN
		fputs( " <->IonNitro()\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);

	/* N+2 => N+ Butler and Dalgarno 1980B
	 * ct with he now from Sun, Sadeghpour, Kirby, Dalgarno, and Lafyatis,
	 * cfa preprint 4208
	 * this agrees with exp Fand&Kwong, ApJ 474, 529 */
	CharExc.CTHeavy[1][1] += (float)(he.hei*0.8e-10);
	/* N+3 => N+2 Butler and Dalgarno 1980B */
	CharExc.CTHeavy[1][2] += (float)(he.hei*1.5e-10);
	/* ce rate from quantal calc of ce, Ap.J. in press,
	 * Feickert, Blint, Surratt, and Watson, (preprint Sep 84).
	 * also Rittby et al J Phys B 17, L677, 1984.
	 * CR = 1.E-9 + 8E-12 * TE10 * SQRTE */
	CharExc.CTHeavy[1][3] += (float)(he.hei*2e-9);

	/* photoionization from 2D of NI */
	if( xIonFracs[1][6] > 0. )
	{
		d5200r = GammaK(OpacPoint.in1[0],OpacPoint.in1[1],OpacPoint.in1[2],1.);
		aeff = d5200r + 1.28e-5;
		cs5200 = 1.32e-4*phycon.te/(TePowers.te10*TePowers.te01);
		/* uses 1 for total pop of atom - so is normalized to unity */
		nitexc.p2nit = (float)(popexc(cs5200,4.,10.,aeff,2.769e4,1.)/aeff);
		/* valence shell photoionization */
		PhotRate.PhotoRate[0][2][0][6] = PhotRate.PhotoRate[0][2][0][6]*
		  (1. - nitexc.p2nit) + nitexc.p2nit*d5200r;
		PhotRate.PhotoRate[1][2][0][6] = PhotRate.PhotoRate[1][2][0][6]*
		  (1. - nitexc.p2nit) + heat.HeatNet*nitexc.p2nit;
	}
	else
	{
		nitexc.p2nit = 0.;
	}

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

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

	/* charge exchange destruction-creation of hydrogen
	 * CTHrec(1) = CTHrec(1) + xIonFracs(7,1) * c1
	 * CTHionAgnt(2) = xIonFracs(7,2) * c2
	 * CTHion(1) = CTHion(1) + CTHionAgnt(2)
	 * chrener = chrener +
	 * (hi*xIonFracs(7,2)*c2-hii*xIonFracs(7,1)*c1)*1.084e4/TE1RYD
	 * */
	if( trace.lgTrace && trace.lgHeavyBug )
	{
		fprintf( ioQQQ, "     IonNitro retun; frac=" );
		for( i=1; i <= 8; i++ )
		{
			fprintf( ioQQQ, "%10.3e", xIonFracs[i][6]/
			  xIonFracs[0][6] );
		}
		fprintf( ioQQQ, "\n" );
	}


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