/*HydroCollid evaluate collision rate for model hydrogen atom */
#include <float.h>
#include "cddefines.h"
#include "physconst.h"
#include "taulines.h"
#include "tepowers.h"
#include "hydrogenic.h"
#include "casbhs.h"
#include "trace.h"
#include "phycon.h"
#include "elecden.h"
#include "hcbr.h"
#include "hcolrt.h"
#include "secondaries.h"
#include "hstat.h"
#include "hdeltat.h"
#include "hlteok.h"
#include "h.h"
#include "zonecnt.h"
#include "charexc.h"
#include "hctonoff.h"
#include "converge.h"
#include "cfit.h"
#include "smeets.h"
#include "atomcwgt.h"
#include "printe82.h"

void HydroCollid(long int ipZ)
{
	long int i, 
	  ipHi, 
	  ipLo, 
	  level, 
	  limit, 
	  lower, 
	  upper;
	double CStemp, 
	  factor ,
	  HCol2s2pp, 
	  ConvLTEPOP;
	static long nzUsed = -1;
	/* too few initializers for aggregate, but it defaults to all zeros */
	static float told[LIMELM] = {0.};

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

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

	/* evaluate collision rate for model hydrogen atom */
	/* ipZ is atomic number minus 1, so H is 0*/


	/***********************************************************************
	 *                                                                     *
	 * get secondary excitation by suprathermal electrons                  *
	 *                                                                     *
	 ***********************************************************************/

	if( ipZ == 0 )
	{
		/* following are for hydrogen ATOMS, ions are in next else */
		/* secondary excitation of Lyman lines */
		Secondaries.Hx12[0][IP2P] = (float)(Secondaries.x12tot*2./3./1.261);
		Secondaries.Hx12[0][IP2S] = (float)(Secondaries.x12tot*1./3./1.261);
		Secondaries.Hx12[0][3] = (float)(Secondaries.x12tot*0.165/1.261);
		Secondaries.Hx12[0][4] = (float)(Secondaries.x12tot*0.059/1.261);

		/* rest are just a wild guess */
		for( i=5; i <= nhlevel; i++ )
		{
			Secondaries.Hx12[0][i] = (float)(Secondaries.Hx12[0][i-1]/2.);
		}
	}
	else if( ipZ==1 )
	{
		/* following for ions, guess from Janev et al. 
		 * only evaluated for He, and there will be a problem if
		 * hydrogenic helium not present, but hydrogenic species
		 * of heavier elements are present.  I can see how this
		 * would ever happen */
		Secondaries.Hx12[1][IP2P] = (float)(Secondaries.x12tot/3.*2./3.);
		Secondaries.Hx12[1][IP2S] = (float)(Secondaries.x12tot/3./3.);
		/* this is just a wild guess */
		for( i=3; i <= nhlevel; i++ )
		{
			Secondaries.Hx12[1][i] = (float)(Secondaries.Hx12[1][i-1]/2.);
		}
	}

	/* skip most of this sub if temperature has not changed by much,
	 * the check on conv.nTotalIoniz is to make sure that we redo this
	 * on the very first call in a grid calc - it is 0 on the first call */
	/* loc1 bug 99 apr 08, tried commenting out first block so always reeval, no joy */
	if( fabs(phycon.te-told[ipZ])/phycon.te < 0.00002 && 
		hydro.hcol[ipZ][IP1S][IP2P] > 0. && conv.nTotalIoniz )
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, 
				"       HydroCollid called ipZ %li - no reval Boltz fac, LTE dens\n",
				ipZ );
		}
		/* (but still eval col below, search for 
		 * following code executed even when temperature did not change" ) */
	}
	else
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, 
				"       HydroCollid called ipZ %li - will reeval Boltz fac, LTE dens\n",
				ipZ );
		}
		told[ipZ] = phycon.te;

		/***********************************************************************
		 *                                                                     *
		 * Boltzmann factors for all levels, ionization and                    *
		 * exciation between levels                                            *
		 *                                                                     *
		 ***********************************************************************/

		/* define boltz facs between adjacent levels, will go from ipHi to one below */
		for( ipHi=IP2S; ipHi <= nhlevel; ipHi++ )
		{
			/* HdeltaTe is array of line temps, upper to lower
			 * their boltzmann factors for this temp,
			 * the array addressing for hlbolt is correct, and backwards from
			 * other similar variables - [up][lo] */
			hydro.hlbolt[ipZ][ipHi-1][ipHi] = 
				exp(-(double)(HydroLines[ipZ][ipHi][ipHi-1].EnergyK/phycon.te));
		}

		/* the highest level, up to the continuum */
		hydro.hcbolt[ipZ][nhlevel] = exp(-(double)(HdeltaT.HCionTe[ipZ][nhlevel]/
		  phycon.te));

		/* define remaining line Boltzmann factors in terms of existing factors between levels */
		for( ipLo=IP1S; ipLo <= (nhlevel - 2); ipLo++ )
		{
			for( ipHi=ipLo + 2; ipHi <= nhlevel; ipHi++ )
			{
				hydro.hlbolt[ipZ][ipLo][ipHi] = 
					hydro.hlbolt[ipZ][ipLo][ipHi-1]*
				  hydro.hlbolt[ipZ][ipHi-1][ipHi];
			}
		}

		for( i=IP1S; i <= (nhlevel - 1); i++ )
		{
			hydro.hcbolt[ipZ][i] = 
				hydro.hlbolt[ipZ][i][nhlevel]*
			  hydro.hcbolt[ipZ][nhlevel];
			/* this is to prevent close-to-overflow problems in model atom */
			if( hydro.hcbolt[ipZ][i] < SMALLDOUBLE )
			{
				hydro.hcbolt[ipZ][i] = 0.;
			}
		}

		/***********************************************************************
		 *                                                                     *
		 * LTE abundances for all levels, ionization and                       *
		 * exciation between levels                                            *
		 *                                                                     *
		 ***********************************************************************/

		/* following factor is actually 4.1412957e-16 (old=4.14158E-16), 
		 * but e- stat weight is in */
		/* ConvLTEPOP = 2.0708e-16/TePowers.te32;*/
		/* >>chng 99 Jun 02, use codata and infinite mass nuc */
		factor = HION_LTE_POP*AtomcWgt.AtomicWeight[ipZ]/
			(AtomcWgt.AtomicWeight[ipZ]+ELECTRON_MASS/ATOMIC_MASS_UNIT) ;

		ConvLTEPOP = pow(factor,1.5)/2./TePowers.te32;

		for( i=IP1S; i <= nhlevel; i++ )
		{
			/* >>>chng 99 mar 16, had gone to gt 0 previous week, but problems
			 * with very very small hcbolt, change back to 1e-100 */
			if( hydro.hcbolt[ipZ][i] > 1e-100 )
			/*if( hydro.hcbolt[ipZ][i] > 0. )*/
			{
				hcolrt.hlte[ipZ][i] = 
					hstat.HStatWght[i]/hydro.hcbolt[ipZ][i]* ConvLTEPOP;
			}
			else
			{
				hcolrt.hlte[ipZ][i] = 0.;
			}
		}

		/* now check for any zeros - if present then matrix cannot be used */
		HlteOk.lgHlteOk[ipZ] = TRUE;
		for( i=IP1S; i <= nhlevel; i++ )
		{
			if( hcolrt.hlte[ipZ][i] <= 0. )
			{
				HlteOk.lgHlteOk[ipZ] = FALSE;
				break;
			}
		}

		/***********************************************************************
		 *                                                                     *
		 * collisional ionization, and three body recombination                *
		 * Vriens and Smeets collisional ionization data set                   *
		 * use for hydrogen itselt                                             *
		 *                                                                     *
		 ***********************************************************************/
		if( ipZ == 0 )
		{
			/* atomic hydrogen
			 * this is total to 2s and 2p, will be split below
			 * HColOn is set to 0 with hydrogen collisions off command, usually 1 */
			hydro.hcolnc[ipZ][IP2P] = (float)(smeets(2)*hydro.HColIonOn);
			/* now do rest of levels, up to second highest */
			for( i=3; i < nhlevel; i++ )
			{
				/* HColOn is option to turn off collisions, "hydrogen collis off" comnd */
				hydro.hcolnc[ipZ][i] = (float)(smeets(i)*hydro.HColIonOn);
			}
			/* always leave highest level coupled to continuum - HColIonOn not there*/
			hydro.hcolnc[ipZ][nhlevel] = (float)(smeets(nhlevel));
		}
		else
		{
			/* hydrogenic ions */
			/* this is total to 2s and 2p */
			hydro.hcolnc[ipZ][IP2P] = (float)(HydColIon(2,ipZ)*hydro.HColIonOn);
			for( i=3; i < nhlevel; i++ )
			{
				/* HColOn is option to turn off collisions, "hydrogen collis off" comnd */
				hydro.hcolnc[ipZ][i] = (float)(HydColIon(i,ipZ)*hydro.HColIonOn);
			}
			/* always leave highest level coupled to continuum - HColIonOn not there*/
			hydro.hcolnc[ipZ][nhlevel] = (float)(HydColIon(nhlevel,ipZ));
		}

		/* total 2s and 2p was set in previous branch, now split into 2s and 2p */
		hydro.hcolnc[ipZ][IP2P] *= 0.75;
		hydro.hcolnc[ipZ][IP2S] = (float)(hydro.hcolnc[ipZ][IP2P]/3.);

		/* collisional ionization from ground */
		cfit(ipZ+1,1,phycon.te,&hydro.hcolnc[ipZ][IP1S]);

		/* option to kill collisional ionization with Hydro collis off */
		hydro.hcolnc[ipZ][IP1S] *= hydro.HColIonOn;
		hydro.hcolnc[ipZ][IP2S] *= hydro.HColIonOn;
		hydro.hcolnc[ipZ][IP2P] *= hydro.HColIonOn;

		/* hummer and storey option for case b kills coll ioniz from 1 and 2 */
		if( casbhs.lgCasBhs )
		{
			hydro.hcolnc[ipZ][IP1S] = 0.;
			hydro.hcolnc[ipZ][IP2P] = 0.;
			hydro.hcolnc[ipZ][IP2S] = 0.;
		}

		/***********************************************************************
		 * 
		 * collisional strengths for all hydrogen lines
		 * 
		 ************************************************************************/

		/* evaluate collisional deexcitation rates,
		 * Hydcs123 evaluates downward collision strength, which is stored in
		 * hydrogen lines vector
		 * Hydcs123, upon call 1=1s, 2=2s, 3=2p, 4=3 */
		for( ipLo=IP1S; ipLo <= (nhlevel - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi <= nhlevel; ipHi++ )
			{
				if( ipHi == 1 )
				{
					/* 2s - 1s */
					CStemp = Hydcs123(1,2,ipZ,'e');
				}
				else if( ipHi == 2 )
				{
					if( ipLo == 0 )
					{
						/* 2p - 1s */
						CStemp = Hydcs123(1,3,ipZ,'e');
					}
					else if( ipLo == 1 )
					{
						/* option to kill 2s2p collisions */
						if( hydro.lgC2psOn )
						{
							/* 2s - 2p for electrons ONLY */
							CStemp = Hydcs123(2,3,ipZ,'e');
							/* 2s - 2p for protons ONLY */
							HCol2s2pp = Hydcs123(2,3,ipZ,'p');
							/* trick, since will later be multiplied by electron,
							 * not proton density */
							/*CStemp += HCol2s2pp*h.hii/phycon.eden;*/
							/* NB - all collision rates will be multiplied by this number */
							CStemp += HCol2s2pp*h.hii/ElecDen.EdenHCorr;
						}
						else
						{
							CStemp = 0.;
						}
					}
					else
					{
						/* this is not possible */
						CStemp = -DBL_MAX;
					}
				}

				else if( ipHi == 3 )
				{
					if( ipLo == 0 )
					{
						/* 3 - 1 */
						CStemp = Hydcs123(1,4,ipZ,'e');
					}
					else if( ipLo == 1 )
					{
						/* 3 - 2s */
						CStemp = Hydcs123(2,4,ipZ,'e');
					}
					else if( ipLo == 2 )
					{
						/* 3 - 2p */
						CStemp = Hydcs123(3,4,ipZ,'e');
					}
					else
					{
						/* this is not possible */
						CStemp = -DBL_MAX;
					}
				}
				else
				{
					/* highly excited levels */
					CStemp = HydColDwn(ipHi,ipLo,ipZ+1);
				}
				HydroLines[ipZ][ipHi][ipLo].cs = (float)CStemp;
			}
		}

		/* save 2s-2p collisions since we always wan to include, even when other
		 * collisions are turned off, since 2s-2p act as one level.  there is a
		 * separate no 2s2p command, which is implemented below */
		CStemp = HydroLines[ipZ][IP2P][IP2S].cs;

		/* kill ALL collisions if this set with hydrogen collisional off set */
		if( hydro.HColExcOn == 0. )
		{
			for( ipLo=IP1S; ipLo <= (nhlevel-1); ipLo++ )
			{
				for( ipHi=ipLo + 1; ipHi <= nhlevel; ipHi++ )
				{
					HydroLines[ipZ][ipHi][ipLo].cs = 0.;
				}
			}
		}

		/* case b hummer&storey does not have any collisions from n=1 or 2 
		 * just kill those collisions here */
		if( casbhs.lgCasBhs )
		{
			for( ipLo=IP1S; ipLo <= IP2P; ipLo++ )
			{
				for( ipHi=3; ipHi <= nhlevel; ipHi++ )
				{
					HydroLines[ipZ][ipHi][ipLo].cs = 0.;
				}
			}
		}

		/* reset the 2s-2p collsions unless turned off with hydrogen collisions 2s2p
		 * command, since we usually want to include these even when line excit
		 * is turned off */
		if( hydro.lgC2psOn )
		{
			HydroLines[ipZ][IP2P][IP2S].cs = (float)CStemp;
		}

		/***********************************************************************
		 *                                                                     *
		 * collisional deexcitation for all lines                              *
		 *                                                                     *
		 ***********************************************************************/

		for( ipLo=IP1S; ipLo <= (nhlevel - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi <= nhlevel; ipHi++ )
			{
				/*This is downward collision rate */
				hydro.hcol[ipZ][ipLo][ipHi] = 
					(float)(HydroLines[ipZ][ipHi][ipLo].cs/
				  TePowers.sqrte*8.629e-6/hstat.HStatWght[ipHi]);
			}
		}

		if( (trace.lgTrace && trace.lgHBugFull) && (ipZ == trace.ipZTrace) )
		{
			fprintf( ioQQQ, "       HydroCollid rate n->n-1:" );
			for( i=IP1S; i <= (nhlevel - 1); i++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", hydro.hcol[ipZ][i][i+1] ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid col ion cof:" );
			for( i=IP1S; i <= nhlevel; i++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e",  hydro.hcolnc[ipZ][i] ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dn2 1s:" );
			for( i=IP1S; i <= nhlevel; i++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", hydro.hcol[ipZ][IP1S][i] ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dn2 2s:" );
			for( i=IP1S; i <= nhlevel; i++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", hydro.hcol[ipZ][IP2S][i] ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dn2 2p:" );
			for( i=IP1S; i <= nhlevel; i++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", hydro.hcol[ipZ][IP2P][i] ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dwn2 3:" );
			for( i=IP1S; i <= nhlevel; i++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", hydro.hcol[ipZ][3][i] ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dwn2 4:" );
			for( i=IP1S; i <= nhlevel; i++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", hydro.hcol[ipZ][4][i] ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dwn2 5:" );
			for( i=IP1S; i <= nhlevel; i++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", hydro.hcol[ipZ][5][i]) );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dwn2 6:" );
			for( i=IP1S; i <= nhlevel; i++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", hydro.hcol[ipZ][6][i] ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dwn2 7:" );
			for( i=IP1S; i <= nhlevel; i++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", hydro.hcol[ipZ][7][i]) );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid boltz fac c:" );
			for( i=IP1S; i <= nhlevel; i++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", hydro.hcbolt[ipZ][i] ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid HLTE(i,Z):  " );
			for( i=IP1S; i <= nhlevel; i++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", hcolrt.hlte[ipZ][i] ));
			}
			fprintf( ioQQQ, "\n" );
		}
	}

	/* following code executed even when temperature did not change*/

	for( level=IP1S; level <= nhlevel; level++ )
	{
		hcbr.hcoldn[ipZ][level] = 0.;
		for( lower=IP1S; lower <= (level - 1); lower++ )
		{
			hcbr.hcoldn[ipZ][level] += hydro.hcol[ipZ][lower][level];
		}

		for( upper=level + 1; upper <= nhlevel; upper++ )
		{
			hcbr.hcoldn[ipZ][level] += 
				(float)(hydro.hcol[ipZ][level][upper]*
				hstat.HStatWght[upper]/hstat.HStatWght[level]*
				hydro.hlbolt[ipZ][level][upper] );
		}
	}

	if( (trace.lgTrace && trace.lgHBugFull) && (ipZ == trace.ipZTrace) )
	{
		limit = MIN2(9,nhlevel);
		fprintf( ioQQQ, "       HydroCollid Z=%2ld", ipZ );
		fprintf( ioQQQ, "hcoldn  " );
		for( i=IP1S; i <= limit; i++ )
		{
			fprintf( ioQQQ,PrintEfmt("%10.3e", hcbr.hcoldn[ipZ][i] ));
		}
		fprintf( ioQQQ, "\n" );

		if( nhlevel > 9 )
		{
			fprintf( ioQQQ, "       HydroCollid     hcoldn  " );
			for( i=10; i <= nhlevel; i++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", hcbr.hcoldn[ipZ][i]) );
			}
			fprintf( ioQQQ, "\n" );
		}
		}

	if( ipZ == 0 )
	{
		/* hydrogen charge transfer charge exchange
		 * following only for hydrogen, only other array element is
		 * 2, all set zero in zero,
		 * HCTOnOff.HCTOn set to 0 with "no charge transfer" command */
		if( ZoneCnt.nzone <= 1 || 
			CharExc.CTHrec[0][1]*CharExc.CTHion[0][1] == 0. )
		{
			if( !conv.nTotalIoniz )
			{
				/* this is zone 0, the first time hydro is called
				 * ct rates have not been evaluated yet by the heavies */
				CharExc.CTHrec[0][1] = 0.;
				CharExc.CTHion[0][1] = 0.;
				CharExc.CTHrec[0][0] = 0.;
				CharExc.CTHion[0][0] = 0.;
			}
			else
			{
				CharExc.CTHrec[0][1] = CharExc.CTHrec[0][0]*HCTOnOff.HCTOn;
				CharExc.CTHion[0][1] = CharExc.CTHion[0][0]*HCTOnOff.HCTOn;
			}
		}
		else if( ZoneCnt.nzone != nzUsed )
		{
			CharExc.CTHrec[0][1] = (float)(CharExc.CTHrec[0][1]*0.75 + 0.25*CharExc.CTHrec[0][0]*
			  HCTOnOff.HCTOn);
			CharExc.CTHion[0][1] = (float)(CharExc.CTHion[0][1]*0.75 + 0.25*CharExc.CTHion[0][0]*
			  HCTOnOff.HCTOn);
		}
		nzUsed = ZoneCnt.nzone;
	}


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