/* 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 */
/*mole_H2_LTE sets Boltzmann factors and LTE unit population of large H2 molecular */
/*H2_init - called to initialize things from cdInit */
/*H2_Zero zero out vars in the large H2 molecule, called from zero 
 * before any commands are parsed */
/* H2_zero_pops_too_low - zero out some H2 variables if we decide not to compute
 * the full sim, called by H2_LevelPops*/
/*H2_Solomon_rate find rates between H2s and H2g and other levels,
 * for eventual use in the chemistry */
/*H2_gs_rates evaluate rates between ground and star states of H2 for use in chemistry */
#include "cddefines.h" 
#include "phycon.h" 
#include "mole.h" 
#include "hmi.h" 
#include "taulines.h" 
#include "h2_priv.h" 
#include "h2.h" 

/*H2_Solomon_rate find rates between H2s and H2g and other levels,
 * for eventual use in the chemistry */
void H2_Solomon_rate( void )
{

	long int iElecLo , iElecHi , iVibLo , iVibHi , iRotLo , iRotHi,
		iplo , ip;

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

	/* eleclo will always be X in this routine */
	iElecLo = 0;

	/* find rate (s-1) h2 dissoc into X continuum by Solomon process and
	 * assign to the TH85 g and s states 
	 * these will go back into the chemistry network */

	/* rates [s-1] for dissociation from s or g, into elec excited states  
	 * followed by dissociation */
	hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.;
	hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0.;

	/* these are used in a print statement - are they needed? */
	hmi.H2_Solomon_elec_decay_H2g = 0.;
	hmi.H2_Solomon_elec_decay_H2s = 0.;

	/* at this point we have already evaluated the sum of the radiative rates out
	 * of the electronic excited states - this is H2_rad_rate_out[elec][vib][rot] 
	 * and this includes both decays into the continuum and bound states of X */

	/* sum over all electronic states, finding dissociation rate */
	for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
	{
		for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
		{
			for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
			{

				/* loop over all lower levels within ground of X to find decay rates from H2g to continuum 
				 * nEner_H2_ground is number of levels in H2g as determined by ENERGY_H2_STAR */
				for( iplo=0; iplo < nEner_H2_ground; ++iplo )
				{
					ip = H2_ipX_ener_sort[iplo];
					iRotLo = ipRot_H2_energy_sort[ip];
					iVibLo = ipVib_H2_energy_sort[ip];
					if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo] )
					{
						/* this is the rate {cm-3 s-1] that mole goes from lower level into elec excited
						 * states then into continuum */
						double rate_up_cont = 
							H2_populations[iElecLo][iVibLo][iRotLo]*
							H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump*
							H2_dissprob[iElecHi][iVibHi][iRotHi] / H2_rad_rate_out[iElecHi][iVibHi][iRotHi];

						/* this is H2g up to exec excit then to continuum - cm-3 s-1 at this point */
						hmi.H2_Solomon_dissoc_rate_BigH2_H2g += rate_up_cont;

						/* rate elec state decays into H2g */
						hmi.H2_Solomon_elec_decay_H2g += 
							H2_populations[iElecHi][iVibHi][iRotHi]*
							H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul*(
							H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pesc+
							H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pdest+
							H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pelec_esc);
					}
				}
				/* loop over levels in H2s */
				for( iplo=nEner_H2_ground; iplo < nLevels_per_elec[0]; ++iplo )
				{
					ip = H2_ipX_ener_sort[iplo];
					iRotLo = ipRot_H2_energy_sort[ip];
					iVibLo = ipVib_H2_energy_sort[ip];
					if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibLo][iRotLo] )
					{
						/* this is the rate {cm-3 s-1] that mole goes from lower level into elec excited
							* states then into continuum */
						double rate_up_cont = 
							H2_populations[iElecLo][iVibLo][iRotLo]*
							H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].pump*
							H2_dissprob[iElecHi][iVibHi][iRotHi] / H2_rad_rate_out[iElecHi][iVibHi][iRotHi];

						/* this is H2s up to exec excit then to continuum - cm-3 s-1 at this point */
						hmi.H2_Solomon_dissoc_rate_BigH2_H2s += rate_up_cont;

						/* rate elec state decays into H2s */
						hmi.H2_Solomon_elec_decay_H2s += 
							H2_populations[iElecHi][iVibHi][iRotHi]*
							H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul*(
							H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pesc+
							H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pdest+
							H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Pelec_esc);
					}
				}
			}/* end iRotHi */
		}/* end iVibHi */
	}/* end iElecHi */
	/* at this point units of hmi.H2_Solomon_elec_decay_H2g, H2s are cm-3 s-1 
	 * since H2_populations are included -
	 * div by pops to get actual dissocation rate, s-1 */
	if( hmi.H2_total > SMALLFLOAT )
	{
		hmi.H2_Solomon_elec_decay_H2g /= SDIV( H2_sum_excit_elec_den );
		hmi.H2_Solomon_elec_decay_H2s /= SDIV( H2_sum_excit_elec_den );

		/* will be used for H2s-> H + H */
		hmi.H2_Solomon_dissoc_rate_BigH2_H2s = hmi.H2_Solomon_dissoc_rate_BigH2_H2s / SDIV(H2_den_s);
		
		/* will be used for H2g-> H + H */
		hmi.H2_Solomon_dissoc_rate_BigH2_H2g = hmi.H2_Solomon_dissoc_rate_BigH2_H2g / SDIV(H2_den_g);
		
	}
	else
	{
		hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0; 
		hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0;	

	}
	/*fprintf(ioQQQ,"DEBUG H2 new %.2e %.2e %.2e %.2e \n",
		hmi.H2_Solomon_elec_decay_H2g ,
		hmi.H2_Solomon_elec_decay_H2s ,
		hmi.H2_Solomon_dissoc_rate_BigH2_H2s,
		hmi.H2_Solomon_dissoc_rate_BigH2_H2g );*/

	/* rate g goes to s */
	hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;


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

	return;
}

/*H2_gs_rates evaluate rates between ground and star states of H2 for use in chemistry */
void H2_gs_rates( void )
{
	long int ipLoX , ip , iRotLoX , iVibLoX ,
		iElecHi , iVibHi , ipOther , iRotOther,
		iRotHi , iVibOther;

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

	/* rate g goes to s */
	hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;

	/* loop over all levels in H2g */
	for( ipLoX=0; ipLoX < nEner_H2_ground; ++ipLoX )
	{
		ip = H2_ipX_ener_sort[ipLoX];
		iRotLoX = ipRot_H2_energy_sort[ip];
		iVibLoX = ipVib_H2_energy_sort[ip];
		/* now find all pumps up to elec excited states */
		/* sum over all electronic states, finding dissociation rate */
		for( iElecHi=1; iElecHi<mole.n_h2_elec_states; ++iElecHi )
		{
			for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
			{
				for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
				{
					if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibLoX][iRotLoX] )
					{
						/* this is the rate {cm-3 s-1] that mole goes from lower level into elec excited
						 * states then into continuum */
						double rate_up_cont = 
							H2_populations[0][iVibLoX][iRotLoX]*
							H2Lines[iElecHi][iVibHi][iRotHi][0][iVibLoX][iRotLoX].pump;

						double decay_star = H2_rad_rate_out[iElecHi][iVibHi][iRotHi] - H2_dissprob[iElecHi][iVibHi][iRotHi];
						/* loop over all other levels in H2g, subtracting their rate - remainder is rate into star 
						 * this is usually only a few levels */
						for( ipOther=0; ipOther < nEner_H2_ground; ++ipOther )
						{
							ip = H2_ipX_ener_sort[ipOther];
							iRotOther = ipRot_H2_energy_sort[ip];
							iVibOther = ipVib_H2_energy_sort[ip];
							if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther] )
							{
								decay_star -= 
									H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Aul*(
									H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Pesc+
									H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Pdest+
									H2Lines[iElecHi][iVibHi][iRotHi][0][iVibOther][iRotOther].Pelec_esc);
							}
						}
						/* MAX because may underflow to negative numbers is rates very large 
						 * this is fraction that returns to H2s */
						decay_star = MAX(0., decay_star )/H2_rad_rate_out[iElecHi][iVibHi][iRotHi];
						hmi.H2_H2g_to_H2s_rate_BigH2 += rate_up_cont*decay_star;

					}/* end if line exists */
				}/* end loop rot elec excit */
			}/* end loop vib elec excit */
		}/* end loop elec elec excit */
	}

	/* at this point units are cm-3 s-1 - convert to rate s-1 */
	hmi.H2_H2g_to_H2s_rate_BigH2 /= SDIV( H2_den_g );

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

	return;
}

/* H2_zero_pops_too_low - zero out some H2 variables if we decide not to compute
 * the full sim, called by H2_LevelPops*/
void H2_zero_pops_too_low( void )
{

	long int iElec, iElecHi, iElecLo, iVib , iVibLo , iVibHi ,
		iRot , iRotHi , iRotLo;

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

	/* >>chng 05 jan 26, add this block to set populations to lte value */
	for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
	{
		for( iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
		{
			for( iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
			{
				/* LTE populations are for unit H2 density, so need to multiply
					* by total H2 density */
				H2_old_populations[iElec][iVib][iRot] = 
					(float)H2_populations_LTE[iElec][iVib][iRot]*hmi.H2_total;
				H2_populations[iElec][iVib][iRot] = H2_old_populations[iElec][iVib][iRot];
			}
		}
	}
	/* zero everything out - loop over all possible lines */
	/* >>chng 05 jan 26, set to lte values, since we still need to accumulate Lyman line
		* optical depths to have correct self-shielding when large h2 does come on */
	for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
	{
		pops_per_elec[iElecHi] = 0.;
		for( iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
		{
			pops_per_vib[iElecHi][iVibHi] = 0.;
			for( iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
			{
				long int lim_elec_lo = 0;
				/* now the lower levels */
				/* NB - iElecLo the lower electronic level is only X 
					* we don't consider excited elec to excit elec trans here, since we are only 
					* concerned with excited electronic levels as a photodissociation process
					* code exists to relax this assumption - simply change following to iElecHi */
				for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
				{
					/* want to include all vib states in lower level if different elec level,
					* but only lower vib levels if same elec level */
					long int nv = nVib_hi[iElecLo];
					if( iElecLo==iElecHi )
						nv = iVibHi;
					for( iVibLo=0; iVibLo<=nv; ++iVibLo )
					{
						long nr = nRot_hi[iElecLo][iVibLo];
						if( iElecLo==iElecHi && iVibHi==iVibLo )
							nr = iRotHi-1;

						for( iRotLo=Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
						{
							if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
							{

								/* population of lower level */
								/* >>chng 05 jan 26, set to lte value */
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo = 
									H2_populations[iElecLo][iVibLo][iRotLo];

								/* population of upper level */
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi = 
									H2_populations[iElecHi][iVibHi][iRotHi];

								/* population of lower level with correction for stim emission */
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopOpc = 
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopLo - 
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].PopHi*
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gLo / 
									H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].gHi;

								/* following two heat exchange excitation, deexcitation */
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].cool = 0.;
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].heat = 0.;

								/* intensity of line */
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].xIntensity = 0.;

								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].phots = 0.;
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ots = 0.;
								H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ColOvTot = 0.;
							}
						}
					}
				}
			}
		}
	}
	hmi.H2_photodissoc_BigH2_H2s = 0.;
	hmi.H2_photodissoc_BigH2_H2g = 0.;
	hmi.HeatH2Dish_BigH2 = 0.;
	hmi.HeatH2Dexc_BigH2 = 0.;
	hmi.deriv_HeatH2Dexc_BigH2 = 0.;
	hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.;
	hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0.;
	hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;

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

	return;
}

/*mole_H2_LTE sets Boltzmann factors and LTE unit population of large H2 molecular */
void mole_H2_LTE( void )
{
	/* used to recall the temperature used for last set of Boltzmann factors */
	static double TeUsedBoltz = -1.;
	double part_fun;
	long int iElec , iVib , iRot;

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

	/* do we need to update the Boltzmann factors and unit LTE populations? */
	/*lint -e777 float test equality */
	if( phycon.te != TeUsedBoltz )
	/*lint +e777 float test equality */
	{
		part_fun = 0.;
		TeUsedBoltz = phycon.te;
		/* loop over all levels setting H2_Boltzmann and deriving partition function */
		for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
		{
			for( iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
			{
				for( iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
				{
					H2_Boltzmann[iElec][iVib][iRot] = 
						/* energy is relative to lowest level in the molecule, v=0, J=0,
						 * so Boltzmann factor is relative to this level */
						sexp( energy_wn[iElec][iVib][iRot] / phycon.te_wn );
					/* sum the partition funciton - Boltzmann factor times statistical weight */
					part_fun += H2_Boltzmann[iElec][iVib][iRot] * H2_stat[iElec][iVib][iRot];
					ASSERT( part_fun > 0 );
				}
			}
		}
		/* have partition function, set H2_populations_LTE (populations for unit H2 density) */
		for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
		{
			for( iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
			{
				for( iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
				{
					/* these are the H2 LTE populations for a unit H2 density -
					 * these populations will sum up to unity */
					H2_populations_LTE[iElec][iVib][iRot] = 
						H2_Boltzmann[iElec][iVib][iRot] * 
						H2_stat[iElec][iVib][iRot] / part_fun;
					/*if( iElec==0 && iVib < 2)
						fprintf(ioQQQ,"DEBUG lte pop\t%i\t%i\t%e\n",
						iVib,iRot,H2_populations_LTE[iElec][iVib][iRot]*hmi.H2_total);*/
				}
			}
		}
		if( mole.lgH2_TRACE >= mole.lgH2_trace_full ) 
			fprintf(ioQQQ,
			"mole_H2_LTE set H2_Boltzmann factors, T=%.2f, partition function is %.2f\n",
			phycon.te,
			part_fun);
	}


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

	return;

}

/*H2_init - called to initialize things from cdInit */
void H2_Init(void)
{

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

	/* the number of electronic quantum states to include.
	 * To do both Lyman and Werner bands want nelec = 3,
	 * default is to do all bands included */
	mole.n_h2_elec_states = N_H2_ELEC;
	h2.nCallH2_this_zone = 0;

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

	return;
}


/*H2_Reset called by IterRestart to reset variables that are needed after an iteration */
void H2_Reset( void )
{

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

	if(mole.lgH2_TRACE) 
		fprintf(ioQQQ,
		"\n***************H2_Reset called, resetting nCallH2_this_iteration, zone %.2f iteration %li\n", 
		fnzone,
		iteration );

	/* number of times large molecules evaluated in this iteration,
	 * is FALSE if never evaluated, on next evaluation will start with lte populations */
	nCallH2_this_iteration = 0;

	/* these remember the largest and smallest factors needed to
	 * renormalize the H2 chemistry */
	h2.renorm_max = 1.;
	h2.renorm_min = 1.;

	/* counters used by H2_itrzn to find number of calls of h2 per zone */
	nH2_pops  = 0;
	nH2_zone = 0;
	/* this is used to establish zone number for evaluation of number of levels in matrix */
	nzone_nlevel_set = 0;

	nzoneAsEval = -1;
	iterationAsEval = -1; 

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

	return;

}

/*H2_Zero zero out vars in the large H2 molecule, called from zero 
 * before any commands are parsed 
 * NB - this routine is called before space allocated - must not zero out
 * any allocated arrays */
void H2_Zero( void )
{

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

	/* this is the smallest ratio of H2/H where we will bother with the large H2 molecule
	 * this value was chosen so that large mole used at very start of TH85 standard pdr,
	 * NB - this appears in headinfo and must be updated there if changed here */
	/* >>chng 03 jun 02, from 1e-6 to 1e-8 - in orion veil large H2 turned on half way
	 * across, and solomon process was very fast since all lines optically thin.  correct
	 * result has some shielding, so reset to lower value so that H2 comes on sooner. */
	mole.H2_to_H_limit = 1e-8;

	/* these are used to set trace levels of output by options on atom h2 trace command 
	 * first is minimum level of trace, keyword is FINAL */
	mole.lgH2_trace_final = 1;
	/* intermediate level of trace, info per iteration, key ITERATION */
	mole.lgH2_trace_iterations = 2;
	/* full trace, keyword is FULL */
	mole.lgH2_trace_full = 3;
	/* print matrices used in solving X */
	mole.lgH2_trace_matrix = 4;

	h2.lgH2ON = FALSE;
	mole.lgH2_TRACE = FALSE;
	/* flag to force using LTE level populations */
	mole.lgH2_LTE = FALSE;

	/* this is number of electronic levels to include in the print and punch output 
	 * changed with the PRINT LINE H2 ELECTRONIC and PUNCH H2 commands 
	 * by default only include X */
	h2.nElecLevelOutput = 1;

	/* counters used by H2_itrzn to find number of calls of h2 per zone */
	nH2_pops  = 0;
	nH2_zone = 0;
	/* this is used to establish zone number for evaluation of number of levels in matrix */
	nzone_nlevel_set = 0;

	/* option to scramble collision data */
	mole.lgH2_NOISE = FALSE;
	mole.lgH2_NOISECOSMIC = FALSE; 

	/* option to turn off or on gbar collisions of the collision rate,
	 * default is to have it on */
	/* turn mole.lgColl_gbar on/off with atom h2 gbar on off */
	mole.lgColl_gbar = TRUE;

	/* include collision rates that come from real calculations,
	 * off with atom h2 collisions off command */
	mole.lgColl_deexec_Calc = TRUE;

	mole.lgColl_dissoc_coll = TRUE;

	/* option to turn off H2 - grain collision & deexcitation,
	 * atom h2 grain collision on/off */
	mole.lgH2_grain_deexcitation = FALSE;

	/* option to turn off ortho-para collisions, command ATOM H2 COLLISIONS ORTHO PARA OFF */
	mole.lgH2_ortho_para_coll_on = TRUE;

	/* which set of H2 - He collision data to use?  default is Le Bourlot Meudon set,
	 * set to Stnacil with command
	 * atom H2 He collisions new */
	mole.lgH2_He_Meudon = TRUE;
	mole.lgH2_He_Stancil = FALSE;

	/* these remember the largest and smallest factors needed to
	 * renormalize the H2 chemistry */
	h2.renorm_max = 1.;
	h2.renorm_min = 1.;

	nCallH2_this_iteration = 0;
	h2.ortho_density = 0.;
	h2.para_density = 0.;

	hmi.H2_Solomon_dissoc_rate_BigH2_H2s = 0.;
	hmi.H2_Solomon_dissoc_rate_BigH2_H2g = 0.;
	hmi.H2_H2g_to_H2s_rate_BigH2 = 0.;
	hmi.H2_photodissoc_BigH2_H2s = 0.;
	hmi.H2_photodissoc_BigH2_H2g = 0.;

	/* say that H2 has never been computed */
	hmi.lgBigH2_evaluated = FALSE;

	hmi.lgH2_Thermal_BigH2 = TRUE;
	hmi.lgH2_Chemistry_BigH2 = TRUE;

	if( !lgH2_READ_DATA )
	{
		/* the number of electronic levels in the H2 molecule,
		 * to just do the Lyman and Werner bands set to 3 -
		 * reset with atom h2 levels command,
		 * default is all levels with data */
		mole.n_h2_elec_states = N_H2_ELEC;
	}

	/* the number of levels used in the matrix solution
	 * of the level H2_populations - set with atom h2 matrix nlevel,
	 * >>chng 04 oct 05, make default 30 levels 
	 * >>chng 04 dec 23, make default 70 levels */
	nXLevelsMatrix = 70;

	/* this is used to establish zone number for evaluation of number of levels in matrix */
	nzone_nlevel_set = -1;

	nzoneAsEval = -1;
	iterationAsEval = -1; 

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

	return;
}

#define SizeOf_line 300
#define Num_fit_par 8   
#define nu 302
#define nl 302   
static double fit_par[nu][nl][Num_fit_par];   
static int defn[nu][nl];   

/* This code was written by Terry Yun, 2005 */ 
/* H2_He_coll_init receives the name of the file that contrains the fitting     
 * coefficeints of all transitions and read into 3d vectors. 
 * It retuns magic number*/
long int H2_He_coll_init(char FILE_NAME_IN[] )
{

	int i, j;   
	long int magic;
	double par[Num_fit_par];
	char line[SizeOf_line];
	int h2_i, h2_f, he_i, he_f;/* target, projectil initial and final indices */
	char quality, space;
	/*'space' varaible is for reading spaces between chracters */
	/* scanf can skip spaces(or tap) when it reads numbers, but not for chracters. */
	double error;

	FILE *ifp;

#	if 0
	FILE *ioQQQ;
	ioQQQ = fopen(FILE_NAME_OUT, "w");
	if( ioQQQ == NULL )
	{
		printf("Can't open %s\n", FILE_NAME_OUT);
		exit(1);
	}
#	endif

#if 0
	/*open the input file to read and count the number of lines*/
	ifp = fopen(FILE_NAME_IN, "r");
	if( ifp == NULL )
	{
		printf("Can't open %s\n", FILE_NAME_IN);
		exit(1);
	}
	while( fgets(line, (int)sizeof(line), ifp) != NULL )
	{
		Num_lines++;
	}
	printf("The number of lines are %i.\n", Num_lines);

	fclose(ifp);
#endif

	/* set a flag: initially everything is '0' */
	for( i=0; i<nu; i++ )
	{
		for( j=0; j<nl; j++ )
		{
			defn[i][j] = 0;
		}
	}

	/*open the input file and put into 3d arrays*/
	ifp = fopen(FILE_NAME_IN, "r");
	if( ifp == NULL )
	{
		printf("Can't open %s\n", FILE_NAME_IN);
		exit(1);
	}

	/*read the magic number*/
	fgets(line, (int)sizeof(line), ifp);
	sscanf(line, "%li", &magic); 

	/* read the file until the end of line */
	while( fgets(line, (int)sizeof(line), ifp) != NULL )
	{
		/* skip any line that starts with '#' */
		if( line[0]!='#' )
		{
			sscanf(line, "%i%i%i%i%c%c%c%c%lf%lf%lf%lf%lf%lf%lf%lf%lf", &h2_i, &h2_f, &he_i, &he_f,&space, &space, &space, &quality, &error, &par[0], &par[1], &par[2], &par[3], &par[4], &par[5], &par[6], &par[7]); 

			/* Set a flag to '1' when the indices are defined */
			defn[h2_i][h2_f] = 1;
			for( i=0; i<Num_fit_par; i++ )
				/* assigning the parameters to 3d array */
				fit_par[h2_i][h2_f][i] = par[i];
		}
	}
	fclose(ifp);  
	/* output fitting coeffs for all defined transitions for testing */
#	if 0
	for( i=0; i<nu; i++ )
	{
		for( j=0; j<nl; j++ )
		{
			if( defn[i][j] == 1 )
			{
				fprintf(ioQQQ, "%i %i %e %e %e %e %e %e %e %e\n",  i, j, fit_par[i][j][0], fit_par[i][j][1], fit_par[i][j][2], fit_par[i][j][3], fit_par[i][j][4], fit_par[i][j][5], fit_par[i][j][6], fit_par[i][j][7]);
			}
		}
	}

	fclose(ofp);
#	endif
	return magic;
}

/* This code was written by Terry Yun, 2005 */ 
/* H2_He_coll Interpolate the rate coefficeints 
 * It receives initial and final transition and temperature
 * The range of the temperature is between 2K - 1e8K */
double H2_He_coll(int init, int final, double temp)
{

	double k, t, b2, c2, f2, t2;
	/*return 0;*/
	/* Invalid entries returns '-1':the initial indices are smaller than the final indices */
	if( temp<2 || temp > 1e8 )
	{
		k = -1;
	}
	else if( init <= final )
	{
		k = -1;
	}
	/* Invalid returns '-1': the indices are greater than 302 or smaller than 0 */
	else if( init < 0 || init >302 || final < 0 || final > 302 )
	{
		k = -1;
	}
	/* Undefined indices returns '0' */
	else if( defn[init][final] == 0 )
	{
		k = -1;
	}
	/* defined indices */
	else if( defn[init][final] == 1 )
	{
		double sum1 , sum2;
		/* t2 = T*10^-3 */
		t2 = temp/1e3; 
		/* t = T*10^-3 + 1*/
		t = t2+1; 
		b2 = fit_par[init][final][1]/(fit_par[init][final][3]*t2+1);
		c2 = fit_par[init][final][2]/(t*t);
		{
			/*@-redef@*/
			enum {DEBUG_LOC=FALSE};
			/*@+redef@*/
			if( DEBUG_LOC )
			{
				fprintf(ioQQQ,"bug H2 He coll\t%i %i %.3e %.3e %.3e \n",
					init,final,
					fit_par[init][final][5],
					fit_par[init][final][6]*t2+1, 
					fit_par[init][final][7]
					/*pow(fit_par[init][final][6]*t2+1,fit_par[init][final][7])*/
					);
			}
		}
		sum1 =  fit_par[init][final][7] * log10( fit_par[init][final][6]*t2+1 );
		if( fabs(sum1)< 38. )
			f2 = fit_par[init][final][5]/pow(fit_par[init][final][6]*t2+1,fit_par[init][final][7]);
		else
			f2= 0.;
		{
			/*@-redef@*/
			enum {DEBUG_LOC=FALSE };
			/*@+redef@*/
			if( DEBUG_LOC )
			{
				fprintf(ioQQQ,"bug H2 He coll\t%i %i %.3e %.3e %.3e %.3e %.3e sum %.3e %.3e \n",
					init,final,
					fit_par[init][final][0],
					b2, 
					c2,
					fit_par[init][final][4],
					f2  ,
					fit_par[init][final][0]+b2+c2,
					fit_par[init][final][4]+f2);
			}
		}

		/*k = pow(10.,fit_par[init][final][0]+b2+c2) + pow(10.,fit_par[init][final][4]+f2);*/
		sum1 = fit_par[init][final][0]+b2+c2;
		sum2 = fit_par[init][final][4]+f2;
		k = 0.;
		if( fabs(sum1) < 38. )
			k += pow(10., sum1 );
		if( fabs(sum2) < 38. )
			k += pow(10., sum2 );
		{
			/*@-redef@*/
			enum {DEBUG_LOC=FALSE};
			/*@+redef@*/
			if( DEBUG_LOC )
			{
				fprintf(ioQQQ,"bug H2 He rate %.3e \n",
					k);
			}
		}
	}
	/*unknown invalid entries returns '99'*/
	else
	{
		k = 99;
	}

	return k;
}

