/*lint -e817 pclint generates tons of bogus nelem < 0 warnings*/
/*MakeRecomb generate recombination coefficients for any species */
#include "cddefines.h"
#include "diffheav.h"
#include "showme.h"
#include "recomrate.h"
#include "nsbig.h"
#include "rcota.h"
#include "tepowers.h"
#include "phycon.h"
#include "cludon.h"
#include "gionrc.h"
#include "supprs.h"
#include "vrrfit.h"
#include "ionrange.h"
#include "ioreco.h"
#include "elementnames.h"
#include "destcrt.h"
#include "drfe.h"
#include "makerecomb.h"
#include "hcaseb.h"
#include "sexp.h"
#include "rrfit.h"

void MakeRecomb(double rec[], 
  double pl[], 
  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 NeBound, 
	  i, 
	  iupper, 
	  limit;
	double SaveRate[3][LIMELM], 
	  fac2, 
	  factor,
	  t4m1, 
	  tefac;
	float radrec;
	static double cludge[4]={3e-13,3e-12,1.5e-11,2.5e-11};

#	ifdef DEBUG_FUN
	fputs( "<+>MakeRecomb()\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( IonRange.IonLow[nelem] > 0 );
	assert( IonRange.IonLow[nelem] <= nelem+2 );

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

	/* zero out ionization, rec, heating save arrays */
	for( i=0; i < nelem+1; i++ )
	{
		RecomRate.RecombinRate[i][nelem] = 0.;
		SaveRate[0][i] = 0.;
		SaveRate[1][i] = 0.;
		SaveRate[2][i] = 0.;
	}

	/* H-like stages use results from hydrogenic series done by hydro... routines */
	/* ionHigh is on physical scale, for stripped carbon is 7, nelem would be 5*/
	if( IonRange.IonHigh[nelem] == (nelem+2) )
	{
		RecomRate.RecombinRate[IonRange.IonHigh[nelem]-2][nelem] = 
			(float)(DestCrt.create[nelem][nelem]);

		/* this was evaluated in hydro routines, and must not be changed here
		 * since we want the bidiag routine to get the same answer */
		SaveRate[0][IonRange.IonHigh[nelem]-2] = hcaseb.HRecRad[nelem] ;

		/* highest stage that need be considered in following loop */
		iupper = IonRange.IonHigh[nelem] - 2;
	}
	else
	{
		iupper = IonRange.IonHigh[nelem] - 1;
	}

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

	/* ionLow is on physical scale, so is 1 for atom */
	for( i=IonRange.IonLow[nelem]-1; i < iupper; i++ )
	{
		/* number of bound electrons of the ion */
		NeBound = nelem - i;
		if( vrrfit.lgVrrFit )
		{
			rrfit(nelem+1,nelem-i+1,phycon.te,&radrec);
		}
		else
		{
			if( (phycon.te > 1e3 || i > 0) || tlow[0] == 0. )
			{
				radrec = (float)(rec[i]*pow(phycon.te,pl[i]));
			}
			else
			{
				/* use special Pequignot and Aldrovandi fit for low temp */
				radrec = (float)(tlow[0]*pow(phycon.te,tlow[1]));
			}
		}
		RecomRate.RecombinRate[i][nelem] = (float)(phycon.eden*(radrec + 
		  rcota.CotaRate[i]));

		/* this is supposed to be the ground term of the recombination
		 * >>chng 97 feb 17, following two did not have eden,
		 * exposed bug when super low density models did not conserve energy */
		DiffHeav.DiffusHeavy[i][nelem] = (float)(phycon.eden*radrec/10.);
		DiffHeav.xLyaHeavy[i][nelem] = (float)(phycon.eden*radrec/10.);

		/* recombinatinon on grain surface */
		RecomRate.RecombinRate[i][nelem] += gionrc.GrnIonRec;

		/* rate coefficient for radiative plus three body, and
		 * recombination on grain surfaces */
		SaveRate[0][i] = radrec + rcota.CotaRate[i] + gionrc.GrnIonRec/phycon.eden;

		/* Burgess dielectronic recombination */
		if( phycon.te > ditcrt[i]*.1 )
		{
			SaveRate[1][i] = supprs.DielSupprs[0][i]/TePowers.te32*
			  DICOEF(0,i)*exp(-DITE(0,i)/phycon.te)*(1. + DICOEF(1,i)*
			  sexp(DITE(1,i)/phycon.te));
			RecomRate.RecombinRate[i][nelem] += (float)(SaveRate[1][i]*
			  phycon.eden);
		}

		/* Nussbaumer and Storey dielectronic recombination
		 * do not ever do it for rec from
		 * a closed shell, NeBound is number of bound electrons in ion */
		if( (NeBound != 2 && NeBound != 10) && NeBound != 18 )
		{
			tefac = ff[i]*t4m1;
			/* fac2 = 1e-14 * sqrte */
			if( ff[i] != 0. && fabs(tefac) < 30. )
			{
				factor = (((aa[i]*t4m1+bb[i])*t4m1+cc[i])*t4m1+dd[i])* exp(-tefac);
				SaveRate[2][i] = supprs.DielSupprs[1][i]*fac2*
					MAX2(0.,factor );
			}
			else if( (phycon.te > 100. && ff[i] == 0.) && i <= 3 )
			{
				SaveRate[2][i] = supprs.DielSupprs[1][i]*cludge[i]*
				  cludon.GuessDiel;
			}
			RecomRate.RecombinRate[i][nelem] += (float)(SaveRate[2][i]*
			  phycon.eden);
		}
	}

	/* for case of iron, add on Dima's fits to diel part, above was zero */
	if( nelem==25 )
	{
		float DieRate;
		/* add on special dima diel coef */
		for( i=IonRange.IonLow[nelem]-1; i < (IonRange.IonHigh[nelem]-1); i++ )
		{
			drfe(i+1,phycon.te,&DieRate);
			SaveRate[2][i] = DieRate;
			/*fprintf(ioQQQ,"%li %g %g\n", i, RecomRate.RecombinRate[i][NELEM-1], DieRate*phycon.eden );*/
			RecomRate.RecombinRate[i][nelem] += (float)(DieRate*phycon.eden);
		}
	}

	/* option to punch rec coef */
	if( lgioRecom )
	{
		/* print name of element */
		fprintf( ioRecom, 
			" %s recombination coefficients \n", 
			elementnames.chElementName[nelem] );

		limit = MIN2(11,IonRange.IonHigh[nelem]-1);
		for( i=0; i < limit; i++ )
		{
			fprintf( ioRecom, "%10.2e", SaveRate[0][i] );
		}
		fprintf( ioRecom, " radiative + 3body vs Z\n" );

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

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

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

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

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

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

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

	/*begin sanity check */
	limit = MIN2(nelem,IonRange.IonHigh[nelem]-1);
	for( i=IonRange.IonLow[nelem]-1; i < limit; i++ )
	{
		if( DiffHeav.DiffusHeavy[i][nelem] <= 0. )
		{
			fprintf( ioQQQ, 
				" MakeRecomb finds insane DiffusHeavy rate, nelem=%3ld rates follow.\n", 
			  nelem );
			fprintf(ioQQQ," current ion range is %li to %li \n",
				IonRange.IonLow[nelem] , IonRange.IonHigh[nelem] );
			fprintf(ioQQQ," specific value was i=%ld val= %g \n",
				i, DiffHeav.DiffusHeavy[i][nelem] );
			for( iupper=0; iupper < IonRange.IonHigh[nelem]; iupper++ )
			{
				fprintf( ioQQQ, "%9.1e", DiffHeav.DiffusHeavy[iupper][nelem] );
			}
			fprintf( ioQQQ, "\n" );
			ShowMe();
			puts( "[Stop in makerecomb]" );
			exit(1);
		}

		if( DiffHeav.xLyaHeavy[i][nelem] <= 0. )
		{
			fprintf( ioQQQ, " MakeRecomb finds insane xLyaHeavy rate, nelem=%3ld rates follow.\n", 
			  nelem );
			for( iupper=0; iupper < IonRange.IonHigh[nelem]; iupper++ )
			{
				fprintf( ioQQQ, "%9.1e", DiffHeav.xLyaHeavy[iupper][nelem] );
			}
			fprintf( ioQQQ, "\n" );
			ShowMe();
			puts( "[Stop in makerecomb]" );
			exit(1);
		}
		if( RecomRate.RecombinRate[i][nelem] <= 0. )
		{
			fprintf( ioQQQ, " MakeRecomb finds insane recombination rate, nelem=%3ld rates follow.\n", 
			  nelem );
			for( iupper=0; iupper < IonRange.IonHigh[nelem]; iupper++ )
			{
				fprintf( ioQQQ, "%9.1e", RecomRate.RecombinRate[iupper][nelem] );
			}
			fprintf( ioQQQ, "\n" );
			ShowMe();
			puts( "[Stop in makerecomb]" );
			exit(1);
		}
	}
	/*end sanity check */

#	ifdef DEBUG_FUN
	fputs( " <->MakeRecomb()\n", debug_fp );
#	endif
	return;
#undef	DITE
#undef	DICOEF
}
/*lint +e817 pclint generates tons of bogus nelem < 0 warnings*/