/*Fe2Lev16 compute populations and cooling due to 16 level Fe II ion */
#include "cddefines.h"
#ifdef NLFE2
#	undef NLFE2
#endif
#define	NLFE2	16
#include "showme.h"
#include "physconst.h"
#include "ionfracs.h"
#include "tcool.h"
#include "phycon.h"
#include "fe2cool.h"
#include "leveln.h"
#include "coladd.h"
#include "fe2lev16.h"

void Fe2Lev16(void)
{
	int lgNegPop;
	long int i, 
	  ipdest[NLFE2][NLFE2], 
	  j, 
	  nUsed, 
	  _r;

	double dfe2dt;
	double data[NLFE2][NLFE2], 
	  dest[NLFE2][NLFE2], 
	  popfe2[NLFE2], 
	  pump[NLFE2][NLFE2];
	static double Fe2A[NLFE2][NLFE2], 
	  Fe2CS[NLFE2][NLFE2];
	static double gfe2[NLFE2]={10.,8.,6.,4.,2.,10.,8.,6.,4.,8.,6.,4.,
	  2.,6.,4.,2.};
	static double ex[NLFE2]={0.,553.646,960.680,1241.151,1405.810,2694.300,
	  3496.489,4083.318,4485.486,11446.297,12074.546,12489.670,12728.966,
	  19387.344,19673.344,20006.633};
	static int _aini = 1;

	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp0[] = {0.00e00,2.13e-03,1.00e-16,1.00e-16,
		  1.00e-16,2.80e-04,1.00e-04,1.00e-16,1.00e-16,4.74e-03,1.00e-16,
		  1.00e-16,1.00e-16,1.00e-16,1.00e-16,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[0][j-1] = _itmp0[_r++];
			}
		}
		{ static double _itmp1[] = {0.00e00,0.00e00,1.57e-03,1.00e-16,
		  1.00e-16,1.00e-16,1.80e-04,6.90e-04,1.00e-16,1.31e-03,4.80e-04,
		  1.00e-16,1.00e-16,6.64e-03,1.00e-16,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[1][j-1] = _itmp1[_r++];
			}
		}
		{ static double _itmp2[] = {0.00e00,0.00e00,0.00e00,7.18e-04,
		  1.00e-16,1.00e-16,3.40e-05,1.00e-04,2.80e-05,9.70e-04,1.98e-03,
		  1.00e-16,1.00e-16,1.00e-16,6.81e-03,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[2][j-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  1.88e-04,1.00e-16,1.00e-16,2.50e-05,4.00e-05,1.00e-16,1.17e-03,
		  2.45e-03,7.80e-04,1.00e-16,1.00e-16,6.23e-03};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[3][j-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,1.00e-16,1.00e-16,1.00e-16,9.90e-06,1.00e-16,1.00e-16,
		  1.08e-03,3.32e-03,1.00e-16,2.20e-04,1.93e-03};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[4][j-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,5.84e-03,1.00e-16,1.00e-16,5.98e-03,3.12e-03,
		  1.00e-16,1.00e-16,3.56e-02,1.00e-16,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[5][j-1] = _itmp5[_r++];
			}
		}
		{ static double _itmp6[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,3.92e-03,1.00e-16,1.32e-03,2.49e-03,
		  4.18e-03,1.00e-16,8.83e-03,2.21e-02,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[6][j-1] = _itmp6[_r++];
			}
		}
		{ static double _itmp7[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,1.41e-03,1.00e-16,1.82e-03,
		  1.18e-03,4.75e-03,1.68e-03,1.29e-02,1.61e-02};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[7][j-1] = _itmp7[_r++];
			}
		}
		{ static double _itmp8[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,1.00e-16,1.00e-16,
		  2.12e-03,2.47e-03,1.00e-16,3.65e-03,2.13e-02};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[8][j-1] = _itmp8[_r++];
			}
		}
		{ static double _itmp9[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,2.56e-03,
		  1.00e-16,1.00e-16,2.23e-03,2.25e-03,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[9][j-1] = _itmp9[_r++];
			}
		}
		{ static double _itmp10[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,
		  1.36e-03,1.00e-16,2.70e-04,1.00e-16,2.76e-03};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[10][j-1] = _itmp10[_r++];
			}
		}
		{ static double _itmp11[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,3.71e-04,7.70e-05,1.60e-04,2.50e-04};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[11][j-1] = _itmp11[_r++];
			}
		}
		{ static double _itmp12[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,1.00e-16,1.50e-04,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[12][j-1] = _itmp12[_r++];
			}
		}
		{ static double _itmp13[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,1.90e-04,1.00e-16};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[13][j-1] = _itmp13[_r++];
			}
		}
		{ static double _itmp14[] = {0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,0.00e00,
		  0.00e00,0.00e00,0.00e00,0.00e00,5.50e-04};
		for( j=1, _r = 0; j <= NLFE2; j++ )
		{
			Fe2A[14][j-1] = _itmp14[_r++];
			}
		}
		for( j=1; j <= NLFE2; j++ )
		{
			Fe2A[15][j-1] = 0.00e00;
			}
		for( i=1; i <= NLFE2; i++ )
		{
			Fe2CS[0][i-1] = 0.000e00;
			}
		{ static double _itmp15[] = {5.521e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[1][i-1] = _itmp15[_r++];
			}
		}
		{ static double _itmp16[] = {1.491e00,5.442e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[2][i-1] = _itmp16[_r++];
			}
		}
		{ static double _itmp17[] = {6.747e-01,1.197e00,4.319e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[3][i-1] = _itmp17[_r++];
			}
		}
		{ static double _itmp18[] = {2.839e-01,4.191e-01,7.255e-01,
		  2.640e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[4][i-1] = _itmp18[_r++];
			}
		}
		{ static double _itmp19[] = {3.595e00,1.633e00,6.373e-01,2.108e-01,
		  5.804e-02,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[5][i-1] = _itmp19[_r++];
			}
		}
		{ static double _itmp20[] = {1.507e00,1.484e00,1.113e00,6.083e-01,
		  2.231e-01,8.365e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[6][i-1] = _itmp20[_r++];
			}
		}
		{ static double _itmp21[] = {4.969e-01,1.026e00,9.581e-01,7.646e-01,
		  4.466e-01,2.227e00,7.386e00,0.000e00,0.000e00,0.000e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[7][i-1] = _itmp21[_r++];
			}
		}
		{ static double _itmp22[] = {1.371e-01,4.457e-01,7.373e-01,
		  7.140e-01,4.239e-01,7.005e-01,2.069e00,5.627e00,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[8][i-1] = _itmp22[_r++];
			}
		}
		{ static double _itmp23[] = {1.098e01,2.479e00,5.033e-01,1.468e-01,
		  5.434e-02,2.106e00,8.906e-01,2.935e-01,6.481e-02,0.000e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[9][i-1] = _itmp23[_r++];
			}
		}
		{ static double _itmp24[] = {5.599e-01,6.540e00,2.847e00,6.264e-01,
		  1.177e-01,7.248e-01,8.511e-01,6.858e-01,3.146e-01,2.230e00,
		  0.000e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[10][i-1] = _itmp24[_r++];
			}
		}
		{ static double _itmp25[] = {1.909e-01,3.552e-01,3.583e00,2.445e00,
		  5.540e-01,1.477e-01,5.556e-01,4.996e-01,5.181e-01,8.177e-01,
		  1.752e00,0.000e00,0.000e00,0.000e00,0.000e00,0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[11][i-1] = _itmp25[_r++];
			}
		}
		{ static double _itmp26[] = {6.104e-02,1.227e-01,1.956e-01,
		  1.534e00,1.649e00,2.123e-02,1.471e-01,3.672e-01,3.274e-01,
		  3.689e-01,3.784e-01,1.015e00,0.000e00,0.000e00,0.000e00,
		  0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[12][i-1] = _itmp26[_r++];
			}
		}
		{ static double _itmp27[] = {9.480e-01,6.813e-01,3.850e-01,
		  1.707e-01,5.327e-02,1.232e00,6.882e-01,3.779e-01,1.971e-01,
		  1.008e00,5.870e-01,2.563e-01,7.694e-02,0.000e00,0.000e00,
		  0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[13][i-1] = _itmp27[_r++];
			}
		}
		{ static double _itmp28[] = {5.024e-01,2.276e-01,2.681e-01,
		  2.971e-01,1.967e-01,4.114e-01,5.541e-01,4.492e-01,2.809e-01,
		  4.690e-01,2.905e-01,3.180e-01,2.262e-01,1.280e00,0.000e00,
		  0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[14][i-1] = _itmp28[_r++];
			}
		}
		{ static double _itmp29[] = {3.278e-02,2.871e-01,2.451e-01,
		  1.319e-01,4.768e-02,1.404e-01,1.990e-01,2.603e-01,2.482e-01,
		  7.626e-02,2.877e-01,2.033e-01,8.596e-02,4.837e-01,6.294e-01,
		  0.000e00};
		for( i=1, _r = 0; i <= NLFE2; i++ )
		{
			Fe2CS[15][i-1] = _itmp29[_r++];
			}
		}
		_aini = 0;
	}

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

	/* number of levels in model ion */

	/* statistical weights */
	/* excitation energies in Kelvin
	 * excitation energies in Kelvin */

	/* 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)
	 * A's are saved into arrays as data(up,lo) */

	/* collision strengths from 
	 * >>refer Zhang, H.L., & Pradhan, A., 1995, A&A, 293, 953
	 * order is cs(low,up), for 10,000 K */

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

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

	/* LevelN itself is well-protected against zero abundances, and low temperatures
	 * >>chng 96 oct 16, reduce number of levels for low temp since
	 * can still be important coolant when highest level not populated */

	if( phycon.te < 500. )
	{
		nUsed = 5;
	}

	else if( phycon.te < 1e3 )
	{
		nUsed = 9;
	}

	else
	{
		nUsed = NLFE2;
	}

	/*void LevelN(long,long,float,double[],double[],double[],double*,
	double*,double*,long*,float*,float*,STRING,int*);*/
	LevelN(nUsed,
		NLFE2,
		xIonFracs[2][25],
		gfe2,/*double*/
		ex,/*double*/
		popfe2,/*double*/
		(double*)data,
	  (double*)dest,
	  (double*)pump,
	  (long*)ipdest,
	  &fe2cool.Fe2L16Tot,
	  &dfe2dt,
	  "FeII",
	  &lgNegPop);

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

	coladd("Fe 2",2,fe2cool.Fe2L16Tot);

	tcool.dCooldT += dfe2dt;

	/* all 10 transitions within ground term */
	fe2cool.fe21308 = popfe2[12]*Fe2A[7][12]*(ex[12] - ex[7])*BOLTZMANN;
	fe2cool.fe21207 = popfe2[11]*Fe2A[6][11]*(ex[11] - ex[6])*BOLTZMANN;
	fe2cool.fe21106 = popfe2[10]*Fe2A[5][10]*(ex[10] - ex[5])*BOLTZMANN;
	fe2cool.fe21006 = popfe2[9]*Fe2A[5][9]*(ex[9] - ex[5])*BOLTZMANN;
	fe2cool.fe21204 = popfe2[11]*Fe2A[3][11]*(ex[11] - ex[3])*BOLTZMANN;
	fe2cool.fe21103 = popfe2[10]*Fe2A[2][10]*(ex[10] - ex[2])*BOLTZMANN;
	fe2cool.fe21104 = popfe2[10]*Fe2A[3][10]*(ex[10] - ex[3])*BOLTZMANN;
	fe2cool.fe21001 = popfe2[9]*Fe2A[0][9]*(ex[9] - ex[0])*BOLTZMANN;
	fe2cool.fe21002 = popfe2[9]*Fe2A[1][9]*(ex[9] - ex[1])*BOLTZMANN;
	fe2cool.fe20201 = popfe2[1]*Fe2A[0][1]*(ex[1] - ex[0])*BOLTZMANN;
	fe2cool.fe20302 = popfe2[2]*Fe2A[1][2]*(ex[2] - ex[1])*BOLTZMANN;
	fe2cool.fe20706 = popfe2[6]*Fe2A[5][6]*(ex[6] - ex[5])*BOLTZMANN;
	fe2cool.fe20807 = popfe2[7]*Fe2A[6][7]*(ex[7] - ex[6])*BOLTZMANN;
	fe2cool.fe20908 = popfe2[8]*Fe2A[7][8]*(ex[8] - ex[7])*BOLTZMANN;
	fe2cool.fe21007 = popfe2[9]*Fe2A[6][9]*(ex[9] - ex[6])*BOLTZMANN;
	fe2cool.fe21107 = popfe2[10]*Fe2A[6][10]*(ex[10] - ex[6])*BOLTZMANN;
	fe2cool.fe21108 = popfe2[10]*Fe2A[7][10]*(ex[10] - ex[7])*BOLTZMANN;
	fe2cool.fe21110 = popfe2[10]*Fe2A[9][10]*(ex[10] - ex[9])*BOLTZMANN;
	fe2cool.fe21208 = popfe2[11]*Fe2A[7][11]*(ex[11] - ex[7])*BOLTZMANN;
	fe2cool.fe21209 = popfe2[11]*Fe2A[8][11]*(ex[11] - ex[8])*BOLTZMANN;
	fe2cool.fe21211 = popfe2[11]*Fe2A[10][11]*(ex[11] - ex[10])*BOLTZMANN;

	/* these are the highest lines in the 16 level atom, and so most uncertain */
	fe2cool.fe21406 = popfe2[13]*Fe2A[5][13]*(ex[13] - ex[5])*BOLTZMANN;
	fe2cool.fe21507 = popfe2[14]*Fe2A[6][14]*(ex[14] - ex[6])*BOLTZMANN;
	fe2cool.fe21508 = popfe2[14]*Fe2A[7][14]*(ex[14] - ex[7])*BOLTZMANN;
	fe2cool.fe21609 = popfe2[15]*Fe2A[8][15]*(ex[15] - ex[8])*BOLTZMANN;

	/* these are only predicted by the very large feii atom */
	fe2cool.fe25to6 = 0.;
	fe2cool.fe27to7 = 0.;
	fe2cool.fe28to8 = 0.;
	fe2cool.fe29to9 = 0.;
 	fe2cool.fe32to6 = 0.;
 	fe2cool.fe33to7 = 0.;
 	fe2cool.fe37to7 = 0.;
 	fe2cool.fe39to8 = 0.;
 	fe2cool.fe40to9 = 0.;
 	fe2cool.fe37to6 = 0.;
 	fe2cool.fe39to7 = 0.;
 	fe2cool.fe40to8 = 0.;
 	fe2cool.fe41to9 = 0.;
 	fe2cool.fe39to6 = 0.;
 	fe2cool.fe40to7 = 0.;
 	fe2cool.fe41to8 = 0.;
	fe2cool.fe42to6 = 0.;
	fe2cool.fe43to7 = 0.;
	fe2cool.fe42to7 = 0.;
	fe2cool.fe80to28 =0.;
	fe2cool.fe36to2 = 0.;
	fe2cool.fe36to3 = 0.;
	fe2cool.fe32to1 = 0.;
	fe2cool.fe33to2 = 0.;
	fe2cool.fe36to5 = 0.;
	fe2cool.fe32to2 = 0.;
	fe2cool.fe33to3 = 0.;
	fe2cool.fe30to3 = 0.;
	fe2cool.fe33to6 = 0.;
	fe2cool.fe24to2 = 0.;
	fe2cool.fe32to7 = 0.;
	fe2cool.fe35to8 = 0.;
	fe2cool.fe34to8 = 0.;
	fe2cool.fe27to6 = 0.;
	fe2cool.fe28to7 = 0.;
	fe2cool.fe30to8 = 0.;
	fe2cool.fe24to6 = 0.;
	fe2cool.fe29to8 = 0.;
	fe2cool.fe24to7 = 0.;
	fe2cool.fe22to7 = 0.;
	fe2cool.fe38to11 =0.;
	fe2cool.fe19to8 = 0.;
	fe2cool.fe17to6 = 0.;
	fe2cool.fe18to7 = 0.;
	fe2cool.fe18to8 = 0.;
	fe2cool.fe17to7 = 0.;

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