/* 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 */
/*CO_step fills in matrix for heavy elements molecular routines */
#include "cddefines.h"
#include "mole.h"
#include "mole_priv.h"
#include "save.h"
#include "dense.h"
#include "atmdat.h"
/* Nick Abel between July and October of 2003 assisted Dr. Ferland in improving the heavy element 
 * molecular network in Cloudy. Before this routine would predict negative abundances if 
 * the fraction of carbon in the form of molecules came close to 100%. A reorganizing of 
 * the reaction network detected several bugs.  Treatment of "coupled reactions",
 * in which both densities in the reaction rate were being predicted by Cloudy, were also 
 * added.  Due to these improvements, Cloudy can now perform calculations
 * where 100% of the carbon is in the form of CO without predicting negative abundances
 *
 * Additional changes were made in November of 2003 so that our reaction 
 * network would include all reactions from the TH85 paper.  This involved 
 * adding silicon to the chemical network.  Also the reaction rates were
 * labeled to make identification with the reaction easier and the matrix 
 * elements of atomic C, O, and Si are now done in a loop, which makes 
 * the addition of future chemical species (like N or S) easy.
 * */
/* Robin Williams in August 2006 onwards reorganized the coding to cut down repetitions.  
 * This isolated several further bugs, and allows a sigificant number of lines of
 * code to be eliminated.  The balance of S2/S2+ amd ClO/ClO+ seems highly sensitive
 * (with small log scale results varying significantly if the order of arithmetic
 * operations is changed) -- I suspect this may imply a bug somewhere.
 * */
/*lint -e778 constant expression evaluatess to 0 in operation '-' */
/*=================================================================*/

STATIC bool lgNucleiConserved(const multi_arr<double,2> &c);

void mole_eval_balance(long int num_total, double *b, bool lgJac, multi_arr<double,2> &c)
{
	long int i, j;
	mole_reaction *rate;
	double rate_tot, rate_deriv[MAXREACTANTS], rk;
	molecule *sp;

	DEBUG_ENTRY("mole_eval_balance()");
	/* zero out array used for formation rates */
	for( i=0; i < num_total; i++ )
	{
		b[i] = 0.;
	}
	if (lgJac)
		c.zero();

	for(mole_reaction_i p=mole_priv::reactab.begin(); 
			p != mole_priv::reactab.end(); ++p) 
	{
		rate = &(*p->second);
		rk = mole.reaction_rks[ rate->index ];

		rate_tot = rk;
		for(i=0;i<rate->nreactants;i++)
		{
			rate_tot *= mole.species[ rate->reactants[i]->index ].den;
		}		
		
		for(i=0;i<rate->nreactants;i++)
		{	
			sp = rate->reactants[i];
			if (rate->rvector[i] == NULL)
			{
				b[sp->index] -= rate_tot;
			}
		}
		
		for(i=0;i<rate->nproducts;i++)
		{
			sp = rate->products[i];
			if (rate->pvector[i] == NULL)
			{
				b[sp->index] += rate_tot;
			}
		}
		
		if (lgJac)
		{
			for(i=0;i<rate->nreactants;i++)
			{
				rate_deriv[i] = rk;
				for(j=0;j<rate->nreactants;j++)
				{
					if(i!=j)
					{
						rate_deriv[i] *= mole.species[ rate->reactants[j]->index ].den;
					}
				}
			}
			for(j=0;j<rate->nreactants;j++)
			{
				sp = rate->reactants[j];
				const double rated = rate_deriv[j];
				for(i=0;i<rate->nreactants;i++)
				{
					if (rate->rvector[i] == NULL)
						c[sp->index][rate->reactants[i]->index] -= rated;
				}
				for(i=0;i<rate->nproducts;i++)
				{
					if (rate->pvector[i] == NULL)
						c[sp->index][rate->products[i]->index] += rated;
				}
			}
		}
	}

	if (lgJac)
	{
		ASSERT( lgNucleiConserved(c) );
	}

	//mole_dominant_rates(findspecies("H+"),ioQQQ);
	return;
}

void mole_eval_sources(long int num_total)
{
	long int i, j, nelem, ion, ion2;
	mole_reaction *rate;
	double rate_tot, rate_deriv[MAXREACTANTS], rk;
	molecule *sp;

	DEBUG_ENTRY("mole_eval_sources()");
	/* zero out array used for formation rates */
	for( i=0; i < num_total; i++ )
	{
		mole.species[i].src = mole.species[i].snk = 0.;
	}

	for( nelem=0; nelem< LIMELM; ++nelem )
	{
		/* these have one more ion than above */
		for( ion=0; ion<nelem+2; ++ion )
		{
			/* zero out the transfer array */
			for( ion2=0; ion2<nelem+2; ++ion2 )
			{
				mole.xMoleChTrRate[nelem][ion][ion2] = 0.;
			}
		}
	}

	for(mole_reaction_i p=mole_priv::reactab.begin(); 
			p != mole_priv::reactab.end(); ++p) 
	{
		rate = &(*p->second);
		rk = mole.reaction_rks[ rate->index ];
		
		for(i=0;i<rate->nreactants;i++)
		{
			rate_deriv[i] = rk;
			for(j=0;j<rate->nreactants;j++)
			{
				if(i!=j)
				{
					rate_deriv[i] *= mole.species[ rate->reactants[j]->index ].den;
				}
			}
		}
		
		rate_tot = rate_deriv[0] * mole.species[ rate->reactants[0]->index ].den;
		
		for(i=0;i<rate->nreactants;i++)
		{	
			sp = rate->reactants[i];
			if (rate->rvector[i] == NULL)
			{
				mole.species[sp->index].snk += rate_deriv[i];
			}
			else
			{
				// exclude charge exchange with isotopes of same element 
				if( rate->nreactants==2 )
				{
					long otherIndex = 1-i;
					molecule *sp2 = rate->reactants[otherIndex];
					if( sp->isMonatomic() && sp2->isMonatomic() && sp->nAtom.begin()->first->el == sp2->nAtom.begin()->first->el )
						continue;
				}

				if ( atmdat.lgCTOn )
				{
					for( ChemAtomList::iterator atom = unresolved_atom_list.begin(); atom != unresolved_atom_list.end(); ++atom)
					{
						nelem = (*atom)->el->Z-1;
						if (sp->nAtom.find(*atom) != sp->nAtom.end() && sp->nAtom[*atom] != 0 && rate->rvector[i]->charge != sp->charge)
						{
							mole.xMoleChTrRate[nelem][sp->charge][rate->rvector[i]->charge] +=
								(realnum) rate_deriv[i];
							break;
						}
					}
				}
			}
		}
		
		for(i=0;i<rate->nproducts;i++)
		{
			sp = rate->products[i];
			if (rate->pvector[i] == NULL)
			{
				mole.species[sp->index].src += rate_tot;
			}
		}
	}

	for (ChemAtomList::iterator atom = unresolved_atom_list.begin();
		  atom != unresolved_atom_list.end(); ++atom)
	{
		const long int nelem=(*atom)->el->Z-1;
		if( !dense.lgElmtOn[nelem] )
			continue;

		for (long int ion=0;ion<nelem+2;ion++) 
		{
			if ((*atom)->ipMl[ion] != -1)
			{
				mole.source[nelem][ion] = mole.species[(*atom)->ipMl[ion]].src;
				mole.sink[nelem][ion] = mole.species[(*atom)->ipMl[ion]].snk;
			}
			else
			{
				mole.source[nelem][ion] = 0.0;
				mole.sink[nelem][ion] = 0.0;
			}
		}
	}

	//mole_dominant_rates(findspecies("H+"),ioQQQ);
	return;
}

STATIC bool lgNucleiConserved(const multi_arr<double,2> &c)
{
	DEBUG_ENTRY ("lgNucleiConserved()");
	bool checkAllOK = true;
	multi_arr<double,2> test(atom_list.size(),mole_global.num_calc),
		tot(atom_list.size(),mole_global.num_calc);
	vector<double> ccache(mole_global.num_calc);
	vector<long> ncache(mole_global.num_calc);
	test.zero();
	tot.zero();

	map<chem_atom*, long> atom_to_index;
	for (unsigned long j=0; j<atom_list.size(); ++j )
	{
		atom_to_index[atom_list[j].get_ptr()] = j;
	}
	for (long j=0;j<mole_global.num_calc;j++) 
	{
		long nc = 0;
		for (long i=0;i<mole_global.num_calc;i++) 
		{
			if (c[i][j] != 0.)
			{
				ccache[nc] = c[i][j];
				ncache[nc] = i;
				++nc;
			}
		}
		if (nc > 0)
		{
			for (molecule::nAtomsMap::const_iterator el = mole_global.list[j]->nAtom.begin(); 
				  el != mole_global.list[j]->nAtom.end(); ++el)
			{
				long natom = atom_to_index[el->first.get_ptr()];
				const int nAtomj = el->second;
				for (long i=0;i<nc;i++) 
				{
					const double term = ccache[i] * nAtomj;
					test[natom][ncache[i]] += term;
					tot[natom][ncache[i]] += fabs(term);
				}
			}
		}
	}

	for( unsigned long natom=0; natom < atom_list.size(); ++natom)
	{
		for (long i=0;i<mole_global.num_calc;i++) 
		{
			const bool checkOK = 
				( fabs(test[natom][i]) <= MAX2(1e-10*tot[natom][i], 1e10*DBL_MIN) );
			if ( UNLIKELY(!checkOK) )
			{
				chem_atom *atom = atom_list[natom].get_ptr();
				fprintf(stdout,"Network conservation error %s %s %g %g %g %g\n",
						  atom->label().c_str(),
						  mole_global.list[i]->label.c_str(),
						  test[natom][i],
						  test[natom][i]/tot[natom][i],
						  mole.species[atom->ipMl[0]].den,
						  mole.species[atom->ipMl[1]].den);
				//fprintf(stdout,"Problem at %s\n",rate->label);
				checkAllOK = false;
			}
		}
	}
	return checkAllOK;
}

void mole_dominant_rates( const molecule *debug_species, FILE *ioOut )
{
	long int i, j;
	mole_reaction *rate, *ratesnk=NULL, *ratesrc=NULL;
	double rate_tot, rate_deriv[MAXREACTANTS], rk;
	molecule *sp;
	double snkx=0.,srcx=0.;

	for(mole_reaction_i p=mole_priv::reactab.begin();
			p != mole_priv::reactab.end(); ++p)
	{
		rate = &(*p->second);
		rk = mole.reaction_rks[ rate->index ];

		for(i=0;i<rate->nreactants;i++)
		{
			rate_deriv[i] = rk;
			for(j=0;j<rate->nreactants;j++)
			{
				if(i!=j)
				{
					rate_deriv[i] *= mole.species[ rate->reactants[j]->index ].den;
				}
			}
		}

		rate_tot = rate_deriv[0] * mole.species[ rate->reactants[0]->index ].den;

		if (debug_species != null_mole)
		{
			for(i=0;i<rate->nproducts;++i)
			{
				sp = rate->products[i];
				if (sp == debug_species && rate->pvector[i] == NULL)
				{
					if (fabs(rate_tot) > srcx)
					{
						srcx = rate_tot;
						ratesrc = rate;
					}
				}
			}
			for(i=0;i<rate->nreactants;++i)
			{
				sp = rate->reactants[i];
				if (sp == debug_species && rate->rvector[i] == NULL)
				{
					if (fabs(rate_deriv[i]) > snkx)
					{
						snkx = rate_deriv[i];
						ratesnk = rate;
					}
				}
			}
		}
	}


	if (debug_species != null_mole)
	{
		if (ratesrc)
		{
			fprintf( ioOut, "%20.20s src %13.7g of %13.7g [",
				ratesrc->label.c_str(),srcx,mole.species[debug_species->index].src);
			for (j=0;j<ratesrc->nreactants;j++)
			{
				if (j)
				{
					fprintf( ioOut, "," );
				}
				fprintf( ioOut, "%-6.6s %13.7g",
					ratesrc->reactants[j]->label.c_str(),
					mole.species[ ratesrc->reactants[j]->index ].den);
			}
			fprintf( ioOut, "]" );
		}
		if (ratesnk)
		{
			fprintf( ioOut, "%20.20s snk %13.7g of %13.7g [",
				ratesnk->label.c_str(), snkx * mole.species[debug_species->index].den,
				mole.species[debug_species->index].snk * mole.species[debug_species->index].den);
			for (j=0;j<ratesnk->nreactants;j++)
			{
				if (j)
				{
					fprintf( ioOut, "," );
				}
				fprintf( ioOut, "%-6.6s %13.7g",
					ratesnk->reactants[j]->label.c_str(),
					mole.species[ ratesnk->reactants[j]->index ].den);
			}
			fprintf( ioOut, "]" );
		}
	}
	fprintf( ioOut, "\n" );

	return;
}