/*oilevl get OI level population with Ly-beta pumping */
#include "cddefines.h"
#include <string.h>
#include "rfield.h"
#include "phycon.h"
#include "edsqte.h"
#include "ipoi.h"
#include "trace.h"
#include "bowoi.h"
#include "popoi.h"
#include "negoi.h"
#include "taulines.h"
#include "tepowers.h"
#include "typmatrx.h"
#include "linpack.h"
#include "putcs.h"
#include "matin1.h"
#include "veclib.h"
#include "oilevl.h"

void oilevl(double abund, 
  double *coloi)
{
	int lgNegPop;

	long int i, 
	  ipiv[6], 
	  j, 
	  job, 
	  ner;

	double a21, 
	  a32, 
	  a41, 
	  a43, 
	  a51, 
	  a53, 
	  a62, 
	  a64, 
	  a65, 
	  c12, 
	  c13, 
	  c14, 
	  c15, 
	  c16, 
	  c21, 
	  c23, 
	  c24, 
	  c25, 
	  c26, 
	  c31, 
	  c32, 
	  c34, 
	  c35, 
	  c36, 
	  c41, 
	  c42, 
	  c43, 
	  c45, 
	  c46, 
	  c51, 
	  c52, 
	  c53, 
	  c54, 
	  c56, 
	  c61, 
	  c62, 
	  c63, 
	  c64, 
	  c65, 
	  cs, 
	  deptoi[6], 
	  e12, 
	  e23, 
	  e34, 
	  e45, 
	  e56, 
	  simple;

	double amat[6][6], 
	  bvec[6], 
	  rcond, 
	  work[6], 
	  zz[7][7];

	static float g[6]={9.,3.,9.,3.,15.,9};

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

	/* following used for linpac matrix inversion */

	/* compute emission from six level OI atom*/

	/* boltzmann factors for ContBoltz since collisions not dominant in UV tran
	 * ipoiex is array lof dEnergy for each level, set in DoPoint */
	e12 = rfield.ContBoltz[ipoi.ipoiex[0]-1];
	e23 = rfield.ContBoltz[ipoi.ipoiex[1]-1];
	e34 = rfield.ContBoltz[ipoi.ipoiex[2]-1];
	e45 = rfield.ContBoltz[ipoi.ipoiex[3]-1];
	e56 = rfield.ContBoltz[ipoi.ipoiex[4]-1];

	/* total rad rates here have dest by background continuum */
	a21 = TauLines[ipT1304].Aul*(TauLines[ipT1304].Pdest + TauLines[ipT1304].Pesc);
	a41 = TauLines[ipT1039].Aul*(TauLines[ipT1039].Pdest + TauLines[ipT1039].Pesc);
	a51 = TauLines[ipTO1025].Aul*(TauLines[ipTO1025].Pdest + TauLines[ipTO1025].Pesc);
	a51 = bowoi.pmpo51;
	a32 = TauLines[ipT8446].Aul*(TauLines[ipT8446].Pdest + TauLines[ipT8446].Pesc);
	a62 = TauLines[ipT4368].Aul*(TauLines[ipT4368].Pdest + TauLines[ipT4368].Pesc);
	a43 = TauLines[ipTOI13].Aul*(TauLines[ipTOI13].Pdest + TauLines[ipTOI13].Pesc);
	a53 = TauLines[ipTOI11].Aul*(TauLines[ipTOI11].Pdest + TauLines[ipTOI11].Pesc);
	a64 = TauLines[ipTOI29].Aul*(TauLines[ipTOI29].Pdest + TauLines[ipTOI29].Pesc);
	a65 = TauLines[ipTOI46].Aul*(TauLines[ipTOI46].Pdest + TauLines[ipTOI46].Pesc);

	/* even at density of 10^17 excited states not in lte due
	 * to fast transitions down - just need to kill a21 to get to unity at 10^17*/

	/* the 2-1 transition is 1302, cs Wang and McConkey '92 Jphys B 25, 5461 */
	cs = 2.151e-5*phycon.te/TePowers.te03;
	PutCS(cs,&TauLines[ipT1304]);

	/* the 5-1 transition is 1027, cs Wang and McConkey '92 Jphys B 25, 5461 */
	cs = 9.25e-7*phycon.te*TePowers.te10/TePowers.te01/TePowers.te01;
	PutCS(cs,&TauLines[ipTO1025]);
	c21 = edsqte.cdsqte*TauLines[ipT1304].cs/g[1];
	c51 = edsqte.cdsqte*TauLines[ipTO1025].cs/g[4];

	/* all following are g-bar approx, g-bar = 0.2 */
	c31 = edsqte.cdsqte*1.0/g[2];
	PutCS(0.27,&TauLines[ipT1039]);
	c41 = edsqte.cdsqte*TauLines[ipT1039].cs/g[3];
	c61 = edsqte.cdsqte*1./g[5];

	c12 = c21*g[1]/g[0]*e12;
	c13 = c31*g[2]/g[0]*e12*e23;
	c14 = c41*g[3]/g[0]*e12*e23*e34;
	c15 = c51*g[4]/g[0]*e12*e23*e34*e45;
	c16 = c61*g[5]/g[0]*e12*e23*e34*e45*e56;

	c32 = edsqte.cdsqte*85./g[2];
	c42 = edsqte.cdsqte*85./g[3];
	c52 = edsqte.cdsqte*85./g[4];
	c62 = edsqte.cdsqte*85./g[5];

	c23 = c32*g[2]/g[1]*e23;
	c24 = c42*g[3]/g[1]*e23*e34;
	c25 = c52*g[4]/g[1]*e23*e34*e45;
	c26 = c62*g[5]/g[1]*e23*e34*e45*e56;

	c43 = edsqte.cdsqte*70./g[3];
	c53 = edsqte.cdsqte*312./g[4];
	c63 = edsqte.cdsqte*1./g[5];

	c34 = c43*g[3]/g[2]*e34;
	c35 = c53*g[4]/g[2]*e34*e45;
	c36 = c63*g[5]/g[2]*e34*e45*e56;

	c54 = edsqte.cdsqte*50./g[4];
	c64 = edsqte.cdsqte*415./g[5];

	c45 = c54*g[4]/g[3]*e45;
	c46 = c64*g[5]/g[3]*e45*e56;

	c65 = edsqte.cdsqte*400./g[5];
	c56 = c65*g[5]/g[4]*e56;

	/* 666 error! this must have all stimulated emission, pump by cont, etc*/

	/* this is check for whether matrix inversion likely to fail */
	simple = (c16 + bowoi.pmpo15)/(c61 + c62 + c64 + a65 + a64 + a62);
	if( simple < 1e-19 )
	{
		popoiCom.popoi[0] = (float)abund;
		for( i=1; i < 6; i++ )
		{
			popoiCom.popoi[i] = 0.;
		}
		*coloi = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->oilevl()\n", debug_fp );
#		endif
		return;
	}

	/*--------------------------------------------------------- */

	for( i=0; i < 6; i++ )
	{
		zz[i][0] = 1.0;
		zz[6][i] = 0.;
	}

	/* first equation is sum of populations */
	zz[6][0] = abund;

	/* level two, 3s 3So */
	zz[0][1] = -c12;
	zz[1][1] = c21 + c23 + c24 + c25 + c26 + a21;
	zz[2][1] = -c32 - a32;
	zz[3][1] = -c42;
	zz[4][1] = -c52;
	zz[5][1] = -c62 - a62;

	/* level three */
	zz[0][2] = -c13;
	zz[1][2] = -c23;
	zz[2][2] = c31 + c32 + c34 + c35 + c36 + a32;
	zz[3][2] = -c43 - a43;
	zz[4][2] = -c53 - a53;
	zz[5][2] = -c63;

	/* level four */
	zz[0][3] = -c14;
	zz[1][3] = -c24;
	zz[2][3] = -c34;
	zz[3][3] = c41 + c42 + c43 + c45 + c46 + a41 + a43;
	zz[4][3] = -c54;
	zz[5][3] = -c64 - a64;

	/* level five */
	zz[0][4] = -c15 - bowoi.pmpo15;
	zz[1][4] = -c25;
	zz[2][4] = -c35;
	zz[3][4] = -c45;
	zz[4][4] = c51 + c52 + c53 + c54 + c56 + a51 + a53;
	zz[5][4] = -c65 - a65;

	/* level six */
	zz[0][5] = -c16;
	zz[1][5] = -c26;
	zz[2][5] = -c36;
	zz[3][5] = -c46;
	zz[4][5] = -c56;
	zz[5][5] = c61 + c62 + c63 + c64 + c65 + a65 + a64 + a62;

	/* which matrix solver? */
	if( strcmp(TypMatrx.chMatrix,"matin1 ") == 0 )
	{
		/*matin1();*/
		ner = 1;
		if( ner != 0 )
		{
			fprintf( ioQQQ, " oi matrix error\n" );
			puts( "[Stop in oilevl]" );
			exit(1);
		}
	}

	else if( strcmp(TypMatrx.chMatrix,"linpack") == 0 )
	{
		/* this one may be more robust */
		for( j=0; j < 6; j++ )
		{
			for( i=0; i < 6; i++ )
			{
				amat[i][j] = zz[i][j];
			}
			bvec[j] = zz[6][j];
		}

		DGETRF(6,6,(double*)amat,6,ipiv,&ner);
		DGETRS('N',6,1,(double*)amat,6,ipiv,bvec,6,&ner);
		if( ner != 0 )
		{
			fprintf( ioQQQ, " oi levl dgetrs error\n" );
			puts( "[Stop in oilevl]" );
			exit(1);
		}

		/* now put results back into z so rest of code treates only
		 * one case - as if matin1 had been used */
		for( i=0; i < 6; i++ )
		{
			zz[6][i] = bvec[i];
		}
	}

	else if( strcmp(TypMatrx.chMatrix,"veclib ") == 0 )
	{
		/* Jason found this one on the Exemplar, distributed source just stops */
		fprintf( ioQQQ, " this has not been checked since H atom conv\n" );
		for( j=0; j < 6; j++ )
		{
			for( i=0; i < 6; i++ )
			{
				amat[i][j] = zz[i][j];
			}
			bvec[j] = zz[6][j];
		}

		job = 0;
		rcond = 0.;
		dgeco((double*)amat,6,6,ipiv,rcond,work);
		dgesl((double*)amat,6,6,ipiv,bvec,job);
		/* now put results back into z so rest of code treates only
		 * one case - as if matin1 had been used */
		for( i=0; i < 6; i++ )
		{
			zz[6][i] = bvec[i];
		}
		puts( "[Stop in oilevl]" );
		exit(1);
	}

	else
	{
		fprintf( ioQQQ, " chMatrix type insane in hydrogen, was%7.7s\n", 
		  TypMatrx.chMatrix );
		puts( "[Stop in oilevl]" );
		exit(1);
	}

	lgNegPop = FALSE;
	for( i=0; i < 6; i++ )
	{
		popoiCom.popoi[i] = (float)zz[6][i];
		if( popoiCom.popoi[i] < 0. )
			lgNegPop = TRUE;
	}

	/* following used to confirm that all dep coef are unity at
	 * density of 1e17, t=10,000, and all A's set to zero */
	if( trace.lgTrace && trace.lgTr8446 )
	{
		deptoi[0] = 1.;
		deptoi[1] = popoiCom.popoi[1]/popoiCom.popoi[0]/(g[1]/g[0]*
		  e12);
		deptoi[2] = popoiCom.popoi[2]/popoiCom.popoi[0]/(g[2]/g[0]*
		  e12*e23);
		deptoi[3] = popoiCom.popoi[3]/popoiCom.popoi[0]/(g[3]/g[0]*
		  e12*e23*e34);
		deptoi[4] = popoiCom.popoi[4]/popoiCom.popoi[0]/(g[4]/g[0]*
		  e12*e23*e34*e45);
		deptoi[5] = popoiCom.popoi[5]/popoiCom.popoi[0]/(g[5]/g[0]*
		  e12*e23*e34*e45*e56);

		fprintf( ioQQQ, " oilevl finds levl pop" );
		for(i=0; i < 6; i++)
			fprintf( ioQQQ, "%11.3e", popoiCom.popoi[i] );
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, " oilevl finds dep coef" );
		for(i=0; i < 6; i++)
			fprintf( ioQQQ, "%11.3e", deptoi[i] );
		fprintf( ioQQQ, "\n" );
	}

	/* this happens due to numerical instability in matrix inversion routine */
	if( lgNegPop )
	{
		negoi.nNegOI += 1;

		fprintf( ioQQQ, " OILEVL finds negative population" );
		for(i=0;i < 6; i++)
			fprintf( ioQQQ, "%10.2e", popoiCom.popoi[i] );
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, " simple 5 =%10.2e\n", simple );

		popoiCom.popoi[5] = 0.;
		popoiCom.popoi[4] = (float)(abund*(c15 + bowoi.pmpo15)/(a51 + a53 + 
		  c51 + c53));
		popoiCom.popoi[3] = 0.;
		popoiCom.popoi[2] = (float)(popoiCom.popoi[4]*(a53 + c53)/(a32 + c32));
		popoiCom.popoi[1] = (float)((popoiCom.popoi[2]*(a32 + c32) + abund*
		  c12)/(a21 + c21));

		popoiCom.popoi[0] = (float)abund;
		/*  write(QQ,'('' OILEVL resets this to simple pop'',1P,6E10.2)')
		 *  1   popoi */
	}

	/* this is total cooling due to model atom, can be neg (heating) */
	*coloi = 
	   (popoiCom.popoi[0]*c12 - popoiCom.popoi[1]*c21)*1.53e-11 + 
	   (popoiCom.popoi[0]*c14 - popoiCom.popoi[3]*c41)*1.92e-11 + 
	   (popoiCom.popoi[0]*c15 - popoiCom.popoi[4]*c51)*1.94e-11 + 
	   (popoiCom.popoi[1]*c23 - popoiCom.popoi[2]*c32)*2.36e-12 + 
	   (popoiCom.popoi[1]*c26 - popoiCom.popoi[5]*c62)*4.55e-12 + 
	   (popoiCom.popoi[2]*c35 - popoiCom.popoi[4]*c53)*1.76e-12 + 
	   (popoiCom.popoi[2]*c34 - popoiCom.popoi[3]*c43)*1.52e-12 + 
	   (popoiCom.popoi[3]*c46 - popoiCom.popoi[5]*c64)*6.86e-13 + 
	   (popoiCom.popoi[4]*c56 - popoiCom.popoi[5]*c65)*4.32e-13;


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