/*he1rad evaluate radiative rates for helium singlets */
#include "cddefines.h"
#include "taulines.h"
#include "nhe1lvl.h"
#include "he1tau.h"
#include "he3tau.h"
#include "trace.h"
#include "otsmin.h"
#include "touton.h"
#include "tepowers.h"
#include "phycon.h"
#include "phe1lv.h"
#include "he1rsp.h"
#include "logte.h"
#include "he1str.h"
#include "he1esc.h"
#include "he1nxt.h"
#include "he1as.h"
#include "he1cbr.h"
#include "heopfr.h"
#include "he1rb.h"
#include "nhe1.h"
#include "he1net.h"
#include "he1stat.h"
#include "converge.h"
#include "caseb.h"
#include "sexp.h"
#include "receff.h"
#include "he1.h"

void he1rad(void)
{
	long int i, 
	  j, 
	  n;

	double a, 
	  b, 
	  fac;
	/* used to save temperature, to reevaluate rec coef */
	static float tused=-1.f;

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

#	define RECF(a,b,c)	((float)((a) + (b)*POW2(logte.alogte - (c))))

	/* evaluate rec coef for H, HE1REC(I,2) is recombination efficiency
	 * OTSMIN normally equal to zero */
	for( n=0; n < NHE1LVL; n++ )
	{
		phe1lv.he1rec[2][n] = (float)receff(nhe1Com.nhe1[n]);
	}

	if( caseb.lgCaseb )
	{
		phe1lv.he1rec[2][0] = 1e-10f;
	}

	for( n=0; n < NHE1LVL; n++ )
	{
		phe1lv.he1rec[2][n] = (float)MAX2(phe1lv.he1rec[2][n],otsminCom.otsmin);
		/* last factor accounts for metal absorption of HeII diffuse
		 * fields, computed in AddOpac, this is counted as recombination */
		phe1lv.he1rec[1][n] = (float)MIN2(1.,phe1lv.he1rec[2][n]+
			(1.-phe1lv.he1rec[2][n])*  heopfr.ophe1f[n]);
	}

	phe1lv.he1rec[1][NHE1LVL] = phe1lv.he1rec[1][1];
	phe1lv.he1rec[2][NHE1LVL] = phe1lv.he1rec[2][1];

	if( trace.lgTrace && trace.lgHe1Bug )
	{
		fprintf( ioQQQ, "       HE1RAD geo recomb effic" );
		for( i=1; i <= 9; i++ )
		{
			fprintf( ioQQQ, "%10.2e", phe1lv.he1rec[2][i-1] );
		}
		fprintf( ioQQQ, "\n" );
		fprintf( ioQQQ, "       HE1RAD net recomb effic" );
		for( i=1; i <= 9; i++ )
		{
			fprintf( ioQQQ, "%10.2e", phe1lv.he1rec[1][i-1] );
		}
		fprintf( ioQQQ, "\n" );
	}

	/* don't reevaluate coefficients if temp unchanged, but always do so
	 * on very first call in multi model grid  -  NOT test, tests whether
	 * this is very first call since cdInit was called */
	if( fabs(tused-phycon.te)/tused > 0.002 || !conv.nTotalIoniz )
	{
		tused = phycon.te;
		phe1lv.he1rec[0][0] = (float)(1.6e-11/TePowers.sqrte);

		he1rsp.he1r2s = (float)(3.84e-13/TePowers.sqrte*TePowers.te03*TePowers.te01);
		he1rsp.he1r2p = (float)(7.768e-12/TePowers.te70*TePowers.te01);
		phe1lv.he1rec[0][1] = he1rsp.he1r2p;
		phe1lv.he1rec[0][NHE1LVL] = he1rsp.he1r2s;
		fac = 4.0*phycon.te;
		if( phycon.te <= 1e5 )
		{
			phe1lv.he1rec[0][2] = (float)(1.256e-11/(TePowers.te70*TePowers.te03*
			  TePowers.te03/TePowers.te01));
		}
		else
		{
			phe1lv.he1rec[0][2] = (float)(1.1*pow(10.,-RECF(9.2150,.12898,
			  4.9933))/fac);
		}
		fac *= 1.1;
		phe1lv.he1rec[0][3] = (float)(pow(10.,-RECF(9.4583,.12236,4.7567))/
		  fac);
		phe1lv.he1rec[0][4] = (float)(pow(10.,-RECF(9.6516,.11481,4.5614))/
		  fac);
		phe1lv.he1rec[0][5] = (float)(pow(10.,-RECF(9.8112,.10747,4.3913))/
		  fac);
		phe1lv.he1rec[0][6] = (float)(pow(10.,-RECF(9.4942,.093811,4.0473))/
		  fac);
		phe1lv.he1rec[0][7] = (float)(pow(10.,-RECF(9.5858,.070572,3.2727))/
		  fac);
		phe1lv.he1rec[0][8] = (float)(pow(10.,-RECF(9.5028,.045593,1.4592))/
		  fac);

		/* >>chng 96 june 5, big jump in rec coef at level7 (old cota level)
		 * causes maser in 7 -6 - don't let this happen, add extro to top level
		 * 666 error! change this to true levels, like Jason's hydrogen
		 * factor is to account for lowering  of rec coef with inc n */
		b = 0.75;
		for( i=6; i < (NHE1LVL - 1); i++ )
		{
			/* amount we will subtract from here, then add to highest level */
			a = phe1lv.he1rec[0][i] - phe1lv.he1rec[0][i-1]*b;
			phe1lv.he1rec[0][i] -= (float)a;
			phe1lv.he1rec[0][NHE1LVL-1] += (float)a;
		}

		/* dielectronic recombination */
		if( phycon.te > 5e4 )
		{
			phe1lv.he1rec[0][NHE1LVL-1] += (float)(1.9e-3/TePowers.te32*sexp(4.7e5/
			  phycon.te)*(1. + 0.3*sexp(9.4e4/phycon.te)));
			/* >>chng 96 oct 22, first coef had been 4.5 below, change to 4.7 as
			 *in aldrovandi and pequignot
			 *  1    1.9e-3/te32*sexp(4.5e5/te)*(1.+0.3*sexp(9.4e4/te))*0. */
		}

		he1rb.he1rtt = phe1lv.he1rec[0][0]*phe1lv.he1rec[1][0];
		/* He1RecB is the case B singlet recombination rate coef */
		he1rb.He1RecB = 0.;
		for( i=1; i < NHE1LVL; i++ )
		{
			he1rb.he1rtt += phe1lv.he1rec[0][i]*phe1lv.he1rec[1][i];
			he1rb.He1RecB += phe1lv.he1rec[0][i];
		}

		if( trace.lgTrace && trace.lgHe1Bug )
		{
			fprintf( ioQQQ, "       HE1RAD eval rec cof" );
			for( i=0; i < 9; i++ )
			{
				fprintf( ioQQQ, "%10.2e", phe1lv.he1rec[0][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HE1RAD total effective recombination coef%10.2e Case B, A=%10.2e%10.2e\n", 
			  he1rb.he1rtt, he1rb.He1RecB, he1rb.He1RecB + phe1lv.he1rec[0][0] );
		}
	}
	else
	{
		he1rb.he1rtt = phe1lv.he1rec[0][0]*phe1lv.he1rec[1][0];
		he1rb.He1RecB = 0.;

		for( i=1; i < NHE1LVL; i++ )
		{
			he1rb.he1rtt += phe1lv.he1rec[0][i]*phe1lv.he1rec[1][i];
			he1rb.He1RecB += phe1lv.he1rec[0][i];
		}

		if( trace.lgTrace && trace.lgHe1Bug )
		{
			fprintf( ioQQQ, "       HE1RAD not reevel rec cof total effective recombination coef is %10.2e\n", 
			  he1rb.he1rtt );
		}
	}

	/* following has MIN(1, removed may 12 96, is it in hydro and he2? */
	for( j=0; j < (NHE1LVL - 1); j++ )
	{
		for( i=j + 1; i < NHE1LVL; i++ )
		{
			/* >>chng 96 may 13 cap of 1 prevented maser emission
			 * he1mis(i,j) = He1EinA(i,j)*MIN(1.,(he1esc(i,j)+peshe1(i,j))) */
			phe1lv.he1mis[j][i] = he1as.He1EinA[j][i]*(he1escCOM.he1esc[j][i] + 
			  he1str.peshe1[j][i]);
			he1netCOM.he1net[j][i] = phe1lv.he1mis[j][i] + phe1lv.he1dst[j][i];
		}
	}

	/* these are used in HLEVEL - total decay rates */
	he1cbr.he1radn[1] = he1netCOM.he1net[0][1] + t206.Pesc*
	  t206.Aul;

	/* this is what the two photon integrates up to */
	he1cbr.he1radn[9] = 8.23f;

	he1cbr.he1radn[2] = he1netCOM.he1net[0][2] + he1netCOM.he1net[1][2];

	he1cbr.he1radn[3] = he1netCOM.he1net[0][3] + he1netCOM.he1net[1][3] + 
	  he1netCOM.he1net[2][3];

	he1cbr.he1radn[4] = he1netCOM.he1net[0][4] + he1netCOM.he1net[1][4] + 
	  he1netCOM.he1net[2][4] + he1netCOM.he1net[3][4];

	he1cbr.he1radn[5] = he1netCOM.he1net[0][5] + he1netCOM.he1net[1][5] + 
	  he1netCOM.he1net[2][5] + he1netCOM.he1net[3][5] + he1netCOM.he1net[4][5];

	he1cbr.he1radn[6] = he1netCOM.he1net[0][6] + he1netCOM.he1net[1][6] + 
	  he1netCOM.he1net[2][6] + he1netCOM.he1net[3][6] + he1netCOM.he1net[4][6] + 
	  he1netCOM.he1net[5][6];

	he1cbr.he1radn[7] = he1netCOM.he1net[0][7] + he1netCOM.he1net[1][7] + 
	  he1netCOM.he1net[2][7] + he1netCOM.he1net[3][7] + he1netCOM.he1net[4][7] + 
	  he1netCOM.he1net[5][7] + he1netCOM.he1net[6][7];

	he1cbr.he1radn[8] = he1netCOM.he1net[0][8] + he1netCOM.he1net[1][8] + 
	  he1netCOM.he1net[2][8] + he1netCOM.he1net[3][8] + he1netCOM.he1net[4][8] + 
	  he1netCOM.he1net[5][8] + he1netCOM.he1net[6][8] + he1netCOM.he1net[7][8];

	/* total rates out of level due to induced processes */
	he1cbr.he1bul[0] = (he1tauCOM.he1con[0][1]*he1statCOM.he1stat[1] + 
	  he1tauCOM.he1con[0][2]*he1statCOM.he1stat[2] + he1tauCOM.he1con[0][3]*
	  he1statCOM.he1stat[3] + he1tauCOM.he1con[0][4]*he1statCOM.he1stat[4] + 
	  he1tauCOM.he1con[0][5]*he1statCOM.he1stat[5] + he1tauCOM.he1con[0][6]*
	  he1statCOM.he1stat[6] + he1tauCOM.he1con[0][7]*he1statCOM.he1stat[7] + 
	  he1tauCOM.he1con[0][8]*he1statCOM.he1stat[8])/he1statCOM.he1stat[0];

	he1cbr.he1bul[1] = (he1tauCOM.he1con[1][2]*he1statCOM.he1stat[2] + 
	  he1tauCOM.he1con[1][3]*he1statCOM.he1stat[3] + he1tauCOM.he1con[1][4]*
	  he1statCOM.he1stat[4] + he1tauCOM.he1con[1][5]*he1statCOM.he1stat[5] + 
	  he1tauCOM.he1con[1][6]*he1statCOM.he1stat[6] + he1tauCOM.he1con[1][7]*
	  he1statCOM.he1stat[7] + he1tauCOM.he1con[1][8]*he1statCOM.he1stat[8])/
	  (he1statCOM.he1stat[1] + he1statCOM.he1stat[9]);

	he1cbr.he1bul[2] = he1tauCOM.he1con[0][2] + he1tauCOM.he1con[1][2] + 
	  (he1tauCOM.he1con[2][3]*he1statCOM.he1stat[3] + he1tauCOM.he1con[2][4]*
	  he1statCOM.he1stat[4] + he1tauCOM.he1con[2][5]*he1statCOM.he1stat[5] + 
	  he1tauCOM.he1con[2][6]*he1statCOM.he1stat[6] + he1tauCOM.he1con[2][7]*
	  he1statCOM.he1stat[7] + he1tauCOM.he1con[2][8]*he1statCOM.he1stat[8])/
	  he1statCOM.he1stat[2];

	he1cbr.he1bul[3] = he1tauCOM.he1con[0][3] + he1tauCOM.he1con[1][3] + 
	  he1tauCOM.he1con[2][3] + (he1tauCOM.he1con[3][4]*he1statCOM.he1stat[4] + 
	  he1tauCOM.he1con[3][5]*he1statCOM.he1stat[5] + he1tauCOM.he1con[3][6]*
	  he1statCOM.he1stat[6] + he1tauCOM.he1con[3][7]*he1statCOM.he1stat[7] + 
	  he1tauCOM.he1con[3][8]*he1statCOM.he1stat[8])/he1statCOM.he1stat[3];

	he1cbr.he1bul[4] = he1tauCOM.he1con[0][4] + he1tauCOM.he1con[1][4] + 
	  he1tauCOM.he1con[2][4] + he1tauCOM.he1con[3][4] + (he1tauCOM.he1con[4][5]*
	  he1statCOM.he1stat[5] + he1tauCOM.he1con[4][6]*he1statCOM.he1stat[6] + 
	  he1tauCOM.he1con[4][7]*he1statCOM.he1stat[7] + he1tauCOM.he1con[4][8]*
	  he1statCOM.he1stat[8])/he1statCOM.he1stat[4];

	he1cbr.he1bul[5] = he1tauCOM.he1con[0][5] + he1tauCOM.he1con[1][5] + 
	  he1tauCOM.he1con[2][5] + he1tauCOM.he1con[3][5] + he1tauCOM.he1con[4][5] + 
	  (he1tauCOM.he1con[5][6]*he1statCOM.he1stat[6] + he1tauCOM.he1con[5][7]*
	  he1statCOM.he1stat[7] + he1tauCOM.he1con[5][8]*he1statCOM.he1stat[8])/
	  he1statCOM.he1stat[5];

	he1cbr.he1bul[6] = he1tauCOM.he1con[0][6] + he1tauCOM.he1con[1][6] + 
	  he1tauCOM.he1con[2][6] + he1tauCOM.he1con[3][6] + he1tauCOM.he1con[4][6] + 
	  he1tauCOM.he1con[5][6] + (he1tauCOM.he1con[6][7]*he1statCOM.he1stat[7] + 
	  he1tauCOM.he1con[6][8]*he1statCOM.he1stat[8])/he1statCOM.he1stat[6];

	he1cbr.he1bul[7] = he1tauCOM.he1con[0][7] + he1tauCOM.he1con[1][7] + 
	  he1tauCOM.he1con[2][7] + he1tauCOM.he1con[3][7] + he1tauCOM.he1con[4][7] + 
	  he1tauCOM.he1con[5][7] + he1tauCOM.he1con[6][7] + he1tauCOM.he1con[7][8]*
	  he1statCOM.he1stat[8]/he1statCOM.he1stat[7];

	he1cbr.he1bul[8] = he1tauCOM.he1con[0][8] + he1tauCOM.he1con[1][8] + 
	  he1tauCOM.he1con[2][8] + he1tauCOM.he1con[3][8] + he1tauCOM.he1con[4][8] + 
	  he1tauCOM.he1con[5][8] + he1tauCOM.he1con[6][8] + he1tauCOM.he1con[7][8];

	if( trace.lgTrace && trace.lgHe1Bug )
	{
		fprintf( ioQQQ, "       HE1RAD finds arrays, with optical depths defined? %1c\n", 
		  TorF(touton.lgTauOutOn) );

		for( i=1; i < NHE1LVL; i++ )
		{
			fprintf( ioQQQ, "      up:%2ld", i );

			for( j=0; j < i; j++ )
			{
				fprintf( ioQQQ, "%10ld", j );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "         A*esc" );
			for( j=0; j < i; j++ )
			{
				fprintf( ioQQQ, "%10.2e", phe1lv.he1mis[j][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "         A-net" );
			for( j=0; j < i; j++ )
			{
				fprintf( ioQQQ, "%10.2e", he1netCOM.he1net[j][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "         tauin" );
			for( j=0; j < i; j++ )
			{
				fprintf( ioQQQ, "%10.2e", he1nxtCOM.he1nxt[j][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "         t tot" );
			for( j=0; j < i; j++ )
			{
				fprintf( ioQQQ, "%10.2e", he1tauCOM.he1lim[j][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "         Esc  " );
			for( j=0; j < i; j++ )
			{
				fprintf( ioQQQ, "%10.2e", he1escCOM.he1esc[j][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "         A*dst" );
			for( j=0; j < i; j++ )
			{
				fprintf( ioQQQ, "%10.2e", phe1lv.he1dst[j][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "         StrkE" );
			for( j=0; j < i; j++ )
			{
				fprintf( ioQQQ, "%10.2e", he1str.peshe1[j][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "         B(ul)" );
			for( j=0; j < i; j++ )
			{
				fprintf( ioQQQ, "%10.2e", he1tauCOM.he1con[j][i] );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "         C(ul)" );
			for( j=0; j < i; j++ )
			{
				fprintf( ioQQQ, "%10.2e", phe1lv.he1cll[j][i]*
				  phycon.eden );
			}
			fprintf( ioQQQ, "\n" );
		}

		fprintf( ioQQQ, "              " );
		for( i=0; i < 10; i++ )
		{
			fprintf( ioQQQ, "%9ld", i );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "       HE1RAD  rad" );
		for( i=0; i < 10; i++ )
		{
			fprintf( ioQQQ, "%10.2e", he1cbr.he1radn[i] );
		}
		fprintf( ioQQQ, "\n" );
	}


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