/* 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 */
/*CO_gasgrain defines matrix elements for gas-grain chemistry of heavy element molecules */
#include "cddefines.h"
#include "physconst.h"
#include "dense.h"
#include "secondaries.h"
#include "grainvar.h"
#include "phycon.h"
#include "mole.h"
/*lint -e778 constant expression evaluatess to 0 in operation '-' */
/*=================================================================*/


/*CO_step fills in matrix for heavy elements molecular routines */
void CO_gasgrain(void)
{

	/* this routine works with the heavy molecules, of which there are NUM_HEAVY_MOLEC,
	 * plus C and O, making a total of NUM_COMOLE_CALC*/

	
	double number_of_sites, 
			binding_energy_CO, 
			binding_energy_H2O, 
			binding_energy_OH, 
			exponent_CO, 
			exponent_H2O, 
			exponent_OH, 
			vib_freq_CO, 
			vib_freq_H2O, 
			vib_freq_OH, 
			kevap_CO, 
			kevap_H2O, 
			kevap_OH, 
			den_times_area,
			rate;
	int nd;

	/******************************** Gas-Grain Chemistry**********************************/

	/*  The Gas-Grain Chemistry rates are taken from: 
	>>refer	Hasegawa, T. I. & Herbst, E. 1993, MNRAS, 261, 83

	So far only CO depletion onto grains is considered, however, this code can be generalized 
	if desired to other molecules, using the data in the above paper.  There are three important reactions 
	to determine the abundance of a molecule on grain surfaces deep in molecular clouds.  The
	rate of accretion of molecule X onto grains is

	R(accretion) = PI*[grain_radius]^2*[thermal velocity]*[density_of_grains]

	Two processes remove molecule X from the grain surface, and that is thermal evaporation, due
	to the heat of the grain, and cosmic ray deabsorption.  The first of these rates come from the 
	above paper, and depends primarily on the dust temperature.  The cosmic ray rate is a constant,
	calculated in Hasegawa and Herbst.  

	For each molecule desired, I have created a new species which is the density of that molecule
	on the grain surface */

	
	/* evaporation rate of molecule on grain is:

	k(evap) = [vibrational absorption frequency]*exp[-binding_energy/dust_temperature]

	The binding energies come from Hasegawa and Herbst, Table 4.  The vibrational frequency comes from
	equation 3 of >>refer Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJSS, 82, 167

	[vibrational absorption frequency] = 
	SQRT[[2*number_of_sites_for_grain_absorption*binding_energy]/[PI^2*mass_of_molecule]]

	**********************************************************************************************/

	/* number of sites on grain */
	number_of_sites = 1.5e15f;
	/* binding energy of CO */
	binding_energy_CO = 1210.0f;
	binding_energy_H2O = 1860.0f;
	binding_energy_OH = 1260.0f;
	/* exponential factor in evaporation rate */
	exponent_CO = 0.0f;
	exponent_H2O = 0.0f;
	exponent_OH = 0.0f;
	/* grain density times grain area */
	den_times_area = 0.0f;

	/* calculate the rates that are dependent on grain physics.  This includes grain density, 
	cross sectional area, and dust temperature of each constituent.  Note that 
	
	gv.bin[nd]->IntVol/gv.bin[nd]->AvVol*gv.bin[nd]->dstfactor*dense.gas_phase[ipHYDROGEN]
	
	is the grain density, according to abundances.c */

	for( nd=0; nd < gv.nBin; nd++ )
	{
		exponent_CO += (double)exp(-binding_energy_CO/gv.bin[nd]->tedust);
		exponent_H2O += (double)exp(-binding_energy_H2O/gv.bin[nd]->tedust);
		exponent_OH += (double)exp(-binding_energy_OH/gv.bin[nd]->tedust);
		den_times_area += (double)((gv.bin[nd]->IntVol/gv.bin[nd]->AvVol)*
			pow(gv.bin[nd]->AvRadius,2)*3.14159f*gv.bin[nd]->dstfactor*
			dense.gas_phase[ipHYDROGEN]);


	}


	vib_freq_CO = (float)(pow(2*number_of_sites*binding_energy_CO/(3.14159f*3.14159f*28*ATOMIC_MASS_UNIT),0.5));
	vib_freq_H2O = (float)(pow(2*number_of_sites*binding_energy_H2O/(3.14159f*3.14159f*18*ATOMIC_MASS_UNIT),0.5));
	vib_freq_OH = (float)(pow(2*number_of_sites*binding_energy_OH/(3.14159f*3.14159f*17*ATOMIC_MASS_UNIT),0.5));
	kevap_CO = vib_freq_CO*exponent_CO;
	kevap_H2O = vib_freq_H2O*exponent_H2O;
	kevap_OH = vib_freq_OH*exponent_OH;


	/*	CO_grain	evaporation =>	CO	+	grain	*/

	mole.c[ipCOGR][ipCOGR] -= kevap_CO;
	mole.c[ipCOGR][ipCO] += kevap_CO;

	/*	CO_grain	+	cosmic_ray	=>	CO	+	grain	*/

	/* >>chng  05 nov 21 NPA.  The cosmic ray desorption rate should be scaled
	 *  by the cosmic ray rate used by the code, as in mole_co_photo.c.  The scale
	 *  factor here assumes a standard cosmic ray ionization rate of 5e-17 */

	mole.c[ipCOGR][ipCOGR] -= 196.*secondaries.csupra[ipHYDROGEN][0];
	mole.c[ipCOGR][ipCO] += 196.*secondaries.csupra[ipHYDROGEN][0];

	/*	CO	+	grain	=>	CO_grain	*/

	mole.c[ipCO][ipCO] -= den_times_area*pow((BOLTZMANN*phycon.te/(28*ATOMIC_MASS_UNIT)),0.5);
	mole.c[ipCO][ipCOGR] += den_times_area*pow((BOLTZMANN*phycon.te/(28*ATOMIC_MASS_UNIT)),0.5);

	/*	H2O_grain	evaporation =>	H2O	+	grain	*/

	mole.c[ipH2OGR][ipH2OGR] -= kevap_H2O;
	mole.c[ipH2OGR][ipH2O] += kevap_H2O;

	/*	H2O_grain	+	cosmic_ray	=>	H2O	+	grain	*/

	mole.c[ipH2OGR][ipH2OGR] -= 0.028*secondaries.csupra[ipHYDROGEN][0];
	mole.c[ipH2OGR][ipH2O] += 0.028*secondaries.csupra[ipHYDROGEN][0];

	/*	H2O	+	grain	=>	H2O_grain	*/

	mole.c[ipH2O][ipH2O] -= den_times_area*pow((BOLTZMANN*phycon.te/(18*ATOMIC_MASS_UNIT)),0.5);
	mole.c[ipH2O][ipH2OGR] += den_times_area*pow((BOLTZMANN*phycon.te/(18*ATOMIC_MASS_UNIT)),0.5);

	/*	OH_grain	evaporation =>	OH	+	grain	*/

	mole.c[ipOHGR][ipOHGR] -= kevap_OH;
	mole.c[ipOHGR][ipOH] += kevap_OH;

	/*	OH_grain	+	cosmic_ray	=>	OH	+	grain	*/

	mole.c[ipOHGR][ipOHGR] -= 126.*secondaries.csupra[ipHYDROGEN][0];
	mole.c[ipOHGR][ipOH] += 126.*secondaries.csupra[ipHYDROGEN][0];

	/*	OH	+	grain	=>	OH_grain	*/

	mole.c[ipOH][ipOH] -= den_times_area*pow((BOLTZMANN*phycon.te/(17*ATOMIC_MASS_UNIT)),0.5);
	mole.c[ipOH][ipOHGR] += den_times_area*pow((BOLTZMANN*phycon.te/(17*ATOMIC_MASS_UNIT)),0.5);

	/*>>chng 06 jan 11, NPA - In some H+ regions the grain temperatures are so low
	that molecular freeze out occurs.  This should not happen, because the ices
	should sublimate in such a harsh environment.  Therefore, we introduce an
	artificial sublimation rate to destroy ices.  THIS IS NOT A PHYSICAL RATE!!!!
	only a rate that gets the desired, realistic result */
	/*TODO	0	find physical theory for this process */
	/* Rate comes from Table curve and assumes that rate is high (~1) in H+
	 * region and negligible ( < 1e-30) in molecular cloud - designed to stop
	 * freeze-out above 100 K */
	rate = sexp( 555.89/phycon.sqrte - 5.55 );

	/* CO_grain + Hot_Ionized_gas => CO + grain */
	mole.c[ipCOGR][ipCOGR] -= rate;
	mole.c[ipCOGR][ipCO] += rate;

	/* OH_grain + Hot_Ionized_gas => COH + grain */
	mole.c[ipOHGR][ipOHGR] -= rate;
	mole.c[ipOHGR][ipOH] += rate;

	/* H2O_grain + Hot_Ionized_gas => CH2O + grain */
	mole.c[ipH2OGR][ipH2OGR] -= rate;
	mole.c[ipH2OGR][ipH2O] += rate;

	if(!mole.lgGrain_mole_deplete || !gv.lgDustOn )
	{
		mole.c[ipCOGR][ipCOGR] = -1;
		mole.c[ipCOGR][ipCO] = 1;

		mole.c[ipH2OGR][ipH2OGR] = -1;
		mole.c[ipH2OGR][ipH2O] = 1;

		mole.c[ipOHGR][ipOHGR] = -1;
		mole.c[ipOHGR][ipOH] = 1;

	}
}

