/*IonIron ionization balance for iron */
#include "cddefines.h"
#ifdef NDIM
#	undef NDIM
#endif
#ifdef NELEM
#	undef NELEM
#endif
#define	NDIM	27
#define	NELEM	26
#include "elmton.h"
#include "fekems.h"
#include "gammas.h"
#include "trace.h"
#include "timed.h"
#include "photrate.h"
#include "ionrange.h"
#include "ionfracs.h"
#include "ionzer.h"
#include "collidionize.h"
#include "makerecomb.h"
#include "theavy.h"
#include "photoionize.h"
#include "bidiag.h"
#include "ionheavy.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 rec[NDIM - 2]={5.2e-10,2.40e-9,3.20e-9,4.17e-9,9.44e-9,
	  2.14e-8,4.48e-8,8.75e-8,1.41e-7,2.26e-7,2.73e-7,3.56e-7,4.42e-7,
	  5.52e-7,7.00e-7,8.59e-7,8.55e-7,8.93e-7,9.64e-7,1.06e-6,1.19e-6,
	  1.22e-6,2.94e-6,5.00e-6,1.97e-6};
	static double pl[NDIM - 2]={-0.891,-0.843,-0.746,-0.682,-0.699,
	  -0.728,-0.759,-0.790,-0.810,-0.829,-0.828,-0.834,-0.836,-0.840,
	  -0.846,-0.850,-0.836,-0.824,-0.816,-0.811,-0.808,-0.800,-0.852,
	  -0.875,-0.787};
	static double tlow[2]={1.28e-10,-0.686};
	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.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( !elmton.lgElmtOn[25] )
	{
		fekems.fekcld = 0.;
		fekems.fekhot = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->IonIron()\n", debug_fp );
#		endif
		return;
	}

	ionzer(NELEM);

	PhotoIonize(NELEM-1,FALSE);
	/* debugging printout for shell photo rates */
	/*GammaPrtRate( ioQQQ , 0 , NELEM-1 );*/
	{
		enum {DEBUG=FALSE};
		if( DEBUG )
		{
#			include "opacpoint.h"
			long int iplow , iphi , ipop , ns , nion , nelem;
			double rate;
			ns = 6;
			nion = 0;
			nelem = 25;
			/* option to redo the rates only on occasion */
			iplow = OpacPoint.ipElement[0][ns][nion][nelem];
			iphi = OpacPoint.ipElement[1][ns][nion][nelem];
			ipop = OpacPoint.ipElement[2][ns][nion][nelem];
			rate = PhotRate.PhotoRate[0][ns][nion][nelem];
			GammaPrt( iplow , iphi , ipop , ioQQQ, rate , rate*0.1 );
		}
	}

	/* 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);

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

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

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

	for( i=IonRange.IonLow[NELEM-1]-1; i < limit; i++ )
	{
		fekems.fekcld += 
			(float)(PhotRate.PhotoRate[0][0][i][NELEM-1]*xIonFracs[i+1][NELEM-1]*
		  fyield[i]);
	}

	/* same sum for hot iron */
	fekems.fekhot = 0.;
	limit = MAX2(19,IonRange.IonLow[NELEM-1]);
	limit2 = MIN2(NELEM,IonRange.IonHigh[NELEM-1]-1);
	for( i=limit-1; i < limit2; i++ )
	{
		fekems.fekhot += 
			(float)(PhotRate.PhotoRate[0][0][i][NELEM-1]*xIonFracs[i+1][NELEM-1]*
		  fyield[i]);
	}

	if( trace.lgTrace && trace.lgHeavyBug )
	{
		fprintf( ioQQQ, "     Fe" );
		for( i=1; i <= 15; i++ )
		{
			fprintf( ioQQQ, "%8.1e", xIonFracs[i][NELEM-1]/
			  xIonFracs[0][NELEM-1] );
		}
		fprintf( ioQQQ, "\n" );
	}

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

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

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