/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*Hydrogenic main routine to call HydroLevel and determine model hydrogen atom level balance */
/* HydroRenorm - renormalize H so that it agrees with the chemistry */
#include "cddefines.h"
#include "taulines.h"
#include "iso.h"
#include "conv.h"
#include "trace.h"
#include "secondaries.h"
#include "hmi.h"
#include "h2.h"
#include "atmdat.h"
#include "dense.h"
#include "ionbal.h"
#include "hydrogenic.h"

void Hydrogenic(void)
{
	long int ipLo ,
		ipHi , 
		nelem,
		ipISO=ipH_LIKE,
		mol;
	double  coltot, 
	  gamtot;
	double sum , error;
	static int lgFinitePop[LIMELM];
	static int lgMustInit=TRUE;
	int lgH_chem_conv;
	int loop_H_chem;
	double solomon_assumed;
	double *PumpSave=NULL;

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

	if( lgMustInit )
	{
		for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
			lgFinitePop[nelem] = TRUE;
		lgMustInit = FALSE;
	}

	/* 
	 * This is an option to account for intrinsic absorption or emission line the lyman 
	 * lines.  This is done without changing the incident coarse continuum.  
	 * Check if continuum pumping of H lyman lines to be multiplied by scale factor
	 * hydro.xLymanPumpingScaleFactor is set with atom h-like lyman pumping scale command 
	 * Lya pump rate is always updated - this is simplest way to finese intrinsic absorption
	 * or emission 
	 */
	if( hydro.xLymanPumpingScaleFactor!=1.f )
	{
		if( (PumpSave = (double *)MALLOC( (unsigned)iso.numLevels_max[ipH_LIKE][ipHYDROGEN]*sizeof(double) ) ) ==NULL )
			BadMalloc();
		ipLo = 0;
		/* do not include the higest levels since pumping to them could create
			* artificial ionizaitons */
		for( ipHi=2; ipHi < iso.numLevels_max[ipH_LIKE][ipHYDROGEN]; ++ipHi )
		{
			PumpSave[ipHi] = EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].pump;
			EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].pump *= hydro.xLymanPumpingScaleFactor;
		}
	}

	for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
	{
		/* do not consider elements that have been turned off */
		if( dense.lgElmtOn[nelem] )
		{
			/* note that nelem scale is totally on c not physical scale, so 0 is h */
			/* evaluate hydrogenic balance if ionization reaches this high */
			if( (dense.IonHigh[nelem] == nelem + 1)  )
			{
				/* truncate atom if physical conditions limit the maximum principal quantum number of a
				 * bound electron to a number less than the malloc'd size */
				/* iso_continuum_lower( ipH_LIKE, nelem ); */

				/* evaluate collisional rates */
				HydroCollid(nelem);

				/* evaluate photoionization rates  */
				iso_photo( ipH_LIKE , nelem );

				/* evaluate recombination rates */
				HydroRecom(nelem);

				/* evaluate state specific creation and destruction processes,
				 * also define iso.xIonSimple */
				iso_ionize_recombine( ipH_LIKE , nelem );

				/* solve for the level populations */
				HydroLevel(nelem);

				/* say that we have set the populations */
				lgFinitePop[nelem] = TRUE;

				if( trace.lgTrace )
				{
					fprintf( ioQQQ, "       Hydrogenic Z=%2ld H2OVH1=",nelem);
					fprintf( ioQQQ,PrintEfmt("%10.3e", iso.pop_ion_ov_neut[ipH_LIKE][nelem]));
					fprintf( ioQQQ, " simple=");
					fprintf( ioQQQ,PrintEfmt("%10.3e", iso.pop_ion_ov_neut[ipH_LIKE][nelem]));
					/* fprintf( ioQQQ, " H-ovH1:");
						 fprintf( ioQQQ,PrintEfmt("%10.2e", hmi.hmovh1 ) ); */
					fprintf( ioQQQ, "\n"); 
				}

			}
			else if( lgFinitePop[nelem] )
			{
				/* this branch, pops were set previously, but now are zero,
				 * so must zero them out this time */
				lgFinitePop[nelem] = FALSE;

				iso.pop_ion_ov_neut[ipH_LIKE][nelem] = 0.;

				/* zero it out since no population*/
				for( ipHi=ipH2s; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
				{
					for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
					{
						/* population of lower level rel to ion */
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopLo = 0.;

						/* population of upper level rel to ion */
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi =  0.;

						/* population of lower level rel to ion, corrected for stim em */
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc =  0.;
					}
				}
			}
			/* if ionization has been set, hose what we have done and redefine */
			if( dense.lgSetIoniz[nelem] )
			{
				dense.xIonDense[nelem][nelem+1-ipISO] = dense.SetIoniz[nelem][nelem+1-ipISO]*dense.gas_phase[nelem];
				dense.xIonDense[nelem][nelem-ipISO] = dense.SetIoniz[nelem][nelem-ipISO]*dense.gas_phase[nelem];
				if( dense.SetIoniz[nelem][nelem+1-ipISO] > SMALLFLOAT )
				{
					iso.Pop2Ion[ipISO][nelem][ipH1s] = dense.SetIoniz[nelem][nelem-ipISO] / dense.SetIoniz[nelem][nelem+1-ipISO];
					EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].PopLo = dense.SetIoniz[nelem][nelem-ipISO] / dense.SetIoniz[nelem][nelem+1-ipISO];
				}
				else
				{
					iso.Pop2Ion[ipISO][nelem][ipH1s] = 0.;
					EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].PopLo = 0.;
				}
			}
		}
	}

	/* ============================================================================== */
	/* rest is for hydrogen only */

	/* this block appears redundant and could be removed? */
	/* >> 02 nov 21 rjrw -- xIonDense is used in hmole but only set in bidiag, 
	 * so we need this at least for the first iteration. */

	/* do molecular balance 
	 * hmovh1 will be ratio of molecular to atomic hydrogen
	 * HIonFrac is fraction of H that is ionized, ratio of ion to atom */

	lgH_chem_conv = FALSE;
	loop_H_chem = 0;
	while( loop_H_chem < 5 && !lgH_chem_conv )
	{
		/* save solomon rate that was assumed in the above - want to make sure that assumed
		 * solomon rate is stable, and does not change after call to large H2 molecule */
		solomon_assumed = hmi.H2_Solomon_dissoc_rate_used_H2g/SDIV(hmi.H2_rate_destroy);

		/* now do chem, this will reset hmi.H2_Solomon_dissoc_rate_used */
		hmole();
		/* now do level populations for H2 */
		H2_LevelPops();
		/* >>chng 05 feb 12 - add logic checking that solmon from big molecule is equal to
		 * rate used in H chem */
		lgH_chem_conv = TRUE;
		/* check whether solomon rate has changed */
		if( h2.lgH2ON  && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
		{
			/* >>chng 05 mar 11, go from "total" to H2g for solomon rate */
			if( fabs( solomon_assumed - hmi.H2_Solomon_dissoc_rate_BigH2_H2g/SDIV(hmi.H2_rate_destroy) ) > 
				conv.EdenErrorAllowed/5.)
			{
				lgH_chem_conv = FALSE;
			}
		}
		++loop_H_chem;
	}

	{
		/*@-redef@*/
		/* often the H- route is the most efficient formation mechanism for H2,
		 * will be through rate called ratach
		 * this debug print statement is to trace h2 oscillations */
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if(DEBUG_LOC )
		{
			fprintf(ioQQQ,"DEBUG\t%.2f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
				fnzone,
				hmi.H2_total ,
				dense.xIonDense[ipHYDROGEN][0],
				dense.xIonDense[ipHYDROGEN][1],
				hmi.H2_Solomon_dissoc_rate_used_H2g,
				hmi.H2_Solomon_dissoc_rate_BD96_H2g,
				hmi.H2_Solomon_dissoc_rate_TH85_H2g);
		}
	}
	/* >>chng 01 may 09, add option to force abundance, with element name ioniz cmmnd */
	if( dense.lgSetIoniz[ipHYDROGEN] )
	{
		dense.xIonDense[ipHYDROGEN][1] = dense.SetIoniz[ipHYDROGEN][1]*dense.gas_phase[ipHYDROGEN];
		dense.xIonDense[ipHYDROGEN][0] = dense.SetIoniz[ipHYDROGEN][0]*dense.gas_phase[ipHYDROGEN];

		/* >>chng 04 dec 03, add these too for this limit */
		iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s] = dense.SetIoniz[ipHYDROGEN][0] / SDIV( dense.SetIoniz[ipHYDROGEN][1]);
		EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopLo = dense.SetIoniz[ipHYDROGEN][0] / SDIV( dense.SetIoniz[ipHYDROGEN][1]);
	}
	else
	{
		/* 
		 * >> chng 03 jan 15 rjrw:- terms are now in ion_solver, to allow for
		 * molecular sources and sinks of H and H+.  ion_solver renormalizes
		 * to keep the total H abundance correct -- only the molecular
		 * network is allowed to change this. 
		 */

		ion_solver( ipHYDROGEN , FALSE );

	}

	/* confirm that species still add up correctly */
	/* this exposes a "leak" that occurs somewhere, almost certainly in hmole
	 * if DEBUG_LOC is set TRUE in the following there will be comments printing
	 * due to a loss of some protons */

	sum = dense.xIonDense[ipHYDROGEN][0] + dense.xIonDense[ipHYDROGEN][1];
	for (mol=0;mol<N_H_MOLEC;mol++) 
	{
		sum += hmi.Hmolec[mol]*hmi.nProton[mol];
	}
	/* do not double count H0 and H+ */
	sum -=  hmi.Hmolec[ipMH]+hmi.Hmolec[ipMHp];

	error = ( dense.gas_phase[ipHYDROGEN] - sum )/dense.gas_phase[ipHYDROGEN];
	{
		/*@-redef@*/
		/* often the H- route is the most efficient formation mechanism for H2,
		 * will be through rate called ratach
		 * this debug print statement is to trace h2 oscillations */
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if(DEBUG_LOC && (fabs(error) > 1e-4) )
			fprintf(ioQQQ,"PROBLEM hydrogenic zone %li hden %.4e, sum %.4e (h-s)/h %.3e \n", nzone, dense.gas_phase[ipHYDROGEN] , sum , 
				error );
	}

	/* >>hcng 05 mar 24,
	 * renormalize the populations and emission of H atom to agree with chemistry */
	HydroRenorm();

	/* this is test whether collisions are important first define ratio of exits
	 * from ground due to coll, rel to everthing else, then flag is large */
	if( iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH1s] > 1e-30 )
	{
		coltot = 
			EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ColUL*
			iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH2p]/iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH1s]*
		  dense.eden*(iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH2p] + iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH2p] + 
		  EmisLines[ipH_LIKE][ipHYDROGEN][3][ipH2p].ColUL*iso.PopLTE[ipH_LIKE][ipHYDROGEN][3]/iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH2p])/
		  (EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ColUL*dense.eden + 
		  EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Aul*
		  (EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pesc + 
		  EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pelec_esc +
		  EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest) );

		/* caution that collisions are important (and so only small changes in
		 * temperature should happen) if more than 20% of the total */
		if( coltot > 0.2 )
		{
			hydro.lgHColionImp = TRUE;
		}
	}
	else
	{
		hydro.lgHColionImp = FALSE;
	}

	/* remember the ratio of pops of 2p to 1s for possible printout in prtComment
	 * and to obtain Lya excitation temps.  the pop of ground is not defined if
	 * NO ionization at all since these pops are relative to ion */
	/* >>chng 99 jun 03, added MAX2 to protect against totally neutral gas */
	if( iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2p]/MAX2(SMALLDOUBLE,iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]) > 0.1 &&
		iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s] > SMALLDOUBLE )
	{
		hydro.lgHiPop2 = TRUE;
		hydro.pop2mx = (float)MAX2(iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2p]/iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s],
		  hydro.pop2mx);
	}

	gamtot = iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s] + secondaries.csupra[ipHYDROGEN][0];

	coltot = iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH1s] + 
	  EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ColUL*4.*
	  iso.Boltzmann[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s];

	/* if ground state destruction rate is significant, recall different dest procceses */
	if( iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s] > SMALLFLOAT )
	{
		hydro.H_ion_frac_photo = 
			(float)(iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]/iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s] );

		/* fraction of ionizations of H from ground, due to thermal collisions */
		hydro.H_ion_frac_collis = 
			(float)(iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.eden/iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s]);

		/* this flag is used in ConvBase to decide whether we
		 * really need to converge the secondary ionization rates */
		secondaries.sec2total = 
			(float)(secondaries.csupra[ipHYDROGEN][0] / iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s]) ;

		/* frac of ionizations due to ct */
		atmdat.HIonFrac = atmdat.HCharExcIonTotal / iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s];
	}
	else
	{
		hydro.H_ion_frac_collis = 0.;
		hydro.H_ion_frac_photo = 0.;
		secondaries.sec2total = 0. ;
		atmdat.HIonFrac = 0.;
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "       Hydrogenic retrn %.2f ",fnzone);
		fprintf(ioQQQ,"H0:%.3e ", dense.xIonDense[ipHYDROGEN][0]);
		fprintf(ioQQQ,"H+:%.3e ", dense.xIonDense[ipHYDROGEN][1]);
		fprintf(ioQQQ,"H2:%.3e ", hmi.H2_total);
		fprintf(ioQQQ,"H-:%.3e ", hmi.Hmolec[ipMHm]);
		fprintf(ioQQQ,"ne:%.3e ", dense.eden);
		fprintf( ioQQQ, " REC, COL, GAMT= ");
		/* recomb rate coef, cm^3 s-1 */
		fprintf(ioQQQ,"%.2e ", iso.RadRec_effec[ipH_LIKE][ipHYDROGEN] );
		fprintf(ioQQQ,"%.2e ", coltot);
		fprintf(ioQQQ,"%.2e ", gamtot);
		fprintf( ioQQQ, " CSUP=");
		PrintE82( ioQQQ, secondaries.csupra[ipHYDROGEN][0]);
		fprintf( ioQQQ, "\n");
	}

	if( hydro.xLymanPumpingScaleFactor!=1.f )
	{
		ipLo = 0;
		/* do not include the higest levels since pumping to them could create
		 * artificial ionizaitons */
		for( ipHi=2; ipHi < iso.numLevels_max[ipH_LIKE][ipHYDROGEN]; ++ipHi )
		{
			EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].pump = PumpSave[ipHi];
		}
		free(PumpSave);
	}

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

/* HydroRenorm - renormalize H so that it agrees with the chemistry */
void HydroRenorm( void )
{
	double sum_atom_iso , renorm;
	long int level,
		nelem,
		ipHi ,
		ipLo;

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

	/*>>chng 04 mar 23, add this renorm */
	/* renormalize the state specific populations, so that they
	 * add up to the results that came from ion_solver 
	 * units at first is sum div by H+ density, since Pop2Ion defn this way */
	sum_atom_iso = 0.;
	for( level=ipH1s; level < iso.numLevels_max[ipH_LIKE][ipHYDROGEN]; level++ )
	{
		sum_atom_iso += iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][level];
	}
	/* now convert to cm-3 - this is total population in iso solved model */
	sum_atom_iso *= dense.xIonDense[ipHYDROGEN][1];

	/* >>chng 04 may 25, humunculus sum_atom_iso is zero */
	if( sum_atom_iso > SMALLFLOAT )
	{
		renorm = dense.xIonDense[ipHYDROGEN][0] / sum_atom_iso;
	}
	else
	{
		renorm = 0.;
	}
	ASSERT( renorm < BIGFLOAT );
	/*fprintf(ioQQQ,"DEBUG renorm\t%.2f\t%.3e\n",fnzone, renorm);*/
	/* renormalize populations from iso-model atom so that they agree with ion solver & chemistry */
	iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s] *= renorm;
	/*fprintf(ioQQQ,"DEBUG h \t%.3e hydrogenic renorm %.3e\n", 
		iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s] ,
		iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]/renorm );*/

	nelem = ipHYDROGEN;
	for( ipHi=ipH2s; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
	{
		iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipHi] *= renorm;
		/*>>chng 04 dec 06, renorm these pops too */
		for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
		{
			/* population of lower level rel to ion */
			EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopLo *= renorm;

			/* population of upper level rel to ion */
			EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi *= renorm;

			/* population of lower level rel to ion, corrected for stim em */
			EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc *= renorm;

			/*ASSERT( EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc < BIGFLOAT );*/
		}
	}

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

	return;
}
