/* 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 */
/*IonIron ionization balance for iron */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#define	NDIM	27
#include "dense.h"
#include "fe.h"
#include "gammas.h"
#include "opacity.h"
#include "trace.h"
#include "grainvar.h"
#include "ionbal.h"
#include "mole.h"

void IonIron(void)
{
	long int i, 
	  limit, 
	  limit2, 
	  _r;

	static double dicoef[2][NDIM - 1], 
	  dite[2][NDIM - 1];
	static double fyield[27]={.34,.34,.35,.35,.36,.37,.37,.38,.39,.40,
	  .41,.42,.43,.44,.45,.46,.47,.47,.48,.48,.49,.49,.11,.75,0.,0.,
	  0.};
	static double ditcrt[NDIM - 1]={1e11,1e11,1e11,1e11,1e11,1e11,1e11,
	  1e11,1e11,1e11,1e11,1e11,1e11,1e11,1e11,1e11,1e11,1e11,1e11,
	  1e11,1e11,1e11,1e11,1e11,1e11,1e11};
	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.,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.,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.,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.,0.,0.,0.,0.,0.,0.};
	static double ff[NDIM - 1]={0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
	  0.,0.,0.,0.,0.,0.,0.};
	/* >>chng 03 aug 15, zero out this array.  having non-zero values in
	 * this array had the effect of using 0 for the guess to the low-T dr
	 * in makerecomb */
	/*static double ff[NDIM - 1]={0.01,0.01,0.01,0.01,0.01,0.01,0.01,
	  0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,
	  0.01,0.01,0.01,0.01,0.01,0.01,0.01};*/
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp2[] = {1.58e-3,8.38e-3,1.54e-2,3.75e-2,
		  0.117,0.254,0.291,0.150,0.140,0.100,0.200,0.240,0.260,0.190,
		  0.120,0.350,0.066,0.10,0.13,0.23,0.14,0.11,0.041,0.747,0.519,
		  0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[0][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {.456,.323,.310,.411,.359,.0975,
		  .229,4.20,3.30,5.30,1.50,0.700,.600,.5,1.,0.,7.8,6.3,5.5,
		  3.6,4.9,1.6,4.2,.284,.279,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dicoef[1][i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {6.00e4,1.94e5,3.31e5,4.32e5,6.28e5,
		  7.50e5,7.73e5,2.62e5,2.50e5,2.57e5,2.84e5,8.69e5,4.21e5,
		  4.57e5,2.85e5,8.18e5,1.51e6,1.30e6,1.19e6,1.09e6,9.62e5,
		  7.23e5,4.23e5,5.84e7,6.01e7,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[0][i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {8.97e4,1.71e5,2.73e5,3.49e5,5.29e5,
		  4.69e5,6.54e5,1.32e6,1.33e6,1.41e6,1.52e6,1.51e6,1.82e6,
		  1.84e6,2.31e6,0.,9.98e6,9.98e6,1.00e7,1.10e7,8.34e6,1.01e7,
		  1.07e7,1.17e7,9.97e6,0.};
		for( i=1, _r = 0; i <= (NDIM - 1); i++ )
		{
			dite[1][i-1] = _itmp5[_r++];
			}
		}
		_aini = 0;
	}

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

	/* iron nelem=26
	 *
	 * Fe rates from Woods et al. Ap.J. 249, 399.
	 * rec from +23, 24 25 from Arnauld et al 85 */

	/* Pequignot and Aldrovandi Ast Ap 161, 169. */

	/* free statement numbers >= 7
	 *
	 * iron, atomic number 26 */
	if( !dense.lgElmtOn[25] )
	{
		fe.fekcld = 0.;
		fe.fekhot = 0.;
		fe.fegrain = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->IonIron()\n", debug_fp );
#		endif
		return;
	}

	ion_zero(ipIRON);

	ion_photo(ipIRON,FALSE);
	/* debugging printout for shell photo rates */
	/*GammaPrtRate( ioQQQ , 0 , ipIRON , TRUE);*/
	/*GammaPrtShells( ipIRON , 13 );
	GammaPrtShells( ipIRON , 12 );
	GammaPrtShells( ipIRON , 0 );*/
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			long int iplow , iphi , ipop , ns , ion ;
			double rate;
			ns = 6;
			ion = 0;
			/* show what some of the ionization sources are */
			iplow = opac.ipElement[ipIRON][ion][ns][0];
			iphi = opac.ipElement[ipIRON][ion][ns][1];
			ipop = opac.ipElement[ipIRON][ion][ns][2];
			rate = ionbal.PhotoRate_Shell[ipIRON][ion][ns][0];
			GammaPrt( iplow , iphi , ipop , ioQQQ, rate , rate*0.1 );
		}
	}

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

	/* get recombination coefficients */
	ion_recomb(FALSE,(double*)dicoef,(double*)dite,ditcrt,aa,bb,cc,
	  dd,ff,ipIRON);

	/* >>chng 06 Feb 28 -- NPA.  Add in charge transfer ionization of Mg with S+, Si+, and C+ */
	/* only include this if molecular network is enabled - otherwise no feedback onto
	 * Si+, S+, and C+ soln */
	if( !co.lgNoCOMole )
	{
		ionbal.PhotoRate_Shell[ipIRON][0][6][0] += 
			rate_COmole.SP_FE_S_FEP*dense.xIonDense[ipSULPHUR][1] +
			rate_COmole.SIP_FE_SI_FEP*dense.xIonDense[ipSILICON][1] + 
			rate_COmole.CP_FE_C_FEP*dense.xIonDense[ipCARBON][1];
	}

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

	/* now find total Auger yield of K-alphas
	 * "cold" iron has M-shell electrons, up to Fe 18 */
	fe.fekcld = 0.;
	limit = MIN2(18,dense.IonHigh[ipIRON]);

	for( i=dense.IonLow[ipIRON]; i < limit; i++ )
	{
		fe.fekcld += 
			(float)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*dense.xIonDense[ipIRON][i]*
		  fyield[i]);
	}

	/* same sum for hot iron */
	fe.fekhot = 0.;
	limit = MAX2(18,dense.IonLow[ipIRON]);

	limit2 = MIN2(ipIRON+1,dense.IonHigh[ipIRON]);
	for( i=limit; i < limit2; i++ )
	{
		fe.fekhot += 
			(float)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*dense.xIonDense[ipIRON][i]*
		  fyield[i]);
	}

	/* Fe Ka from grains - Fe in grains assumed to be atomic
	 * gv.elmSumAbund[ipIRON] is number density of iron added over all grain species */
	i = 0;
	/* fyield is 0.34 for atomic fe */
	fe.fegrain = (float)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*fyield[i]*
				 gv.elmSumAbund[ipIRON]);

	if( trace.lgTrace && trace.lgHeavyBug )
	{
		/* print densities of each stage of ionization */
		fprintf( ioQQQ, "     Fe" );
		for( i=1; i <= 15; i++ )
		{
			fprintf( ioQQQ, "\t%.1e", dense.xIonDense[ipIRON][i-1]/dense.gas_phase[ipIRON] );
		}
		fprintf( ioQQQ, "\n" );
	}

	if( trace.lgTrace && trace.lgFeBug )
	{
		fprintf( ioQQQ, " IonIron-Abund:" );
		for( i=1; i <= 27; i++ )
		{
			fprintf( ioQQQ, "%10.2e", dense.xIonDense[ipIRON][i-1] );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, " IonIron - Ka(Auger)%10.2e\n", fe.fekcld + 
		  fe.fekhot );
	}

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

