/* This file is part of Cloudy and is copyright (C) 1978-2004 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 */
#include "cddefines.h"
#include "taulines.h"
#include "iso.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,
		mol;
	double  coltot, 
	  gamtot;
	double sum , error;
	double sum_excit , renorm;
	long int level;
	static int lgFinitePop[LIMELM];
	static int lgMustInit=TRUE;

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

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

	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)  )
			{
				/* 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[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.;
					}
				}
			}
		}
	}

	/* ============================================================================== */
	/* 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. */

#if 0
	/* update neutral and ion densities for current pop ratio and molecules */
	fprintf(ioQQQ,"<<%g %g %g\n",
					(float)(dense.gas_phase[ipHYDROGEN]/(1. + iso.pop_ion_ov_neut[ipH_LIKE][ipHYDROGEN])),
					dense.gas_phase[ipHYDROGEN]*iso.pop_ion_ov_neut[ipH_LIKE][ipHYDROGEN]/
					(1. + iso.pop_ion_ov_neut[ipH_LIKE][ipHYDROGEN]),
					hmi.hmovh1);
	if(1 || conv.nPres2Ioniz > 0 )
	{
		dense.xIonDense[ipHYDROGEN][1] = (float)(dense.gas_phase[ipHYDROGEN]*iso.pop_ion_ov_neut[ipH_LIKE][ipHYDROGEN]/
			(1. + iso.pop_ion_ov_neut[ipH_LIKE][ipHYDROGEN] + hmi.hmovh1));

		dense.xIonDense[ipHYDROGEN][0] = (float)(dense.gas_phase[ipHYDROGEN]/(1. + iso.pop_ion_ov_neut[ipH_LIKE][ipHYDROGEN] + 
			hmi.hmovh1));
	}
	fprintf(ioQQQ,">>%g %g\n",dense.xIonDense[ipHYDROGEN][0],dense.xIonDense[ipHYDROGEN][1]);
#endif

	/* 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 */
	
	if (0) 
	{
		fprintf(ioQQQ,"Ions into Hmole %g %g\n",dense.xIonDense[ipHYDROGEN][0],dense.xIonDense[ipHYDROGEN][1]); 
	}

	hmole();
	{
		/*@-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,
				hmi.H2_Solomon_dissoc_rate_BD96,
				hmi.H2_Solomon_dissoc_rate_TH85);
		}
	}
	/* >>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];
	}
	else
	{
		/* 
		 * >> chng 03 jan 15 rjrw:- terms are now in bidiag, to allow for
		 * molecular sources and sinks of H and H+.  Bidiag 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 );
	}

	/*>>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 */
	sum_excit = 0.;
	for( level=ipH1s; level < iso.numLevels[ipH_LIKE][ipHYDROGEN]; level++ )
	{
		sum_excit += iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][level];
	}
	sum_excit *= dense.xIonDense[ipHYDROGEN][1];

	/* >>chng 04 may 25, humunculus sum_excit is zero */
	if( sum_excit > SMALLFLOAT )
	{
		renorm = dense.xIonDense[ipHYDROGEN][0] / sum_excit;
	}
	else
	{
		renorm = 0.;
	}
	/*fprintf(ioQQQ,"DEBUG renorm\t%.2f\t%.3e\n",fnzone, renorm);*/
	for( level=ipH1s; level < iso.numLevels[ipH_LIKE][ipHYDROGEN]; level++ )
	{
		iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][level] *= renorm;
	}

	/* now do level populations for H2 */
	H2_LevelPops();

	/* 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");
	}

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