/* 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 */
/*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 "heavy.h"
#include "hmi.h"
#include "grainvar.h"
#include "dense.h"
#include "conv.h"
#include "thermal.h"
#include "iso.h"
#include "abund.h"
#include "punch.h"
#include "elementnames.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,
	  DielRecomRateCoef_HiT[LIMELM] ,
	  DielRecomRateCoef_LowT[LIMELM] ,
	  ChargeTransfer[LIMELM] ,
	  t4m1, 
	  tefac;

	/* these are used for adding noise to rec coef */
	static double RecNoise[LIMELM][LIMELM];
	static int lgNoiseNeedEval=TRUE;

#	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;

	/* option to put noise into rec coef -
	 * one time initialization of noise  - this is set with command
	 * set dielectroninc recombiantion kludge noise */
	if( lgNoiseNeedEval )
	{
		int n;
		if( ionbal.guess_noise !=0. )
		{
			for( n=ipHYDROGEN; n<LIMELM; ++n )
			{
				for( ion=0; ion<=n; ++ion )
				{
					/* log normal noise with dispersion entered on command line */
					/* NB the seed for rand was set when the command was parsed */
					RecNoise[n][ion] = pow(10., RandGauss( 0. , ionbal.guess_noise ) );
				}
			}
		}
		else
		{
			for( n=ipHYDROGEN; n<LIMELM; ++n )
			{
				for( ion=0; ion<=n; ++ion )
				{
					RecNoise[n][ion] = 1.;
				}
			}
		}
		lgNoiseNeedEval = FALSE;
	}

	/* 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.;
		ChargeTransfer[ion] = 0.;
		DielRecomRateCoef_LowT[ion] = 0.;
		DielRecomRateCoef_HiT[ion] = 0.;
		ionbal.RR_rate_coef_used[nelem][ion] = 0.;
		ionbal.DR_rate_coef_used[nelem][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 */
		ChargeTransfer[ion] = -FLT_MAX;
		DielRecomRateCoef_LowT[ion] = -FLT_MAX;
		DielRecomRateCoef_HiT[ion] = -FLT_MAX;
	}

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

	DielRecomRateCoef_LowT[nelem] = 0. ;
	DielRecomRateCoef_LowT[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 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] = 
			/* He0 + ion charge transfer recombination */
			atmdat.HeCharExcRecTo[nelem][ion]*
			/* following product is density [cm-3] of ground state of He0 */
			iso.Pop2Ion[ipHE_LIKE][ipHELIUM][ipHe1s1S]*dense.xIonDense[ipHELIUM][1] + 
			/* H0 + ion charge transfer recombination */
			atmdat.HCharExcRecTo[nelem][ion]*
			/* following product is density [cm-3] of ground state of H0 */
			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];

		/* >>chng 06 feb 01, add option to use Badnell RR data rather than Verner */
		if( ionbal.lgRR_recom_Badnell_use && ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] )
		{
			ionbal.RR_rate_coef_used[nelem][ion] = ionbal.RR_Badnell_rate_coef[nelem][ion];
		}
		else
		{
			ionbal.RR_rate_coef_used[nelem][ion] = ionbal.RR_Verner_rate_coef[nelem][ion];
		}

		/* Burgess or high-T dielectronic recombination */
		DielRecomRateCoef_HiT[ion] = 0.;
		/* >>chng 03 oct 29, do fe here rather than after loop */
		if( nelem==ipIRON )
		{
			/* implements dn > 0 DR from Arnaud & Raymond 1992 */
			DielRecomRateCoef_HiT[ion] = atmdat_dielrec_fe(ion+1,phycon.te);
		}
		else if( phycon.te > (ditcrt[ion]*0.1) )
		{
			DielRecomRateCoef_HiT[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));
		}

		/* begin 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 */
		DielRecomRateCoef_LowT[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;
			/* do not do iron here since all dn = 0 DR either Badnell or a guess */
			if( ff[ion] != 0. && nelem!=ipIRON )
			{
				/* >>chng 06 feb 14, O+3 ff[ion] is very negative, as a result exp goes to +inf
				 * at very low tempertatures.  This is a error in the Nussbaumer & Storey fits to DR.
				 * do not use them is tefac = ff[ion] / t4 is very negative */
				if( tefac > -5. )
				{
					factor = (((aa[ion]*t4m1+bb[ion])*t4m1+cc[ion])*t4m1+dd[ion])* sexp(tefac);
					DielRecomRateCoef_LowT[ion] = ionbal.DielSupprs[1][ion]*fac2*
						MAX2(0.,factor );
				}
				else
				{
					DielRecomRateCoef_LowT[ion] = 0.;
				}
			}
			else if( ionbal.lg_guess_coef )
			{
				/* first guesses are those based on Nussbaumer & Storey and are from
				 * >>refer	ion	DR	Ali, B., Blum, R. D., Bumgardner, T. E., Cranmer, S. R., 
				 * >>refercon	Ferland, G. J., Haefner, R. I., & Tiede, G. P. 1991, PASP, 103, 1182*/
				if( (ff[ion] == 0.) && (ion <= 3)  )
				{
					/* 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};

					/* >>chng 01 jun 30, make GuessDiel array */
					DielRecomRateCoef_LowT[ion] = ionbal.DielSupprs[1][ion]*cludge[ion]*
						/* this is optional scale factor on set diel rec command 
						 * >>chng 06 feb 07, add Boltzmann factor to induce behavior
						 * like in Nussbaumer & Storey */
						ionbal.GuessDiel[ion] * sexp( 1e3 / phycon.te );
				}
				/* this implements the low-T kludge dielec command */
				else
				{
					/* assume the recom rate is the coefficient scanned off the steve option, times the charge
					 * before recombination raised to a power */
					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} 
					};

					/* these implement the guesses in 
					 * Kraemer, S.B., Ferland, G.J., & Gabel, J.R. 2004, ApJ, 604, 556  */
					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 */
						DielRecomRateCoef_LowT[ion] = fitcoef[nshell][0] +
							fitcoef[nshell][1]*(ion+1) + 
							fitcoef[nshell][2]*POW2(ion+1.) ;
						DielRecomRateCoef_LowT[ion] = 1e-10*pow(10.,DielRecomRateCoef_LowT[ion]);

					}
					else
					{
						/* this is guess for all other elements presented in 
						 * >>refer	all	DR	Kraemer, S.B., Ferland, G.J., & Gabel, J.R. 2004, ApJ, 604, 556 */
						DielRecomRateCoef_LowT[ion] = 3e-12*pow(10.,(double)(ion+1)*0.1);
					}
					/* >>chng 06 feb 02, add option to use mean of Badnell DR in place
					 * of these hacks */
					if( ionbal.lg_use_DR_Badnell_rate_coef_mean_ion )
						DielRecomRateCoef_LowT[ion] = ionbal.DR_Badnell_rate_coef_mean_ion[ion];
				}
				/* include optional noise here 
				 * >>chng 06 feb 07, move noise down to here so that use for both
				 * guesses of DR rates */
				DielRecomRateCoef_LowT[ion] *= RecNoise[nelem][ion];
			}
		}
		/* >>chng 05 dec 19, add option to use Badnell numbers */
		/* this is total old DR rates - may not use it */
		ionbal.DR_old_rate_coef[nelem][ion] = DielRecomRateCoef_HiT[ion] + DielRecomRateCoef_LowT[ion];

		/* set total DR rate - either Badnell if it exists or on with set diel recom badnell */
		if( ionbal.lgDR_recom_Badnell_use && ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] )
		{
			ionbal.DR_rate_coef_used[nelem][ion] = ionbal.DR_Badnell_rate_coef[nelem][ion];
		}
		else
		{
			ionbal.DR_rate_coef_used[nelem][ion] = ionbal.DR_old_rate_coef[nelem][ion];
		}

		/* sum of recombination rates [units s-1] for radiative, three body, charge transfer */
		ionbal.RateRecomTot[nelem][ion] = 
		  dense.eden* (
		  ionbal.RR_rate_coef_used[nelem][ion] + 
		  ionbal.DR_rate_coef_used[nelem][ion] +
		  ionbal.CotaRate[ion] ) + 
		  ChargeTransfer[ion];

		/* >>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.RR_rate_coef_used[nelem][ion]+DielRecomRateCoef_LowT[ion]+DielRecomRateCoef_HiT[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]);*/
		/* >>chng 05 sep 24, just print one long line - need info */
		limit = dense.IonHigh[nelem];
		for( i=0; i < limit; i++ )
		{
			fprintf( ioOut, "%10.2e", ionbal.RR_rate_coef_used[nelem][i] );
		}
		fprintf( ioOut, " radiative vs Z\n" );

		for( i=0; i < limit; i++ )
		{
			/* >>chng 06 jan 19, from div by eden to div by H0 - want units of cm3 s-1 but
			 * no single collider does this so not possible to get rate coefficient easily
			 * H0 is more appropriate than ne */
			fprintf( ioOut, "%10.2e", ChargeTransfer[i]/SDIV(dense.xIonDense[ipHYDROGEN][0]) );
		}
		fprintf( ioOut, " CT/n(H0)\n" );

		for( i=0; i < limit; i++ )
		{
			fprintf( ioOut, "%10.2e", ionbal.CotaRate[ion] );
		}
		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", gv.GrainChTrRate[nelem][i+1][i]/dense.eden );
		}
		fprintf( ioOut, " Grain vs Z /ne\n" );

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

		for( i=0; i < limit; i++ )
		{
			fprintf( ioOut, "%10.2e", DielRecomRateCoef_LowT[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", dense.eden*ionbal.CotaRate[ion] );
			}
			fprintf( ioOut, "\n" );

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

			fprintf( ioOut, " NS" );
			for( i=11; i < limit; i++ )
			{
				fprintf( ioOut, "%10.2e", DielRecomRateCoef_LowT[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	*/
	double 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 = 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 = te1[i];
					tfidle(FALSE);
					dense.IonLow[nelem] = 0;
					dense.IonHigh[nelem] = nelem+1;
					if( ConvBase(0) )
						fprintf(ioQQQ,"PROBLEM ConvBase returned error.\n");
					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 = 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 = te2[i];
					tfidle(FALSE);
					dense.IonLow[nelem] = 0;
					dense.IonHigh[nelem] = nelem+1;
					if( ConvBase(0) )
						fprintf(ioQQQ,"PROBLEM ConvBase returned error.\n");
					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*/

