/*he3lev compute ionization balance level populations for helium triplets */
#include <string.h>
#include "cddefines.h"
#include "taulines.h"
#include "nhe1lvl.h"
#include "nhe3lvl.h"
#include "showme.h"
#include "opacpoint.h"
#include "herec.h"
#include "he3lte.h"
#include "he3c.h"
#include "phe1lv.h"
#include "he3bn.h"
#include "che2sn.h"
#include "he3lines.h"
#include "he3as.h"
#include "ionfracs.h"
#include "d10830.h"
#include "he3rr.h"
#include "he3tau.h"
#include "phycon.h"
#include "rfield.h"
#include "he.h"
#include "heinwd.h"
#include "he3n.h"
#include "he3gams.h"
#include "r2s3bd.h"
#include "trace.h"
#include "he3pht.h"
#include "he3hsv.h"
#include "cionhe.h"
#include "he3rate.h"
#include "zonecnt.h"
#include "typmatrx.h"
#include "linpack.h"
#include "printe82.h"
#include "matin1.h"
#include "veclib.h"
#include "csphot.h"
#include "he3.h"

void he3lev(double *he2ov3)
{
	int lgNegPop;
	long int i, 
	  ipiv[5], 
	  j, 
	  job, 
	  nerror;
	double amat[5][5], 
	  bvec[5], 
	  dest, 
	  out2, 
	  ratio, 
	  rcond, 
	  sum ,
	  work[5], 
	  zz[6][6];

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

	/* following used for linpac matrix inversion */
	/* free statement labels >= 9
	 *
	 * NOTATION;
	 * 1-2S  2-2P  3-3S  4-3P  5-3D
	 * */
	if( trace.lgTrace && trace.lgHe3Bug )
	{
		fprintf( ioQQQ, " HE3LEV called.  HeI, HeII, HeIII=%10.2e%10.2e%10.2e\n", 
		  he.hei, he.heii, he.heiii );
	}

	/* last number was 1e-15 before 95 nov 14, changed to get ldl.in to work
	 * several bugs present in following code, he3n was pop cm^-3, not rel
	 * to ion */

	/* estimate simple he triplet population */
	if( phycon.eden > 0. )
	{
		/* estimate of HeII/He3I
		 * >>chng 96 oct 27, to exact expression from hazy
		 *  1  ( 0.365/(te*te10*te10) ) */
		he3rate.she2ov3 = (float)((1. + 3110.*pow(phycon.te/1e4,-0.51)/phycon.eden)/
		  (5.79e-6*pow(phycon.te/1e4,-1.18)));
	}
	else
	{
		he3rate.she2ov3 = 0.;
	}

	/* quit if little abundance or turned off */
	if( xIonFracs[2][1]/xIonFracs[0][1] < 1e-20 )
	{
		for( i=0; i < NHE3LVL; i++ )
		{
			he3nCom.he3n[i] = 0.;
			/* >>chng 96 july 16, was set to 0, change to 1 */
			he3bnCom.he3bn[i] = 1.;
		}

		*he2ov3 = he3rate.she2ov3;
		he.hei3 = (float)(he.heii/ *he2ov3);
		he3pht.he3local = 0.;
		he3pht.he3relax = 1.;

		/* >>chng 96 july 12, following ratio was inverted, fixed 90.01
		 * he3n(1) = he2ov3 */
		if( *he2ov3 > 1e-25 )
		{
			/* this is ratio of 23S to ion */
			he3nCom.he3n[0] = (float)(1./ *he2ov3);
		}
		else
		{
			he3nCom.he3n[0] = 0.;
		}

		he3hsvCom.he3hsv[0] = 0.;
		he3hsvCom.he3hsv[1] = 0.;
		cionhe.He3IonCool = 0.;
		he3lines.p2s = he3nCom.he3n[0]*he.heii;
		he3lines.p2p = 0.;
		he3lines.p3s = 0.;
		he3lines.p3p = 0.;
		he3lines.p3d = 0.;
		he3lines.p10830 = 0.;
		heinwd.he10in = 0.;
		he3lines.p5876 = 0.;
		he3lines.p7065 = 0.;
		he3lines.p3889 = 0.;
		he3lines.he3clx = 0.;

		if( trace.lgTrace && trace.lgHe3Bug )
		{
			fprintf( ioQQQ, " HE3LEV return, all neutral, hei3, he3n(1)=%10.2e%10.2e\n", 
			  he.hei3, he3nCom.he3n[0] );
		}
		
#	ifdef DEBUG_FUN
		fputs( " <->he3lev()\n", debug_fp );
#	endif
		return;
	}

	/* following are recombination to various levels */
	zz[NHE3LVL][0] = (he3rr.r2s + r2s3bdCom.r2s3bd)*phycon.eden + 
	  phe1lv.he1n[0]*che2sn.csg2s + phe1lv.he12s*che2sn.cs2s2s + phe1lv.he12p*
	  che2sn.cs2p2s;

	zz[NHE3LVL][1] = (he3rr.r2p + r2s3bdCom.r2p3bd)*phycon.eden + 
	  phe1lv.he1n[0]*che2sn.csg2p + phe1lv.he12s*che2sn.cs2s2p + phe1lv.he12p*
	  che2sn.cs2p2p;

	zz[NHE3LVL][2] = he3rr.r3s*phycon.eden;
	zz[NHE3LVL][3] = he3rr.r3p*phycon.eden;
	zz[NHE3LVL][4] = he3rr.r3d*phycon.eden;

	/* level 1= 2s balance equation */
	zz[0][0] = he3c.c2s2p + he3c.c2s3s + he3c.c2s3p + he3c.c2s3d + 
	  he3gams.he3gam[0] + che2sn.c2sg + he3as.a2sg + he3c.c2sion + 
	  che2sn.c2ss2s + che2sn.c2ss2p;

	out2 = he3gams.he3gam[0] + che2sn.c2sg + he3as.a2sg + he3c.c2sion + 
	  che2sn.c2ss2s + che2sn.c2ss2p;

	/* fix simple estimate by this factor
	 * >>chng 96 oct 21, c2sion had been in he3relax and not in he3local */
	he3pht.he3local = (float)((he3gams.he3gam[0] + he3c.c2sion)/out2);
	he3pht.he3relax = (float)((che2sn.c2sg + he3as.a2sg + che2sn.c2ss2s + 
	  che2sn.c2ss2p)/out2);

	he3rate.she2ov3 /= he3pht.he3relax;
	he3rate.collhe3 = (float)(out2 - he3gams.he3gam[0]);

	/* if local photoionization rate is greater than stored max value, 
	 * then update both stored max value as well as fraction due to Lya */
	if( he3pht.he3local > he3pht.he3photo )
	{
		/* remember max ratio of photionzation and collisional ionization relative
		 * to all rates out */
		he3pht.he3photo = he3pht.he3local;
		/* remember the zone where it happened */
		he3pht.nzone = ZoneCnt.nzone ;
		/* remember fraction of exits due to Lya photoionization at this peak*/
		he3pht.he3lya = 
			(float)(
			csphot(
			HydroLines[0][IP2P][IP1S].ipCont , he.nhei3 , OpacPoint.ioptri)*
		  rfield.otslin[HydroLines[0][IP2P][IP1S].ipCont-1]/
		  (he3gams.he3gam[0] + he3c.c2sion));
	}

	zz[1][0] = -(d10830.r10830 + he3c.c2p2s);
	zz[2][0] = -he3c.c3s2s;
	zz[3][0] = -(he3as.a3889 + he3c.c3p2s);
	zz[4][0] = -he3c.c3d2s;

	/* level 2=2p balance equation */
	zz[1][1] = he3c.c2p2s + d10830.r10830 + he3as.a591 + he3gams.he3gam[1] + 
	  he3c.c2p3s + he3c.c2p3d + he3c.c2p3p + che2sn.c2pg + he3c.c2pion + 
	  che2sn.c2ps2s + che2sn.c2ps2p;

	zz[0][1] = -he3c.c2s2p;
	zz[2][1] = -(he3as.a7065 + he3c.c3s2p);
	zz[3][1] = -he3c.c3p2p;
	zz[4][1] = -(he3as.a5876 + he3c.c3d2p);

	/* level 3=3s balance equation */
	zz[2][2] = he3as.a7065 + he3c.c3s2p + he3c.c3s2s;
	zz[0][2] = -he3c.c2s3s;
	zz[1][2] = -he3c.c2p3s;
	zz[3][2] = -he3as.a4;
	zz[4][2] = 0.;

	/* LEVEL 4=3P BALANCE EQUATION */
	zz[3][3] = he3as.a4 + he3as.a3889 + he3c.c3p2s + he3c.c3p2p;
	zz[0][3] = -he3c.c2s3p;
	zz[1][3] = -he3c.c2p3p;
	zz[2][3] = 0.;
	zz[4][3] = -he3as.a19;

	/* LEVEL 5=3D BALANCE EQUATION */
	zz[4][4] = he3c.c3d2s + he3c.c3d2p + he3as.a5876 + he3as.a19;
	zz[0][4] = -he3c.c2s3d;
	zz[1][4] = -he3c.c2p3d;
	zz[2][4] = 0.;
	zz[3][4] = 0.;

	if( trace.lgTrace && trace.lgHe3Bug)
	{
		fprintf( ioQQQ, "     He triplet matrix\n" );
		for( j=0; j<NHE3LVL; ++j )
		{
			for( i=0; i<=NHE3LVL; ++i )
			{
				fprintf(ioQQQ,PrintEfmt("%9.2e", zz[i][j] ) );
			}
			fprintf( ioQQQ,"\n");
		}
	}


	for( j=0; j < NHE3LVL; j++ )
	{
		ratio = 1./zz[j][j];
		for( i=0; i < (NHE3LVL + 1); i++ )
		{
			zz[i][j] *= ratio;
		}
	}


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

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

		DGETRF(5,5,(double*)amat,5,ipiv,&nerror);
		DGETRS('N',5,1,(double*)amat,5,ipiv,bvec,5,&nerror);
		if( nerror != 0 )
		{
			fprintf( ioQQQ, " fe2 ir dgetrs error\n" );
			puts( "[Stop in he3lev]" );
			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 < 5; i++ )
		{
			zz[5][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 < 5; j++ )
		{
			for( i=0; i < 5; i++ )
			{
				amat[i][j] = zz[i][j];
			}
			bvec[j] = zz[5][j];
		}

		job = 0;
		rcond = 0.;
		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 < 5; i++ )
		{
			zz[5][i] = bvec[i];
		}

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

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

	/* HE3N(i) must always be population/HeII */
	he3nCom.he3n[0] = (float)zz[5][0];
	he3lines.p2s = he3nCom.he3n[0]*he.heii;

	he3nCom.he3n[1] = (float)zz[5][1];
	he3lines.p2p = he3nCom.he3n[1]*he.heii;

	he3nCom.he3n[2] = (float)zz[5][2];
	he3lines.p3s = he3nCom.he3n[2]*he.heii;

	he3nCom.he3n[3] = (float)zz[5][3];
	he3lines.p3p = he3nCom.he3n[3]*he.heii;

	he3nCom.he3n[4] = (float)zz[5][4];
	he3lines.p3d = he3nCom.he3n[4]*he.heii;

	sum = 0.;
	lgNegPop = FALSE;
	for( i=0; i < NHE3LVL; i++ )
	{
		if( he3nCom.he3n[i] <= 0. )
			lgNegPop = TRUE;
		sum += he3nCom.he3n[i];
		if( he3lteCom.he3lte[i] > 0. )
		{
			he3bnCom.he3bn[i] = he3nCom.he3n[i]/he3lteCom.he3lte[i]/phycon.eden;
		}
		else
		{
			he3bnCom.he3bn[i] = 1.;
		}
	}

	/* make sure all pops are positive */
	if( lgNegPop )
	{
		fprintf( ioQQQ, " negative helium triplet populations, zone%4ld\n", 
		  ZoneCnt.nzone );
		fprintf( ioQQQ, " pops are; " );
		for( i=0; i<NHE3LVL ; i++ )
		{
			fprintf( ioQQQ, " %f", he3nCom.he3n[i] );
		}
		fprintf( ioQQQ, "\n" );
		ShowMe();
	}

	if( sum > 0. )
	{
		*he2ov3 = 1./sum;
	}
	else
	{
		*he2ov3 = 0.;
	}

	he3lines.p10830 = (float)(he3lines.p2p*he3as.a10830*1.84e-12);
	heinwd.he10in = he3lines.p10830*he3tau[IPT10830-1].FracInwd;
	he3lines.p5876 = (float)(he3lines.p3d*he3as.a5876*3.39e-12);
	he3lines.p7065 = (float)(he3lines.p3s*he3as.a7065*2.82e-12);
	he3lines.p3889 = (float)(he3lines.p3p*he3as.a3889*5.12e-12);

	/* total collisional excitation cooling of triplet lines */
	he3lines.he3clx = (float)((he3lines.p2s*he3c.c2s2p - he3lines.p2p*he3c.c2p2s)*
	  1.84e-12 + (he3lines.p2s*(he3c.c2s3s + he3c.c2s3p + he3c.c2s3d) - 
	  he3lines.p3s*he3c.c3s2s - he3lines.p3p*he3c.c3p2s - he3lines.p3d*
	  he3c.c3d2s)*5e-12 + (he3lines.p2p*(he3c.c2p3s + he3c.c2p3d) - 
	  he3lines.p3s*he3c.c3s2p - he3lines.p3d*he3c.c3d2p)*3e-12);

	/* collisional ionization cooling
	 * c2sion, c2pion, both already have eden included
	 * >>chng 96 july 15, added cor for 3 body heating */
	cionhe.He3IonCool = (float)(he3lines.p2s*he3c.c2sion*(1. - 1./he3bnCom.he3bn[0])*
	  7.735e-12 + he3lines.p2p*he3c.c2pion*(1. - 1./he3bnCom.he3bn[1])*
	  5.8018e-12);

	if( trace.lgTrace && trace.lgHeBug )
	{
		fprintf( ioQQQ, 
			"     HeI tri gam%10.2e rec%10.2e 3bd%10.2e cion2S%10.2e col 23S-sin%10.2e\n", 
		  he3gams.he3gam[0], herecCom.rectri, r2s3bdCom.r2s3bd + r2s3bdCom.r2p3bd, 
		  he3c.c2sion, (che2sn.c2ss2s + che2sn.c2ss2p)/phycon.eden );

		if( trace.lgHe3Bug )
		{
			fprintf( ioQQQ, "     HeI tri lev" );
			for(i=0; i < NHE3LVL; i++)
			{
				fprintf( ioQQQ, "%9.1e", he3nCom.he3n[i] );
			}
			fprintf( ioQQQ, "\n" );

			for( i=0; i < 5; i++ )
			{
				fprintf( ioQQQ, "%9.1e", he3nCom.he3n[i]/he3nCom.he3n[0] );
			}
			fprintf( ioQQQ, "\n" );
		}

		dest = he3c.c2sion + he3gams.he3gam[0] + che2sn.c2ss2s + che2sn.c2ss2p + 
		  he3as.a2sg;

		if( dest > 0. )
		{
			fprintf( ioQQQ, "     HeI3 dest frc; col-ion%10.2e H Ly-a%10.2e C+A to sing%10.2e\n", 
			  he3c.c2sion/dest, 
			  rfield.otslin[HydroLines[0][IP2P][IP1S].ipCont-1]*1.97e-18/dest, 
			  (che2sn.c2ss2s + che2sn.c2ss2p + he3as.a2sg)/dest );
		}
		else
		{
			fprintf( ioQQQ, "     HeI3 no dest rates\n" );
		}
	}

	if( *he2ov3 <= 0. )
	{
		fprintf( ioQQQ, " HE3LEV finds relative population <=0;%10.3e\n", 
		  *he2ov3 );
		puts( "[Stop in he3lev]" );
		exit(1);
	}

	if( trace.lgTrace && trace.lgHeBug )
	{
		fprintf( ioQQQ, 
			"     HE3LEV returns; HE2OV3=%10.2e simple=%10.2e relx=%10.2e", 
		  *he2ov3, he3rate.she2ov3,he3pht.he3relax );
		fprintf( ioQQQ, " pops" );
		for(i=0; i < NHE3LVL; i++)
		{
			fprintf( ioQQQ, "%9.1e", he3nCom.he3n[i] );
		}
		fprintf( ioQQQ, "\n" );
	}

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