/* 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 */
/*HeRecom - do recomb coef for He, called by HeLike */
/*cross_section - calculates the photoionization cross_section for a given level and photon energy*/
/*radrecomb - calculates radiative recombination coefficients. */
/*RecomInt - Integral in milne relation.  Called by qg32.	*/
/*He_cross_section returns cross section (cm^-2), 
 * given EgammaRyd, the photon energy in Ryd,
 * ipLevel, the index of the level, 0 is ground, 3 within 2 3P,
 * nelem is charge, equal to 1 for Helium 
 * this is a wrapper for cross_section */
/*He_RRCoef_Te evaluated radiative recombination coef at some temperature */
/*HelikeCheckRecomb - called by SanityCheck to confirm that recombination coef are ok,
 * return value is relative error between new calculation of recom, and interp value */
/*Recomb_Seaton59 - find recombination for given n,using Seaton 59 approximation.
 * The following three are needed by Recomb_Seaton59:
 *     ExponentialInt
 *     X1Int
 *     X2Int	*/

#include "cddefines.h" 
#include "trace.h"
#include "rt.h"
#include "path.h"
#include "physconst.h" 
#include "hydro_bauman.h"
#include "iso.h"
#include "helike.h"
#include "helike_pcs.h"
#include "helike_recom.h"
#include "interpolate.h"
#include "dense.h"
#include "phycon.h"
#include "tfidle.h"
#include "opacity.h"
#include "ionbal.h"
#include "atmdat.h"

/*lint -e662 creation of  out of bound pointer */
/*lint -e661 creation of  out of bound pointer */

/* this will save log of radiative recombination rate coefficients at N_HE_TE_RECOMB temperatures.
 * there will be NumHeLevRecomb[nelem] levels saved in RRCoef[nelem][level][temp] */
static double ***RRCoef/*[LIMELM][NumHeLevRecomb[nelem]][N_HE_TE_RECOMB]*/;

static double **TotalRecomb;	/*[nelem][i]*/

/* the designation of the levels, chLevel[n][string]
static char **chLevel; */

/* the array of logs of temperatures at which RRCoef was defined */
static double TeRRCoef[N_HE_TE_RECOMB];

static double TempInterp( double* TempArray , double* ValueArray, long NumElements );

static double He_RRCoef_Te( long nelem , long n );

/* This computes that amount of the total recombination into an element (at a
 * specific temperature) that goes into levels having n greater than that otherwise
 * accounted for.	The reference is Seaton 59, main equation is 29. */
static double Recomb_Seaton59( long nelem, double temp, long n );

/* The three of these are used in the computation of Recomb_Seaton59	*/
static double ExponentialInt( double v );
static double X1Int( double u );
static double X2Int( double u );

static double cross_section(double EgammaRyd);
	
static double radrecomb(double temp, long nelem, long ipLo);

static double RecomInt(double EE);

static double kTRyd,EthRyd,Xn_S59; 
static long ipLev,globalZ;

/*He_cross_section returns cross section (cm^-2), 
 * given EgammaRyd, the photon energy in Ryd,
 * ipLevel, the index of the level, 0 is ground, 3 within 2 3P,
 * nelem is charge, equal to 1 for Helium 
 * this is a wrapper for cross_section */
double He_cross_section( double EgammaRyd , long ipLevel , long nelem )
{
	double cs;
	static double ThreshCS; 
	static long ipCurrentLevel=-1;

#define NUM_PREF_TCS 31

	/* These corrections ensure helium cross-sections agree with the highly accurate 
	 * values given by Hummer and Storey (1998).  The index is [ipLev - 1]	*/
	/* >>refer He	PCS	Hummer, D.G., & Storey, P.J. 1998, MNRAS 297, 1073	*/  
	double PreferredThres[NUM_PREF_TCS] = {
		/* n = 1 */
		7.394,
		/* n = 2 */
		5.480,9.253,15.98,15.98,15.98,13.48,
		/* n = 3 */
		8.025,14.49,28.52,18.48,18.13,26.99,
		/* n = 4 */
		10.77,20.38,41.59,36.70,35.75,-1.00,-1.00,41.75,
		/* n = 5, these are calculated using the Hummer and Storey 
		 * oscillator strength extrapolation procedure */
		13.54,26.58,55.71,-1.00,-1.00,-1.00,-1.00,-1.00,-1.00,58.07};

	globalZ = nelem;

	ipLev = ipLevel;

	EthRyd = iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLev];
	
	/* make sure this routine not called within collapsed high levels */
	ASSERT( ipLevel < iso.numLevels[ipHE_LIKE][nelem] - iso.nCollapsed[ipHE_LIKE][nelem] );

	/* this is a resolved level */
	cs = (1.e-18)*cross_section( EgammaRyd );

	/* Scale to preferred threshold values */
	if( nelem==ipHELIUM && ipLevel<NUM_PREF_TCS && !helike.lgSetBenjamin )
	{
		if( ipCurrentLevel!=ipLevel )
		{
			ThreshCS = cross_section( EthRyd );
			ipCurrentLevel = ipLevel;
		}

		if( PreferredThres[ipLevel] > 0. )
			cs *= PreferredThres[ipLevel]/ThreshCS;
	}

	return cs;
}

/*cross_section calculates the photoionization cross_section for a given level and photon energy*/
static double cross_section(double EgammaRyd)
{
	double E0[29] = {
	1.36E+01,2.01E+01,1.76E+01,3.34E+01,4.62E+01,6.94E+01,8.71E+01,1.13E+02,1.59E+02,2.27E+02,
	2.04E+02,2.74E+02,2.75E+02,3.38E+02,4.39E+02,4.17E+02,4.47E+02,5.18E+02,6.30E+02,6.27E+02,
	8.66E+02,7.67E+02,9.70E+02,9.66E+02,1.06E+03,1.25E+03,1.35E+03,1.43E+03,1.56E+03};
	double sigma[29] = {
	9.49E+02,3.20E+02,5.46E+02,2.85E+02,2.34E+02,1.52E+02,1.33E+02,1.04E+02,6.70E+01,4.00E+01,
	6.14E+01,4.04E+01,4.75E+01,3.65E+01,2.45E+01,3.14E+01,3.11E+01,2.59E+01,1.94E+01,2.18E+01,
	1.23E+01,1.76E+01,1.19E+01,1.31E+01,1.20E+01,9.05E+00,8.38E+00,8.06E+00,7.17E+00};
	double y_a[29] = {
	1.47E+00,7.39E+00,1.72E+01,2.16E+01,2.18E+01,2.63E+01,2.54E+01,2.66E+01,3.35E+01,5.32E+01,
	2.78E+01,3.57E+01,2.85E+01,3.25E+01,4.41E+01,3.16E+01,3.04E+01,3.28E+01,3.92E+01,3.45E+01,
	5.89E+01,3.88E+01,5.35E+01,4.83E+01,5.77E+01,6.79E+01,7.43E+01,7.91E+01,9.10E+01};
	double P[29] = {
	3.19E+00,2.92E+00,3.16E+00,2.62E+00,2.58E+00,2.32E+00,2.34E+00,2.26E+00,2.00E+00,1.68E+00,
	2.16E+00,1.92E+00,2.14E+00,2.00E+00,1.77E+00,2.04E+00,2.09E+00,2.02E+00,1.86E+00,2.00E+00,
	1.62E+00,1.93E+00,1.70E+00,1.79E+00,1.72E+00,1.61E+00,1.59E+00,1.58E+00,1.54E+00};
	double y_w[29] = 
	{2.039,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
	double yzero[29] =
	{0.4434,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
	double yone[29] =
	{2.136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
	
	/* These numbers are for the calculation of S cross-sections for n = 2 to n = 10. */
	/* indices are [n-2][s][a,b]	*/
	double heliumSparameters[9][2][2] = {
		{{0.0242024555,0.9702079},	{0.095921966,0.7337243}},

		{{0.0141220876,3.5176515},	{0.062991353,3.3995626}},

		{{0.0091014916,8.4545141},	{0.044113776,9.5657406}},
		
		{{0.0059539267,16.705745},	{0.032224774,20.933772}},
	
		{{0.0039485621,28.762086},	{0.024926753,38.464649}},
	
		{{0.0026053536,45.352333},	{0.020215452,63.043223}},
		
		{{0.0016632439,67.296132},	{0.016490176,96.296457}},
		
		{{0.0009867019,95.189826},	{0.013689704,138.82733}},
		
		{{0.0005155132,129.46774},	{0.011708145,191.73963}}};

	/* These are for the calculation of P cross-sections ...
	 * for n = 2 and n = 3 of triplets... */
	/* indices are [n-2][a,b]	*/
	double heTripP_params[2][2] = {
		{0.00489878046,3.0995085}, {0.01284809745,14.191429}};
	/* and for n=2 through n=9 of singlets. */
	/* indices are [n-2][a,b,c]	*/
	double heSingP_params[8][3] = {
		{-2.23352104,	-0.25979099,	-3.81078671},
		{-3.34344798,	-0.22891612,	-3.51168670},
		{-4.06915276,	-0.18554441,	-3.32255666},
		{-4.64178338,	-0.15390749,	-3.19447260},
		{-5.09796105,	-0.13059073,	-3.09491264},
		{-5.47775263,	-0.11232976,	-3.01326975},
		{-5.83522711,	-0.10848632,	-2.98479642},
		{-5.98242671,	-0.06820050,	-2.79942712}};

	double pcs,Egamma,lnE,a,b,c,y,F,x;
	long nelem = globalZ;

	Egamma = EgammaRyd * EVRYD;

	/* >>refer Helike	PCS	Verner, D.A., Ferland, G.J., Korista, K.T., & Yakovlev, D.G.
	 * >>refercon	1996a, ApJ 465,487	*/
	/* All ground state cross-sections calculated here.	*/
	if (ipLev == ipHe1s1S)
	{
		x = Egamma/E0[nelem-1] - yzero[nelem-1];
		y = sqrt(x*x + yone[nelem-1]*yone[nelem-1]);
		F = ((x-1)*(x-1)+y_w[nelem-1]*y_w[nelem-1])
			* pow(y,0.5*P[nelem-1]-5.5) * pow((1+sqrt(y/y_a[nelem-1])),-P[nelem-1]);
		pcs = sigma[nelem-1]*F;
		ASSERT( pcs > 0. );
	}
	/* Due to the highly non-hydrogenic nature of S orbitals, both singlet and triplet,
	 * a special formula needs to be applied to them.  The hydrogenic formula is 
	 * unsatisfactory since its errors get larger with increasing n.	*/
	
	/* Fix a better one for n = 2 as well?	*/
	else if ( ( L_(ipLev) == 0 ) && ( N_(ipLev) > 2 ) && ( nelem >= ipCARBON ) )
	{
		lnE = log(EgammaRyd * POW2(1.667*N_(ipLev)/nelem));
		pcs = POW2(5./(double)nelem) * exp(1.790166293 - 0.127439504 * POW2(lnE) - 1.660569112*lnE);
		pcs *= pow(N_(ipLev)/3.,1.365);
		ASSERT( pcs > 0. );
	}
	/* In addition to the above, the S orbitals of helium need to be handled separately.	*/
	else if ( ( nelem == ipHELIUM ) && ( L_(ipLev) == 0 ) )
	{
		if ( N_(ipLev)<=10 )
		{
			/* These are fits to TOPbase data	*/
			a = heliumSparameters[N_(ipLev)-2][S_(ipLev)][0];
			b = heliumSparameters[N_(ipLev)-2][S_(ipLev)][1];
		}
		/* These are fits to the parameters of the above fits!	*/
		else if ( S_(ipLev) == 0 )
		{
			a = 4.94 / pow((double)N_(ipLev),3.9);
			b = 0.147 * pow((double)N_(ipLev), 2.947);
		}
		else
		{
			a = 0.364/ pow((double)N_(ipLev), 1.49);
			b = 0.138 * pow((double)N_(ipLev), 3.146);
		}
			
		/* Put E^-3 tail on helium at energies beyond first resonance.	*/
		if ( EgammaRyd > 2.75 + EthRyd)
		{
			/* The value of the above fit at the point of the first resonance is approximately 
			 * equal to 1/(2.75^2 * b).  An E^-3 tail has a value of about 1/(2.75^3).  
			 * (EthRyd is much smaller than 2.75 even at n=2, just as parameter a is dwarved by b.)
			 * So the scale factor is 2.75/b.	*/
			 pcs = 2.75 * POW3(1./EgammaRyd) / b;
		}
		else
		{
			pcs = 1./(a+b*POW2(EgammaRyd));
		}

		/* TODO check that these are okay */
		if( N_(ipLev)>=6 )
		{
			double tweakle_dee[2][20] = {
			{0.9540,0.9485,0.9463,0.9449,0.9426,0.9505,0.9533,0.9562,0.9590,0.9617,
				0.9644,0.9670,0.9694,0.9717,0.9740,0.9761,0.9782,0.9802,0.9821,0.9839},
			{1.0418,1.0400,1.0378,1.0322,1.0265,1.0383,1.0439,1.0502,1.0569,1.0638,
				1.0709,1.0779,1.0850,1.0919,1.0987,1.1054,1.1120,1.1185,1.1248,1.1309}};

			if( N_(ipLev) <=25 )
			{
				pcs *= tweakle_dee[S_(ipLev)][N_(ipLev)-6];
			}
			else
			{
				pcs *= tweakle_dee[S_(ipLev)][19];
			}
		}

		ASSERT( pcs > 0. );
	}
	else if ( ( nelem == ipHELIUM ) && ( L_(ipLev) == 1 ) && (S_(ipLev) == 1) && ( N_(ipLev) <= 10 ) )
	{
		/* Fit to TOPbase data for triplet P...only n=2 and 3 are available.	*/
		if( N_(ipLev) <= 3 )
		{
			a = heTripP_params[N_(ipLev)-2][0];
			b = heTripP_params[N_(ipLev)-2][1];
			pcs = 1./(a+b*POW3(EgammaRyd));
		}
		/* Interpolate on a set of values obtained from Rob Bauman's Peach implementation */
		else
		{
			if( EgammaRyd - EthRyd < 1.0 )
			{
				pcs = pwlinear( PeachE, TripP[N_(ipLev)-4], 460, EgammaRyd - EthRyd );
				/* must convert from atomic units	*/
				pcs /= 0.035710595746;
			}
			else
			{
				/* need a tail */
				pcs = TripP[N_(ipLev)-4][458]/0.035710595746*pow(EgammaRyd - EthRyd, -2.6)*pow(0.98807,2.6);
			}
		}

		ASSERT( pcs > 0. );
	}
	/* D levels: interpolate on a set of values obtained from Rob Bauman's Peach implementation */
	else if ( ( nelem == ipHELIUM ) && ( L_(ipLev) == 2 ) && ( N_(ipLev) <= 10 ) )
	{
		double EelecRyd = EgammaRyd - EthRyd;
		EelecRyd = MAX2( 0., EelecRyd );

		if( S_(ipLev) == 1 )
		{
			if( EelecRyd < 1.0 )
			{
				pcs = pwlinear( PeachE, TripD[N_(ipLev)-3], 460, EelecRyd );
				/* must convert from atomic units	*/
				pcs /= 0.035710595746;
			}
			else
			{
				/* need a tail */
				pcs = TripD[N_(ipLev)-3][458]/0.035710595746*pow(EelecRyd, -3.65)*pow(0.98807,-3.65);
			}
		}
		else
		{
			if( EelecRyd < 1.0 )
			{
				pcs = pwlinear( PeachE, SingD[N_(ipLev)-3], 460, EelecRyd );
				/* must convert from atomic units	*/
				pcs /= 0.035710595746;
			}
			else
			{
				/* need a tail */
				pcs = SingD[N_(ipLev)-3][458]/0.035710595746*pow(EelecRyd, -3.75)*pow(0.98807,-3.75);
			}
		}

		ASSERT( pcs > 0. );
	}
	/* Fit to TOPbase data for singlet P.	*/
	else if ( ( nelem == ipHELIUM ) && ( L_(ipLev) == 1 ) && (S_(ipLev) == 0) && ( N_(ipLev) <= 7 ) )
	{	
		a = heSingP_params[N_(ipLev)-2][0];
		b = heSingP_params[N_(ipLev)-2][1];
		c = heSingP_params[N_(ipLev)-2][2];
		pcs = exp( a ) * pow( EgammaRyd, b*log(EgammaRyd)+c );

		ASSERT( pcs > 0. );
	}
	/* To everything else we apply a hydrogenic routine.	*/
	else 
	{
		double rel_photon_energy;
		
		/* >>chng 02 apr 24, more protection against calling with too small an energy */
		/* evaluating H-like photo cs at He energies, may be below threshold -
		 * prevent this from happening */
		rel_photon_energy = EgammaRyd / EthRyd;
		rel_photon_energy = MAX2( rel_photon_energy , 1. + FLT_EPSILON*2. );

		/* >>chng 02 apr 25, only one routine exposed, chooses which to evaluate */
		pcs = (1.e18)*H_photo_cs(rel_photon_energy , N_(ipLev), L_(ipLev), nelem);

		/* Rescale by ratio of real energy to hydrogenic energy.*/	
		pcs *= EthRyd * POW2( (double)N_(ipLev)/(double)nelem );

		/* This scale factor does a good job of making the SINGLET p threshold values smooth
		 * with increasing n */
		if( ( L_(ipLev) == 1 ) && (S_(ipLev) == 0) )
			pcs *= 1.06;

		/* This scale factor does a good job of making the TRIPLET p threshold values smooth
		 * with increasing n */
		/* TODO is this good enough? Need some justification! */
		if( ( L_(ipLev) == 1 ) && (S_(ipLev) == 1) )
			pcs *= 1.07;
			

		ASSERT( pcs > 0. );

	}

	return pcs;
}

/*radrecomb calculates radiative recombination coefficients. */
static double radrecomb(double temp, long nelem, long ipLo)
{
	double alpha,RecomIntegral=0.,b,E1,E2,step,OldRecomIntegral,TotChangeLastFive;
	double change[5] = {0.,0.,0.,0.,0.};
	
	/* Factors outside integral in Milne relation.	*/
	 b = MILNE_CONST * (2.*L_(ipLo)+1.) * (2.*S_(ipLo)+1.) * pow(temp,-1.5) / 4;
	/* kT in Rydbergs.	*/
	kTRyd = temp / TE1RYD;
	globalZ = nelem;
	ipLev = ipLo;

	/* Begin integration.	*/
	/* First define characteristic step as the minimum of kTRyd and 3EthRyd.	*/
	E1 = EthRyd;
	step = MIN( 0.25*kTRyd, 0.5*E1 );
	E2 = E1 + step;
	/* Perform initial integration, from threshold to threshold + step.	*/
	RecomIntegral = qg32( E1, E2, RecomInt);
	/* Repeat the integration, adding each new result to the total, 
	 * except that the step size is doubled every time, since values away from 
	 * threshold tend to fall off more slowly.	*/
	do
	{
		OldRecomIntegral = RecomIntegral;
		E1 = E2;
		step *= 1.25;
		E2 = E1 + step;
		RecomIntegral += qg32( E1, E2, RecomInt);
		change[4] = change[3];
		change[3] = change[2];
		change[2] = change[1];
		change[1] = change[0];
		change[0] = (RecomIntegral - OldRecomIntegral)/RecomIntegral;
		TotChangeLastFive = change[0] + change[1] + change[2] + change[3] + change[4];
	/* Continue integration until the upper limit exceeds 100kTRyd, an arbitrary
	 * point at which the integrand component exp(electron energy/kT) is very small,
	 * making neglible cross-sections at photon energies beyond that point,
	 * OR when the last five steps resulted in less than a 1 percent change.	*/
	} while ( ((E2-EthRyd) < 100.*kTRyd) && ( TotChangeLastFive > 0.0001) );

	/* Calculate recombination coefficient.	*/
	alpha = b * RecomIntegral;

	if ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) )
	{
		/* Split 2 trip P term.	*/
		alpha *= (2.*(ipLo-3L)+1.) / 9.;
	}

	if ( alpha < SMALLFLOAT )
		alpha = SMALLFLOAT;

	return alpha;
}

/*RecomInt, used in comput milne relation for he rec - the energy is photon Rydbergs.	*/
static double RecomInt(double ERyd)
{
	double x1, temp;

	/* Milne relation integrand	*/
	x1 = ERyd * ERyd * exp(-1.0 * ( ERyd - EthRyd ) / kTRyd);
	temp = He_cross_section( ERyd , ipLev, globalZ );
	x1 *= temp;
	
	return x1;
}

/* >>refer	He-like	RR	Seaton, M.J. 1959, MNRAS 119, 81S */
static double Recomb_Seaton59( long nelem, double temp, long n)
{
	double lambda = TE1RYD * nelem * nelem / temp;
	/* smallest x ever used here should be lowest Z, highest T, highest n...
	 * using helium, logt = 10., and n = 100, we get xmin = 1.5789E-9.	*/
	double x = lambda / n / n;
	double AlphaN;
	double SzeroOfX = 0.;
	double SoneOfX = 0.;
	double StwoOfX = 0.;
	double SnOfLambda = 0.;
	double lowerlimit, upperlimit, step;

	Xn_S59 = x;

	/* Equation 12	*/
	lowerlimit = x;
	step = 3. * x;
	upperlimit = lowerlimit + step;
	SzeroOfX = qg32( lowerlimit, upperlimit, ExponentialInt);

	do
	{
		lowerlimit = upperlimit;
		step *= 2;
		upperlimit = lowerlimit + step;
		SzeroOfX += qg32( lowerlimit, upperlimit, ExponentialInt);
	} while ( upperlimit < 20. );

	/* This must be placed inside integral...too big to be 
	 * handled separately.	
	SzeroOfX *= exp( x );	*/

	/* Equations 13 and 14 */
	lowerlimit = 0.;
	step = 0.5;
	upperlimit = lowerlimit + step;
	SoneOfX = qg32( lowerlimit, upperlimit, X1Int);
	StwoOfX = qg32( lowerlimit, upperlimit, X2Int);

	do
	{
		lowerlimit = upperlimit;
		step *= 2;
		upperlimit = lowerlimit + step;
		SoneOfX += qg32( lowerlimit, upperlimit, X1Int);
		SoneOfX += qg32( lowerlimit, upperlimit, X2Int);
	} while ( upperlimit < 200. );

	SoneOfX *= 0.1728 * pow( x, 1./3. );
	StwoOfX *= -0.0496 * pow( x, 2./3. );

	/* Equation 11	*/
	SnOfLambda = SzeroOfX + pow(1./lambda, 1./3.)*SoneOfX + pow(1./lambda, 2./3.)*StwoOfX;

	AlphaN = 5.197E-14 * nelem * pow(x, 1.5) * SnOfLambda;

	return AlphaN;

}

static double ExponentialInt( double v )
{
	double Integrand;

	Integrand = exp( -1. * v + Xn_S59) / v;

	return Integrand;
}

static double X1Int( double u )
{

	double Integrand;

	Integrand = pow(1./(u + 1.), 5./3.) * (u - 1.) * exp( -1. * Xn_S59 * u );

	return Integrand;
}

static double X2Int( double u )
{

	double Integrand;

	Integrand = pow(1./(u + 1.), 7./3.) * (u*u + 4./3.*u + 1.) * exp( -1. * Xn_S59 * u );

	return Integrand;
}

#undef DEBUG_LOC 


void HelikeRecombSetup( void )
{

	double RadRecombReturn;
	long int i, i1, i2, i3, i4, i5;
	long int ipLo, nelem, *NumHeLevRecomb;

#	define chLine_LENGTH 1000
	char chLine[chLine_LENGTH] , 
		/* this must be longer than chDataPath, set in path.h */
		chFilename[FILENAME_PATH_LENGTH_2];

	FILE *ioDATA;
	int lgEOL;

	/*double ***RRCoef[LIMELM][MaxLevels][N_HE_TE_RECOMB]*/
	if( (RRCoef = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	/*double **TotalRecomb[LIMELM][N_HE_TE_RECOMB]*/
	if( (TotalRecomb = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	/* The number of recombination coefficients to be read from file for each element.	*/
	if( (NumHeLevRecomb = (long*)MALLOC(sizeof(long)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
	{

		long int MaxLevels;

		if( nelem == ipHELIUM )
			NumHeLevRecomb[nelem]= ( ( 1 + HE_RREC_MAXN ) * HE_RREC_MAXN + 1 );
		else
			NumHeLevRecomb[nelem]= ( ( 1 + HE_LIKE_RREC_MAXN ) * HE_LIKE_RREC_MAXN + 1 );

		/* must always have at least NumHeLevRecomb[nelem] levels since that is number 
		 * that will be read in from he rec data file, but possible to require more */
		MaxLevels = MAX2( NumHeLevRecomb[nelem] , iso.numLevels[ipHE_LIKE][nelem] );

		/* always define this */
		/* >>chng 02 jan 24, RRCoef will be iso.numLevels[ipHE_LIKE][nelem] levels, no iso.numLevels,
		 * code will stop if more than this is requested */
		if( (RRCoef[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MaxLevels) ))==NULL )
			BadMalloc();

		/*double ***TotalRecomb[LIMELM][N_HE_TE_RECOMB]*/
		if( (TotalRecomb[nelem] = (double*)MALLOC(sizeof(double)*(unsigned)N_HE_TE_RECOMB ))==NULL )
			BadMalloc();
		
		for( ipLo=ipHe1s1S; ipLo < MaxLevels ;++ipLo )
		{
			if( (RRCoef[nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)N_HE_TE_RECOMB ))==NULL )
				BadMalloc();
		}
	}

	/******************************************************************/
	/**  Establish radiative recombination rate coefficients - RRC	***/
	/******************************************************************/
	for (i = 0; i < N_HE_TE_RECOMB; i++)
	{
		/* this is the vector of temperatures */
		TeRRCoef[i] = 0.5*(i);
	}

	/* This flag says we are not compiling the data file	*/
	if ( !helike.lgCompileRecomb)
	{
		/* check on path if file not here and path set */
		/* path was parsed in getset */
		if( lgDataPathSet == TRUE )
		{
			/*path set, so look only there */
			strcpy( chFilename , chDataPath );
			strcat( chFilename , "he_iso_recomb.dat" );
		}
		else
		{
			/* path not set, check local space only */
			strcpy( chFilename , "he_iso_recomb.dat" );
		}

		if( trace.lgTrace )
			fprintf( ioQQQ," HeCreate opening he_iso_recomb.dat:");

		/* Now try to read from file...*/
		if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
		{
			fprintf( ioQQQ, " HeCreate could not open he_iso_recomb.dat\n" );
			if( lgDataPathSet == TRUE )
				fprintf( ioQQQ, " even tried path\n" );

			if( lgDataPathSet == TRUE )
			{
				fprintf( ioQQQ, " HeCreate could not open he_iso_recomb.dat\n");
				fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
				fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
				fprintf( ioQQQ, " Defaulting to on-the-fly computation, ");
				fprintf( ioQQQ, " but the code runs much faster if you compile he_iso_recomb.dat!\n");
			}
			
			/* Do on the fly computation of R.R. Coef's instead.	*/
			for( nelem = ipHELIUM; nelem < LIMELM; nelem++ )
			{
				if( dense.lgElmtOn[nelem] )
				{
					/* Zero out the recombination sum array.	*/
					for(i = 0; i < N_HE_TE_RECOMB; i++)
					{
						TotalRecomb[nelem][i] = 0.;
					}
					
					/* NumHeLevRecomb[nelem] corresponds to n = 40 for He and 10 for ions, at present	*/
					/* There is no need to fill in values for collapsed levels, because we do not need to
					 * interpolate for a given temperature, just calculate it directly with a hydrogenic routine.	*/
					for( ipLo=ipHe1s1S; ipLo < iso.numLevels[ipHE_LIKE][nelem]-iso.nCollapsed[ipHE_LIKE][nelem]; ipLo++ )
					{
						EthRyd = iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo];

						/* loop over temperatures to produce array of recombination coefficients	*/
						for (i = 0; i < N_HE_TE_RECOMB; i++)
						{
							/* Store log of recombination coefficients, in N_HE_TE_RECOMB half dec steps */
							RadRecombReturn = radrecomb( pow( 10.,TeRRCoef[i] ) ,nelem,ipLo);
							TotalRecomb[nelem][i] += RadRecombReturn;
							RRCoef[nelem][ipLo][i] = log10(RadRecombReturn);
						}
					}
					for (i = 0; i < N_HE_TE_RECOMB; i++)
					{
						for( i1 = ( (nelem == ipHELIUM) ? (HE_RREC_MAXN + 1) : (HE_LIKE_RREC_MAXN + 1) ); i1<=SumUpToThisN; i1++ )
						{
							TotalRecomb[nelem][i] += Recomb_Seaton59( nelem, pow(10.,TeRRCoef[i]), i1 );
						}
						TotalRecomb[nelem][i] = log10( TotalRecomb[nelem][i] );
					}
				}
			}
		}
		/* Data file is present and readable...begin read.	*/
		else 
		{
			/* check that magic number is ok */
			if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			{
				fprintf( ioQQQ, " HeCreate could not read first line of he_iso_recomb.dat.\n");
				puts( "[Stop in HeCreate]" );
				cdEXIT(EXIT_FAILURE);
			}
			i = 1;
			i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			i4 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( i1 !=RECOMBMAGIC || i2 !=NumHeLevRecomb[ipHELIUM] || i3 !=NumHeLevRecomb[ipLITHIUM] || i4 !=N_HE_TE_RECOMB )
			{
				fprintf( ioQQQ, 
					" HeCreate: the version of he_iso_recomb.dat is not the current version.\n" );
				fprintf( ioQQQ, 
					" HeCreate: I expected to find the numbers  %i %li %li %i and got %li %li %li %li instead.\n" ,
					RECOMBMAGIC ,
					NumHeLevRecomb[ipHELIUM],
					NumHeLevRecomb[ipLITHIUM],
					N_HE_TE_RECOMB,
					i1 , i2 , i3, i4 );
				fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
				fprintf( ioQQQ, 
					" HeCreate: please recompile the data file with the COMPILE HE-LIKE command.\n" );
				puts( "[Stop in HeCreate]" );
				cdEXIT(EXIT_FAILURE);
			}

			i5 = 1;
			/* now read in the data */
			for( nelem = ipHELIUM; nelem < LIMELM; nelem++ )
			{
				for( ipLo=ipHe1s1S; ipLo <= NumHeLevRecomb[nelem]; ipLo++ )
				{
					i5++;
					/* get next line image */
					if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
					{
						fprintf( ioQQQ, " HeCreate could not read line %li of he_iso_recomb.dat.\n", i5);
						puts( "[Stop in HeCreate]" );
						cdEXIT(EXIT_FAILURE);
					}
					/* each line starts with element and level number */
					i3 = 1;
					i1 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
					i2 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
					/* check that these number are correct */
					if( i1!=nelem || i2!=ipLo )
					{
						fprintf( ioQQQ, " HeCreate detected insanity in he_iso_recomb.dat.\n");
						fprintf( ioQQQ, 
							" HeCreate: please recompile the data file with the COMPILE HE-LIKE command.\n" );
						puts( "[Stop in HeCreate]" );
						cdEXIT(EXIT_FAILURE);
					}

					/* loop over temperatures to produce array of recombination coefficients	*/
					for (i = 0; i < N_HE_TE_RECOMB; i++)
					{
						/* The last line for each element is the total recombination for each temp.	*/
						if( ipLo == NumHeLevRecomb[nelem] )
						{
							TotalRecomb[nelem][i] = FFmtRead(chLine,&i3,chLine_LENGTH,&lgEOL);
						}
						else
							RRCoef[nelem][ipLo][i] = FFmtRead(chLine,&i3,chLine_LENGTH,&lgEOL);

						if( lgEOL )
						{
							fprintf( ioQQQ, " HeCreate detected insanity in he_iso_recomb.dat.\n");
							fprintf( ioQQQ, 
								" HeCreate: please recompile the data file with the COMPILE HE-LIKE command.\n" );
							puts( "[Stop in HeCreate]" );
							cdEXIT(EXIT_FAILURE);
						}
					}
				}

				/* following loop only executed if we need more levels than are
				 * stored in the recom coef data set
				 * do not do collapsed levels since will use H-like recom there */
				if( nelem == ipHELIUM || dense.lgElmtOn[nelem] ) 
				{
					for( ipLo=NumHeLevRecomb[nelem]; ipLo < iso.numLevels[ipHE_LIKE][nelem]-iso.nCollapsed[ipHE_LIKE][nelem]; ipLo++ )
				    {
					     EthRyd = iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo];
						 
						 for (i = 0; i < N_HE_TE_RECOMB; i++)
					     {
						     /* Store log of recombination coefficients, in N_HE_TE_RECOMB half dec steps */
						     RRCoef[nelem][ipLo][i] = log10(radrecomb( pow( 10.,TeRRCoef[i] ) ,nelem,ipLo));
					     }
				     }
			     }
			}

			/* check that ending magic number is ok */
			if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			{
				fprintf( ioQQQ, " HeCreate could not read last line of he_iso_recomb.dat.\n");
				puts( "[Stop in HeCreate]" );
				cdEXIT(EXIT_FAILURE);
			}
			i = 1;
			i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			i4 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			
			if( i1 !=RECOMBMAGIC || i2 !=NumHeLevRecomb[ipHELIUM] || i3 !=NumHeLevRecomb[ipLITHIUM] || i4 !=N_HE_TE_RECOMB )
			{
				fprintf( ioQQQ, 
					" HeCreate: the version of he_iso_recomb.dat is not the current version.\n" );
				fprintf( ioQQQ, 
					" HeCreate: I expected to find the numbers  %i %li %li %i and got %li %li %li %li instead.\n" ,
					RECOMBMAGIC ,
					NumHeLevRecomb[ipHELIUM],
					NumHeLevRecomb[ipLITHIUM],
					N_HE_TE_RECOMB,
					i1 , i2 , i3, i4 );
				fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
				fprintf( ioQQQ, 
					" HeCreate: please recompile the data file with the COMPILE HE-LIKE command.\n" );
				puts( "[Stop in HeCreate]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* close the data file */
			fclose( ioDATA );
		}
	}
	/* We are compiling the he_iso_recomb.dat data file.	*/
	else
	{
		/* option to create table of recombination coefficients,
		 * executed with the compile he-like command */
		FILE *ioRECOMB;
		
		/*RECOMBMAGIC the magic number for the table of recombination coefficients */
		/*NumHeLevRecomb[nelem] the number of levels that will be done */
		/* create recombination coefficients  */
		if( ( ioRECOMB = fopen( "he_iso_recomb.dat" , "w" ) ) == NULL )
		{
			fprintf( ioQQQ, " HeCreate could not open he_iso_recomb.dat for writing.\n" );
			puts( "[Stop in HeCreate]" );
			cdEXIT(EXIT_FAILURE);
		}
		fprintf(ioRECOMB,"%i\t%li\t%li\t%i\tHe-like isoelectronic sequence recomb data, created by COMPILE HE-LIKE command, with %li He levels, %li ion levels, and %i temperatures.\n",
			RECOMBMAGIC ,
			NumHeLevRecomb[ipHELIUM],
			NumHeLevRecomb[ipLITHIUM],
			N_HE_TE_RECOMB,
			NumHeLevRecomb[ipHELIUM],
			NumHeLevRecomb[ipLITHIUM],
			N_HE_TE_RECOMB );

		for( nelem = ipHELIUM; nelem < LIMELM; nelem++ )
		{
			/* this must pass since compile he-like command reset numLevels to the macro */
			ASSERT( NumHeLevRecomb[nelem] <= iso.numLevels[ipHE_LIKE][nelem] );

			/* Zero out the recombination sum array.	*/
			for(i = 0; i < N_HE_TE_RECOMB; i++)
			{
				TotalRecomb[nelem][i] = 0.;
			}
			
			for( ipLo=ipHe1s1S; ipLo < NumHeLevRecomb[nelem] ; ipLo++ )
			{
				EthRyd = iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo];
				
				fprintf(ioRECOMB, "%li\t%li", nelem, ipLo );
				/* loop over temperatures to produce array of recombination coefficients	*/
				for (i = 0; i < N_HE_TE_RECOMB; i++)
				{
					/* Store log of recombination coefficients, in N_HE_TE_RECOMB half dec steps */
					RadRecombReturn = radrecomb( pow( 10.,TeRRCoef[i] ) ,nelem,ipLo);
					TotalRecomb[nelem][i] += RadRecombReturn;
					RRCoef[nelem][ipLo][i] = log10(RadRecombReturn);
					fprintf(ioRECOMB, "\t%f", RRCoef[nelem][ipLo][i] );
				}
				fprintf(ioRECOMB, "\n" );
			}
			
			/* Store one additional line in he_iso_recomb.dat that gives the total recombination,
			 * as computed by the sum so far, plus levels up to SumUpToThisN using Seaton 59, 
			 * for each element and each temperature.	*/
			fprintf(ioRECOMB, "%li\t%li", nelem, NumHeLevRecomb[nelem] );
			for (i = 0; i < N_HE_TE_RECOMB; i++)
			{
				for( i1 = ( (nelem == ipHELIUM) ? (HE_RREC_MAXN + 1) : (HE_LIKE_RREC_MAXN + 1) ); i1<=SumUpToThisN; i1++ )
				{
					TotalRecomb[nelem][i] += Recomb_Seaton59( nelem, pow(10.,TeRRCoef[i]), i1 );
				}
				fprintf(ioRECOMB, "\t%f", log10( TotalRecomb[nelem][i] ) );
			}
			fprintf(ioRECOMB, "\n" );
		}
		/* end the file with the same information */
		fprintf(ioRECOMB,"%i\t%li\t%li\t%i\tHe-like isoelectronic sequence recomb data, created by COMPILE HE-LIKE command, with %li He levels, %li ion levels, and %i temperatures.\n",
			RECOMBMAGIC ,
			NumHeLevRecomb[ipHELIUM],
			NumHeLevRecomb[ipLITHIUM],
			N_HE_TE_RECOMB,
			NumHeLevRecomb[ipHELIUM],
			NumHeLevRecomb[ipLITHIUM],
			N_HE_TE_RECOMB );

		fclose( ioRECOMB );

		fprintf( ioQQQ, "HeCreate: compilation complete, he_iso_recomb.dat created.\n" );
		puts( "[Stop in HeCreate]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* End of recombination setup.	*/

	free( NumHeLevRecomb );

	return;
}

/*===================================================================================*/
/* HeRecom do recomb coef for He, called by HeLike */
void HeRecom( 
			/* nelem on the c scale, He is 1 */
			long nelem )
{
	static double ditcrt[LIMELM]= {0.,5e4,1.20E+04,1.10E+04,4.40E+05,7.00E+05,1.00E+06,1.50E+06,
	  1.50E+06,3.80E+06,3.40E+04,4.00E+06,5.50E+06,6.30E+06,5.00E+06,9.00E+06,1.40E+07,
	  1.40E+07,2.00E+05,2.00E+05,2.00E+05,1.00E+05,7.00E+04,4.00E+04,6.00E+06,1.00E+11,
	  1.00E+20,1.00E+20,1.00E+20,1.00E+20};

	static double dicoef0[LIMELM]= {0.,1.9e-3,6.15E-03,1.62E-03,4.78E-02,3.22E-02,4.61E-02,6.23E-02,
	  2.52E-01,1.44E-01,1.80E-01,1.44E-01,3.40E+00,1.89E-01,4.02E-01,2.41E-01,1.10E+01,2.97E-01,
	  5.49E-01,3.55E-01,2.30E-01,1.40E-01,1.10E-01,4.10E-02,7.47E-01,5.19E-01,5.75E-01,5.75E-01,
	  5.75E-01,5.75E-01};

	static double dite0[LIMELM]= {0.,4.7e5,1.41E+05,8.19E+04,3.44E+06,4.06E+06,5.44E+06,7.01E+06,1.40E+07,
	  1.50E+07,1.30E+07,1.50E+07,1.70E+07,1.99E+07,2.41E+07,2.54E+07,2.80E+07,3.13E+07,3.65E+07,
	  3.78E+07,1.09E+06,9.62E+05,7.23E+05,4.23E+05,5.84E+07,6.01E+07,6.81E+07,6.81E+07,6.81E+07,
	  6.81E+07};

	static double dicoef1[LIMELM]= {0.,0.3,5.88E-02,3.43E-01,3.62E-01,3.15E-01,3.09E-01,3.04E-01,
	  3.04E-01,2.96E-01,0.00E+00,2.91E-01,0.00E+00,2.86E-01,2.98E-01,2.81E-01,0.00E+00,2.77E-01,
	  2.92E-01,2.75E-01,3.60E+00,4.90E+00,1.60E+00,4.20E+00,2.84E-01,2.79E-01,0.00E+00,2.86E-01,
	  2.86E-01,2.86E-01};

	static double dite1[LIMELM]= {0.,9.4e4,1.41E+05,1.59E+05,5.87E+05,8.31E+05,1.14E+06,1.47E+06,
	  1.80E+06,2.24E+06,0.00E+00,3.09E+06,0.00E+00,4.14E+06,4.67E+06,5.30E+06,0.00E+00,6.54E+06,
	  7.25E+06,7.68E+06,1.10E+07,8.34E+06,1.01E+07,1.07E+07,1.17E+07,9.97E+06,0.00E+00,9.08E+06,
	  9.08E+06,9.08E+06};

	long ipFirstCollapsed, LastN=0L, ThisN=1L, ipLevel; 
	double RecUnitVolTotal, topoff, TotMinusExplicitResolved,
		TotRRThisN=0., TotRRLastN=0., DielRec, t;
	static double RecExplictLevels[LIMELM], TotalRadRecomb[LIMELM], RecCollapsed[LIMELM];

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

	/* evaluate recombination escape probability for all levels */

	/* define radiative recombination rates for all levels */ 
	/* this will be the sum of all levels explicitly included in the model atom */
	RecExplictLevels[nelem] = 0.;
	
	/* this is the number of resolved levels, so first collapsed level is [ipFirstCollapsed] */
	ipFirstCollapsed = iso.numLevels[ipHE_LIKE][nelem]-iso.nCollapsed[ipHE_LIKE][nelem]; 

	if( helike.lgSetBenjamin && (nelem == ipHELIUM) )
        /* total topoff to Benjamin numbers */
		TotalRadRecomb[nelem] = 4.27E-13 * pow(10000./phycon.te, 0.678);
	else
	{
		TotalRadRecomb[nelem] = He_RRCoef_Te( nelem , ipFirstCollapsed );
		if( helike.lgRandErrGen )
		{
			TotalRadRecomb[nelem] *= 
				helike.ErrorFactor[nelem][iso.numLevels[ipHE_LIKE][nelem]][iso.numLevels[ipHE_LIKE][nelem]][ipRad];
		}
	}

	for( ipLevel=0; ipLevel<ipFirstCollapsed; ++ipLevel )
	{
		double BenDirectRadRecHe[31] = {
		1.54E-13,1.49E-14,5.55E-15,6.23E-15,1.87E-14,3.12E-14,1.26E-14,3.72E-15,1.62E-15,
		1.95E-14,1.33E-14,4.31E-15,5.01E-15,1.50E-15,7.00E-16,9.04E-15,8.29E-15,2.67E-15,
		4.15E-15,1.38E-15,2.43E-15,7.51E-16,3.66E-16,4.84E-15,4.99E-15,1.60E-15,3.65E-15,
		1.22E-15,1.34E-15,4.46E-16,1.35E-15};

		double BenDirectRadExpHe[31] = {
		0.486,0.381,0.451,0.639,0.639,0.639,0.695,0.344,0.444,0.632,0.868,0.872,0.700,0.328,0.445,
		0.632,0.867,0.871,1.046,1.046,0.708,0.327,0.454,0.636,0.870,0.873,1.048,1.048,1.183,1.183,0.718};

		double BenIndirectRadRecHe[31] = {
		0.      ,7.07E-15,2.25E-15,1.12E-15,3.36E-15,5.61E-15,3.11E-15,1.51E-15,6.09E-16,
		3.63E-15,8.19E-15,2.64E-15,1.30E-15,5.16E-16,2.44E-16,1.62E-15,5.19E-15,1.64E-15,
		6.94E-15,2.31E-15,6.50E-16,2.09E-16,1.16E-16,8.63E-16,3.21E-15,9.86E-16,5.83E-15,
		1.94E-15,8.39E-15,2.80E-15,3.82E-16};

		double BenIndirectRadExpHe[31] = {
		0.000,0.675,0.766,0.901,0.901,0.901,0.879,0.696,0.782,0.873,1.097,1.114,0.876,0.701,0.783,
		0.845,1.078,1.106,1.271,1.271,0.854,0.711,0.786,0.813,1.055,1.096,1.267,1.267,1.375,1.375,0.832};

		double BenIndirectDensExpHe[31] = {
		0.000,-0.001,0.000,0.009,0.009,0.009,0.008,-0.001,0.000,0.010,0.002,0.003,0.008,
		-0.001,0.000,0.010,0.002,0.003,-0.001,-0.001,0.008,-0.001,0.000,0.010,0.002,0.003,
		-0.001,-0.001,-0.007,-0.007,0.008};

		if( helike.lgSetBenjamin && (nelem == ipHELIUM) && ( ipLevel <= 30 ) 
			&& ( phycon.te < 1000000.1f ) && ( phycon.te > 99.9f ) )
		{
			iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecRad] = 
				BenDirectRadRecHe[ipLevel] * pow( 10000./phycon.te, BenDirectRadExpHe[ipLevel] ) +
				BenIndirectRadRecHe[ipLevel] * pow( 10000./phycon.te, BenIndirectRadExpHe[ipLevel] ) * 
				pow( dense.eden/100., BenIndirectDensExpHe[ipLevel] );
		}
		else
		{
			/* this is radiative recombination rate coefficient */
			double RadRec = He_RRCoef_Te( nelem , ipLevel );
			iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecRad] = RadRec;

			
			if( helike.lgRandErrGen )
			{
				iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecRad] *=
					helike.ErrorFactor[nelem][ iso.numLevels[ipHE_LIKE][nelem] ][ipLevel][ipRad];
			}
		}

		RecExplictLevels[nelem] += iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecRad];
		ASSERT( iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecRad] > 0. );
	}

	/* do not include ground in this sum	*/
	RecExplictLevels[nelem] -= iso.RadRecomb[ipHE_LIKE][nelem][ipHe1s1S][ipRecRad];	

	/**************************************************/
	/***  Add on recombination to collapsed levels  ***/
	/**************************************************/
	RecCollapsed[nelem] = 0.;
	for( ipLevel=ipFirstCollapsed; ipLevel<iso.numLevels[ipHE_LIKE][nelem]; ++ipLevel )
	{
		double RadRec = atmdat_H_rad_rec(nelem, N_(ipLevel), phycon.te);

		/* this is radiative recombination rate coefficient */
		iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecRad] = RadRec;

		/* This kills recombination into the collapsed level so that the forced
		 * statistical weighting can be bypassed */
		if( helike.lgTopoff==FALSE )
		{
			iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecRad] *= 1E-10;
		}

		if( helike.lgRandErrGen )
		{
			iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecRad] *= 
				helike.ErrorFactor[nelem][ iso.numLevels[ipHE_LIKE][nelem] ][ipLevel][ipRad];
		}

		RecCollapsed[nelem] += iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecRad];

		ASSERT( iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecRad] > 0. );
	}

	for( ipLevel = ipHe1s1S; ipLevel<iso.numLevels[ipHE_LIKE][nelem]; ipLevel++ )
	{

		if( N_(ipLevel) == ThisN )
		{
			TotRRThisN += iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecRad];
		}
		else
		{
			ASSERT( (TotRRThisN<TotRRLastN) || (ThisN<=2L) || (phycon.te>3E7) || (nelem!=ipHELIUM) );
			LastN = ThisN;
			ThisN = N_(ipLevel);
			TotRRLastN = TotRRThisN;
			TotRRThisN = iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecRad];

			{
				/* Print the sum of recombination coefficients per n at current temp.	*/
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				static long RUNONCE = FALSE;
				
				if( !RUNONCE && DEBUG_LOC )
				{
					static long FIRSTTIME = TRUE;

					if( FIRSTTIME )
					{
						fprintf( ioQQQ,"Sum of Radiative Recombination at current temp, T=%.2f\n", phycon.te);
						FIRSTTIME= FALSE;
					}
					
					fprintf( ioQQQ,"%li\t%.2e\n",LastN,TotRRLastN );
				}
				RUNONCE = TRUE;
			}
		}
	}


	/* this is case B recombination, sum without the ground included */
	iso.RadRec_caseB[ipHE_LIKE][nelem] = TotalRadRecomb[nelem] - iso.RadRecomb[ipHE_LIKE][nelem][ipHe1s1S][ipRecRad];
	ASSERT( iso.RadRec_caseB[ipHE_LIKE][nelem] > 0.);

	/* this is used for radiative recombination cooling in coolr, this does not include
	 * optical depth effects - may be reevaluated in level solver with opt dep included */
	ionbal.RadRecomRateCoef[nelem][nelem-1] = TotalRadRecomb[nelem];

	/* >>chng 02 jan 07, do not add three body, complete model atom will automatically include it
	 * problem is that most elements will not go to high n */
	/* add three-body recombination, convert to per unit vol 
	RecUnitVolTotal = dense.eden*(TotalRadRecomb[nelem] + ionbal.CotaRate[nelem-1]);*/
	/* convert to per unit vol */
	RecUnitVolTotal = dense.eden*TotalRadRecomb[nelem];

	/**********************************************************************/
	/***  Add topoff (excess) recombination	to top level.  This is only	***/
	/***  done if not using Benjamin command and atom is not full size. ***/
	/**********************************************************************/
	if( !helike.lgSetBenjamin && !helike.lgFullSize[nelem] )
	{
		/* at this point we have RecExplictLevels[nelem], the sum of radiative recombination 
		 * to all levels explicitly included in the model atom, and RecUnitVolTotal, the total 
		 * recombination rate.  The difference is the amount of "topoff" we will need to do */
		TotMinusExplicitResolved = RecUnitVolTotal/dense.eden - 
			iso.RadRecomb[ipHE_LIKE][nelem][0][ipRecRad] - RecExplictLevels[nelem] ;
		
		topoff = TotMinusExplicitResolved - RecCollapsed[nelem];

		/* the atmdat_rad_rec fits are too small at high temperatures, so this atom is
		 * better than the topoff.  Only a problem for helium itself, at high temperatures,
		 * when physics disabled with the atom he-like benjamin command.
		 * complain if the negative topoff is not for this case */
		/* >>chng 02 jul 30, do not print if benjamin is set - this causes several known problems */
		if( topoff < 0. && (nelem!=ipHELIUM || phycon.te < 1e5 ) )
		{
			fprintf(ioQQQ," PROBLEM  negative topoff for %li, rel error was %.2e\n", 
				nelem, topoff/(RecUnitVolTotal/dense.eden) );
		}
		
		if( helike.lgTopoff==FALSE )
			topoff *= 1E-20;
		
		/* >>chng 02 feb 24...now always have at least one collapsed level...so
		 * this bit greatly simplifies.	*/
		iso.RadRecomb[ipHE_LIKE][nelem][iso.numLevels[ipHE_LIKE][nelem]-1][ipRecRad] += 
			MAX2(0., topoff );

		/*for( ipLevel=QuantumNumbers2Index[nelem][16][0][1]; ipLevel<ipFirstCollapsed; ipLevel++ )
		{
			iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecRad] += 
				(2.*L_(ipLevel)+1)*(2.*S_(ipLevel)+1)*topoff/1024;
		}*/
	}
		

	/***********************************************/
	/***  Add on Dielectronic recombination	- DR ***/
	/***********************************************/
	DielRec = 0.;
	if( helike.lgDielRecom > 0 )
	{
		if( nelem==ipIRON )
		{
			/* for case of iron, add on Dima's fits to diel part, above was zero,
			 * this requests He-like Fe */
			DielRec = atmdat_dielrec_fe(25,phycon.te);
		}
		else if( phycon.te > (float)(ditcrt[nelem]*0.1) )
		{
			/* the suppression factor is on a charge scale, so rec to neutral
			 * is always 0, hence the nelem-1 */
			if( nelem == ipHELIUM )
			{
				/* this is the dr rate that was in the original He singlets coding, and has been
				 * part of cloudy since the beginning - based on Aldrovandi & Pequignot */
				/* >>refer He	DR	Aldrovandi, S.M.V., & Pequignot, D. 1973, A&A 25, 137	*/
				if( phycon.te > 5e4 )
				{
					t = (float)(1.9e-3/phycon.te32*sexp(4.7e5/phycon.te)*(1. + 0.3*sexp(9.4e4/phycon.te)));
				}
				else
					t = 0;
				DielRec = t;
			}
			else
			{
				if( helike.lgDielRecom == 1 && nelem >= 12 )
				{
					long i;
					double EoverT;
					double ParamsE[28][3] ={
						{5.966E+01,	7.165E+01,	7.686E+01}, {1.066E+02,	1.289E+02,	1.398E+02},
						{1.672E+02,	2.032E+02,	2.223E+02},	{2.415E+02,	2.948E+02,	3.248E+02},
						{3.296E+02,	4.038E+02,	4.476E+02},	{4.314E+02,	5.303E+02,	5.909E+02},
						{5.471E+02,	6.744E+02,	7.548E+02},	{6.767E+02,	8.362E+02,	9.397E+02},
						{8.202E+02,	1.016E+03,	1.146E+03},	{9.749E+02,	1.200E+03,	1.353E+03},
						{1.149E+03,	1.428E+03,	1.622E+03},	{1.335E+03,	1.662E+03,	1.892E+03},
						{1.533E+03,	1.913E+03,	2.184E+03},	{1.749E+03,	2.195E+03,	2.520E+03},
						{1.974E+03,	2.470E+03,	2.833E+03},	{2.218E+03,	2.799E+03,	3.228E+03},
						{2.470E+03,	3.099E+03,	3.570E+03},	{2.741E+03,	3.469E+03,	4.014E+03},
						{3.023E+03,	3.801E+03,	4.396E+03},	{3.320E+03,	4.180E+03,	4.843E+03},
						{3.631E+03,	4.577E+03,	5.312E+03},	{3.957E+03,	4.992E+03,	5.803E+03},
						{4.297E+03,	5.426E+03,	6.317E+03},	{4.641E+03,	5.797E+03,	6.740E+03},
						{5.018E+03,	6.348E+03,	7.414E+03},	{5.399E+03,	6.845E+03,	8.000E+03},
						{5.796E+03,	7.345E+03,	8.601E+03},	{6.207E+03,	7.871E+03,	9.230E+03}	};
					double ParamsC[28][3] ={
						{6.774E-13,	2.089E-12,	1.116E-12}, {6.774E-13,	2.089E-12,	1.116E-12},
						{6.774E-13,	2.089E-12,	1.116E-12}, {6.774E-13,	2.089E-12,	1.116E-12},
						{6.774E-13,	2.089E-12,	1.116E-12}, {6.774E-13,	2.089E-12,	1.116E-12},
						{6.774E-13,	2.089E-12,	1.116E-12}, {6.774E-13,	2.089E-12,	1.116E-12},
						{6.774E-13,	2.089E-12,	1.116E-12}, {6.774E-13,	2.089E-12,	1.116E-12},
						{7.301E-13,	2.039E-12,	9.005E-13}, {7.827E-13,	1.988E-12,	6.850E-13},
						{8.197E-13,	1.867E-12,	5.788E-13}, {8.567E-13,	1.746E-12,	4.726E-13},
						{8.762E-13,	1.609E-12,	4.190E-13}, {8.956E-13,	1.471E-12,	3.653E-13},
						{8.942E-13,	1.337E-12,	3.374E-13}, {8.927E-13,	1.202E-12,	3.095E-13},
						{8.584E-13,	1.100E-12,	3.059E-13}, {8.241E-13,	9.970E-13,	3.022E-13},
						{7.899E-13,	8.945E-13,	2.986E-13}, {7.556E-13,	7.920E-13,	2.950E-13},
						{7.213E-13,	6.895E-13,	2.913E-13}, {6.870E-13,	5.870E-13,	2.877E-13},
						{6.548E-13,	5.686E-13,	2.348E-13}, {6.226E-13,	5.502E-13,	1.818E-13},
						{6.226E-13,	5.502E-13,	1.818E-13}, {6.226E-13,	5.502E-13,	1.818E-13}  };

					/* Use Gu data	*/
					/* TODO	2	insert refer here.	*/
					DielRec = 0.;
					for( i = 1; i <= 3; i++ )
					{
						ASSERT( nelem <= 29 );
						EoverT = ParamsE[nelem-2][i-1]*TE1RYD/phycon.te/EVRYD;
						DielRec += ParamsC[nelem-2][i-1]*pow(EoverT,1.5)*sexp(EoverT);
					}
				}
				else
				{
					/* Use old treatment	*/
					/* TODO	2	put Burgess reference here	*/
					DielRec = ionbal.DielSupprs[0][nelem-1]/phycon.te32*dicoef0[nelem]*
						exp(-dite0[nelem]/phycon.te)*(1. + dicoef1[nelem]*sexp(dite1[nelem]/phycon.te));
				}
			}
		}
	}

	iso.RadRecomb[ipHE_LIKE][nelem][iso.numLevels[ipHE_LIKE][nelem]-1][ipRecRad] += 
			MAX2(0., DielRec );

	RecUnitVolTotal += DielRec*dense.eden;
	/* fprintf( ioQQQ,"nelem\t%li\DielRec\t%.3e\n",nelem,DielRec); */

	/* option for case b conditions, kill ground state recombination */
	if( opac.lgCaseB )
	{
		iso.RadRecomb[ipHE_LIKE][nelem][ipHe1s1S][ipRecEsc] = 1e-10;
		iso.RadRecomb[ipHE_LIKE][nelem][ipHe1s1S][ipRecNetEsc] = 1e-10;
	}
	else
	{
		double temp = RT_recom_effic(iso.ipIsoLevNIonCon[ipHE_LIKE][nelem][ipHe1s1S]);
		iso.RadRecomb[ipHE_LIKE][nelem][ipHe1s1S][ipRecEsc] = 
			MAX2( temp , opac.otsmin);
		iso.RadRecomb[ipHE_LIKE][nelem][ipHe1s1S][ipRecNetEsc] = 
			MIN2(1.,iso.RadRecomb[ipHE_LIKE][nelem][ipHe1s1S][ipRecEsc]+
			(1.-iso.RadRecomb[ipHE_LIKE][nelem][ipHe1s1S][ipRecEsc])*
			iso.ConOpacRatio[ipHE_LIKE][nelem][ipHe1s1S]);
	}
	/* total continuum effective escape prob for all levels  */
	/* RadRec_effec is total effective radiative recombination */
	iso.RadRec_effec[ipHE_LIKE][nelem] = iso.RadRecomb[ipHE_LIKE][nelem][ipHe1s1S][ipRecRad]*
		  iso.RadRecomb[ipHE_LIKE][nelem][ipHe1s1S][ipRecNetEsc];

	/**************************************************************/
	/***  Stuff escape probabilities, and effective rad recomb  ***/
	/**************************************************************/
	for( ipLevel=ipHe1s1S+1; ipLevel<iso.numLevels[ipHE_LIKE][nelem]; ++ipLevel )
	{
		/* this is escape probability */
		iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecEsc] = 
			RT_recom_effic(iso.ipIsoLevNIonCon[ipHE_LIKE][nelem][ipLevel]);

		/* correct above value for possible otsmin value
		 * otsmin set to zero in zerologic, set to 1 with 
		 * the NO ON THE SPOT command */
		iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecEsc] = 
			MAX2(iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecEsc], opac.otsmin);

		/* net escape prob includes dest by background opacity */
		iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecNetEsc] = 
			MIN2(1.,iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecEsc]+
		  (1.-iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecEsc])*iso.ConOpacRatio[ipHE_LIKE][nelem][ipLevel]);

		ASSERT( iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecEsc] >= 0. );
		ASSERT( iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecNetEsc] >= 0. );

		/* sum of all effective rad rec */
		iso.RadRec_effec[ipHE_LIKE][nelem] += iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecRad]*
		  iso.RadRecomb[ipHE_LIKE][nelem][ipLevel][ipRecNetEsc];
	}

	return;

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

/* TempInterp - interpolate on an array
 * TODO use a canned interpolation routine, no need for special one here */
static double TempInterp( double* TempArray , double* ValueArray, long NumElements )
{
	static long int ipTe=-1;
	double rate = 0.;
	long i;

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

	ASSERT( fabs( 1. - (double)phycon.alogte/log10((double)phycon.te) ) < 0.0001 );

	if( ipTe < 0 )
	{
		/* te totally unknown */
		if( ( phycon.alogte < (float)TempArray[0] ) || 
			( phycon.alogte > (float)TempArray[NumElements-1] ) )
		{
			fprintf(ioQQQ," TempInterp called with te out of bounds \n");
			puts( "[Stop in TempInterp]" );
			cdEXIT(EXIT_FAILURE);
		}
		/* now search for temperature */
		for( i=0; i<NumElements-1; ++i )
		{
			if( ( phycon.alogte > (float)TempArray[i] ) && 
				( phycon.alogte <= (float)TempArray[i+1] ) )
			{
				/* found the temperature - use it */
				ipTe = i;
				break;
			}
		}
		ASSERT( (ipTe >=0) && (ipTe < NumElements-1)  );
			
	}
	else if( phycon.alogte < (float)TempArray[ipTe] )
	{
		/* temp is too low, must also lower ipTe */
		ASSERT( phycon.alogte > (float)TempArray[0] );
		/* decrement ipTe until it is correct */
		while( ( phycon.alogte < (float)TempArray[ipTe] ) && ipTe > 0)
			--ipTe;
	}
	else if( phycon.alogte > (float)TempArray[ipTe+1] )
	{
		/* temp is too high */
		ASSERT( phycon.alogte <= (float)TempArray[NumElements-1] );
		/* increment ipTe until it is correct */
		while( ( phycon.alogte > (float)TempArray[ipTe+1] ) && ipTe < NumElements-1)
			++ipTe;
	}

	ASSERT( (ipTe >=0) && (ipTe < NumElements-1)  );

	/* ipTe should now be valid */
	ASSERT( ( phycon.alogte >= (float)TempArray[ipTe] )
		&& ( phycon.alogte <= (float)TempArray[ipTe+1] ) && ( ipTe < NumElements-1 ) );

	rate = ((double)phycon.alogte - TempArray[ipTe]) / (TempArray[ipTe+1] - TempArray[ipTe]) *
		(ValueArray[ipTe+1] - ValueArray[ipTe]) + ValueArray[ipTe] ;

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

	return rate;
}

/*He_RRCoef_Te evaluated radiative recombination coef at some temperature */
static double He_RRCoef_Te( long nelem , long n )
{
	double rate = 0.;

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

	/* if n is equal to the number of levels, return the total recomb, else recomb for given level.	*/
	if( n == iso.numLevels[ipHE_LIKE][nelem] - iso.nCollapsed[ipHE_LIKE][nelem] )
	{
		rate = TempInterp( TeRRCoef, TotalRecomb[nelem], N_HE_TE_RECOMB );
	}
	else
	{
		rate = TempInterp( TeRRCoef, RRCoef[nelem][n], N_HE_TE_RECOMB );
	}

	/* that was the log, now make linear */
	rate = pow( 10. , rate );

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

	return rate;
}

/*HelikeCheckRecomb called by SanityCheck to confirm that recombination coef are ok
 * return value is relative error between new calculation of recom, and interp value */
double HelikeCheckRecomb(
	/* the chemical element, 1 for He */
	long nelem ,
	/* the level, 0 for ground */
	long level ,
	/* the temperature to be used */
	double temperature )
{
	double RecombRelError ,
		RecombInterp,
		RecombCalc,
		SaveTemp;

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

	/* this global variable is needed by radrecomb */
	EthRyd = iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][level];

	/* first set temp todesired value */
	SaveTemp = phycon.te;
	phycon.te = (float)temperature;
	tfidle(TRUE);

	/* actually calculate the recomb coef from the Milne relation,
	 * normally only done due to compile he-like command */
	RecombCalc = radrecomb( temperature , nelem , level );

	/* interpolate the recomb coef, this is the usual method */
	RecombInterp = He_RRCoef_Te( nelem , level );

	/* reset temp */
	phycon.te = (float)SaveTemp;
	tfidle(TRUE);

	RecombRelError = ( RecombInterp - RecombCalc ) / MAX2( RecombInterp , RecombCalc );

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

	return RecombRelError ;
}


/*lint +e662 creation of  out of bound pointer */

#undef PARALLEL
/*lint +e662 creation of  out of bound pointer */
/*lint +e661 creation of  out of bound pointer */









