/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
#include "cddefines.h"
#include "physconst.h"
#include "h2.h"
#include "h2_priv.h"
#include "hmi.h"
#include "rfield.h"
#include "thirdparty.h"
#include "ipoint.h"
#include "mole.h"
	
STATIC void GetDissociationRateCoeff( diss_tran& tran );

void diatomics::Read_Mol_Diss_cross_sections()
{
	DEBUG_ENTRY( "diatomics::Read_Mol_Diss_cross_sections()" );

	FILE *ioDATA;

	char chPath[FILENAME_PATH_LENGTH_2],
		chDirectory[FILENAME_PATH_LENGTH_2],
		chLine[FILENAME_PATH_LENGTH_2];

	/* these can be generalized later for other molecules */
	const int ipNUM_FILES = 1;

	char chFileNames[ipNUM_FILES][17] =
		{"cont_diss.dat"} ;

	// these cross sections are not meaningful without big H2 enabled. 
	ASSERT( lgEnabled );

	// for the summed-over-final-state rates, allocate space for all states in model. 
	// Will be zero-filled and overwritten with whatever data is available.
	Cont_Dissoc_Rate.reserve( n_elec_states );
	for( int iElecLo=0; iElecLo<n_elec_states; ++iElecLo )
	{
		Cont_Dissoc_Rate.reserve( iElecLo, nVib_hi[iElecLo]+1 );
		for( int iVibLo=0; iVibLo<=nVib_hi[iElecLo]; ++iVibLo )
		{
			Cont_Dissoc_Rate.reserve( iElecLo, iVibLo, nRot_hi[iElecLo][iVibLo]+1 );
		}
	}
	Cont_Dissoc_Rate.alloc();

	/* drill down to the directory where H2_cont_diss.dat lives */
#	ifdef _WIN32
	strcpy( chDirectory, "h2\\" );
#	else
	strcpy( chDirectory, "h2/" );
#	endif

	/* now open the data file */
	strcpy( chPath, chDirectory );
	strcat( chPath, chFileNames[0] );
	ioDATA = open_data( chPath, "r" );

	Diss_Trans.clear();

	/* track number of datasets in file */
	while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
	{
		static bool skipData = false;
		long ini=0, inf, iv, ij;

		/* ignore all lines that begin with a '#', unless they are telling us the next quantum numbers.  If so,
		 * then read the quantum numbers and start a new dataset */
		if(sizeof(chLine) >= 2 && chLine[0] == '#' && chLine[1] == '!' && chLine[4] == 'n'  && chLine[5] =='e' )
		{
			sscanf(chLine,"#!  nei=%li, nef=%li, vi= %li, ji= %li",&ini, &inf, &iv, &ij);
			// skip data that is for levels beyond our model
			if( ini > n_elec_states )
				skipData = true;
			else if( iv > nVib_hi[ini] )
				skipData = true;
			else if( ij > nRot_hi[ini][iv] )
				skipData = true;
			else
			{
				skipData = false;
				diss_level Leveli = {ini, iv, ij};
				diss_level Levelf = {inf, -1, -1};
				diss_tran thisTran( Leveli, Levelf );
				Diss_Trans.push_back( thisTran );
			}
		}
		
		if( skipData )
			continue;
	
		/* if line does not begin with a '#', then they are numbers in the datset.  Read in energy in 
		 * wavenumbers and cross section in Angstroms ^2 */
		if(chLine[0] != '#')
		{
			realnum energy, xsection;
			const realnum AngSquared = 1e-16f;
			sscanf(chLine,"%f,%f",&energy, &xsection);

			/* convert energy in wavenumbers to Ryd */
			Diss_Trans.back().energies.push_back( energy*WAVNRYD );

			/* convert cross section from Angstrom^2 to cm^2 by multiplying by 1e-16 */
			Diss_Trans.back().xsections.push_back( xsection*AngSquared );
		}
	}

	fclose(ioDATA);
	return;
}

double diatomics::MolDissocOpacity( const diss_tran& tran, const double& Mol_Ene )
{
	DEBUG_ENTRY( "diatomics::MolDissocOpacity()" );

	double cross_section = MolDissocCrossSection( tran, Mol_Ene ) * 
		states[ ipEnergySort[ tran.initial.n ][ tran.initial.v ][ tran.initial.j ] ].Pop();
	return cross_section;
}

double MolDissocCrossSection( const diss_tran& tran,  const double& Mol_Ene )
{
	DEBUG_ENTRY( "MolDissocCrossSection()" );

	double photodiss_cs = 0.;

	/* if energy is less than threshold, then set cross section is zero */
	if (Mol_Ene < tran.energies[0])
		photodiss_cs = 0.;
	/* if energy is greater than the highest energy defined in Philip's data, set 
	 * cross section to zero for now.  This may need to be changed later */
	else if(Mol_Ene > tran.energies.back())
		photodiss_cs = tran.xsections.back()/sqrt(powi(Mol_Ene/tran.energies.back(),7));
	/* If energy is in between high and low energy limits, do linear interpolation
	 * on Philip's data to compute cross section */
	else
	{
		ASSERT( Mol_Ene > tran.energies[0] && Mol_Ene < tran.energies.back() );
		photodiss_cs = linint(&tran.energies[0],
			&tran.xsections[0],
			tran.xsections.size(),
			Mol_Ene);
	}

	return photodiss_cs;
}

void diatomics::Mol_Photo_Diss_Rates( void )
{
	DEBUG_ENTRY( "diatomics::Mol_Photo_Diss_Rates()" );

	/* Start Calculation of H2 continuum dissociation rate here */

	// these cross sections are not meaningful without big H2 enabled. 
	ASSERT( lgEnabled && mole_global.lgStancil );

	/* Zero out dissociation rates */
	Cont_Dissoc_Rate.zero();
	Cont_Dissoc_Rate_H2g = 0.;
	Cont_Dissoc_Rate_H2s = 0.;
	
	for_each( Diss_Trans.begin(), Diss_Trans.end(), GetDissociationRateCoeff );
	for( vector< diss_tran >::const_iterator dt = Diss_Trans.begin(); dt != Diss_Trans.end(); ++dt )
	{
		double rate = (*this).GetDissociationRate( *dt );
		// this is summed over final electron state, tran.rate_coeff is not (because it's used to conserve energy)	
		Cont_Dissoc_Rate[dt->initial.n][dt->initial.v][dt->initial.j] += dt->rate_coeff;
		if( states[ ipEnergySort[dt->initial.n][dt->initial.v][dt->initial.j] ].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
		{
			/* Add up total rate due to levels which are above arbitrarily defined 
			 * H2* threshold.  Call this dissociation rate of H2* */
			Cont_Dissoc_Rate_H2s += rate;
		}
		else
		{
			/* Add up total rate due to levels which are below arbitrarily defined 
			 * H2* threshold.  Call this dissociation rate of H2g (H2 ground) */
			Cont_Dissoc_Rate_H2g += rate;
		}
	}

	/* This has current units of cm-3 s-1.  We want just the total rate coefficient for
	 * insertion into h_step.  Therefore, divide rate by density of H2* and H2g to get rate 
	 * coefficient for reaction */
	Cont_Dissoc_Rate_H2g = Cont_Dissoc_Rate_H2g / SDIV(H2_den_g);
	Cont_Dissoc_Rate_H2s = Cont_Dissoc_Rate_H2s / SDIV(H2_den_s);

	return;
}

STATIC void GetDissociationRateCoeff( diss_tran& tran )
{
	DEBUG_ENTRY( "GetDissociationRateCoeff()" );
	
	long index_min = ipoint( tran.energies[0] );
	long index_max = ipoint( tran.energies.back() ); 
	index_max = MIN2( index_max, rfield.nflux-1 );

	tran.rate_coeff = 0.;
	for(long i = index_min; i <= index_max; ++i)
	{
		/* Find cross section at given n, v, j, and photon energy */
		double photodiss_cs = MolDissocCrossSection( tran, rfield.anu[i] );
				
		/* Compute integral of cross section and photon energy -- rate coefficient */
		tran.rate_coeff += photodiss_cs*( rfield.flux[0][i] + rfield.ConInterOut[i]+
			rfield.outlin[0][i]+ rfield.outlin_noplot[i]  );
	}
}

double diatomics::GetDissociationRate( const diss_tran& tran )
{
	DEBUG_ENTRY( "diatomics::GetDissociationRate()" );

	long iElecLo = tran.initial.n;
	long iVibLo = tran.initial.v;
	long iRotLo = tran.initial.j;
					
	/* Compute rate, product of rate coefficient times density */
	return states[ ipEnergySort[iElecLo][iVibLo][iRotLo] ].Pop() * tran.rate_coeff;
}

double diatomics::Cont_Diss_Heat_Rate( void )
{
	DEBUG_ENTRY( "diatomics::Cont_Diss_Heat_Rate()" );

	/* 29 Jul 10 - NPA - this tells code to not compute detailed H2 dissociation heating rate
	 * if Stancil rate is not used.  However, we also do not want to compute this rate if Stancil 
	 * rate is on, but small H2 molecule is used.  Insert logic hear to not use rate if 
	 * small H2 atom is on */

	if( !mole_global.lgStancil || !lgEnabled )
		return 0.;
	
	Mol_Photo_Diss_Rates();

	double Cont_Dissoc_Heat_Rate = 0.0;
	for( vector< diss_tran >::const_iterator dt = Diss_Trans.begin(); dt != Diss_Trans.end(); ++dt )
		Cont_Dissoc_Heat_Rate += (*this).GetHeatRate( *dt );

	return Cont_Dissoc_Heat_Rate;
}

double diatomics::GetHeatRate( const diss_tran& tran )
{
	DEBUG_ENTRY( "diatomics::GetHeatRate()" );
					
	long index_min = ipoint( tran.energies[0] );
	long index_max = ipoint( tran.energies.back() ); 
	index_max = MIN2( index_max, rfield.nflux-1 );
 
	long iElecLo = tran.initial.n;
	long iVibLo = tran.initial.v;
	long iRotLo = tran.initial.j;
	double rate = 0.;

	for( long i = index_min; i<= index_max; ++i )
	{
		/* The energy of H2(v,J) above ground */
		double E_H2 = states[ ipEnergySort[iElecLo][iVibLo][iRotLo] ].energy().Ryd();

		/* The energy of the reactants, E[H2(v,J)] + E(h(nu)) */
		double E_Reactants = rfield.anu[i] + E_H2;

		/* The energy of the products  E[H(1s)] + E[H(2p)]
		 * In Ryd, this is just 0.75 Ryd */
		double E_Products = 0.75;

		/* The density of the current H2(v,J) level */
		double density = states[ ipEnergySort[iElecLo][iVibLo][iRotLo] ].Pop();

		/* The dissociation rate coefficient of the H2(v,J) level at the current energy
		 * computed in Mol_Photo_Diss_Rates().  This is the integral of the cross section
		 * times the energy over which the cross section is defined */
		double Rate_Coeff = tran.rate_coeff;

		/* The heating rate.  This is the product of the H2(v,J) density, the dissociation rate,
		 * and the energy difference between the reactants and products.  Any excess energy goes
		 * into heating the gas.  Has the units of erg cm-3 s-1 by multiplying rate by the converstion
		 * factor going from erg to Ryd.  Sum over all possible H2 levels and dissociation energies. */
		rate += EN1RYD*(E_Reactants - E_Products) * Rate_Coeff * density;
	}

	return rate;
}