/*he1col evaluate collisional rates for helium singlets */
#include "cddefines.h"
#include "nhe1lvl.h"
#include "tepowers.h"
#include "phe1lv.h"
#include "heclon.h"
#include "casbhs.h"
#include "trace.h"
#include "phycon.h"
#include "he1blt.h"
#include "he1crt.h"
#include "he1cbr.h"
#include "he1stat.h"
#include "he1deltat.h"
#include "cfit.h"
#include "converge.h"
#include "hdexct.h"
#include "he1.h"

void he1col(void)
{
	long int i, 
	  ihi, 
	  iup, 
	  j, 
	  lo, 
	  n;
	double cs12p, 
	  cs12s, 
	  f, 
	  hd, 
	  hsfac, 
	  sum;

	static double rcion[9]={1.04e-10,8.24e-9,1.08e-7,6.61e-7,3.45e-6,
	  8.06e-6,9.47e-5,2.38e-3,3.36e-1};
	static double r1n[9]={0.,0.,0.,1.05e-8,3.45e-9,1.54e-9,3.23e-10,
	  2.14e-11,7.09e-14};
	static double r2n[9]={0.,0.,7.88e-8,1.14e-8,3.10e-9,1.10e-9,2.0e-10,
	  1.0e-11,2.8e-14};
	static double told = 0.;

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

	/* define boltz facs between adjacent levels */
	for( iup=1; iup < NHE1LVL; iup++ )
	{
		/* He1deltaTe is array of line temps, upper to lower
		 * their boltzmann factors for this temp */
		he1blt.he1bol[iup][iup-1] = exp(-(double)(He1deltaT.He1deltaTe[iup-1][iup]/
		  phycon.te));
	}

	/* this is 1s-2s */
	he1blt.he1bol[NHE1LVL][0] = he1blt.he1bol[1][0];

	/* the highest level, up to the continuum */
	he1blt.he1cbt[NHE1LVL-1] = exp(-(double)(He1deltaT.He1CionTe[NHE1LVL-1]/
	  phycon.te));

	/* define remaining Boltzmann factors */
	for( j=0; j < (NHE1LVL - 2); j++ )
	{
		for( i=j + 2; i < NHE1LVL; i++ )
		{
			he1blt.he1bol[i][j] = he1blt.he1bol[i-1][j]*he1blt.he1bol[i][i-1];
		}
	}

	for( i=0; i < (NHE1LVL - 1); i++ )
	{
		he1blt.he1cbt[i] = he1blt.he1bol[NHE1LVL-1][i]*he1blt.he1cbt[NHE1LVL-1];
	}
	he1blt.he1cbt[NHE1LVL] = he1blt.he1cbt[1];

	/*-------------------------------------------------------------
	 *
	 * define HE1LTE(I) = lte population (cm^3) / Ne Np
	 * only used in HE1LEV if Te>1000K
	 * extra div by 2 due to He+ partition function */
	f = 2.0708e-16/TePowers.te32/2.;
	for( i=0; i < NHE1LVL; i++ )
	{
		if( (float)(he1blt.he1cbt[i]) > 1e-70 )
		{
			/* n.b.; on following two, order is VERY imp fror 32 bit machines */
			he1crt.he1lte[i] = (f/he1blt.he1cbt[i])*he1statCOM.he1stat[i];
		}
		else
		{
			he1crt.he1lte[i] = 0.;
		}
	}

	he1crt.he1lte[NHE1LVL] = he1crt.he1lte[1]*he1statCOM.he1stat[NHE1LVL]/
	  he1statCOM.he1stat[1];

	/*-----------------------------------------------------------------------
	 * evaluate all needed collisional and rad rates
	 * collisions turned off, HD=0
	 *----------------------------------------------------------------------
	 * */
	if( heclon.lgHeClON )
	{
		hd = 1.;
	}
	else
	{
		hd = 0.;
	}

	/* collisional ionization
	 * he1clc(1) = rcion(1) * he1cbt(1)*sqrte*hd
	 * call colfit(2,2,0,te,he1clc(1) )
	 * >>chng 97 mar 19, to new data set for collisional ionization */
	cfit(2,2,phycon.te,&phe1lv.he1clc[0]);

	phe1lv.he1clc[0] *= (float)(hd);

	phe1lv.he1clc[1] = (float)(rcion[1]*he1blt.he1cbt[1]*TePowers.te10*TePowers.te10*
	  TePowers.te10*TePowers.te01*TePowers.te01*hd);

	phe1lv.he1clc[2] = (float)(rcion[2]*he1blt.he1cbt[2]*TePowers.te10*TePowers.te10*
	  hd);

	phe1lv.he1clc[3] = (float)(rcion[3]*he1blt.he1cbt[3]*TePowers.te10*TePowers.te01*
	  hd);

	phe1lv.he1clc[4] = (float)(rcion[4]*he1blt.he1cbt[4]*TePowers.te01*hd);

	phe1lv.he1clc[5] = (float)(rcion[5]*he1blt.he1cbt[5]/TePowers.te01/TePowers.te01*
	  hd);

	phe1lv.he1clc[6] = (float)(rcion[6]*he1blt.he1cbt[6]/(TePowers.te10*TePowers.te03*
	  TePowers.te03)*hd);

	phe1lv.he1clc[7] = (float)(rcion[7]*he1blt.he1cbt[7]/(TePowers.te10*TePowers.te10*
	  TePowers.te10*TePowers.te03)*hd);

	phe1lv.he1clc[8] = (float)(rcion[8]*he1blt.he1cbt[8]/TePowers.sqrte*hd);

	/* don't reevaluate downward coll if dT/T not large
	 * he1cll is zero on very first call since zeroed in block */
	if( fabs(phycon.te-told)/phycon.te < 0.002 && phe1lv.he1cll[0][1] > 0. &&
		conv.nTotalIoniz )
	{
		if( trace.lgTrace && trace.lgHeBug )
		{
			fprintf( ioQQQ, "       HE1COL called-no reval Cul\n" );
		}
		/*  (but still eval col below) */
	}
	else
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, "       HE1COL called-will reeval Cul\n" );
		}
		told = phycon.te;

		/* from Osterbrock's book, fit to 10,000K - 20,000K
		 * >>chng 97 july 25, updated to 
		 *>>refer Sawey, P.M.J., and Berrington, K.A., 1993,
		 *>>refer At. Data Nucl. Data Tables 55, 81 */
		f = MIN2(26.63,5.112e-3*phycon.te/TePowers.te10/TePowers.te01);
		he1crt.he1c2p2s = (float)(f*8.629e-6/TePowers.sqrte/2.);

		/*  collision rates for 1-2S, 1-2P, and 1-3 from
		 *  >>refer Aggarwal, K.M. 1983, MNRAS, 202, 15P
		 *  >>chng 97 july 25, updated to 
		 *  >>refer Sawey, P.M.J., and Berrington, K.A., 1993,
		 *  >>refer At. Data Nucl. Data Tables 55, 81 */
		cs12s = MIN2(0.0425,9.055e-3*TePowers.te10*TePowers.te05);

		/*  f = cs12p */
		cs12p = MIN2(0.0311,4.242e-5*TePowers.sqrte*TePowers.te02*
		  TePowers.te02);

		/*  C1N is deexcitation rate
		 *  following is 2P - 1S only */
		phe1lv.he1cll[0][1] = (float)(cs12p/TePowers.sqrte*8.629e-6/6.*hd);
		/*  following is 2S - 1S only */

		he1crt.he1c2s1s = (float)(cs12s/TePowers.sqrte*8.629e-6/2.*hd);

		/*  option to turn off  collisions frrom ground and n=2 */
		if( casbhs.lgCasBhs || (hd == 0.) )
		{
			hsfac = 0.;
		}
		else
		{
			hsfac = 1.;
		}

		/*  HE1CLL(u,l) is de-excitation rate from u to l
		 *  rate coef from Callaway et al. preprint */
		phe1lv.he1cll[0][2] = (float)(1.4e-9*hsfac);

		/*  rate from i to both 2s and 2p */
		phe1lv.he1cll[1][2] = (float)((3.39*TePowers.te10*TePowers.te10*TePowers.te10/
		  TePowers.te01/TePowers.te01/TePowers.te01)/TePowers.sqrte*
		  8.629e-6/18.*hsfac);

		for( ihi=3; ihi < NHE1LVL; ihi++ )
		{
			phe1lv.he1cll[0][ihi] = (float)(r1n[ihi]/TePowers.sqrte*hsfac);
			/* different temp denpendence is not a bug */
			phe1lv.he1cll[1][ihi] = (float)(r2n[ihi]*hsfac);
			for( lo=2; lo < ihi ; lo++ )
			{
				/* hdexct returns collision strength in C version,
				 * had been deexcitation rate in C90.  Convert following
				 * to deexcitation rate */
				phe1lv.he1cll[lo][ihi] = (float)(hdexct(ihi+1,lo+1)*hd*
					8.629e-6 / TePowers.sqrte / he1statCOM.he1stat[ihi]);
			}
		}


		if( trace.lgTrace && trace.lgHe1Bug )
		{
			fprintf( ioQQQ, "       HE1COL rate n->n-1:" );
			for( i=1; i <= 8; i++ )
			{
				fprintf( ioQQQ, "%10.2e", phe1lv.he1cll[i-1][i] );
			}
			fprintf( ioQQQ, "%10.2e\n", he1crt.he1c2s1s );

			fprintf( ioQQQ, "       HE1COL col ion cof:" );

			for(i=0; i < 9; i++)
				fprintf( ioQQQ, "%10.2e", phe1lv.he1clc[i] );
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HE1COL rate down 1:" );
			for( i=0; i < 9; i++ )
			{
				fprintf( ioQQQ, "%10.2e", phe1lv.he1cll[0][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HE1COL rate down 2:" );
			for( i=0; i < 9; i++ )
			{
				fprintf( ioQQQ, "%10.2e", phe1lv.he1cll[1][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HE1COL rate down 3:" );
			for( i=0; i < 9; i++ )
			{
				fprintf( ioQQQ, "%10.2e", phe1lv.he1cll[2][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HE1COL rate down 4:" );
			for( i=0; i < 9; i++ )
			{
				fprintf( ioQQQ, "%10.2e", phe1lv.he1cll[3][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HE1COL rate down 5:" );
			for( i=0; i < 9; i++ )
			{
				fprintf( ioQQQ, "%10.2e", phe1lv.he1cll[4][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HE1COL boltz fac c:" );
			for(i=0; i < 9; i++)
				fprintf( ioQQQ, "%10.2e", he1blt.he1cbt[i] );
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HE1COL HE1LTE(i):  " );
			for(i=0; i < 9; i++)
				fprintf( ioQQQ, "%10.2e", he1crt.he1lte[i] );
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       he1col       c2p2s=%10.2e\n", 
			  he1crt.he1c2p2s );
		}
	}

	if( he1crt.he1lte[0] > 0. )
	{
		/* following are eden independent, used in hlevel */
		he1cbr.he1coldn[0] = (float)((phe1lv.he1cll[0][1]*he1crt.he1lte[1] + 
		  he1crt.he1c2s1s*he1crt.he1lte[9] + phe1lv.he1cll[0][2]*he1crt.he1lte[2] + 
		  phe1lv.he1cll[0][3]*he1crt.he1lte[3] + phe1lv.he1cll[0][4]*
		  he1crt.he1lte[4] + phe1lv.he1cll[0][5]*he1crt.he1lte[5] + 
		  phe1lv.he1cll[0][6]*he1crt.he1lte[6] + phe1lv.he1cll[0][7]*
		  he1crt.he1lte[7] + phe1lv.he1cll[0][8]*he1crt.he1lte[8])/
		  he1crt.he1lte[0]);

		/* SUM is total rate from high states to both 2s and 2p */
		sum = (phe1lv.he1cll[1][2]*he1crt.he1lte[2] + phe1lv.he1cll[1][3]*
		  he1crt.he1lte[3] + phe1lv.he1cll[1][4]*he1crt.he1lte[4] + 
		  phe1lv.he1cll[1][5]*he1crt.he1lte[5] + phe1lv.he1cll[1][6]*
		  he1crt.he1lte[6] + phe1lv.he1cll[1][7]*he1crt.he1lte[7] + 
		  phe1lv.he1cll[1][8]*he1crt.he1lte[8])/(he1crt.he1lte[1] + 
		  he1crt.he1lte[9]);

		/* HE1CLL(2,1) is 2p - 1s rate; c2s1s is 2s - 1s rate */
		he1cbr.he1coldn[1] = (float)(sum + phe1lv.he1cll[0][1]);
		he1cbr.he1coldn[9] = (float)(sum + he1crt.he1c2s1s);

		he1cbr.he1coldn[2] = (float)(phe1lv.he1cll[0][2] + phe1lv.he1cll[1][2] + 
		  (phe1lv.he1cll[2][3]*he1crt.he1lte[3] + phe1lv.he1cll[2][4]*
		  he1crt.he1lte[4] + phe1lv.he1cll[2][5]*he1crt.he1lte[5] + 
		  phe1lv.he1cll[2][6]*he1crt.he1lte[6] + phe1lv.he1cll[2][7]*
		  he1crt.he1lte[7] + phe1lv.he1cll[2][8]*he1crt.he1lte[8])/
		  he1crt.he1lte[2]);

		he1cbr.he1coldn[3] = (float)(phe1lv.he1cll[0][3] + phe1lv.he1cll[1][3] + 
		  phe1lv.he1cll[2][3] + (phe1lv.he1cll[3][4]*he1crt.he1lte[4] + 
		  phe1lv.he1cll[3][5]*he1crt.he1lte[5] + phe1lv.he1cll[3][6]*
		  he1crt.he1lte[6] + phe1lv.he1cll[3][7]*he1crt.he1lte[7] + 
		  phe1lv.he1cll[3][8]*he1crt.he1lte[8])/he1crt.he1lte[3]);

		he1cbr.he1coldn[4] = (float)((phe1lv.he1cll[4][5]*he1crt.he1lte[5] + 
		  phe1lv.he1cll[4][6]*he1crt.he1lte[6] + phe1lv.he1cll[4][7]*
		  he1crt.he1lte[7] + phe1lv.he1cll[4][8]*he1crt.he1lte[8])/
		  he1crt.he1lte[4] + phe1lv.he1cll[0][4] + phe1lv.he1cll[1][4] + 
		  phe1lv.he1cll[2][4] + phe1lv.he1cll[3][4]);

		he1cbr.he1coldn[5] = (float)(phe1lv.he1cll[0][5] + phe1lv.he1cll[1][5] + 
		  phe1lv.he1cll[2][5] + phe1lv.he1cll[3][5] + phe1lv.he1cll[4][5] + 
		  (phe1lv.he1cll[5][6]*he1crt.he1lte[6] + phe1lv.he1cll[5][7]*
		  he1crt.he1lte[7] + phe1lv.he1cll[5][8]*he1crt.he1lte[8])/
		  he1crt.he1lte[5]);

		he1cbr.he1coldn[6] = (float)(phe1lv.he1cll[0][6] + phe1lv.he1cll[1][6] + 
		  phe1lv.he1cll[2][6] + phe1lv.he1cll[3][6] + phe1lv.he1cll[4][6] + 
		  phe1lv.he1cll[5][6] + (phe1lv.he1cll[6][7]*he1crt.he1lte[7] + 
		  phe1lv.he1cll[6][8]*he1crt.he1lte[8])/he1crt.he1lte[6]);

		he1cbr.he1coldn[7] = (float)(phe1lv.he1cll[0][7] + phe1lv.he1cll[1][7] + 
		  phe1lv.he1cll[2][7] + phe1lv.he1cll[3][7] + phe1lv.he1cll[4][7] + 
		  phe1lv.he1cll[5][7] + phe1lv.he1cll[6][7] + phe1lv.he1cll[7][8]*
		  he1crt.he1lte[8]/he1crt.he1lte[7]);

		he1cbr.he1coldn[8] = phe1lv.he1cll[0][8] + phe1lv.he1cll[1][8] + 
		  phe1lv.he1cll[2][8] + phe1lv.he1cll[3][8] + phe1lv.he1cll[4][8] + 
		  phe1lv.he1cll[5][8] + phe1lv.he1cll[6][8] + phe1lv.he1cll[7][8];

		if( trace.lgTrace && trace.lgHe1Bug )
		{
			fprintf( ioQQQ, "       HE1CLL he1coldn(i)      " );
			for( i=0; i < 10; i++ )
			{
				fprintf( ioQQQ, "%10.2e", he1cbr.he1coldn[i] );
			}
			fprintf( ioQQQ, "\n" );
		}
	}
	else
	{
		for( n=0; n < NHE1LVL; n++ )
		{
			he1cbr.he1coldn[n] = 0.;
		}
	}


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