/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*HydroLevelPop solve for ionization balance level populations of model hydrogen atom */
#include "cddefines.h"
#include "physconst.h"
#include "taulines.h"
#include "iso.h"
#include "secondaries.h"
#include "elementnames.h"
#include "atmdat.h"
#include "rfield.h"
#include "dense.h"
#include "trace.h"
#include "lapack.h"
#include "hydrogenic.h"
#include "dynamics.h"
/*lint -e662 Possible access of out-of-bounds pointer*/

static double **z/*[iso.numLevels[ipH_LIKE][nelem]+2][iso.numLevels[ipH_LIKE]+2]*/;

static void PrtHydroTrace2( long int nelem )
{
	long ipHi , j , n;
	if( nelem == ipHYDROGEN )
	{
		double sum=0 , HRateDestGnd[10];
		for(ipHi=1; ipHi<iso.numLevels[ipH_LIKE][nelem]; ++ipHi )
		{
			sum += EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].pump ;
		}
		fprintf(ioQQQ,"DEBUG pump\t%.2f", fnzone);
		for(ipHi=1; ipHi<iso.numLevels[ipH_LIKE][nelem]; ++ipHi )
		{
			fprintf(ioQQQ,"\t%.3f",EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].pump/SDIV(sum)) ;
		}
		fprintf(ioQQQ,"\n");
		ipHi = iso.numLevels[ipH_LIKE][nelem]-2;
		fprintf(ioQQQ,"DEBUG pumphi\t%.2f\t%.3e\t%.3e\t%.3e\n",
			fnzone,
			EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].pump,
			rfield.OccNumbIncidCont[EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].ipCont-1],
			rfield.OccNumbDiffCont[EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].ipCont-1]);
		/* remember ground state destruction rate */
		HRateDestGnd[0] = iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];
		HRateDestGnd[1] = iso.gamnc[ipH_LIKE][nelem][ipH1s]/
			iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];
		HRateDestGnd[2] = iso.ColIoniz[ipH_LIKE][nelem][ipH1s]*
			dense.EdenHCorr/iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];
		HRateDestGnd[3] = secondaries.csupra[nelem][nelem]/iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];
		HRateDestGnd[4] = secondaries.Hx12[MIN2(nelem,1)][ipH2p]*
			9./iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];
		HRateDestGnd[5] = atmdat.HCharExcIonTotal/
			iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];
		HRateDestGnd[6] = atmdat.HCharExcIonTotal;
		HRateDestGnd[7] = sum/iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];
		HRateDestGnd[8] = sum;

		fprintf(ioQQQ," grnd dest fracs\t%.2f",fnzone);
		for( j=ipH1s; j < 9; j++ )
		{
			fprintf( ioQQQ,"\t");
			fprintf( ioQQQ,PrintEfmt("%8.1e", HRateDestGnd[j] ));
		}
		fprintf(ioQQQ,"\thcoldc\t%e\n", dense.EdenHCorr );
		fprintf( ioQQQ, "\n" );
	}

	fprintf( ioQQQ, "  pop level     others => (HydroLevelPop)\n" );
	for( n=ipH1s; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
	{
		fprintf( ioQQQ, "       HII%2ld", n );
		for( j=ipH1s; j < iso.numLevels[ipH_LIKE][nelem]; j++ )
		{
			fprintf( ioQQQ," ");
			/*fprintf( ioQQQ,PrintEfmt("%8.1e", z[j][n] ) );*/
			fprintf( ioQQQ,"%.9e\t", z[j][n] );
		}
		fprintf( ioQQQ, "\n" );
	}
}

void HydroLevelPop(long int nelem)
{
	long int n, 
	  ipHi, 
	  ipLo, 
	  i,
	  j, 
	  level, 
	  level_error;
	double *amat;

	double BigError ;

	int32 nerror;
	int32 *ipiv ; /* malloc out to [iso.numLevels[ipH_LIKE][nelem]+1] */

	double 
		*creation,
		*Save_creation,
		**SaveZ/*[iso.numLevels[ipH_LIKE][nelem]+2][iso.numLevels[ipH_LIKE]+2]*/, 
		*work/*[iso.numLevels[ipH_LIKE][nelem]+2]*/;
	double *error;/*[iso.numLevels[ipH_LIKE][nelem]+2]*/

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

	/* check that we were called with valid charge */
	ASSERT( nelem >= 0 );
	ASSERT( nelem < LIMELM );

	if( (ipiv = (int32 *)MALLOC(sizeof(int32)*(unsigned)(iso.numLevels[ipH_LIKE][nelem]) ) ) == NULL )
		BadMalloc();
	if( (creation = (double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels[ipH_LIKE][nelem]) ) ) == NULL )
		BadMalloc();
	if( (Save_creation = (double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels[ipH_LIKE][nelem]) ) ) == NULL )
		BadMalloc();
	if( (work = (double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels[ipH_LIKE][nelem]) ) ) == NULL )
		BadMalloc();
	if( (error = (double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels[ipH_LIKE][nelem]) ) ) == NULL )
		BadMalloc();

	/* now do the 2D arrays */
	if( (SaveZ = (double **)MALLOC(sizeof(double *)*(unsigned)(iso.numLevels[ipH_LIKE][nelem]) ) ) ==NULL )
		BadMalloc();

	if(  (z = (double **)MALLOC(sizeof(double *)*(unsigned)(iso.numLevels[ipH_LIKE][nelem]) ))==NULL  )
		BadMalloc();

	/* now do the second dimension */
	for( i=0; i<(iso.numLevels[ipH_LIKE][nelem]); ++i )
	{
		if( (SaveZ[i] = (double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels[ipH_LIKE][nelem]) ) ) ==NULL )
			BadMalloc();

		if( (z[i] = (double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels[ipH_LIKE][nelem]) ) ) ==NULL )
			BadMalloc();
	}


	/* fill in recombination vector - values were set in iso_ionize_recombine.c */
	for( level=ipH1s; level < iso.numLevels[ipH_LIKE][nelem]; ++level )
	{
		/* units are s-1 */
		creation[level] = iso.RateCont2Level[ipH_LIKE][nelem][level];
	}

	/* master balance equation */
	for( level=ipH1s; level < iso.numLevels[ipH_LIKE][nelem]; level++ )
	{
		/* all process depopulating level */
		/* this includes grain charge transfer ionization */
		z[level][level] = iso.RateLevel2Cont[ipH_LIKE][nelem][level];

		/* all processes populating level from below */
		for( ipLo=ipH1s; ipLo < level; ipLo++ )
		{
			double CollRate = EmisLines[ipH_LIKE][nelem][level][ipLo].ColUL* dense.EdenHCorr;

			z[ipLo][level] = 
				-CollRate * 
				(double)iso.stat[ipH_LIKE][nelem][level]/(double)iso.stat[ipH_LIKE][nelem][ipLo]*
				iso.Boltzmann[ipH_LIKE][nelem][level][ipLo]  -
				EmisLines[ipH_LIKE][nelem][level][ipLo].pump *
				KillIfBelowPlasma(EmisLines[ipH_LIKE][nelem][level][ipLo].EnergyWN/WAVNRYD);

			/* pumping out of here to lower level */
			z[level][level] += (double)EmisLines[ipH_LIKE][nelem][level][ipLo].pump *
				(double)iso.stat[ipH_LIKE][nelem][ipLo]/(double)iso.stat[ipH_LIKE][nelem][level] *
				KillIfBelowPlasma(EmisLines[ipH_LIKE][nelem][level][ipLo].EnergyWN/WAVNRYD);

			/* collisions out of here to lower level */
			z[level][level] += CollRate;

			/* radiative decays out of here to lower level */
			z[level][level] += 
				EmisLines[ipH_LIKE][nelem][level][ipLo].Aul*
				(EmisLines[ipH_LIKE][nelem][level][ipLo].Pesc + 
				EmisLines[ipH_LIKE][nelem][level][ipLo].Pelec_esc +
				EmisLines[ipH_LIKE][nelem][level][ipLo].Pdest) *
				KillIfBelowPlasma(EmisLines[ipH_LIKE][nelem][level][ipLo].EnergyWN/WAVNRYD);
		}

		/* all processes populating level from above */
		for( ipHi=level + 1; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
		{
			double RadDecay , CollRate ;
			RadDecay =
				EmisLines[ipH_LIKE][nelem][ipHi][level].Aul*
				(EmisLines[ipH_LIKE][nelem][ipHi][level].Pesc + 
				EmisLines[ipH_LIKE][nelem][ipHi][level].Pelec_esc +
				EmisLines[ipH_LIKE][nelem][ipHi][level].Pdest);

			/* >>chng 02 feb 06, define this and use in each direction below -
			 * attempt at solving roundoff problems on alphas */
			CollRate = EmisLines[ipH_LIKE][nelem][ipHi][level].ColUL *dense.EdenHCorr;

			z[ipHi][level] = 
				-( RadDecay + 
				(double)EmisLines[ipH_LIKE][nelem][ipHi][level].pump *
				(double)iso.stat[ipH_LIKE][nelem][level]/(double)iso.stat[ipH_LIKE][nelem][ipHi] ) *
				KillIfBelowPlasma(EmisLines[ipH_LIKE][nelem][ipHi][level].EnergyWN/WAVNRYD) - 
				CollRate;

			/* pumping out of here to upper level */
			z[level][level] += EmisLines[ipH_LIKE][nelem][ipHi][level].pump *
				KillIfBelowPlasma(EmisLines[ipH_LIKE][nelem][ipHi][level].EnergyWN/WAVNRYD);

			/* collisions out of here to upper level */
			z[level][level] += (double)iso.stat[ipH_LIKE][nelem][ipHi] / (double)iso.stat[ipH_LIKE][nelem][level] *
				iso.Boltzmann[ipH_LIKE][nelem][ipHi][level]*CollRate;
		}
	}

	/* >>chng 02 jul 22, add induced 2-nu to level populations */
	/* induced two photon emission - special because upward and downward are
	 * not related by ratio of statistical weights */
	/* iso.lgInd2nu_On is controlled with SET IND2 ON/OFF command */
	/*TODO	0	 this does not look right - check logic */
	z[ipH2s][ipH1s] -= iso.TwoNu_induc_dn[ipH_LIKE][nelem]*iso.lgInd2nu_On ;
	z[ipH1s][ipH2s] -= iso.TwoNu_induc_up[ipH_LIKE][nelem]*iso.lgInd2nu_On ;

	/* rates out of 1s, and out of 2s */
	z[ipH1s][ipH1s] += iso.TwoNu_induc_up[ipH_LIKE][nelem]*iso.lgInd2nu_On ;
	z[ipH2s][ipH2s] += iso.TwoNu_induc_dn[ipH_LIKE][nelem]*iso.lgInd2nu_On ;

	if( secondaries.Hx12[0][1] > 0. )
	{
		/* now add on supra thermal excitation */
		for( level=ipH2s; level < iso.numLevels[ipH_LIKE][nelem]; level++ )
		{
			double RateUp , RateDown;

			RateUp = secondaries.Hx12[MIN2(nelem,1)][level];

			RateDown = RateUp * (double)iso.stat[ipH_LIKE][nelem][ipH1s] /
				(double)iso.stat[ipH_LIKE][nelem][level];

			/* stuff in min after Hx12 evaluates to 0 for atom, or 1 for ion */
			/* total rate out of lower level */
			z[ipH1s][ipH1s] += RateUp;

			/* rate from the upper level to ground */
			z[level][ipH1s] -= RateDown ;

			/* rate from ground to upper level */
			z[ipH1s][level] -= RateUp ;

			z[level][level] += RateDown;  
		}
	}


	/* debug oscillations in h=2 population */
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && nelem==ipHYDROGEN && nzone > 600 )
		{
			fprintf(ioQQQ, "DEBUG Hn=2\t%.2f\t%.3e\t%.3e\t%.3e\n",
			fnzone,
			z[ipH1s][ipH1s],
			z[ipH2p][ipH2p],
			EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pdest );
		}
	}

#if FOO
	{ 
		float mytot,err;
		int myeq;

		for (myeq=ipH1s; myeq<iso.numLevels[ipH_LIKE][nelem]; myeq++ ) {
			mytot = 0.;
			for (level=ipH1s; level<iso.numLevels[ipH_LIKE][nelem]; level++ ) {
				mytot += z[myeq][level];
			}
			err = fabs((mytot-iso.RateLevel2Cont[ipH_LIKE][nelem][myeq])/z[myeq][myeq]);
			if (err > 1e-6)
				printf("BalChk: %d %d e %g\n",nelem,myeq,err);
		}
	}
#endif

	/* >>chng 02 Sep 06 rjrw -- all elements have these terms */
	/*>>>chng 02 oct 01, only include if lgAdvection is set */
	if( dynamics.lgAdvection )
	{
		/* add in advection - these terms normally zero */
		/* assume for present that all advection is into ground state */
		creation[ipH1s] += 
			dynamics.Source[nelem][nelem]/SDIV(dense.xIonDense[nelem][nelem+1])*
			dynamics.lgISO[ipH_LIKE];
		/* >>chng 02 Sep 06 rjrw -- advective term not recombination */
		/* can sink from all components (must do, for conservation) */
		for( ipLo=ipH1s; ipLo<iso.numLevels[ipH_LIKE][nelem]; ++ipLo )
		{
			z[ipLo][ipLo] += dynamics.Rate*dynamics.lgISO[ipH_LIKE];
		}
	}

	/* =================================================================== 
	 *
	 * at this point all matrix elements have been established 
	 *
	 * ==================================================================== */

	if( (trace.lgTrace && trace.lgIsoTraceFull[ipH_LIKE]) && (nelem == trace.ipIsoTrace[ipH_LIKE]) )
		PrtHydroTrace2(nelem);

	/* save matrix */
	for( j=ipH1s; j < iso.numLevels[ipH_LIKE][nelem]; j++ )
	{
		for( n=ipH1s; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
		{
			SaveZ[n][j] = z[n][j];
		}
		Save_creation[j] = creation[j];
	}

	/*double amat[LMHLVL+2][LMHLVL+2];*/
	/* iso.numLevels[ipH_LIKE][nelem]+1 is right dimension of matrix */
#	ifdef AMAT
#		undef AMAT
#	endif
#	define AMAT(I_,J_)	(*(amat+(I_)*iso.numLevels[ipH_LIKE][nelem]+(J_)))
	/* MALLOC space for the  1-d array */
	if( (amat=(double*)MALLOC( (sizeof(double)*(unsigned)(iso.numLevels[ipH_LIKE][nelem]*(iso.numLevels[ipH_LIKE][nelem])) ))) == NULL )
		BadMalloc();

	/* this one may be more robust */
	for( j=ipH1s; j < iso.numLevels[ipH_LIKE][nelem]; j++ )
	{
		for( n=ipH1s; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
		{
			/*amat[n][j] = z[n][j];*/
			AMAT(n,j) = z[n][j];
		}
	}

	/*DGETRF(iso.numLevels[ipH_LIKE][nelem],iso.numLevels[ipH_LIKE][nelem],
		amat,(iso.numLevels[ipH_LIKE][nelem]),ipiv,&nerror1);

	DGETRS('N',iso.numLevels[ipH_LIKE][nelem],1,amat,(iso.numLevels[ipH_LIKE][nelem]),ipiv,
		bvec,(iso.numLevels[ipH_LIKE][nelem]),&nerror2);*/

	nerror = 0;

	getrf_wrapper(iso.numLevels[ipH_LIKE][nelem],iso.numLevels[ipH_LIKE][nelem],
		      amat,iso.numLevels[ipH_LIKE][nelem],ipiv,&nerror);


	getrs_wrapper('N',iso.numLevels[ipH_LIKE][nelem],1,amat,iso.numLevels[ipH_LIKE][nelem],ipiv,
		      creation,iso.numLevels[ipH_LIKE][nelem],&nerror);

	/*DGETRF(iso.numLevels[ipH_LIKE][nelem],iso.numLevels[ipH_LIKE][nelem],
		amat,iso.numLevels[ipH_LIKE][nelem],ipiv,&nerror1);

	DGETRS('N',iso.numLevels[ipH_LIKE][nelem],1,amat,iso.numLevels[ipH_LIKE][nelem],ipiv,
		creation,iso.numLevels[ipH_LIKE][nelem],&nerror2);*/

	if( nerror != 0 )
	{
		fprintf( ioQQQ, " HydroLevelPop: dgetrs finds singular or ill-conditioned matrix\n" );
		puts( "[Stop in HydroLevelPop]" );
		cdEXIT(EXIT_FAILURE);
	}

	free(amat);

	/* check whether solution is valid */
	for( level=ipH1s; level < iso.numLevels[ipH_LIKE][nelem]; level++ )
	{
		double BigSoln= 0.;
		error[level] = 0.;
		for( n=ipH1s; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
		{
			/* remember the largest value of the soln matrix to div by below */
			if( fabs(SaveZ[n][level] ) > BigSoln )
				BigSoln = fabs(SaveZ[n][level]);

			error[level] += SaveZ[n][level]*creation[n];
		}

		if( BigSoln > 0. )
		{
			error[level] = (error[level] - Save_creation[level])/ BigSoln;
		}
		else
		{
			error[level] = 0.;
		}
	}

	/* remember largest residual in matrix inversion */
	BigError = -1.;
	level_error = -1;
	for( level=ipH1s; level < iso.numLevels[ipH_LIKE][nelem]; level++ )
	{
		double abserror;
		abserror = fabs( error[level]) ;
		/* this will be the largest residual in the matrix inversion */
		if( abserror > BigError )
		{
			BigError = abserror ;
			level_error = level;
		}
	}

	/* matrix inversion should be nearly as good as the accuracy of a double,
	 * but demand that it is better than epsilon for a float */
	if( BigError > FLT_EPSILON ) 
	{
		fprintf(ioQQQ,
			" PROBLEM  HydroLevelPop zone %.2f - largest residual in hydrogenic %s nelem=%li matrix inversion is %g "
			"level was %li \n", 
			fnzone,
			elementnames.chElementName[nelem],
			nelem , 
			BigError , 
			level_error);
		/* >>chng 04 jul 18, do not abort, do not even do the show me
		 * the code had done showme, but the cdexit was commented out for some time */
		/*ShowMe();
		puts( "[Stop in HydroLevelPop]" );*/
		/* cdEXIT(EXIT_FAILURE); */
	}

	/* convert from departure coef in Z(I) to level pop rel to HII */

	/* put departure coefficients and level populations into master array */
	for( ipLo=ipH1s; ipLo < iso.numLevels[ipH_LIKE][nelem]; ipLo++ )
	{
		iso.Pop2Ion[ipH_LIKE][nelem][ipLo] = creation[ipLo];
		if( iso.Pop2Ion[ipH_LIKE][nelem][ipLo] <= 0 )
		{
			fprintf(ioQQQ," non-positive level pop for H iso, nelem = %li = %s, level=%li val=%.3e\n",
				nelem , 
				elementnames.chElementSym[nelem],
				ipLo,
				iso.Pop2Ion[ipH_LIKE][nelem][ipLo] );
		}

		if( iso.PopLTE[ipH_LIKE][nelem][ipLo] > 0. )
		{
			iso.DepartCoef[ipH_LIKE][nelem][ipLo] = 
				(iso.Pop2Ion[ipH_LIKE][nelem][ipLo]/
				(iso.PopLTE[ipH_LIKE][nelem][ipLo]* dense.eden) );
		}
		else
		{
			iso.DepartCoef[ipH_LIKE][nelem][ipLo] = 0.;
		}
	}

	free( ipiv );
	free( work );
	free( creation );
	free( Save_creation );
	free( error );

	/* now free up the 2D arrays */
	for( i=0; i<iso.numLevels[ipH_LIKE][nelem]; ++i )
	{
		free( SaveZ[i] ) ;
		free( z[i] ) ;
	}
	free( SaveZ ) ;
	free( z ) ;

#	ifdef DEBUG_FUN
	fputs( " <->HydroLevelPop()\n", debug_fp );
#	endif
	return;
}
/*lint +e662 Possible access of out-of-bounds pointer*/

