/* 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 */
#include "cddefines.h"
#include "dense.h"
#include "ionbal.h"
#include "converge.h"
#include "grainvar.h"
#include "grains.h"

/*grain_ion_charge_transfer generate grain surface recombination and ionization rates for all ions */
void grain_ion_charge_transfer( void )
{
	int ion,
	  ion_from,
	  ion_to,
	  nelem;
	int lgFirstPass;
	static int lgMustInit=TRUE;
	static long int nzone_used[LIMELM];
	static double abund_old[LIMELM][LIMELM+1];

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

	/* ionbal.lgGrainIonRecom - option to turn this physics off with 
	 * no grain neutralization command;  also do not do if grains, or grain physics,
	 * is not enabled */
	if( !ionbal.lgGrainIonRecom || !gv.lgGrainPhysicsOn || !gv.lgDustOn )
		return;

	/* much of this code is copied from ion_solver and is to stabalize ratio of ion densities,
	 * so that we can keep bidiagonal solver */
	if( lgMustInit )
	{
		/* loop is really over nelem not ion */
		for( nelem=0; nelem < LIMELM; ++nelem )
			nzone_used[nelem] = -1;
		lgMustInit = FALSE;
	}

	for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
	{
		/* skip if element is not turned on */
		if( !dense.lgElmtOn[nelem] )
			continue;

		/* set flag saying whether we have already saved previous
		 * ion abundances for this zone */
		lgFirstPass = FALSE;
		if( nzone_used[nelem] != nzone )
		{
			lgFirstPass = TRUE;
			nzone_used[nelem] = nzone;
		}
		/* in this case some entries in abund_old might be invalid */
		if( conv.lgIonStageTrimed )
			lgFirstPass = TRUE;

		/* find grain destruction rates for all ions that exist */
		for( ion_from=dense.IonLow[nelem]; ion_from < dense.IonHigh[nelem]; ion_from++ )
		{
			/* the total efective grain destruction rate is just the sum of all processes
			 * that change this ion into something else */
			ionbal.GrainDestr[nelem][ion_from] = 0.;

			/*>>chng 04 feb 16, -> PvH treatment of recombinatinon on grain surface
			 *>>chng 04 apr 10, complete upgrade to formal treatment of all ion rec processes */
			for( ion_to=ion_from+1; ion_to <= nelem+1; ++ion_to )
			{
				/* do not do ion onto itself */
				if( ion_from != ion_to )
				{
					/* this is a destruction process for ion, representing changes to other stages
					 * of ionization, either above or below current ion
					 * in routine ion_recomb this will be added to total recombination rate 
					 * units of GrainDestr is s-1 - used to multiply ion density to get destruction
					 * per unit time & vol
					 * this is addded to ionbal.RateIonizTot[nelem][ion] in ion_solver */
					ionbal.GrainDestr[nelem][ion_from] += gv.GrainChTrRate[nelem][ion_from][ion_to];
				}
			}
		}

		/* find grain creation rates for all ions that exist */
		for( ion_to=dense.IonLow[nelem]; ion_to < dense.IonHigh[nelem]; ion_to++ )
		{
			ionbal.GrainCreat[nelem][ion_to] = 0.;

			/*>>chng 04 feb 16, -> PvH treatment of recombinatinon on grain surface
			 *>>chng 04 apr 10, complete upgrade to formal treatment of all ion rec processes */
			for( ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from )
			{
				/* double ratio_ion_density;
				do not do ion onto itself */
				if( ion_from == ion_to )
					continue;

				/* to avoid off diagonal matrix elements do trick to keep bi-diagonal,
				 * as explained towards ened of section on heavy elements in part 2 of hazy */
				/* units are cm3 s-1, will be multiplied by ion above ion to get rate is s-1 */
				if( ion_from == ion_to+1 )
				{
					/* this is "classical" recombination, from one stage above to this ion,
					 * this is what the solver will expect to see, so no trick necessary */
					ionbal.GrainCreat[nelem][ion_to] += gv.GrainChTrRate[nelem][ion_from][ion_to];
				}
#				if 0
				else
				{
					/* will div by this number - make sure it is safe */
					if( dense.xIonDense[nelem][ion_to+1] > 1e-30 )
					{
						if( /*nzone < 3 ||*/ lgFirstPass || abund_old[nelem][ion_to+1] < 1e-30 )
						{
							/* ion_to points to current stage of ionization */
							ratio_ion_density = (double)(dense.xIonDense[nelem][ion_from])/
								(double)(dense.xIonDense[nelem][ion_to+1]);
						}
						else
						{
							double frac_old = 0.5;
							/* ion_to points to current stage of ionization */
							ratio_ion_density = 
								frac_old*abund_old[nelem][ion_from]/abund_old[nelem][ion_to+1] +
								(1.-frac_old)*(double)dense.xIonDense[nelem][ion_from]/
								(double)dense.xIonDense[nelem][ion_to+1];
						}
					}
					else
					{
						ratio_ion_density = 1.;
					}
					ionbal.GrainCreat[nelem][ion_to] += gv.GrainChTrRate[nelem][ion_from][ion_to]*
						(float)ratio_ion_density;
				}
#				endif
			}/* loop over incoming ion */
		}/* loop over outgoing ion */

		/*if( nelem==ipZINC )
		{
			fprintf(ioQQQ,"DEBUG zn grn\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
				gv.GrainChTrRate[nelem][1][0],
				gv.GrainChTrRate[nelem][2][0],
				gv.GrainChTrRate[nelem][2][1],
				ionbal.GrainCreat[nelem][0],
				ionbal.GrainCreat[nelem][1],
				ionbal.GrainCreat[nelem][2],
				ionbal.GrainDestr[nelem][0],
				ionbal.GrainDestr[nelem][1],
				ionbal.GrainDestr[nelem][2]);
		}*/
		/* save ionization balance for this element */
		for( ion=dense.IonLow[nelem]; ion <= dense.IonHigh[nelem]; ion++ )
		{
			abund_old[nelem][ion] = dense.xIonDense[nelem][ion];
		}

	}/* loop over elements */

	/* option to produce some debug information */
	{
		/*@-redef@*/
		enum{DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			nelem = ipCARBON;
			fprintf(ioQQQ,"DEBUG grain ion charge transf\t%.2f\t%li\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
				fnzone,
				dense.IonLow[nelem], 
				dense.IonHigh[nelem],
				ionbal.GrainCreat[nelem][0],
				gv.GrainChTrRate[nelem][1][0],
				gv.GrainChTrRate[nelem][2][0],
				dense.xIonDense[nelem][1]/SDIV(dense.xIonDense[nelem][0]),
				dense.xIonDense[nelem][2]/SDIV(dense.xIonDense[nelem][0]));
		}
	}

	/*fprintf(ioQQQ,"DEBUG grn recion\t%.2f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
		fnzone,
		gv.GrainChTrRate[ipSILICON][0][1],
		gv.GrainChTrRate[ipSILICON][1][2],
		gv.GrainChTrRate[ipSILICON][0][2],
		gv.GrainChTrRate[ipSILICON][1][0],
		gv.GrainChTrRate[ipSILICON][2][1],
		gv.GrainChTrRate[ipSILICON][2][0]
	);*/

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

	return;
}
