/*HydroLevelDep solve for ionization balance level populations of model hydrogen atom */
#include <string.h>
#include "cddefines.h"
#include "taulines.h"
#include "hcolrt.h"
#include "secondaries.h"
#include "hphoto.h"
#include "charexc.h"
#include "trace.h"
#include "phycon.h"
#include "hcbr.h"
#include "hstat.h"
#include "hcaseb.h"
#include "typmatrx.h"
#include "hratedest.h"
#include "elecden.h"
#include "linpack.h"
#include "matin1.h"
#include "veclib.h"
#include "printe82.h"
#include "hydrogenic.h"

void HydroLevelDep(long int ipZ ,
	double **SaveZ/*[nhlevel+2][nhlevel+2]*/, 
	double *bvec/*[nhlevel+2]*/, 
	double *error/*[nhlevel+2]*/, 
	double *work/*[nhlevel+2]*/, 
	double **z/*[nhlevel+2][nhlevel+2]*/ ,
	long int *ipiv , /* malloc out to [nhlevel+1] */
	double *totcap/* malloc out to [nhlevel+1]*/)
{
	char chType[5];
	long int n, 
	  ipHi, 
	  ipLo, 
	  j, 
	  job, 
	  level, 
	  nerror;

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

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

	/* 
	 * following electron density has approximate correction for neutrals
	 * corr of hi*1.7e-4 accounts for col ion by HI; Drawin Zs Phys 225, 483.
	 * used EdenHCorr instead
	 * edhi = eden + hi * 1.7e-4
	 */

	/* this is main solver */
	strcpy( chType, "mtrx" );
	hcaseb.HRecEffec[ipZ] = 0.;
	for( n=IP1S; n <= nhlevel; n++ )
	{
		z[nhlevel+1][n] = ((double)(hydro.hrec[ipZ][n][ipRecRad]*
		  hydro.hrec[ipZ][n][ipRecNetEsc]) + 
		  (double)(HPhoto.rinduc[ipZ][n]))/
		  hcolrt.hlte[ipZ][n] + (double)(hydro.hcolnc[ipZ][n]*
		  ElecDen.EdenHCorr);

		/* hlte(n,ipZ) is only LTE pop when mult by Ne Nh */
		totcap[n] = z[nhlevel+1][n]*hcolrt.hlte[ipZ][n];
		hcaseb.HRecEffec[ipZ] += (float)totcap[n];
	}

	/* add on charge transfer recombination into ground */
	z[nhlevel+1][IP1S] += CharExc.CTHrec[MIN2(ipZ+1,2)-1][1]/
	  hcolrt.hlte[ipZ][IP1S]/phycon.eden;

	totcap[IP1S] = z[nhlevel+1][IP1S]*hcolrt.hlte[ipZ][IP1S];

	/* master balance equation */
	for( level=IP1S; level <= nhlevel; level++ )
	{
		/* 
		 * all process depopulating level
		 * hradn(ipZ,level) is sum of spontaneous decays to ipLo levels
		 * hradn(ipZ,level) is sum of HRadNet, HRadNet is net A down
		 * hbul(ipZ,n) is sum of induced exits from level n, both up and down
		 */
		z[level][level] = HPhoto.hgamnc[ipZ][level] + 
			hcbr.hradn[ipZ][level] + 
		  (hcbr.hcoldn[ipZ][level] + 
		  hydro.hcolnc[ipZ][level])* ElecDen.EdenHCorr + 
		  hcbr.hbul[ipZ][level] + Secondaries.csupra;

		/* all processes populating level from below */
		for( ipLo=IP1S; ipLo <= (level - 1); ipLo++ )
		{
			z[ipLo][level] = -hydro.hcol[ipZ][ipLo][level]*
			  ElecDen.EdenHCorr - 
			  HydroLines[ipZ][level][ipLo].pump*
			  hstat.HStatWght[ipLo]/hstat.HStatWght[level]/
			  hydro.hlbolt[ipZ][ipLo][level];
		}

		/* all processes populating level from above */
		for( ipHi=level + 1; ipHi <= nhlevel; ipHi++ )
		{
			z[ipHi][level] = -(hydro.HRadNet[ipZ][level][ipHi] + 
			  HydroLines[ipZ][ipHi][level].pump*
			  hstat.HStatWght[level]/hstat.HStatWght[ipHi] + 
			  hydro.hcol[ipZ][level][ipHi]*ElecDen.EdenHCorr)*
			  hcolrt.hlte[ipZ][ipHi]/hcolrt.hlte[ipZ][level];
		}
	}
	if( Secondaries.Hx12[0][1] > 0. )
	{
		/* now add on supra thermal excitation */
		for( level=IP2S; level <= nhlevel; level++ )
		{
			/* stuff in min after Hx12 evaluates to 0 for atom, or 1 for atom */
			z[IP1S][IP1S] += Secondaries.Hx12[MIN2(ipZ,1)][level];

			z[level][IP1S] += -Secondaries.Hx12[MIN2(ipZ,1)][level]*hstat.HStatWght[IP1S]/
			  hstat.HStatWght[level]*hcolrt.hlte[ipZ][level]/
			  hcolrt.hlte[ipZ][IP1S];

			z[IP1S][level] += -Secondaries.Hx12[MIN2(ipZ,1)][level]*hstat.HStatWght[IP1S]/
			  hstat.HStatWght[level]/hydro.hlbolt[ipZ][IP1S][level];

			z[level][level] += Secondaries.Hx12[MIN2(ipZ,1)][level]*2./
			  hstat.HStatWght[level];
		}
	}

	/* charge transfer ionization */
	z[IP1S][IP1S] += CharExc.CTHion[MIN2(ipZ,1)][1];

	/* remember ground state destruction rate */
	if( ipZ == 0 )
	{
		HRateDest.HRateDestGnd[0] = (float)z[IP1S][IP1S];
		HRateDest.HRateDestGnd[1] = (float)(HPhoto.hgamnc[ipZ][IP1S]/
		  z[IP1S][IP1S]);
		HRateDest.HRateDestGnd[2] = (float)(hydro.hcolnc[ipZ][IP1S]*
		  ElecDen.EdenHCorr/z[IP1S][IP1S]);
		HRateDest.HRateDestGnd[3] = (float)(hcbr.hcoldn[ipZ][IP1S]*
		  ElecDen.EdenHCorr/z[IP1S][IP1S]);
		HRateDest.HRateDestGnd[4] = (float)(Secondaries.csupra/z[IP1S][IP1S]);
		HRateDest.HRateDestGnd[5] = (float)(Secondaries.Hx12[MIN2(ipZ,1)][IP2P]*
		  9./z[IP1S][IP1S]);
		HRateDest.HRateDestGnd[6] = (float)(hcbr.hbul[ipZ][IP1S]/z[IP1S][IP1S]);
		HRateDest.HRateDestGnd[7] = (float)(CharExc.CTHion[MIN2(ipZ,1)][1]/
		  z[IP1S][IP1S]);
	}
#if 0
	/* test for helium bug */
	if( ipZ == 1 )
	{
		HRateDest.HRateDestGnd[0] = (float)z[3][3];
		HRateDest.HRateDestGnd[1] = (float)(HPhoto.hgamnc[ipZ][3]/
		  z[3][3]);
		HRateDest.HRateDestGnd[2] = (float)(hydro.hcolnc[ipZ][3]*
		  ElecDen.EdenHCorr/z[3][3]);
		HRateDest.HRateDestGnd[3] = (float)(hcbr.hcoldn[ipZ][3]*
		  ElecDen.EdenHCorr/z[3][3]);
		HRateDest.HRateDestGnd[4] = (float)(hcbr.hradn[ipZ][3]/z[3][3]);
		HRateDest.HRateDestGnd[5] = (float)(hcbr.hbul[ipZ][3]/z[3][3]);
		fprintf(ioQQQ,"heliumbugdest");
		for( ipHi=0; ipHi<6 ;++ipHi )
		{
			fprintf(ioQQQ,"%11.2e",HRateDest.HRateDestGnd[ipHi] );
		}
		fprintf(ioQQQ,"\n" );
	}
#	endif

	if( (trace.lgTrace && trace.lgHBugFull) && (ipZ == trace.ipZTrace) )
	{
		if( ipZ == 0 )
		{
			fprintf(ioQQQ," HydroLevelDep grnd dest fracs:");
			for( j=IP1S; j < 8; j++ )
			{
				fprintf( ioQQQ," ");
				fprintf( ioQQQ,PrintEfmt("%8.1e", HRateDest.HRateDestGnd[j] ));
			}
			fprintf( ioQQQ, "\n" );
		}
		fprintf( ioQQQ, "  dep level     others => \n" );
		for( n=IP1S; n <= nhlevel; n++ )
		{
			fprintf( ioQQQ, "       HII%2ld", n );
			for( j=IP1S; j <= (nhlevel + 1); j++ )
			{
				fprintf( ioQQQ," ");
				fprintf( ioQQQ,PrintEfmt("%8.1e", z[j][n] ));
			}
			fprintf( ioQQQ, "\n" );
		}
	}

	/* scale to largest number of order unity */
	for( j=IP1S; j <= nhlevel; j++ )
	{
		for( n=IP1S; n <= (nhlevel + 1); n++ )
		{
			z[n][j] *= hcolrt.hlte[ipZ][j];
		}
	}

	/* save matrix */
	for( j=IP1S; j <= nhlevel; j++ )
	{
		for( n=IP1S; n <= (nhlevel + 1); n++ )
		{
			SaveZ[n][j] = z[n][j];
		}
	}

	/* which matrix solver? */
	if( strcmp(TypMatrx.chMatrix,"matin1 ") == 0 )
	{
		/* this is the usual matrix inversion method, slightly faster */
		/*matin1();*/
		nerror = 0;
		if( nerror != 0 )
		{
			fprintf( ioQQQ, " hydrogen matrix error, stop.\n" );
			puts( "[Stop in HydroLevelDep]" );
			exit(1);
		}
	}
	/* the default matrix inversion routine */
	else if( strcmp(TypMatrx.chMatrix,"linpack") == 0 )
	{
		/*double amat[LMHLVL+2][LMHLVL+2] ;*/
		double *amat;
	  /* nhlevel+1 is right dimension of matrix */
#		ifdef AMAT
#			undef AMAT
#		endif
#		define AMAT(I_,J_)	(*(amat+(I_)*(nhlevel+1)+(J_)))
		/* malloc space for the  1-d array */
		if( (amat=(double*)malloc( (sizeof(double)*(unsigned)((nhlevel+1)*(nhlevel+1)) ))) == NULL )
		{
			fprintf( ioQQQ, " HydroLevelDep malloc amat error\n" );
			puts( "[Stop in HydroLevelDep]" );
			exit(1);
		}

		/* this one may be more robust */
		for( j=IP1S; j <= nhlevel; j++ )
		{
			for( n=IP1S; n <= nhlevel; n++ )
			{
				/*amat[n][j] = z[n][j];*/
				AMAT(n,j) = z[n][j];
			}
			bvec[j] = z[nhlevel+1][j];
		}
		/*DGETRF(nhlevel+1,nhlevel+1,
		  (double*)amat,LMHLVL+2,ipiv,&nerror);
		DGETRS('N',nhlevel+1,1,(double*)amat,LMHLVL+2,ipiv,
		  bvec,LMHLVL+1,&nerror);*/
		DGETRF(nhlevel+1,nhlevel+1,
		  amat,(nhlevel+1),ipiv,&nerror);
		DGETRS('N',nhlevel+1,1,amat,(nhlevel+1),ipiv,
		  bvec,(nhlevel+1),&nerror);
		if( nerror != 0 )
		{
			fprintf( ioQQQ, " HydroLevelDep dgetrs error\n" );
			puts( "[Stop in HydroLevelDep]" );
			exit(1);
		}
		/* now put results back into z so rest of code treates only
		 * one case - as if matin1 had been used */
		for( n=IP1S; n <= nhlevel; n++ )
		{
			z[nhlevel+1][n] = bvec[n];
		}
		free( amat);
	}
	else if( strcmp(TypMatrx.chMatrix,"veclib ") == 0 )
	{
		double /*amat[LMHLVL+2][LMHLVL+2],*/ rcond;
		/* 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=IP1S; j <= nhlevel; j++ )
		{
			/*for( n=IP1S; n <= nhlevel; n++ )
			{
				amat[n][j] = z[n][j];
			}*/
			bvec[j] = z[nhlevel+1][j];
		}
		job = 0;
		rcond = 0.;
		/*dgeco((double*)amat,nhlevel+2,nhlevel+1,ipiv,rcond, work);
		dgesl((double*)amat,nhlevel+2,nhlevel+1,ipiv,bvec, job);*/
		dgeco((double*)z,nhlevel+2,nhlevel+1,ipiv,rcond, work);
		dgesl((double*)z,nhlevel+2,nhlevel+1,ipiv,bvec, job);
		/* now put results back into z so rest of code treates only
		 * one case - as if matin1 had been used */
		for( n=IP1S; n <= nhlevel; n++ )
		{
			z[nhlevel+1][n] = bvec[n];
		}
		puts( "[Stop in HydroLevelDep]" );
		exit(1);
	}
	else
	{
		fprintf( ioQQQ, " chMatrix type insane in HydroLevelDep, was%7.7s\n", 
		  TypMatrx.chMatrix );
		puts( "[Stop in HydroLevelDep]" );
		exit(1);
	}
	/* end of branches for which matrix solution, now check if valid */

	/* check whether solution is valid */
	for( level=IP1S; level <= nhlevel; level++ )
	{
		double BigSoln= 0.;
		error[level] = 0.;
		for( n=IP1S; n <= nhlevel; 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]*z[nhlevel+1][n];
		}

		if( BigSoln > 0. )
		{
			error[level] = (error[level] - SaveZ[nhlevel+1][level])/ BigSoln;
		}
		else
		{
			error[level] = 0.;
		}
	}

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

	/* put departure coefficients and level populations into master array */
	for( ipLo=IP1S; ipLo <= nhlevel; ipLo++ )
	{
		hydro.hbn[ipZ][ipLo] = z[nhlevel+1][ipLo];

		/* hn(n,ipZ) are pop per level rel to HII; n.e., N(n)/N(HII) */
		hydro.hn[ipZ][ipLo] = 
			(float)(hydro.hbn[ipZ][ipLo]*hcolrt.hlte[ipZ][ipLo]*
		  phycon.eden);
	}

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