/*LevelN compute an arbitrary N level atom */
#include <string.h>
#include "cddefines.h"
#include "physconst.h"
#include "poplevls.h"
#include "edsqte.h"
#include "phycon.h"
#include "trace.h"
#include "typmatrx.h"
#include "tcool.h"
#include "sexp.h"
#include "veclib.h"
#include "addotslin.h"
#include "matin1.h"
#include "linpack.h"
#include "leveln.h"

void LevelN(long int nlev, 
  long int ndim, 
  float abund, 
  double g[], 
  double ex[], 
  double p[], 
  double *data, 
  double *dest, 
  double *pump, 
  long int *ipdest, 
  double *cooltl, 
  double *coolder, 
  char *chLabel, 
  int *lgNegPop)
{
#define DATA(I_,J_)	(*(data+(I_)*(ndim)+(J_)))
#define DEST(I_,J_)	(*(dest+(I_)*(ndim)+(J_)))
#define PUMP(I_,J_)	(*(pump+(I_)*(ndim)+(J_)))
#define IPDEST(I_,J_)	(*(ipdest+(I_)*(ndim)+(J_)))

	long int i, 
	  ihi, 
	  ilo, 
	  ipiv[LIMLEVELN - 1], 
	  j, 
	  job, 
	  ner;

	double cool1, 
	  crate[LIMLEVELN][LIMLEVELN], 
	  excit[LIMLEVELN][LIMLEVELN], 
	  ots;

	double amat[LIMLEVELN - 1][LIMLEVELN - 1], 
	  bvec[LIMLEVELN - 1], 
	  rcond, 
	  work[LIMLEVELN - 1], 
	  zz[LIMLEVELN][LIMLEVELN];

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

	/* nlev is the number of levels to compute, ndim is dimension of array */
	/* following used for linpac matrix inversion */

	/* PARAMETER statement:
	 * LIMLEVELN is upper limit to number of levels allowed by present dim
	 *
	 * entering variables
	 * NLEV is number of present dimensions
	 * ABUND is total abundance of species, used for nth equation
	 * G(ndim) is stat weight of levels
	 * EX(ndim) is excitation potential of levels, deg K
	 * 0 for first one, NOT d(ENER), but energy rel to ground
	 * DATA(ihi,ilo) is transit prob from up to low
	 * DATA(ilo,ihi) is col str from up to low
	 * lgNegPop is true is negative populations occurred
	 *
	 * deduced variables
	 * EXCIT(IHI,ILO) excitation Boltzman factor
	 * CRATE(IHI,ILO) collis deexcit rate
	 * CRATE(ILO,IHI) collis excit rate
	 *
	 * free statements >=20
	 *
	 * check if too many levels used
	 * >>chng 96 aug 22, had been .gt., off by one at limit since z(n,n+1) used */

	if( ndim >= LIMLEVELN )
	{
		fprintf( ioQQQ, " LevelN called with too many levels, the limit is%5ld\n", 
		  LIMLEVELN );
		puts( "[Stop in leveln]" );
		exit(1);
	}

	/*begin sanity checks
	 * excitation temperature of lowest level must be zero */
	if( ex[0] != 0. )
	{
		fprintf( ioQQQ, " LevelN called with non-zero ex(1)\n" );
		puts( "[Stop in leveln]" );
		exit(1);
	}

	/* both dest and pump (low, hi) are not used, should be zero */
	for( ihi=1; ihi < nlev; ihi++ )
	{
		for( ilo=0; ilo < ihi; ilo++ )
		{
			/* zero ots destruction rate */
			if( DEST(ihi,ilo) != 0. )
			{
				fprintf( ioQQQ, 
					" LevelN called with non-zero dest(l,u), atom=%s \n", chLabel);
			}

			if( PUMP(ihi,ilo) != 0. )
			{
				fprintf( ioQQQ, 
					" LevelN called with non-zero pump(l,u), atom=%s \n", chLabel);
			}

			if( IPDEST(ihi,ilo) != 0 )
			{
				fprintf( ioQQQ, 
					" LevelN called with non-zero ipdest(l,u), atom=%s \n", chLabel);
			}
		}
	}

	/* data(i,i) should be zero */
	for( ihi=0; ihi < nlev; ihi++ )
	{
		if( DATA(ihi,ihi) != 0 )
		{
			fprintf( ioQQQ, " LevelN called with non-zero data(i,i), atom=%s \n", chLabel);
		}
	}
	/*end sanity checks */

	if( trace.lgTrace && trace.lgTrLevN )
	{
		fprintf( ioQQQ, " LevelN trace printout for atom=%s \n", chLabel);
		fprintf( ioQQQ, " dest\n" );

		fprintf( ioQQQ, "  hi  lo" );

		for( ilo=1; ilo <= nlev; ilo++ )
		{
			fprintf( ioQQQ, "%4ld", ilo );
		}
		fprintf( ioQQQ, "      \n" );

		for( ihi=0; ihi < nlev; ihi++ )
		{
			fprintf( ioQQQ, "%3ld", ihi );
			for( ilo=0; ilo < nlev; ilo++ )
			{
				fprintf( ioQQQ, "%10.2e", DEST(ilo,ihi) );
			}
			fprintf( ioQQQ, "\n" );
		}

		fprintf( ioQQQ, " data\n" );

		fprintf( ioQQQ, "  hi  lo" );
		for( ilo=1; ilo <= nlev; ilo++ )
		{
			fprintf( ioQQQ, "%4ld", ilo );
		}
		fprintf( ioQQQ, "      \n" );

		for( ihi=0; ihi < nlev; ihi++ )
		{
			fprintf( ioQQQ, "%3ld", ihi );
			for( ilo=0; ilo < nlev; ilo++ )
			{
				fprintf( ioQQQ, "%10.2e", DATA(ilo,ihi) );
			}
			fprintf( ioQQQ, "\n" );
		}

		fprintf( ioQQQ, " pump\n" );

		fprintf( ioQQQ, "  hi  lo" );
		for( ilo=1; ilo <= nlev; ilo++ )
		{
			fprintf( ioQQQ, "%4ld", ilo );
		}
		fprintf( ioQQQ, "      \n" );

		for( ihi=0; ihi < nlev; ihi++ )
		{
			fprintf( ioQQQ, "%3ld", ihi );
			for( ilo=0; ilo < nlev; ilo++ )
			{
				fprintf( ioQQQ, "%10.2e", PUMP(ilo,ihi) );
			}
			fprintf( ioQQQ, "\n" );
		}
	}

	/* check for zero abundance and exit if so */
	if( abund <= 0. )
	{
		*cooltl = 0.;
		*coolder = 0.;
		*lgNegPop = FALSE;

		for( i=0; i < nlev; i++ )
		{
			p[i] = 0.;
			PopLevls.PopLevels[i] = 0.;
			PopLevls.DepLTELevels[i] = 0.;
		}

		PopLevls.PopLevels[0] = (float)abund;
		PopLevls.DepLTELevels[0] = 1.;
		/* there are TWO abort returns in this sub */
		
#	ifdef DEBUG_FUN
		fputs( " <->LevelN()\n", debug_fp );
#	endif
		return;
	}

	for( ilo=0; ilo < (nlev - 1); ilo++ )
	{
		for( ihi=ilo + 1; ihi < nlev; ihi++ )
		{
			excit[ilo][ihi] = sexp((ex[ihi]-ex[ilo])/phycon.te);
		}
	}

	if( trace.lgTrace && trace.lgTrLevN )
	{
		fprintf( ioQQQ, " excit, te=%10.2e\n", phycon.te );
		fprintf( ioQQQ, "  hi  lo" );

		for( ilo=1; ilo <= nlev; ilo++ )
		{
			fprintf( ioQQQ, "%4ld", ilo );
		}
		fprintf( ioQQQ, "      \n" );

		for( ihi=0; ihi < nlev; ihi++ )
		{
			fprintf( ioQQQ, "%3ld", ihi );
			for( ilo=0; ilo < nlev; ilo++ )
			{
				fprintf( ioQQQ, "%10.2e", excit[ilo][ihi] );
			}
			fprintf( ioQQQ, "\n" );
		}
	}

	/* punt if total excitation rate is zero */
	if( excit[0][nlev-1] + PUMP(0,nlev-1) < 1e-13 )
	{
		*cooltl = 0.;
		*coolder = 0.;
		*lgNegPop = FALSE;

		for( i=1; i < nlev; i++ )
		{
			p[i] = 0.;
			/* these are off by one - lowest index is zero */
			PopLevls.PopLevels[i] = 0.;
			PopLevls.DepLTELevels[i] = 0.;
		}

		/* everything in ground */
		p[0] = (double)abund;
		PopLevls.PopLevels[0] = (float)abund;
		PopLevls.DepLTELevels[0] = 1.;
		
#		ifdef DEBUG_FUN
		fputs( " <->LevelN()\n", debug_fp );
#		endif
		return;
	}

	/* factor is 8.629E-6 / SQRTE * EDEN */
	for( ilo=0; ilo < (nlev - 1); ilo++ )
	{
		for( ihi=ilo + 1; ihi < nlev; ihi++ )
		{
			/* this is deexcitation rate
			 * factor is 8.629E-6 / SQRTE * EDEN */
			crate[ilo][ihi] = DATA(ihi,ilo)/g[ihi]*edsqte.cdsqte;
			/* this is excitation rate */
			crate[ihi][ilo] = crate[ilo][ihi]*g[ihi]/g[ilo]*
			  excit[ilo][ihi];
		}
	}

	/* rate equations equal zero */
	for( i=0; i < nlev; i++ )
	{
		for( j=1; j < nlev; j++ )
		{
			zz[i][j] = 0.;
		}
		zz[i][0] = 1.0;
	}

	/* following is colum of vector */
	for( i=0; i < nlev; i++ )
	{
		zz[nlev][i] = 0.;
	}

	/* zz(nlev,nlev+1) = abund
	 * solve the problem with everything normalzied to pop of unity
	 * for numerical stbility
	 * zz(nlev,nlev+1) = 1. */
	zz[nlev][0] = 1.;

	/* eqns for destruction of level
	 * data(hi,lo) are A's, crate is coll excit in direction
	 * doi=1,nlev-1 */
	for( i=1; i < nlev; i++ )
	{
		zz[i][i] = 0.;

		/* leaving I to below */
		for( ilo=0; ilo < i; ilo++ )
		{
			zz[i][i] += crate[ilo][i] + DATA(ilo,i) + DEST(ilo,i) + 
			  PUMP(ilo,i)*g[ilo]/g[i];
			/* >>chng 97 jul 31, added pumping down
			 * coming to level I from below */
			zz[ilo][i] = -crate[i][ilo] - PUMP(ilo,i);
		}

		for( ihi=i + 1; ihi < nlev; ihi++ )
		{
			/* leaving I to above */
			zz[i][i] += crate[ihi][i] + PUMP(i,ihi);
			/* coming to level I from above */
			zz[ihi][i] = -crate[i][ihi] - DATA(i,ihi) - DEST(i,ihi) - 
			  PUMP(i,ihi)*g[i]/g[ihi];
			/* >>chng 97 jul 31, added pumping down */
		}
	}

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

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

		DGETRF(nlev,nlev,(double*)amat,LIMLEVELN-1,ipiv,&ner);
		/* usage DGETRS, 'N' = no transpose
		 * order of matrix,
		 * number of cols in bvec, =1
		 * array
		 * leading dim of array */
		DGETRS('N',nlev,1,(double*)amat,LIMLEVELN-1,ipiv,bvec,LIMLEVELN-1,&ner);

		if( ner != 0 )
		{
			fprintf( ioQQQ, " fe2 ir dgetrs error\n" );
			puts( "[Stop in leveln]" );
			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 < nlev; i++ )
		{
			zz[nlev][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 < nlev; j++ )
		{
			for( i=0; i < nlev; i++ )
			{
				amat[i][j] = zz[i][j];
			}
			bvec[j] = zz[nlev][j];
		}

		job = 0;
		rcond = 0.;
		/*666 error! in following one is dimension, second is number of levels */
		dgeco((double*)amat,5,5,ipiv,rcond,work);
		dgesl((double*)amat,5,5,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 < nlev; i++ )
		{
			zz[nlev][i] = bvec[i];
		}

		puts( "[Stop in leveln]" );
		exit(1);
	}

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

	/* set populations */
	for( i=0; i < nlev; i++ )
	{
		/* convert to real abundances */
		p[i] = (double)(zz[nlev][i]*abund);
	}

	/* now find total cooling and its derivative */
	*cooltl = 0.;
	*coolder = 0.;
	for( ihi=1; ihi < nlev; ihi++ )
	{
		for( ilo=0; ilo < ihi; ilo++ )
		{
			/* this is net cooling in oK */
			cool1 = (p[ilo]*crate[ihi][ilo] - p[ihi]*crate[ilo][ihi])*
			  (ex[ihi] - ex[ilo]);
			*cooltl += (float)cool1;
			/* simplistic way to handle line heating - better would
			 * be to go  over to method in level2 and level3 */
			*coolder += (float)(MAX2(0.,cool1)*((ex[ihi] - ex[0])*tcool.tsq1 - 
			  tcool.halfte));
		}
	}
	/* convert to ergs */
	*cooltl *= (float)BOLTZMANN;
	*coolder *= (float)BOLTZMANN;

	/* fill in departure coefficients */
	if( p[0] > 1e-20 && excit[0][nlev-1] > 1e-20 )
	{
		for( ihi=1; ihi < nlev; ihi++ )
		{
			/* these are off by one - lowest index is zero */
			PopLevls.PopLevels[ihi] = (float)p[ihi];
			PopLevls.DepLTELevels[ihi] = (float)((p[ihi]/p[0])*(g[0]/g[ihi])/
			  excit[0][ihi]);
		}
	}

	else
	{
		for( ihi=1; ihi < nlev; ihi++ )
		{
			/* these are off by one - lowest index is zero */
			PopLevls.PopLevels[ihi] = 0.;
			PopLevls.DepLTELevels[ihi] = 0.;
		}
	}

	PopLevls.PopLevels[0] = (float)abund;
	PopLevls.DepLTELevels[0] = 1.;

	/* in case atom has been trimed down */
	for( i=nlev; i < ndim; i++ )
	{
		p[i] = 0.;
		PopLevls.PopLevels[i] = 0.;
		PopLevls.DepLTELevels[i] = 0.;
	}

	/* sanity check for stability of inversion */
	*lgNegPop = FALSE;
	for( i=0; i < nlev; i++ )
	{
		if( p[i] < 0. )
		{
			*lgNegPop = TRUE;
		}
	}

	if( *lgNegPop )
	{
		fprintf( ioQQQ, " LevelN found negative population, atom=%4.4s\n", 
		  chLabel );
		fprintf( ioQQQ, " absolute population=" );

		for( i=0; i < nlev; i++ )
		{
			fprintf( ioQQQ, "%10.2e", p[i] );
		}

		fprintf( ioQQQ, "\n" );
		for( i=0; i < nlev; i++ )
		{
			p[i] = (double)MAX2(0.,p[i]);
		}
	}

	if( trace.lgTrace && trace.lgTrLevN )
	{
		fprintf( ioQQQ, " absolute population=" );
		for( i=0; i < nlev; i++ )
		{
			fprintf( ioQQQ, " %10.2e", p[i] );
		}
		fprintf( ioQQQ, "\n" );
	}


	for( ihi=1; ihi < nlev; ihi++ )
	{
		for( ilo=0; ilo < (ihi - 1); ilo++ )
		{
			/* set ots destruction rate */
			if( IPDEST(ilo,ihi) > 0 )
			{
				ots = p[ihi]*DEST(ilo,ihi);
				AddOTSLin((float)ots,IPDEST(ilo,ihi));
			}
		}
	}


#	ifdef DEBUG_FUN
	fputs( " <->LevelN()\n", debug_fp );
#	endif
	return;
#undef	IPDEST
#undef	PUMP
#undef	DEST
#undef	DATA
}