/* 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 */
/*lint -e817 pclint generates tons of bogus nelem < 0 warnings*/
/*ion_recomb generate recombination coefficients for any species */
/*ion_recombAGN generate recombination coefficients for AGN table */
#include "cddefines.h"
#include "phycon.h"
#include "tfidle.h"
#include "heavy.h"
#include "hmi.h"
#include "dense.h"
#include "iso.h"
#include "abund.h"
#include "punch.h"
#include "elementnames.h"
#include "converge.h"
#include "atmdat.h"
#include "ionbal.h"

/*ion_recomb generate recombination coefficients for any species */
void ion_recomb(
  /* this is debug flag */
  int lgPrintIt,
  double *dicoef, 
  double *dite, 
  double ditcrt[], 
  double aa[], 
  double bb[], 
  double cc[], 
  double dd[], 
  double ff[], 
  /* nelem is the atomic number on the C scale, 0 for H */
  long int nelem/*, 
  double tlow[]*/)
{
#define DICOEF(I_,J_)	(*(dicoef+(I_)*(nelem+1)+(J_)))
#define DITE(I_,J_)	(*(dite+(I_)*(nelem+1)+(J_)))
	long int i, 
	  ion, 
	  limit;
	double 
	  fac2, 
	  factor,
	  DielRecomHiT[LIMELM] ,
	  DielRecomLowT[LIMELM] ,
	  Recom3Bod[LIMELM] ,
	  ChargeTransfer[LIMELM] ,
	  t4m1, 
	  tefac;

	/* these are guessed diel rec rates, taken from 
	  >>refer	all	diel	Ali, B., Blum, R. D., Bumgardner, T. E., Cranmer, S. R., Ferland, G. J., 
	  >>refercon Haefner, R. I., & Tiede, G. P. 1991, PASP, 103, 1182 */
	static double cludge[4]={3e-13,3e-12,1.5e-11,2.5e-11};

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

	/* set up ionization balance matrix, C(I,1)=destruction, 2=creation
	 * heating rates saved in array B(I) in same scrat block
	 * factor is for Aldrovandi+Pequignot fit, FAC2 is for Nuss+Storey
	 * fit for dielectronic recomb
	 * GrnIonRec is rate ions recombine on grain surface, normally zero;
	 * set in hmole, already has factor of hydrogen density
	 * */

	/* routine only used for Li on up */
	ASSERT( nelem < LIMELM);
	ASSERT( nelem > 1 );

	/* check that range of ionization is correct */
	ASSERT( dense.IonLow[nelem] >= 0 );
	ASSERT( dense.IonLow[nelem] <= nelem+1 );

	atmdat.nsbig = MAX2(dense.IonHigh[nelem]+1,atmdat.nsbig);
	t4m1 = 1e4/phycon.te;
	fac2 = 1e-14*phycon.sqrte;

	/* total grain surface ionization and recombination */

	/* this routine only does simple two-level species, 
	 * loop over ions will be <= limit, IonHigh is -1 since very
	 * highest stage of ionization is not recombined into.  
	 * for Li, will do only atom, since ions are H and He like,
	 * so limit is zero */
	limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);

	/* zero-out loop comes before main loop since there are off-diagonal
	 * elements in the main ionization loop, due to multi-electron processes */
	/* >>chng 00 dec 07, limit changed to identical to bidiag */
	for( ion=0; ion <= limit; ion++ )
	{
		ionbal.RateRecomTot[nelem][ion] = 0.;
		ionbal.RadRecomRateCoef[nelem][ion] = 0.;
		Recom3Bod[ion] = 0.;
		DielRecomHiT[ion] = 0.;
		ChargeTransfer[ion] = 0.;
		DielRecomLowT[ion] = 0.;
	}
	for( ion=limit+1; ion < LIMELM; ion++ )
	{
		/* >>chng 01 dec 18, do not set this to FLT_MAX since it clobbers what
		 * had been set in h-like and he-like routines - that would only affect
		 * the printout */
		Recom3Bod[ion] = -FLT_MAX;
		DielRecomHiT[ion] = -FLT_MAX;
		ChargeTransfer[ion] = -FLT_MAX;
		DielRecomLowT[ion] = -FLT_MAX;
	}

	/* H-like and He-like stages use results from hydrogenic series done by hydro... routines */
	/* this was evaluated in hydro routines, and must not be changed here
	 * since we want the bidiag routine to get the same answer */
	Recom3Bod[nelem] = ionbal.RateRecomTot[nelem][nelem]/dense.eden ;
	/* now do he-like */
	Recom3Bod[nelem-1] = ionbal.RateRecomTot[nelem][nelem-1]/dense.eden ;

	DielRecomHiT[nelem] = 0. ;
	DielRecomHiT[nelem-1] = 0. ;

	DielRecomLowT[nelem] = 0. ;
	DielRecomLowT[nelem-1] = 0. ;

	/* these are counted elsewhere and must not be added here */
	Heavy.xLyaHeavy[nelem][nelem] = 0.;
	Heavy.xLyaHeavy[nelem][nelem-1] = 0.;

	/*--------------------------------------------------------------------- */

	/* IonLow is 0 for the atom, limit chosen to NOT include iso sequences  */
	for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
	{
		/* number of bound electrons of the ion after recombination, 
		 * for an atom (ion=0) this is
		 * equal to nelem+1, the element on the physical scale, since nelem is 
		 * on the C scale, being 5 for carbon */
		long n_bnd_elec_after_recomb = nelem+1 - ion;

		/* >>chng 03 nov 03, n_bnd_elec_after_recomb had been on scale, so needed to add 1,
		 * now on physics scale.  Also removed lgVrrFit which was never used */
		/* this branch dima's fits - routine takes arguments on physics scale,
		 * units are cm3 s-1 */
		/* save radiative recombination rate coef since needed for rec cooling,
		 * this should be pure rad rec, and not include density, so rate coef, cm3 s-1 */
		ionbal.RadRecomRateCoef[nelem][ion] = 
			atmdat_rad_rec( nelem+1 , n_bnd_elec_after_recomb , phycon.te );

		/* this is total, has units s-1 */
		ionbal.RateRecomTot[nelem][ion] = dense.eden*(ionbal.RadRecomRateCoef[nelem][ion] + 
		  ionbal.CotaRate[ion]) ;

		/* >>chng 02 nov 06 add charge transfer here rather than in ion_solver */
		/* charge transfer recombination of this species by ionizing hydrogen and helium */
		ChargeTransfer[ion] = 
			/* >>chng 02 may 05, should be just ground state, make it so */
			/*atmdat.HeCharExcRecTo[nelem][ion]*dense.xIonDense[ipHELIUM][0] + */
			atmdat.HeCharExcRecTo[nelem][ion]*iso.Pop2Ion[ipHE_LIKE][ipHELIUM][ipHe1s1S]*dense.xIonDense[ipHELIUM][1] + 
			/*atmdat.HCharExcRecTo[nelem][ion]*dense.xIonDense[ipHYDROGEN][0];*/
			/* >>chng 01 06 01 should be ground state only, change to do so */
			atmdat.HCharExcRecTo[nelem][ion]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.xIonDense[ipHYDROGEN][1];
		/*>>chng 04 feb 20, add this, had always been in for dest for H- */
		/* charge transfer recombination of first ion to neutral, by reaction with H- 
		 * the ion==0 is right, the first array element is the */
		if( ion==0 && nelem>ipHELIUM && atmdat.lgCTOn )
			ChargeTransfer[ion] += hmi.hmin_ct_firstions * hmi.Hmolec[ipMHm];
		/*if(nzone>140&&nelem==5&&ion==0 )
			fprintf(ioQQQ,"DEBUG ct\t%.2f\t%.4e\t%.4e\t%.4\n",
			fnzone,
			dense.eden,
			hmi.Hmolec[ipMHm],
			dense.xIonDense[ipCARBON][1]);*/

		/* now copy this into main array */
		ionbal.RateRecomTot[nelem][ion] += ChargeTransfer[ion];

		/* sum of rate coefficients for radiative, three body, charge transfer, and
		 * recombination on grain surfaces, this is only used for trace
		 * printout below */
		ionbal.RateRecomTot[nelem][ion] += Recom3Bod[ion];

		/* net recombination rate on grain surfaces, calculated in grain_ionrec.c
		ionbal.RateRecomTot[nelem][ion] += ionbal.GrainCreat[nelem][ion]; */

		/* Burgess or high-T dielectronic recombination */
		DielRecomHiT[ion] = 0.;
		/* >>chng 03 oct 29, do fe here rather than after loop */
		if( phycon.te > (float)(ditcrt[ion]*0.1) && nelem!=ipIRON )
		{
			DielRecomHiT[ion] = ionbal.DielSupprs[0][ion]/phycon.te32*
			  DICOEF(0,ion)*exp(-DITE(0,ion)/phycon.te)*
			  (1. + DICOEF(1,ion)*
			  sexp(DITE(1,ion)/phycon.te));
		}
		else if( nelem==ipIRON )
		{
			/* >>chng 03 oct 29, do this here rather than after loop */
			DielRecomHiT[ion] = atmdat_dielrec_fe(ion+1,phycon.te);
		}
		ionbal.RateRecomTot[nelem][ion] += DielRecomHiT[ion]*dense.eden;

		/* dn = 0 dielectronic recombination
		 * do not include it for rec from
		 * a closed shell, n_bnd_elec_after_recomb-1 is number of bound electrons in parent ion */
		DielRecomLowT[ion] = 0.;
		if( ((n_bnd_elec_after_recomb-1) !=  2) && 
			((n_bnd_elec_after_recomb-1) != 10) && 
			((n_bnd_elec_after_recomb-1) != 18) )
		{
			tefac = ff[ion]*t4m1;
			/* fac2 = 1e-14 * sqrte */
			/* >>chng 03 oct 29, do not do iron here since done above */
			if( ff[ion] != 0. && fabs(tefac) < 30. && nelem!=ipIRON )
			{
				factor = (((aa[ion]*t4m1+bb[ion])*t4m1+cc[ion])*t4m1+dd[ion])* exp(-tefac);
				DielRecomLowT[ion] = ionbal.DielSupprs[1][ion]*fac2*
					MAX2(0.,factor );
			}
			/* >>chng 03 oct 29, do not do iron here since done above */
			else if( (phycon.te > 100. && ff[ion] == 0.) && ion <= 3  )
			{
				/* >>chng 01 jun 30, make GuessDiel array */
				DielRecomLowT[ion] = ionbal.DielSupprs[1][ion]*cludge[ion]*
				  ionbal.GuessDiel[ion];
			}
			/* this implements the steve option on the dielec command */
			else if( ionbal.lg_guess_coef && ion>3)
			{
				/* assume the recom rate is the coefficient scanned off the steve option, times the charge
				 * before recombination raised to a power */
				static double RecNoise[LIMELM][LIMELM];
				static int lgNoiseEval=FALSE;
				double fitcoef[3][3] = 
				{
					/* L- shell, Z=17-23 */
					{-52.5073,+5.19385,-0.126099} ,
					/*  M-shell (Z=9-15) */
					{-10.9679,1.66397,-0.0605965} ,
					/* N shell */
					{-3.95599,1.61884,-0.194540} 
				};

				long int nelem_save = nelem , ion_save = ion;

				/* option to put noise into rec coef */
				if( ionbal.guess_noise !=0. && !lgNoiseEval )
				{
					for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
					{
						for( ion=0; ion<=nelem; ++ion )
						{
							/* log normal noise with dispersion entered on command line */
							/* NB the seed for rand was set when the command was parsed */
							RecNoise[nelem][ion] = pow(10., RandGauss( 0. , ionbal.guess_noise ) );
						}
					}
				}
				else if( !lgNoiseEval )
				{
					for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
					{
						for( ion=0; ion<=nelem; ++ion )
						{
							RecNoise[nelem][ion] = 1.;
						}
					}
				}
				lgNoiseEval = TRUE;
				nelem = nelem_save;
				ion = ion_save;
				if( nelem==ipIRON )
				{
					int nshell;
					if( (n_bnd_elec_after_recomb>=4) && (n_bnd_elec_after_recomb<=11) )
						nshell = 0;
					else if( (n_bnd_elec_after_recomb>=12) && (n_bnd_elec_after_recomb<=19 ))
						nshell = 1;
					else
						nshell = 2;
					/* n_bnd_elec_after_recomb is the number of bound electrons */
					DielRecomLowT[ion] = fitcoef[nshell][0] +
						fitcoef[nshell][1]*(ion+1) + 
						fitcoef[nshell][2]*POW2(ion+1.) ;
					DielRecomLowT[ion] = 1e-10*pow(10.,DielRecomLowT[ion]);

				}
				else
				{
					/* this is guess for all other elements */
					DielRecomLowT[ion] = 3e-12*pow(10.,(double)(ion+1)*0.1);
				}
				/* include optional noise here */
				DielRecomLowT[ion] *= RecNoise[nelem][ion];
			}
            ionbal.RateRecomTot[nelem][ion] += DielRecomLowT[ion]* dense.eden;
		}

		/* >>chng 01 jun 30, FRAC_LINE was 0.1, not 1, did not include anything except
		 * radiative recombination, the radrec term */
#		define FRAC_LINE 1.
		/* was 0.1 */
		/*Heavy.xLyaHeavy[nelem][ion] = (float)(dense.eden*radrec*FRAC_LINE );*/
		Heavy.xLyaHeavy[nelem][ion] = (float)(dense.eden*
			(ionbal.RadRecomRateCoef[nelem][ion]+DielRecomLowT[ion]+DielRecomHiT[ion])*FRAC_LINE );
	}

	/* option to punch rec coef */
	if( punch.lgioRecom || lgPrintIt )
	{
		/* >>chng 04 feb 22, make option to print ions for single element */
		FILE *ioOut;
		if( lgPrintIt )
			ioOut = ioQQQ;
		else
			ioOut = punch.ioRecom;

		/* print name of element */
		fprintf( ioOut, 
			" %s recombination coefficients %.2f \tte\t%.4e\tne\t%.4e\n", 
			elementnames.chElementName[nelem] , fnzone , phycon.te , dense.eden );

		limit = MIN2(11,dense.IonHigh[nelem]);
		for( i=0; i < limit; i++ )
		{
			fprintf( ioOut, "%10.2e", ionbal.RadRecomRateCoef[nelem][i] );
		}
		fprintf( ioOut, " radiative vs Z\n" );

		for( i=0; i < limit; i++ )
		{
			fprintf( ioOut, "%10.2e", ChargeTransfer[i]/dense.eden );
		}
		fprintf( ioOut, " CT/ne\n" );

		for( i=0; i < limit; i++ )
		{
			fprintf( ioOut, "%10.2e", Recom3Bod[i]/dense.eden );
		}
		fprintf( ioOut, " 3body vs Z /ne\n" );

		/* note different upper limit - this routine does grain rec for all ions */
		for( i=0; i < dense.IonHigh[nelem]; i++ )
		{
			fprintf( ioOut, "%10.2e", ionbal.GrainCreat[nelem][i]/dense.eden );
		}
		fprintf( ioOut, " Grain vs Z /ne\n" );

		for( i=0; i < limit; i++ )
		{
			fprintf( ioOut, "%10.2e", DielRecomHiT[i] );
		}
		fprintf( ioOut, " Burgess vs Z\n" );

		for( i=0; i < limit; i++ )
		{
			fprintf( ioOut, "%10.2e", DielRecomLowT[i] );
		}
		fprintf( ioOut, " Nussbaumer Storey vs Z\n" );

		/* total recombination rate, with density included - this goes into the matrix */
		for( i=0; i < limit; i++ )
		{
			fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i] );
		}
		fprintf( ioOut, 
			" total rec rate (with density) for %s\n", 
			elementnames.chElementSym[nelem] );
		for( i=0; i < limit; i++ )
		{
			fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i]/dense.eden );
		}
		fprintf( ioOut, 
			" total rec rate / ne for %s\n\n", 
			elementnames.chElementSym[nelem] );

		/* spill over to next line for many stages of ionization */
		if( dense.IonHigh[nelem] > 11 )
		{
			limit = MIN2(29,dense.IonHigh[nelem]);
			fprintf( ioOut, " R " );
			for( i=11; i < limit; i++ )
			{
				fprintf( ioOut, "%10.2e", Recom3Bod[i] );
			}
			fprintf( ioOut, "\n" );

			fprintf( ioOut, " B " );
			for( i=11; i < limit; i++ )
			{
				fprintf( ioOut, "%10.2e", DielRecomHiT[i] );
			}
			fprintf( ioOut, "\n" );

			fprintf( ioOut, " NS" );
			for( i=11; i < limit; i++ )
			{
				fprintf( ioOut, "%10.2e", DielRecomLowT[i] );
			}
			fprintf( ioOut, "\n" );

			fprintf( ioOut, "   " );
			for( i=11; i < limit; i++ )
			{
				fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i] );
			}
			fprintf( ioOut, "\n\n" );
		}
	}

	/* >>chng 02 nov 09, from -2 to -NISO */
	/*limit = MIN2(nelem-2,dense.IonHigh[nelem]-1);*/
	limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
	for( i=dense.IonLow[nelem]; i <= limit; i++ )
	{
		ASSERT( Heavy.xLyaHeavy[nelem][i] > 0. );
		ASSERT( ionbal.RateRecomTot[nelem][i] > 0. );
	}

#	ifdef DEBUG_FUN
	fputs( " <->ion_recomb(FALSE,)\n", debug_fp );
#	endif
	return;
#	undef	DITE
#	undef	DICOEF
}

/*ion_recombAGN generate recombination coefficients for AGN table */
void ion_recombAGN( FILE * io )
{
#	define N1LIM 3
#	define N2LIM 4
	double te1[N1LIM]={ 5000., 10000., 20000.};
	double te2[N2LIM]={ 20000.,50000.,100000.,1e6};
	/* this is boundary between two tables */
	double BreakEnergy = 100./13.0;
	long int nelem, ion , i;
	/* this will hold element symbol + ionization */
	char chString[100],
		chOutput[100];
	/* save temp here	*/
	float TempSave = phycon.te;
	/* save ne here	*/
	double EdenSave = dense.eden;

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

	dense.eden = 1.;
	/*atmdat_readin();*/

	/* first put header on file */
	fprintf(io,"X+i\\Te");
	for( i=0; i<N1LIM; ++i )
	{
		phycon.te = (float)te1[i];
		fprintf(io,"\t%.0f K",phycon.te);
	}
	fprintf(io,"\n");

	/* now do loop over temp, but add elements */
	for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem )
	{
		/* this list of elements included in the AGN tables is defined in zeroabun.c */
		if( abund.lgAGN[nelem] )
		{
			for( ion=0; ion<=nelem; ++ion )
			{
				ASSERT( Heavy.Valence_IP_Ryd[nelem][ion] > 0.05 );

				if( Heavy.Valence_IP_Ryd[nelem][ion] > BreakEnergy )
					break;

				/* print chemical symbol */
				sprintf(chOutput,"%s", 
					elementnames.chElementSym[nelem]);
				/* some elements have only one letter - this avoids leaving a space */
				if( chOutput[1]==' ' )
					chOutput[1] = chOutput[2];
				/* now ionization stage */
				if( ion==0 )
				{
					sprintf(chString,"0 ");
				}
				else if( ion==1 )
				{
					sprintf(chString,"+ ");
				}
				else
				{
					sprintf(chString,"+%li ",ion);
				}
				strcat( chOutput , chString );
				fprintf(io,"%5s",chOutput );

				for( i=0; i<N1LIM; ++i )
				{
					phycon.te = (float)te1[i];
					tfidle(FALSE);
					dense.IonLow[nelem] = 0;
					dense.IonHigh[nelem] = nelem+1;
					ConvBase(0);
					fprintf(io,"\t%.2e",ionbal.RateRecomTot[nelem][ion]);
				}
				fprintf(io,"\n");
			}
			fprintf(io,"\n");
		}
	}

	/* second put header on file */
	fprintf(io,"X+i\\Te");
	for( i=0; i<N2LIM; ++i )
	{
		phycon.te = (float)te2[i];
		tfidle(FALSE);
		fprintf(io,"\t%.0f K",phycon.te);
	}
	fprintf(io,"\n");

	/* now do same loop over temp, but add elements */
	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
	{
		/* this list of elements included in the AGN tables is defined in zeroabun.c */
		if( abund.lgAGN[nelem] )
		{
			for( ion=0; ion<=nelem; ++ion )
			{
				ASSERT( Heavy.Valence_IP_Ryd[nelem][ion] > 0.05 );

				if( Heavy.Valence_IP_Ryd[nelem][ion] <= BreakEnergy )
					continue;

				/* print chemical symbol */
				fprintf(io,"%s", 
					elementnames.chElementSym[nelem]);
				/* now ionization stage */
				if( ion==0 )
				{
					fprintf(io,"0 ");
				}
				else if( ion==1 )
				{
					fprintf(io,"+ ");
				}
				else
				{
					fprintf(io,"+%li",ion);
				}

				for( i=0; i<N2LIM; ++i )
				{
					phycon.te = (float)te2[i];
					tfidle(FALSE);
					dense.IonLow[nelem] = 0;
					dense.IonHigh[nelem] = nelem+1;
					ConvBase(0);
					fprintf(io,"\t%.2e",ionbal.RateRecomTot[nelem][ion]);
				}
				fprintf(io,"\n");
			}
			fprintf(io,"\n");
		}
	}

	phycon.te = TempSave;
	tfidle(TRUE);
	dense.eden = EdenSave;

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

	return;
}
/*lint +e817 pclint generates tons of bogus nelem < 0 warnings*/

