/*beseq compute level populations and emissivity for Be-sequence ions */
#include <string.h>
#include "cddefines.h"
#include "taulines.h"
#include "ionfracs.h"
#include "edsqte.h"
#include "poplevls.h"
#include "phycon.h"
#include "typmatrx.h"
#include "tcool.h"
#include "linpack.h"
#include "chionlbl.h"
#include "addotslin.h"
#include "coladd.h"
#include "matin1.h"
#include "veclib.h"
#include "beseq.h"

void beseq(double cs12, 
  double cs13, 
  double cs23, 
  EmLine * t, 
  double a30)
{

	char chLab[5];

	int lgNegPop;

	long int i, 
	  ipiv[3-(0)+1], 
	  j, 
	  job, 
	  ner, 
	  nerror;

	double AbunxIon, 
	  Enr01, 
	  Enr10, 
	  a20, 
	  aem, 
	  boltz, 
	  c01, 
	  c02, 
	  c03, 
	  c10, 
	  c12, 
	  c13, 
	  c20, 
	  c21, 
	  c23, 
	  c30, 
	  c31, 
	  c32, 
	  coolng, 
	  cs1u, 
	  dest, 
	  excit, 
	  r01, 
	  r02, 
	  r03, 
	  r10, 
	  r12, 
	  r13, 
	  r20, 
	  r21, 
	  r23, 
	  r30, 
	  r31, 
	  r32, 
	  rcond, 
	  ri02, 
	  ri20, 
	  tot20;

	double amat[3-(0)+1][3-(0)+1], 
	  bvec[3-(0)+1], 
	  work[3-(0)+1], 
	  zz[4-(0)+1][4-(0)+1];

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

	/* function returns total emission in both components of line
	 * destruction by background opacity computed and added to otslin field
	 * here,
	 * total cooling computed and added via coladd
	 *
	 * finds population of four level Be-like atom
	 * level 1 = ground, J=0
	 * levels 2,3,4 are J=0,1,2, OF 3P.
	 * levels 3-1 is the fast one, j=1 to ground j=0
	 * cs1u is coll str to all J sub-levels of 3P; C.S. goes as 2J+1*/

	/* beseq(cs12,cs13,cs23,t->t,a30,chLab)
	 * xIonFracs(nelem,i) is density of ith ionization stage (cm^-3) */
	AbunxIon = xIonFracs[t->IonStg][t->nelem -1];

	/* put null terminated line label into chLab using line array*/
	chIonLbl(chLab, t) ;
	boltz = t->EnergyK/phycon.te;

	/* low density cutoff to keep matrix happy */
	if( AbunxIon <= 1e-20 || boltz > 30. )
	{
		/* put in pop since possible just too cool */
		t->PopLo = (float)AbunxIon;
		t->PopOpc = (float)AbunxIon;
		t->PopHi = 0.;
		t->xIntensity = 0.;
		t->cool = 0.;
		t->phots = 0.;
		t->ots = 0.;
		t->AovTot = 0.;
		t->ColOvTot = 0.;
		t->heat = 0.;
		/* level populations */
		PopLevls.PopLevels[0] = (float)AbunxIon;
		PopLevls.PopLevels[1] = 0.;
		PopLevls.PopLevels[2] = 0.;
		PopLevls.PopLevels[3] = 0.;
		PopLevls.DepLTELevels[0] = 1.;
		PopLevls.DepLTELevels[1] = 0.;
		PopLevls.DepLTELevels[2] = 0.;
		PopLevls.DepLTELevels[3] = 0.;
		coladd(chLab, (long)t->WLAng ,0.);
		
#		ifdef DEBUG_FUN
		fputs( " <->beseq()\n", debug_fp );
#		endif
		return;
	}
	excit = exp(-boltz);
	cs1u = t->cs;

	/* these must be the statistical weights */
	assert( t->gLo == 1. );
	assert( t->gHi == 3. );

	/* collision strength must be positive */
	assert( cs1u > 0. );

	/* incuded rates for fastest transition */
	ri02 = t->pump;

	/* back reaction has ratio of stat weights */
	ri20 = t->pump*1./3.;

	/* net rate out of level 3, with destruction */
	a20 = t->Aul*(t->Pesc + t->Pdest);
	tot20 = a20 + ri20;

	/* 1/9 is ratio of level stat weight to term stat weight */
	c10 = cs1u*edsqte.cdsqte/9.;
	c01 = c10*excit;
	r01 = c01;
	r10 = c10;

	/* stat weights canceled out here */
	c20 = c10;
	c02 = c01*3.;
	r02 = c02 + ri02;
	r20 = c20 + tot20;

	c30 = c10;
	c03 = c01*5.;
	r30 = c30 + a30;
	r03 = c03;

	c21 = cs12*edsqte.cdsqte/3.;
	c12 = c21*3.;
	r12 = c12;
	r21 = c21;

	c31 = cs13*edsqte.cdsqte/5.;
	c13 = c31*5.;
	r31 = c31;
	r13 = c13;

	c32 = cs23*edsqte.cdsqte/5.;
	c23 = c32*1.667;
	r32 = c32;
	r23 = c23;

	/* set up matrix */
	for( i=0; i <= 3; i++ )
	{
		/* first equation will be sum to abund */
		zz[i][0] = 1.;
		zz[4][i] = 0.;
	}

	/* zz(0,4) = AbunxIon */
	zz[4][0] = 1.;

	/* ground level 0 is implicit
	 * level 1 balance equation */
	zz[0][1] = -r01;
	zz[1][1] = r10 + r12 + r13;
	zz[2][1] = -r21;
	zz[3][1] = -r31;

	/* level 2 balance equation */
	zz[0][2] = -r02;
	zz[1][2] = -r12;
	zz[2][2] = r20 + r21 + r23;
	zz[3][2] = -r32;

	/* level 3 balance equation */
	zz[0][3] = -r03;
	zz[1][3] = -r13;
	zz[2][3] = -r23;
	zz[3][3] = r30 + r31 + r32;

	/* which matrix solution? */
	if( strcmp(TypMatrx.chMatrix,"matin1 ") == 0 )
	{
		/* this is the usual matrix inversion method, slightly faster */
		/*matin1();*/
		ner = 0;
		if( ner != 0 )
		{
			fprintf( ioQQQ, " beseq matin1 matrix error\n" );
			puts( "[Stop in beseq]" );
			exit(1);
		}
	}

	else if( strcmp(TypMatrx.chMatrix,"linpack") == 0 )
	{
		/* this one may be more robust */
		for( j=0; j <= 3; j++ )
		{
			for( i=0; i <= 3; i++ )
			{
				amat[i][j] = zz[i][j];
			}
			bvec[j] = zz[3+1][j];
		}
		DGETRF(4,4,(double*)amat,4,ipiv,&nerror);
		DGETRS('N',4,1,(double*)amat,4,ipiv,bvec,4,&nerror);
		if( nerror != 0 )
		{
			fprintf( ioQQQ, " beseq dgetrs error\n" );
			puts( "[Stop in beseq]" );
			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 <= 3; i++ )
		{
			zz[3+1][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 <= 3; j++ )
		{
			for( i=0; i <= 3; i++ )
			{
				amat[i][j] = zz[i][j];
			}
			bvec[j] = zz[3+1][j];
		}
		job = 0;
		rcond = 0.;
		dgeco((double*)amat,4,4,ipiv,rcond,work);
		dgesl((double*)amat,4,4,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 <= 3; i++ )
		{
			zz[3+1][i] = bvec[i];
		}
		puts( "[Stop in beseq]" );
		exit(1);
	}

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

	lgNegPop = FALSE;
	for( i=0; i <= 3; i++ )
	{
		PopLevls.PopLevels[i] = (float)(zz[4][i]*AbunxIon);
		if( PopLevls.PopLevels[i] < 0. )
			lgNegPop = TRUE;
	}

	/* check for negative level populations, this would be a major error */
	if( lgNegPop )
	{
		fprintf( ioQQQ, " beseq finds non-positive pop,=" );
		for( i=0; i <= 3; i++ )
		{
			fprintf( ioQQQ, "%g ", PopLevls.PopLevels[i] );
		}
		fprintf( ioQQQ, "%s \n", chLab );
		fprintf( ioQQQ, " te=%g  total abund=%g  boltz=%g \n", 
		  phycon.te, AbunxIon, boltz );
		puts( "[Stop in beseq]" );
		exit(1);
	}

	/* convert level populations over to departure coeficients relative to ground */
	PopLevls.DepLTELevels[0] = 1.;
	PopLevls.DepLTELevels[1] = (float)((PopLevls.PopLevels[1]/PopLevls.PopLevels[0])/
	  excit);
	PopLevls.DepLTELevels[2] = (float)((PopLevls.PopLevels[2]/PopLevls.PopLevels[0])/
	  (excit*3.));
	PopLevls.DepLTELevels[3] = (float)((PopLevls.PopLevels[3]/PopLevls.PopLevels[0])/
	  (excit*5.));

	/* compute ratio Aul/(Aul+Cul) */
	aem = t->Aul*t->Pesc;

	t->AovTot = (float)(aem/(aem + c20));

	t->ColOvTot = (float)(c02/r02);

	t->PopLo = (float)AbunxIon;

	/* >>chng 96 sep 12, ipLnPopu had not been set before */
	t->PopHi = PopLevls.PopLevels[2];

	t->PopOpc = (float)(AbunxIon - PopLevls.PopLevels[2]*1./3.);

	/* this will be escaping part of line
	 * number of escaping line photons, used elsewhere for outward beam */
	t->phots = PopLevls.PopLevels[2] * t->Aul * t->Pesc;

	t->xIntensity = t->phots * t->EnergyErg;

	/* follow is total cooling due to this line, inclues destroyed part
	 *
	 * following are correct heat exchange, it will mix to get correct deriv */
	Enr01 = PopLevls.PopLevels[0]*(c01 + c02 + c03)*t->EnergyErg;
	Enr10 = (PopLevls.PopLevels[1]*c10 + PopLevls.PopLevels[2]*c20 + 
	  PopLevls.PopLevels[3]*c30)*t->EnergyErg;

	/* net cooling due to excit minus part of de-excit */
	t->cool = (float)(Enr01 - Enr10*t->ColOvTot);

	/* net heating is remainder */
	t->heat = (float)(Enr10*(1. - t->ColOvTot));

	/* put line cooling into cooling stack */
	coolng = t->cool;
	coladd(chLab, (long)t->WLAng ,coolng);

	/* derivative of cooling function */
	tcool.dCooldT += coolng*(t->EnergyK*tcool.tsq1 - tcool.halfte);

	/* destroyed part of line */
	dest = PopLevls.PopLevels[2]*t->Aul*t->Pdest;
	t->ots = (float)dest;

	/* now add to ots field */
	AddOTSLin((float)dest , t->ipCont );

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