/*Fe4Lev12 compute populations and cooling due to 12 level Fe IV ion */
#include "cddefines.h"
#include "showme.h"
#define	NLFE4	12
#include "physconst.h"
#include "ionfracs.h"
#include "fe4cool.h"
#include "tcool.h"
#include "leveln.h"
#include "coladd.h"
#include "fe4lev12.h"

void Fe4Lev12(void)
{
	int lgNegPop;
	long int i, 
	  ipdest[NLFE4][NLFE4], 
	  j, 
	  _r;

	double dfe4dt;

	double data[NLFE4][NLFE4], 
	  dest[NLFE4][NLFE4], 
	  popfe4[NLFE4], 
	  pump[NLFE4][NLFE4];

	static double Fe4A[NLFE4][NLFE4], 
	  Fe4CS[NLFE4][NLFE4];

	static double gfe4[NLFE4]={6.,12.,10.,6.,8.,6.,4.,2.,8.,2.,6.,4.};

	static double ex[NLFE4]={0.,46395.8,46464.,46475.6,46483.,50725.,
	  50838.,50945.,55796.,55966.,56021.,56025.};

	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp0[] = {0.,0.,0.,1.e-5,0.,1.368,.89,0.,
		  1.3e-3,1.8e-4,.056,.028};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[0][j-1] = _itmp0[_r++];
			}
		}
		{ static double _itmp1[] = {0.,0.,2.6e-8,0.,0.,0.,0.,0.,1.7e-7,
		  0.,0.,0.};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[1][j-1] = _itmp1[_r++];
			}
		}
		{ static double _itmp2[] = {0.,0.,0.,0.,3.5e-7,6.4e-10,0.,0.,
		  6.315e-4,0.,6.7e-7,0.};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[2][j-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {0.,0.,0.,0.,1.1e-6,6.8e-5,8.6e-6,
		  3.4e-10,7.6e-5,1.e-7,5.8e-4,2.8e-4};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[3][j-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {0.,0.,0.,0.,0.,1.5e-5,1.3e-9,0.,
		  7.6e-4,0.,1.1e-6,6.0e-7};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[4][j-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {0.,0.,0.,0.,0.,0.,1.1e-5,1.2e-13,
		  .038,9.9e-7,.022,.018};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[5][j-1] = _itmp5[_r++];
			}
		}
		{ static double _itmp6[] = {0.,0.,0.,0.,0.,0.,0.,3.7e-5,2.9e-6,
		  .034,3.5e-3,.039};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[6][j-1] = _itmp6[_r++];
			}
		}
		{ static double _itmp7[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,.058,
		  3.1e-6,1.4e-3};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[7][j-1] = _itmp7[_r++];
			}
		}
		{ static double _itmp8[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.3e-4,
		  3.1e-14};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[8][j-1] = _itmp8[_r++];
			}
		}
		{ static double _itmp9[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.9e-19,
		  1.0e-5};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[9][j-1] = _itmp9[_r++];
			}
		}
		{ static double _itmp10[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
		  0.,1.3e-7};
		for( j=1, _r = 0; j <= NLFE4; j++ )
		{
			Fe4A[10][j-1] = _itmp10[_r++];
			}
		}
		for( j=1; j <= NLFE4; j++ )
		{
			Fe4A[11][j-1] = 0.;
			}
		for( i=1; i <= NLFE4; i++ )
		{
			Fe4CS[0][i-1] = 0.;
			}
		{ static double _itmp11[] = {0.98,0.,0.,0.,0.,0.,0.,0.,0.,0.,
		  0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[1][i-1] = _itmp11[_r++];
			}
		}
		{ static double _itmp12[] = {0.8167,3.72,0.,0.,0.,0.,0.,0.,
		  0.,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[2][i-1] = _itmp12[_r++];
			}
		}
		{ static double _itmp13[] = {0.49,0.0475,0.330,0.,0.,0.,0.,
		  0.,0.,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[3][i-1] = _itmp13[_r++];
			}
		}
		{ static double _itmp14[] = {0.6533,0.473,2.26,1.64,0.,0.,0.,
		  0.,0.,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[4][i-1] = _itmp14[_r++];
			}
		}
		{ static double _itmp15[] = {0.45,0.686,0.446,0.106,0.254,0.,
		  0.,0.,0.,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[5][i-1] = _itmp15[_r++];
			}
		}
		{ static double _itmp16[] = {0.30,0.392,0.152,0.269,0.199,0.605,
		  0.,0.,0.,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[6][i-1] = _itmp16[_r++];
			}
		}
		{ static double _itmp17[] = {0.15,0.0207,0.190,0.0857,0.166,
		  0.195,0.327,0.,0.,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[7][i-1] = _itmp17[_r++];
			}
		}
		{ static double _itmp18[] = {0.512,1.23,0.733,0.174,0.398,0.623,
		  0.335,0.102,0.,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[8][i-1] = _itmp18[_r++];
			}
		}
		{ static double _itmp19[] = {0.128,0.0583,0.185,0.200,0.188,
		  0.0835,0.127,0.0498,0.0787,0.,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[9][i-1] = _itmp19[_r++];
			}
		}
		{ static double _itmp20[] = {0.384,0.578,0.534,0.363,0.417,
		  0.396,0.210,0.171,0.810,0.101,0.,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[10][i-1] = _itmp20[_r++];
			}
		}
		{ static double _itmp21[] = {0.256,0.234,0.306,0.318,0.403,
		  0.209,0.195,0.112,0.195,0.458,0.727,0.};
		for( i=1, _r = 0; i <= NLFE4; i++ )
		{
			Fe4CS[11][i-1] = _itmp21[_r++];
			}
		}
		_aini = 0;
	}

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

	/* number of levels in model ion */

	/* these are in wavenumbers
	 * data ex/ 0., 32245.5, 32293., 32301.2, 
	 *  1  32306., 35254., 35333., 35407., 38779., 38897., 38935.,
	 *  2  38938./ 
	 * excitation energies in Kelvin */

	/* A's are from Garstang, R.H., MNRAS 118, 572 (1958).
	 * each set is for a lower level indicated by second element in array,
	 * index runs over upper level
	 * A's are saved into arrays as data(up,lo) */

	/* collision strengths from Berrington and Pelan  Ast Ap S 114, 367.
	 * order is cs(low,up) */

	/* all elements are used, and must be set to zero if zero */
	for( i=0; i < NLFE4; i++ )
	{
		for( j=0; j < NLFE4; j++ )
		{
			data[j][i] = 1e33;
			dest[j][i] = 0.;
			pump[j][i] = 0.;
			ipdest[j][i] = 0;
		}
		data[i][i] = 0.;
		popfe4[i] = 0.;
	}

	/* fill in einstein As and collision strengths */
	for( i=0; i < NLFE4; i++ )
	{
		for( j=i + 1; j < NLFE4; j++ )
		{
			data[i][j] = Fe4A[i][j];
			data[j][i] = Fe4CS[j][i];
		}
	}

	/* leveln itself is well-protected against zero abundances,
	 * low temperatures */

	LevelN(NLFE4,
		NLFE4,
		xIonFracs[4][25],
		gfe4,
		ex,
	  popfe4,
	  (double*)data,
	  (double*)dest,
	  (double*)pump,
	  (long*)ipdest,
	  &fe4cool.Fe4CoolTot,
	  &dfe4dt,
	  "FeIV",
	  &lgNegPop);

	if( lgNegPop )
	{
		fprintf( ioQQQ, " fe4levl2 found negative populations\n" );
		ShowMe();
		puts( "[Stop in fe4levl2]" );
		exit(1);
	}

	coladd("Fe 4",0,fe4cool.Fe4CoolTot);

	tcool.dCooldT += dfe4dt;

	/* three transitions hst can observe
	 * 4 -1 3095.9A and 5 -1 3095.9A */
	fe4cool.fe40401 = (popfe4[3]*Fe4A[0][3] + popfe4[4]*Fe4A[0][4])*
	  6.42e-12;
	fe4cool.fe42836 = popfe4[5]*Fe4A[0][5]*7.02e-12;
	fe4cool.fe42829 = popfe4[6]*Fe4A[0][6]*7.03e-12;
	fe4cool.fe42567 = (popfe4[10]*Fe4A[0][10] + popfe4[11]*Fe4A[0][11])*
	  7.752e-12;
	fe4cool.fe41207 = popfe4[11]*Fe4A[6][11]*(ex[11] - ex[6])*BOLTZMANN;
	fe4cool.fe41206 = popfe4[11]*Fe4A[5][11]*(ex[11] - ex[5])*BOLTZMANN;
	fe4cool.fe41106 = popfe4[10]*Fe4A[5][10]*(ex[10] - ex[5])*BOLTZMANN;
	fe4cool.fe41007 = popfe4[9]*Fe4A[6][9]*(ex[9] - ex[6])*BOLTZMANN;
	fe4cool.fe41008 = popfe4[9]*Fe4A[7][9]*(ex[9] - ex[7])*BOLTZMANN;
	fe4cool.fe40906 = popfe4[8]*Fe4A[5][8]*(ex[8] - ex[5])*BOLTZMANN;


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