/* 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 */
/*iso_ionize_recombine find state specific creation and destruction rates for iso sequences */
/*ChargeTransferUpdate update rate of ct ionization and recombination for H atoms */
#include "cddefines.h"
#include "ionbal.h"
#include "phycon.h"
#include "atmdat.h"
#include "dense.h"
#include "grainvar.h"
#include "hmi.h"
#include "hmrate.h"
#include "co.h"
#include "secondaries.h"
#include "iso.h"
/*lint -e778 constant expression evaluates to zero - in hmrate macro */
/*ChargeTransferUpdate update rate of ct ionization and recombination for H atoms */
static void ChargeTransferUpdate(void)
{
	long int ion;
	long int nelem;
	/* find total rate for charge transfer ionization of hydrogen,
	 * recombination for other element is ionization of hydrogen */
	atmdat.HCharExcIonTotal = 0.;
	for( nelem=ipHELIUM; nelem<LIMELM; ++nelem)
	{
		/* ion on the C scale, 0 is atom, so goes up to nelem+1,
		 * for helium nelem=1, ions go from 1 to 2 */
		for( ion=1; ion<=nelem+1; ++ion )
		{
			double one;
			/* we intentionally skip CT with O+ since this is in hmole */
			if( (nelem == ipOXYGEN) && (ion == 1) ) 
				continue;
			/* charge transfer ionization of H, recombination for other species */
			one = atmdat.HCharExcRecTo[nelem][ion-1]*dense.xIonDense[nelem][ion];
			atmdat.HCharExcIonTotal += one;
		}
	}

	/* >>chng 01 may 07,  add this in */
	/* charge transfer recombination of hydrogen,
	 * which is ionization of the heavy element */
	atmdat.HCharExcRecTotal = 0.;
	for( nelem=ipHELIUM; nelem<LIMELM; ++nelem)
	{
		/* this is ion on the abundances scale, 1 is atom, so goes up to nelem+1,
		 * for helium nelem=1, ion must go up to 3 */
		for( ion=0; ion<=nelem; ++ion )
		{
			/* option to skip Oo => O+ */
			if( (nelem == ipOXYGEN) && (ion == 0) ) 
				continue;
			/* charge transfer ionization of H, recombination for other species */
			atmdat.HCharExcRecTotal += 
				atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[nelem][ion];
		}
	}

	/* >>logic checked 02 may 02, think it's right */
	/* add on charge transfer ionization of helium,
	 * recombination for other element is ionization of helium */
	atmdat.HeCharExcIonTotal = 0.;
	/* loop up from Li, the next heavier element */
	for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem)
	{
		double hold_one = 0.;
		/* check that array bounds not exceeded */
		/* ion on the C scale, 0 is atom, so goes up to nelem+1,
		 * for helium nelem=1, ions go from 1 to 2 */
		for( ion=1; ion<=nelem+1; ++ion )
		{
			/* charge transfer ionization of He, recombination for other species */
			hold_one = atmdat.HeCharExcRecTo[nelem][ion-1]*dense.xIonDense[nelem][ion];
			atmdat.HeCharExcIonTotal += hold_one;
		}
	}

	/* >>chng 04 jul 02,
	 * add on charge transfer reactions of He-H */
	atmdat.HeCharExcIonTotal += atmdat.HCharExcIonOf[ipHELIUM][0]*dense.xIonDense[ipHYDROGEN][1];

	/* charge transfer recombination of helium,
	 * which is ionization of the heavy element */
	atmdat.HeCharExcRecTotal = 0.;
	for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem)
	{
		/* this is ion on the physics scale, 1 is atom, so goes up to nelem+1,
		 * for helium nelem=1, ion must go up to 3 */
		for( ion=0; ion<=nelem; ++ion )
		{
			/* charge transfer recombination of He, ionization for other species */
			atmdat.HeCharExcRecTotal += 
				atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[nelem][ion];
		}
	}
	/* >>chng 04 jul 02
	 * Add on charge transfer reactions of He+ +H0 -> He0 + H+ */
	atmdat.HeCharExcRecTotal += atmdat.HCharExcRecTo[ipHELIUM][0]*dense.xIonDense[ipHYDROGEN][0];

	return;
}

/*iso_ionize_recombine find state specific creation and destruction rates for iso sequences */
void iso_ionize_recombine(
	/* iso sequence, hydrogen or helium for now */
	long ipISO , 
	/* the chemical element, 0 for hydrogen */
	long int nelem )
{
	long int level,
		ion;
	double Recom3Body,
		sum;

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

	/* get total charge transfer ionization rate if this is hydrogen itself,
	 * routine also does CT for helium but all models must have hydrogen, so
	 * he is done at this stage */
	if( ipISO==ipH_LIKE && nelem == ipHYDROGEN )
		ChargeTransferUpdate();

	/* find recombination and ionization elements, will use to get simple estimate
	 * of ionization ratio below */
	for( level=ipH1s; level< iso.numLevels[ipISO][nelem]; ++level)
	{
		/* all process moving level to continuum, units s-1 */
		iso.RateLevel2Cont[ipISO][nelem][level] = iso.gamnc[ipISO][nelem][level] + 
		  iso.ColIoniz[ipISO][nelem][level]* dense.EdenHCorr + 
		  secondaries.csupra[nelem][nelem-ipISO];

		/* all processes from continuum to level n, units s-1 */
		iso.RateCont2Level[ipISO][nelem][level] = (
			/* radiative recombination */
			iso.RadRecomb[ipISO][nelem][level][ipRecRad]*
			iso.RadRecomb[ipISO][nelem][level][ipRecNetEsc] + 

			/* induced recombination */
			iso.RecomInducRate[ipISO][nelem][level]*iso.PopLTE[ipISO][nelem][level] + 

			/* collisional or three body recombination */
			/* PopLTE(level,nelem) is only LTE pop when mult by n_e n_H */
			iso.ColIoniz[ipISO][nelem][level]*dense.EdenHCorr*iso.PopLTE[ipISO][nelem][level]
			) * dense.eden;

	}

	/* all following go into or out of ground state */
	level = 0;

	/* grain charge transfer recombination and ionization, 
	 * assume goes into and comes from ground state */
	/* first grain surface losses - both to lower and higher stages of ionization */
	sum = 0.;
	for(ion=0; ion<nelem+1; ++ion )
		if( ion!=nelem-ipISO )
			sum += gv.GrainChTrRate[nelem][nelem-ipISO][ion];
	/*>>chng 04 sep 08, replaced grain ct rate from old sum with current sum */
	iso.RateLevel2Cont[ipISO][nelem][level] += sum;

	/* >>chng 04 apr 10, add formally correct grain surface ionization and recombination - 
	 * this HAD NOT been included for H-like species */
	/* >>chng 04 sep 08, use explicit term from higher stage of ionization - note that added
	 * GrainCreat[nelem][nelem-ipISO] for case of atomic He incorrectly added two-stage recombination
	 * to existing matrix logic */
	/*iso.RateCont2Level[ipISO][nelem][level] += ionbal.GrainCreat[nelem][nelem-ipISO];*/
	iso.RateCont2Level[ipISO][nelem][level] += gv.GrainChTrRate[nelem][nelem+1-ipISO][nelem-ipISO];

	/* now charge transfer - all into/from ground, two cases, H and not H */
	if( nelem==ipHYDROGEN )
	{
		/* charge transfer, hydrogen onto everything else */
		/* charge exchange ionization */
		iso.RateLevel2Cont[ipISO][nelem][ipH1s] += atmdat.HCharExcIonTotal;
		/* charge exchange recombination */
		iso.RateCont2Level[ipISO][nelem][ipH1s] += atmdat.HCharExcRecTotal;
	}
	else if( nelem==ipHELIUM && ipISO==ipHE_LIKE )
	{
		/* add CO here, only for he itself, 
		 * destruction of CO but recombination for he */
		/* >>chng 04 feb 14, update date to 1.4ee-9 from 1.1e-9,
		 * make var so auto keep parallel with co.c */

		/* The following terms in co.c destroy He+ and thereby need to be included 
		 * here.  Also included in the sum of co.hep_destroy are the contributions
		 * to the recombination due to hmole_step.*/

		co.hep_destroy = 0.;

		co.rate_HeP_CH_CHP_He           = (float)hmrate(0.0000000005,0,0)*co.hevmol[ipCH];
		co.rate_HeP_H2O_H2OP_He         = (float)hmrate(0.0000000000605,0,0)*co.hevmol[ipH2O];
		co.rate_HeP_O2_O2P_He           = (float)hmrate(0.000000000033,0,0)*co.hevmol[ipO2];
		co.rate_HeP_Si_SiP_He           = (float)hmrate(0.0000000033,0,0)*co.hevmol[ipATSI];
		co.rate_HeP_CH_CP_He_H          = (float)hmrate(0.0000000011,0,0)*co.hevmol[ipCH];
		co.rate_HeP_CH2_CHP_He_H        = (float)hmrate(0.00000000075,0,0)*co.hevmol[ipCH2];
		co.rate_HeP_OH_OP_He_H          = (float)hmrate(0.0000000011,0,0)*co.hevmol[ipOH];
		co.rate_HeP_H2O_OHP_He_H        = (float)hmrate(0.000000000286,0,0)*co.hevmol[ipH2O];
		co.rate_HeP_SiH_SiP_He_H        = (float)hmrate(0.0000000018,0,0)*co.hevmol[ipSIH];
		co.rate_HeP_H2O_OH_He_HP        = (float)hmrate(0.000000000204,0,0)*co.hevmol[ipH2O];
		co.rate_HeP_CH2_CP_He_H2        = (float)hmrate(0.00000000075,0,0)*co.hevmol[ipCH2];
		co.rate_HeP_O2_OP_O_He          = (float)hmrate(0.000000001,0,0)*co.hevmol[ipO2];
		co.rate_HeP_SiO_SiP_O_He        = (float)hmrate(0.00000000086,0,0)*co.hevmol[ipSIO];
		co.rate_HeP_SiO_Si_OP_He        = (float)hmrate(0.00000000086,0,0)*co.hevmol[ipSIO];
		co.rate_HeP_CH3_CHP_He          = (float)hmrate(0.0000000018,0,0)*co.hevmol[ipCH3];
		co.rate_HeP_CH4_CHP_He          = (float)hmrate(0.00000000024,0,0)*co.hevmol[ipCH4];
		co.rate_HeP_CH4_CH2P_He         = (float)hmrate(0.00000000095,0,0)*co.hevmol[ipCH4];
		co.rate_HeP_CH4_CH3_He          = (float)hmrate(0.00000000048,0,0)*co.hevmol[ipCH4];
		co.rate_HeP_CH4_CH3P_He         = (float)hmrate(0.000000000085,0,0)*co.hevmol[ipCH4];
		co.rate_HeP_CH4_CH4P_He         = (float)hmrate(0.000000000051,0,0)*co.hevmol[ipCH4];
		co.rate_HeP_NH_NP_He_H		 = (float)hmrate(0.0000000011,0,0)*co.hevmol[ipNH];
		co.rate_HeP_NH2_NP_He_H2    = (float)hmrate(0.0000000008,0,0)*co.hevmol[ipNH2];
		co.rate_HeP_NH2_NHP_He_H    = (float)hmrate(0.0000000008,0,0)*co.hevmol[ipNH2];
		co.rate_HeP_NH3_NHP_He_H2   = (float)hmrate(0.000000000176,0,0)*co.hevmol[ipNH3];
		co.rate_HeP_NH3_NH2P_He_H   = (float)hmrate(0.00000000176,0,0)*co.hevmol[ipNH3];
		co.rate_HeP_CN_N_CP_He      = (float)hmrate(0.00000000088,0,0)*co.hevmol[ipCN];
		co.rate_HeP_CN_NP_C_He      = (float)hmrate(0.00000000088,0,0)*co.hevmol[ipCN];
		co.rate_HeP_HCN_N_CHP_He    = (float)hmrate(0.000000000651,0,0)*co.hevmol[ipHCN];
		co.rate_HeP_HCN_NP_CH_He    = (float)hmrate(0.000000000217,0,0)*co.hevmol[ipHCN];
		co.rate_HeP_HCN_N_CP_He_H   = (float)hmrate(0.000000000775,0,0)*co.hevmol[ipHCN];
		co.rate_HeP_HCN_CNP_He_H    = (float)hmrate(0.00000000146,0,0)*co.hevmol[ipHCN];
		co.rate_HeP_N2_NP_N_He      = (float)hmrate(0.00000000096,0,0)*co.hevmol[ipN2];
		co.rate_HeP_NO_OP_N_He      = (float)hmrate(0.0000000002,0,0)*co.hevmol[ipNO];
		co.rate_HeP_NO_O_NP_He      = (float)hmrate(0.0000000014,0,0)*co.hevmol[ipNO];
		co.rate_HeP_HNO_NOP_He_H    = (float)hmrate(0.000000001,0,0)*co.hevmol[ipHNO];
		co.rate_HeP_HNO_NO_He_HP    = (float)hmrate(0.000000001,0,0)*co.hevmol[ipHNO];
		co.rate_HeP_HS_SP_He_H      = (float)hmrate(0.0000000017,0,0)*co.hevmol[ipHS];
		co.rate_HeP_OCN_CN_OP_He    = (float)hmrate(0.000000003,0,0)*co.hevmol[ipOCN];
		co.rate_HeP_OCN_CNP_O_He    = (float)hmrate(0.000000003,0,0)*co.hevmol[ipOCN];
		co.rate_HeP_SIN_SIP_N_He    = (float)hmrate(0.000000002,0,0)*co.hevmol[ipSIN];
		co.rate_HeP_N2O_N2_OP_He    = (float)hmrate(0.000000000276,0,0)*co.hevmol[ipN2O];
		co.rate_HeP_N2O_N2P_O_He    = (float)hmrate(0.00000000124,0,0)*co.hevmol[ipN2O];
		co.rate_HeP_N2O_NO_NP_He    = (float)hmrate(0.0000000003,0,0)*co.hevmol[ipN2O];
		co.rate_HeP_N2O_NOP_N_He    = (float)hmrate(0.000000000483,0,0)*co.hevmol[ipN2O];
		co.rate_HeP_CS_SP_C_He      = (float)hmrate(0.0000000013,0,0)*co.hevmol[ipCS];
		co.rate_HeP_CS_S_CP_He      = (float)hmrate(0.0000000013,0,0)*co.hevmol[ipCS];
		co.rate_HeP_NS_S_NP_He      = (float)hmrate(0.0000000012,0,0)*co.hevmol[ipNS];
		co.rate_HeP_NS_SP_N_He      = (float)hmrate(0.0000000012,0,0)*co.hevmol[ipNS];
		co.rate_HeP_SO_S_OP_He      = (float)hmrate(0.00000000083,0,0)*co.hevmol[ipSO];
		co.rate_HeP_SO_SP_O_He      = (float)hmrate(0.00000000083,0,0)*co.hevmol[ipSO];
		co.rate_HeP_OCS_SP_CO_He    = (float)hmrate(0.00000000076,0,0)*co.hevmol[ipOCS];
		co.rate_HeP_OCS_S_COP_He    = (float)hmrate(0.00000000076,0,0)*co.hevmol[ipOCS];
		co.rate_HeP_OCS_CSP_O_He    = (float)hmrate(0.00000000076,0,0)*co.hevmol[ipOCS];
		co.rate_HeP_OCS_CS_OP_He    = (float)hmrate(0.00000000076,0,0)*co.hevmol[ipOCS];
		co.rate_HeP_S2_SP_S_He      = (float)hmrate(0.000000002,0,0)*co.hevmol[ipOCS];
		co.rate_HeP_NH3_NH3P_He     = (float)hmrate(0.000000000264,0,0)*co.hevmol[ipNH3];
		co.rate_HeP_N2_N2P_He       = (float)hmrate(0.00000000064,0,0)*co.hevmol[ipN2];

		

		co.hep_destroy = co.rate_HeP_CH_CHP_He    +
						co.rate_HeP_H2O_H2OP_He	     + 		co.rate_HeP_O2_O2P_He    +       
						co.rate_HeP_Si_SiP_He        +		co.rate_HeP_CH_CP_He_H   +    
						co.rate_HeP_CH2_CHP_He_H     +		co.rate_HeP_OH_OP_He_H   +    
						co.rate_HeP_H2O_OHP_He_H     +		co.rate_HeP_SiH_SiP_He_H +   
						co.rate_HeP_H2O_OH_He_HP     +		co.rate_HeP_CH2_CP_He_H2 + 
						co.rate_HeP_O2_OP_O_He       +		co.rate_HeP_SiO_SiP_O_He +
						co.rate_HeP_SiO_Si_OP_He     +		co.rate_HeP_CH3_CHP_He   +
						co.rate_HeP_CH4_CHP_He       +		co.rate_HeP_CH4_CH2P_He  +
						co.rate_HeP_CH4_CH3_He       +		co.rate_HeP_CH4_CH3P_He  +
						co.rate_HeP_CH4_CH4P_He		 +      co.rate_HeP_NH_NP_He_H	 +
						co.rate_HeP_NH2_NP_He_H2     +		co.rate_HeP_NH2_NHP_He_H    +
						co.rate_HeP_NH3_NHP_He_H2   +		co.rate_HeP_NH3_NH2P_He_H   +
						co.rate_HeP_CN_N_CP_He      +		co.rate_HeP_CN_NP_C_He      +
						co.rate_HeP_HCN_N_CHP_He    +		co.rate_HeP_HCN_NP_CH_He    +
						co.rate_HeP_HCN_N_CP_He_H   +		co.rate_HeP_HCN_CNP_He_H    +
						co.rate_HeP_N2_NP_N_He      +		co.rate_HeP_NO_OP_N_He      +
						co.rate_HeP_NO_O_NP_He      +		co.rate_HeP_HNO_NOP_He_H    +
						co.rate_HeP_HNO_NO_He_HP    +		co.rate_HeP_HS_SP_He_H      +
						co.rate_HeP_OCN_CN_OP_He    +		co.rate_HeP_OCN_CNP_O_He    +
						co.rate_HeP_SIN_SIP_N_He    +		co.rate_HeP_N2O_N2_OP_He    +
						co.rate_HeP_N2O_N2P_O_He    +		co.rate_HeP_N2O_NO_NP_He    +
						co.rate_HeP_N2O_NOP_N_He    +		co.rate_HeP_CS_SP_C_He      +
						co.rate_HeP_CS_S_CP_He      +		co.rate_HeP_NS_S_NP_He      +
						co.rate_HeP_NS_SP_N_He      +		co.rate_HeP_SO_S_OP_He      +
						co.rate_HeP_SO_SP_O_He      +		co.rate_HeP_OCS_SP_CO_He    +
						co.rate_HeP_OCS_S_COP_He    +		co.rate_HeP_OCS_CSP_O_He    +
						co.rate_HeP_OCS_CS_OP_He    +		co.rate_HeP_S2_SP_S_He      +
						co.rate_HeP_NH3_NH3P_He     +		co.rate_HeP_N2_N2P_He       +
						co.hevmol[ipCO]*co.rate_co_hep_o_cp_he +  
						hmi.H2_total*hmi.rheph2hpheh +      hmi.heph2heh2p*hmi.H2_total;
			
		iso.RateCont2Level[ipISO][nelem][ipHe1s1S] += co.hep_destroy;

		/* this is ioniz of He due to ct with all other gas constituents */
		iso.RateLevel2Cont[ipISO][nelem][ipHe1s1S] += atmdat.HeCharExcIonTotal;
		/* this is recom of He due to ct with all other gas constituents */
		iso.RateCont2Level[ipISO][nelem][ipHe1s1S] += atmdat.HeCharExcRecTotal;
	}
	else
	{
		iso.RateCont2Level[ipISO][nelem][level] +=
			atmdat.HeCharExcRecTo[nelem][nelem-ipISO]*iso.Pop2Ion[ipHE_LIKE][ipHELIUM][ipHe1s1S]*dense.xIonDense[ipHELIUM][1] + 
			atmdat.HCharExcRecTo[nelem][nelem-ipISO]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.xIonDense[ipHYDROGEN][1];

		iso.RateLevel2Cont[ipISO][nelem][level] +=
			atmdat.HeCharExcIonOf[nelem][nelem-ipISO]*dense.xIonDense[ipHELIUM][1] +
			atmdat.HCharExcIonOf[nelem][nelem-ipISO]*dense.xIonDense[ipHYDROGEN][1];
	}


	/* now create sums of recombination and ionization rates */
	ionbal.RateRecomTot[nelem][nelem-ipISO] = 0.;
	ionbal.RadRecomRateCoef[nelem][nelem-ipISO] = 0.;
	Recom3Body = 0.;
	for( level=ipH1s; level< iso.numLevels[ipISO][nelem]; ++level)
	{

		/* units of ionbal.RateRecomTot are s-1,
		 * equivalent ionization term is done after level populations are known */
		ionbal.RateRecomTot[nelem][nelem-ipISO] += iso.RateCont2Level[ipISO][nelem][level];

		/* just the radiative recombination rate coef, cm3 s-1 */
		ionbal.RadRecomRateCoef[nelem][nelem-ipISO] += iso.RadRecomb[ipISO][nelem][level][ipRecRad]*
			iso.RadRecomb[ipISO][nelem][level][ipRecNetEsc];

		ASSERT( ionbal.RadRecomRateCoef[nelem][nelem-ipISO]> 0. );

		/* this is three-body recombination rate coef by itself - 
		 * need factor of eden to become rate */
		Recom3Body += iso.ColIoniz[ipISO][nelem][level]*dense.EdenHCorr*iso.PopLTE[ipISO][nelem][level];
	}

	/* fraction of total recombs due to three body - when most are due to this
	 * small changes in temperature can produce large changes in rec coef,
	 * and so in ionization */
	iso.RecomCollisFrac[ipISO][nelem] = (float)(Recom3Body* dense.eden / ionbal.RateRecomTot[nelem][nelem-ipISO] );

	/* get simple estimate of level of ionization */
	if( ionbal.RateRecomTot[nelem][nelem-ipISO] > 0. )
	{
		iso.xIonSimple[ipISO][nelem] = iso.RateLevel2Cont[ipISO][nelem][ipH1s]/ionbal.RateRecomTot[nelem][nelem-ipISO];
	}
	else
	{
		iso.xIonSimple[ipISO][nelem] = 0.;
	}

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

	return;
}
/*lint +e778 constant expression evaluates to zero - in hmrate macro */
