/* 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 */
/*ChargTranEval fill in the HCharExcIonOf and Rec arrays with Kingdon's fitted CT with H */
/*ChargTranSumHeat sum net heating due to charge transfer, called in HeatSum */
/*HCTIon H charge transfer ionization*/ 
/*HCTRecom H charge transfer recombination */
/*ChargTranPun, punch charge transfer coef */
/*MakeHCTData holds data for charge tranfer fits */
#include "cddefines.h"
#include "phycon.h"
#include "physconst.h"
#include "abund.h"
#include "dense.h"
#include "iso.h"
#include "thermal.h"
#include "co.h"
#include "elementnames.h"
#include "tfidle.h"
#include "heavy.h"
#include "trace.h"
#include "converge.h"
#include "atmdat.h"
/*lint -e662 possible creation out of bounds pointer */

/*HCTion H charge transfer ionization, H+ + A => H + A+ */
static double HCTIon(long int ion, 
  long int nelem);

/*HCTRecom H charge transfer recombination, H + A+ => H+ + A */
static double HCTRecom(long int ion, 
  long int nelem);

/* the translated block data */
static void MakeHCTData(void);

/* the structures for storing the charge transfer data, these are filled in
 * at the end of this file, in what used to be a block data  */
static double CTIonData[LIMELM][4][8];
static double CTRecombData[LIMELM][4][7];
/* this will be flag for whether or not charge transfer data
 * have been initialized */
static int lgCTDataDefined = FALSE;

/*ChargTranEval update charge trans rate coefficients if temperature has changed */
void ChargTranEval(
		/* return value is H ionization rate (s-1) due to O+ charge transfer */
		double *O_HIonRate )
{
	long int ion, 
	  nelem;
	double Ln_T, a, b, c, a_op, b_op, c_op, d_op, e_op, f_op, a_o, 
			b_o, c_o, d_o, e_o, f_o, g_o; 

	static double TeUsed = -1.;

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

	/* first is to force reevaluation on very first call */
	if( !conv.nTotalIoniz || fabs(phycon.te-TeUsed)/phycon.te > 0.01 )
	{
		/* refresh hydrogen charge transfer arrays */
		/* >>chng 01 apr 25, lower limit had been 2, lithium, changed to 1, helium */
		for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
		{
			for( ion=0; ion <= nelem; ion++ )
			{
				/* >>chng 01 apr 28, add factors  to turn off ct,
				 * had previously been handled downstream */

				/* HCharExcIonOf[nelem][ion]*hii  is ion => ion+1 for nelem */
				/* charge transfer ionization O^0 + H+ -> O+ + H0 
				 * is HCharExcIonOf[ipOXYGEN][0]*dense.xIonDense[ipHYDROGEN][1]
				 * charge transfer recombination of atomic hydrogen is
				 * HCharExcIonOf[ipOXYGEN][0]*dense.xIonDense[ipOXYGEN][0] */
				atmdat.HCharExcIonOf[nelem][ion] = HCTIon(ion,nelem);

				/* HCharExcRecTo[nelem][ion]*hi is ion+1 => ion of nelem */
				/* charge transfer recombination O+ + H0 -> O^0 + H+ is
				 * HCharExcRecTo[ipOXYGEN][0]*dense.xIonDense[ipHYDROGEN][0]
				 * charge transfer ionization of atomic hydrogen is
				 * HCharExcRecTo[ipOXYGEN][0]*dense.xIonDense[ipOXYGEN][1] */
				atmdat.HCharExcRecTo[nelem][ion] = HCTRecom(ion,nelem);
			}

			/* zero out helium charge trasnfer arrays */
			for( ion=0; ion < LIMELM; ion++ )
			{
				atmdat.HeCharExcIonOf[nelem][ion] = 0.;
				atmdat.HeCharExcRecTo[nelem][ion] = 0.;
			}
		}

		/* >>chng 04 apr 27, the above included only the radiative charge transfer from
		 * stancil et al 1998.  must explicitly add on the ct fitted by Kingdon & Ferland,
		 * the code must have originally had the KF fit, but got hosed at some point 
		 * these are two different sets of physics 
		 * He0 + H+ -> He+ + H0 */
		if( phycon.te > 6000. )
			atmdat.HCharExcIonOf[ipHELIUM][0] += 7.47e-15*pow(phycon.te/1e4,2.06)*
			(1.+9.93*sexp(3.89e-4*phycon.te) );

		/* >>chng 04 jun 21 -- NPA.  Put in the charge transfer rate for:
			He+ + H => He + H+ as defined in the UMIST database.  This is only
			used if the "set UMIST rates" command is used */
		if(!co.lgUMISTrates)
		{
			atmdat.HCharExcRecTo[ipHELIUM][0] = 4.85e-15*pow(phycon.te/300, 0.18);
		}
		/* >>chng 04 jun 21 -- NPA.  update the charge transfer rates between hydrogen
		   and oxygen to:  
		   >>refer	O	CT	Stancil et al. 1999, A&AS, 140, 225-234 
		   Instead of using the UMIST rate, the program TableCurve was used
		   to generate a fit to the data listed in Tables 2, 3, and 4 of the 
		   aforementioned reference. The following fitting equations agree 
		   very well with the published data. */

		Ln_T = log(phycon.te);

		/* At or below 10K, just use the value the formula's below give
		   at 10K.*/
		/* do both O+ -> O and O -> O+ for low T limit */
		if (phycon.te <= 10. )
		{
			atmdat.HCharExcRecTo[ipOXYGEN][0] = 3.744e-10;
			atmdat.HCharExcIonOf[ipOXYGEN][0] = 4.749e-20;
		}

		/* this does O+ -> O for all higher temperatures */
		if( phycon.te > 10.)
		{
			a_op = 2.3344302e-10;
			b_op = 2.3651505e-10;
			c_op = -1.3146803e-10;
			d_op = 2.9979994e-11;
			e_op = -2.8577012e-12;
			f_op = 1.1963502e-13;

			/* equation rank 53 of TableCurve */
			atmdat.HCharExcRecTo[ipOXYGEN][0] =  a_op + b_op*Ln_T + c_op*pow(Ln_T,2) + d_op*pow(Ln_T,3)
				+ e_op*pow(Ln_T,4) + f_op*pow(Ln_T, 5);
		}

		/* now do O -> O+
		 * The next two equations were chosen to match up at 200K, so that there
		 *are no discontinuities */
		if ((phycon.te > 10.) && (phycon.te <= 190.))
		{
			a = -21.134531;
			b = -242.06831;
			c = 84.761441;
			
			/* equation rank 2 of TableCurve */
			atmdat.HCharExcIonOf[ipOXYGEN][0] = exp((a + b/SDIV(phycon.te) + c/SDIV(phycon.tesqrd)));
		}

		else if ((phycon.te > 190.) && (phycon.te <= 200.))
		{
		
			/* We are using two fitting formula's for this rate, in order to assure no
			sudden "jumps" in the rate, the rate between 190-200K is made to 
			increase linearly.  The formula below gets the same answer as the equation
			above at 190K, and gets the same answer as the the formula below this one 
			at 200K*/
			atmdat.HCharExcIonOf[ipOXYGEN][0] = 2.18733e-12*(phycon.te-190) + 1.85823e-10;
		}

		else if (phycon.te > 200.)
		{
		
			a_o = -7.6767404e-14;
			b_o = -3.7282001e-13;
			c_o = -1.488594e-12;
			d_o = -3.6606214e-12; 
			e_o = 2.0699463e-12;
			f_o = -2.6139493e-13;
			g_o = 1.1580844e-14;

			/* equation rank 120 of TableCurve */
			atmdat.HCharExcIonOf[ipOXYGEN][0] = a_o + b_o*Ln_T + c_o*pow(Ln_T,2) + d_o*pow(Ln_T,3)
										+ e_o*pow(Ln_T,4) + f_o*pow(Ln_T, 5) + g_o*pow(Ln_T,6);
		}

		/* >>chng 01 may 07, following had all been +=, ok if above was zero.
		 * changed to = and added HCTOn */
		/* >>chng 01 jan 30, add following block of CT reactions */
		/* ionization, added as per Phillip Stancil OK, number comes from 
		 * >>refer	Fe	CT	Tielens, A.G.G.M., & Hollenbach, D., 1985a, ApJ, 294, 722-746 */
		/* >>refer	Fe	CT	Prasad, S.S., & Huntress, W.T., 1980, ApJS, 43, 1-35 */
		/* the actual rate comes from the following paper: */
		/* >>refer	Fe	CT	Pequignot, D., & Aldrovandi, S.M.V., 1986, A&A, 161, 169-176 */
		/* Fe0 + H+ => Fe+ + H0 */
		atmdat.HCharExcIonOf[ipIRON][0] = 3e-9;

		/* all remaining entries are from Pequignot & Aldrovandi*/
		/* >>refer	Al	CT	Pequignot, D., & Aldrovandi, S.M.V., 1986, A&A, 161, 169-176 */
		/* Al0 + H+ => Al+ + H0 */
		atmdat.HCharExcIonOf[ipALUMINIUM][0] = 3e-9;

		/* >>refer	P	CT	Pequignot, D., & Aldrovandi, S.M.V., 1986, A&A, 161, 169-176 */
		/* P0 + H+ => P+ + H0 */
		atmdat.HCharExcIonOf[ipPHOSPHORUS][0] = 1e-9;

		/* >>refer	Cl	CT	Pequignot, D., & Aldrovandi, S.M.V., 1986, A&A, 161, 169-176 */
		/* Cl0 + H+ => Cl+ + H0 */
		atmdat.HCharExcIonOf[ipCHLORINE][0] = 1e-9;

		/* >>refer	Ti	CT	Pequignot, D., & Aldrovandi, S.M.V., 1986, A&A, 161, 169-176 */
		/* Cl0 + H+ => Cl+ + H0 */
		atmdat.HCharExcIonOf[ipTITANIUM][0] = 3e-9;

		/* >>refer	Mn	CT	Pequignot, D., & Aldrovandi, S.M.V., 1986, A&A, 161, 169-176 */
		/* Mn0 + H+ => Mn+ + H0 */
		atmdat.HCharExcIonOf[ipMANGANESE][0] = 3e-9;

		/* >>refer	Ni	CT	Pequignot, D., & Aldrovandi, S.M.V., 1986, A&A, 161, 169-176 */
		/* Ni0 + H+ => Ni+ + H0 */
		atmdat.HCharExcIonOf[ipNICKEL][0] = 3e-9;

		/* >>chng 01 feb 02, add following CT reaction from */
		/* >>refer	Na0	CT	Dutta, C.M., Nordlander, P., Kimura, M., & Dalgarno, A., 2001, Phys REv A, 63, 022709 */
		/* this is roughly their number around 500K - they do not give explicit values, rather
		 * a small figure.  Previous cals were 5 orders of mag smaller at this temp.  
		 * ND this deposites electron into n=2 */
		/* Na0 + H+ => Na+ + H0(n=2) */
		atmdat.HCharExcIonOf[ipSODIUM][0] = 7e-12;

		/* ===================== helium charge transfer ====================*/

		/* atmdat.HeCharExcIonOf is ionization, */
		/* [0] is Atom^0 + He+ => Atom+1 + He0
		 * [n] is Atom^+n + He+ => Atom^+n-1 + He0 */

		/* atmdat.HeCharExcRecTo is recombination */
		/* [0] is Atom^+1 + He0 => Atom^0 + He^+
		 * [n] is Atom^+n+1 + He0 => Atom^+n + He^+ */

		/* Carbon */
		/* recombination */
		/* C+3 + He => C+2 + He+ */
		/* >>refer	C+3	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcRecTo[ipCARBON][2] = 4.6e-19*phycon.tesqrd;

		/* C+4 + He => C+3 + He+ */
		/* >>refer	C+4	CT	Butler, S.E., & Dalgarno, 1980b */
		atmdat.HeCharExcRecTo[ipCARBON][3] = 1e-14;

		/* ionization */
		/* C0 + He+ => C+ + He0 */
		/* unknown reference, from older version of the code */
		/*atmdat.HeCharExcIonOf[ipCARBON][0] = 4.17e-17*(phycon.te/phycon.te03);*/
		
		/* >>chng 04 jun 21 -- update this rate to that given in the UMIST database - NPA */
		atmdat.HeCharExcIonOf[ipCARBON][0] = 6.3e-15*pow((phycon.te/300),0.75);
		
		/* C+1 + He+ => C+2 + He */
		/* >>refer	C+1	CT	Butler, S.E., Heil,T.G., & Dalgarno, A. 1980, ApJ, 241, 442*/
		atmdat.HeCharExcIonOf[ipCARBON][1] = 
			5e-20*phycon.tesqrd*sexp(0.07e-4*phycon.te)*sexp(6.29/phycon.te_eV) ;

		/* nitrogen */
		/* recombination */
		/* N+2 => N+ Butler and Dalgarno 1980B
		 * ct with update
		 * >>refer	N+2	CT	Sun, Sadeghpour, Kirby, Dalgarno, and Lafyatis, cfa preprint 4208
		 * this agrees with exp 
		 * >>refer	N+2	CT	Fand&Kwong, ApJ 474, 529 */
		atmdat.HeCharExcRecTo[ipNITROGEN][1] = 0.8e-10;

		/* N+3 => N+2 */
		/* >>refer	N+3	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcRecTo[ipNITROGEN][2] = 1.5e-10;

		/* ce rate from quantal calc of ce,
		 * >>refer	N+4	CT	Feickert, Blint, Surratt, and Watson, (preprint Sep 84). Ap.J. in press,
		 * >>refer	N+4	CT	Rittby et al J Phys B 17, L677, 1984.
		 * CR = 1.E-9 + 8E-12 * TE10 * SQRTE */
		atmdat.HeCharExcRecTo[ipNITROGEN][3] = 2e-9;

		/* ionization */
		/* N+1 + He+ => N+2 + He */
		/* >>refer	N+1	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcIonOf[ipNITROGEN][1] = 
			3.7e-20*phycon.tesqrd*sexp(0.063e-4*phycon.te)*sexp(1.44/phycon.te_eV) ;

		/* oxygen */
		/* recombination */
		/* O+2 + He  => O+1 + He+ */
		/* >>refer	O+2	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcRecTo[ipOXYGEN][1] = 3.2e-14*phycon.te/phycon.te05;
		/* O+3 + He => O+2 + He+ */
		/* >>refer	O+3	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcRecTo[ipOXYGEN][2] = 1e-9;
		/* O+4 + He  => O+3 + He+ */
		/* >>refer	O+4	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcRecTo[ipOXYGEN][3] = 6e-10;

		/* neon */
		/* recombination */
		/* Ne+2 + He  => Ne+1 + He+ */
		/* >>refer	Ne+2	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcRecTo[ipNEON][1] = 1e-14;
		/* Ne+3 + He => Ne+2 + He+ */
		/* >>refer	Ne+3	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcRecTo[ipNEON][2] = 1e-16*phycon.sqrte;
		/* Ne+4 + He  => Ne+3 + He+ */
		/* >>refer	Ne+4	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcRecTo[ipNEON][3] = 1.7e-11*phycon.sqrte;

		/* magnesium */
		/* recombination */
		/* Mg+3 + Heo => Mg+2 */
		/* >>refer	Mg+3	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcRecTo[ipMAGNESIUM][2] = 7.5e-10;
		/* Mg+4 + Heo => Mg+3 */
		/* >>refer	Mg+4	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcRecTo[ipMAGNESIUM][3] = 1.4e-10*phycon.te30;


		/* silicon */
		/* recombination */
		/* Si+3 +He => Si+2 + He+ */
		/* >>refer	Si+3	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 
		 * scale by 1.3 to bring into agreement with
		 * >>refer	Si+3	CT	Fang, Z., & Kwong, H.S. 1997 ApJ 483, 527 */
		atmdat.HeCharExcRecTo[ipSILICON][2] += phycon.sqrte*phycon.te10*phycon.te10*
		  1.3*1.5e-12;

		/* Si+4 + Heo => Si+3
		 * >>refer	Si+4	CT	Opradolce et al., 1985, A&A, 148, 229 */
		atmdat.HeCharExcRecTo[ipSILICON][3] = 2.54e-11*phycon.sqrte/phycon.te03/
		  phycon.te01/phycon.te01;

		/* ionization */
		/* Si0 + He+ => Si+ + He0 */
		/* >>refer	Si0	CT	Prasad, S.S., & Huntress, W.T., 1980, ApJS, 43, 1-35 */
		atmdat.HeCharExcIonOf[ipSILICON][0] = 3.3e-9;

		/* Si+1 + He+ => Si+2 + He */
		/* >>refer	Si+1	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcIonOf[ipSILICON][1] = 
			1.5e-11*phycon.te20*phycon.te05*sexp(6.91/phycon.te_eV) ;

		/* Si+2 + He+ => Si+3 + He */
		/* >>refer	Si+2	CT	Gargaud, M., McCarroll, R., & Valiron, P. 1982, A&ASup, 45, 603 */
		atmdat.HeCharExcIonOf[ipSILICON][2] = 
			1.15e-11*phycon.sqrte*sexp(8.88/phycon.te_eV) ;

		/* sulphur */
		/* recombination */
		/* S+3 + Heo => S+2 */
		/* >>refer	S+3	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcRecTo[ipSULPHUR][2] = phycon.sqrte*1.1e-11;

		/* S+4 + Heo => S+3 */
		/* >>refer	S+4	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		/* >>chng 04 jul 01, from [ipSULPHUR][2] to [3] - bug */
		atmdat.HeCharExcRecTo[ipSULPHUR][3] = 4.8e-14*phycon.te30;

		/* ionization */
		/* S+1 + He+ => S+2 + He */
		/* >>refer	S+1	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcIonOf[ipSULPHUR][1] = 
			4.4e-16*phycon.te*phycon.te20*sexp(0.036e-4*phycon.te)*sexp(9.2/phycon.te_eV) ;

		/* S+2 + He+ => S+3 + He */
		/* >>refer	S+2	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcIonOf[ipSULPHUR][2] = 
			5.5e-18*phycon.te*phycon.sqrte*phycon.te10*sexp(0.046e-4*phycon.te)*sexp(10.5/phycon.te_eV) ;

		/* Argon */
		/* recombination */
		/* >>refer	Ar+2	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcRecTo[ipARGON][1] = 1.3e-10;

		/* >>refer	Ar+3	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcRecTo[ipARGON][2] = 1.e-14;

		/* >>refer	Ar+4	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcRecTo[ipARGON][3] = 1.6e-8/phycon.te30;

		/* ionization */
		/* Ar+1 + He+ => Ar+2 + He0 */
		/* >>refer	Ar+1	CT	Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */
		atmdat.HeCharExcIonOf[ipARGON][1] = 1.1e-10;

		TeUsed = phycon.te;

		if (!co.lgUMISTrates)
		{
			/* Set all charge transfer rates equal to zero that do not appear
			   in the UMIST database.  This if statement is only performed
			   if the "set UMIST rates" command is used */

			atmdat.HCharExcIonOf[ipHELIUM][0] = 0;
			atmdat.HCharExcIonOf[ipCARBON][0] = 0;
			atmdat.HCharExcRecTo[ipCARBON][0] = 0;

			atmdat.HeCharExcRecTo[ipCARBON][0] = 0;
			atmdat.HeCharExcIonOf[ipOXYGEN][0] = 0;
			atmdat.HeCharExcRecTo[ipOXYGEN][0] = 0;
		}
			

		/* this is set FALSE with the no charge transfer command */
		if( !atmdat.lgCTOn )
		{
			for( nelem=0; nelem< LIMELM; ++nelem )
			{
				for( ion=0; ion<LIMELM; ++ion )
				{
					atmdat.HCharExcIonOf[nelem][ion] = 0.;
					atmdat.HCharExcRecTo[nelem][ion] = 0.;
					atmdat.HeCharExcIonOf[nelem][ion] = 0.;
					atmdat.HeCharExcRecTo[nelem][ion] = 0.;
				}
			}
		}
	}

	/* return value is H ionization rate (s-1) due to O+ charge transfer */
	*O_HIonRate = atmdat.HCharExcRecTo[ipOXYGEN][0]*dense.xIonDense[ipOXYGEN][1];

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

/*================================================================================*
 *================================================================================*/
double ChargTranSumHeat(void)
{
	long int ion, 
	  nelem;
	double SumCTHeat_v;

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

	/* second dimension is ionization stage,
	 * 1=+1 for parent, etc
	 * third dimension is atomic weight of atom */

	/* make sure data are initialized */
	ASSERT( lgCTDataDefined );

	SumCTHeat_v = 0.;
	/* >>chng 01 apr 25, lower limit had been 0 should have been 1 (helium) */
	for( nelem=1; nelem < LIMELM; nelem++ )
	{
		/* >>chng >>01 apr 25, loops had been to LIMELM, which may have done no harm
		 * since extra elements were set to zero, but is incorrect since the physical
		 * limit is the number of stages of ionization */
		int limit = MIN2(4, nelem);
		/* this first group of lower stages of ionization have exact rate coefficients */
		for( ion=0; ion < limit; ion++ )
		{
			/* CTIonData[nelem][ion][7] and CTRecombData[nelem][ion][6] are the energy deficits in eV,
			 * atmdat.HCharExcIonOf[nelem][ion] and atmdat.HCharExcIonOf[nelem][ion] 
			 * save the rate coefficients 
			 * this is sum of heat exchange in eV s^-1 cm^-3 */
			SumCTHeat_v += 

				/* heating due to ionization of heavy element, recombination of hydrogen */
				atmdat.HCharExcIonOf[nelem][ion]*CTIonData[nelem][ion][7]*
				(double)dense.xIonDense[ipHYDROGEN][1]*(double)dense.xIonDense[nelem][ion] + 

				/* heating due to recombination of heavy element, ionization of hydrogen */
				atmdat.HCharExcRecTo[nelem][ion]*CTRecombData[nelem][ion][6]*
				iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.xIonDense[ipHYDROGEN][1]*
				(double)dense.xIonDense[nelem][ion+1];
		}

		/* >>chng >>01 apr 25, following loop had been to LIMELM, change to nelem */
		/* following do not have exact energies, so use 2.86*(Z-1) */
		for( ion=4; ion < nelem; ion++ )
		{
			SumCTHeat_v += 
				atmdat.HCharExcRecTo[nelem][ion]* 2.86*(double)ion *
				(double)dense.xIonDense[ipHYDROGEN][0]*(double)dense.xIonDense[nelem][ion+1];
		}
	}

	/* convert from eV to ergs, HCharHeatOn usually 1, set to 0 with no CTHeat,  
	 * EN1EV is ergs in 1 eV, 1.602176e-012*/
	SumCTHeat_v *= EN1EV * atmdat.HCharHeatOn;

	if( thermal.htot > 1e-35 )
	{
		/* remember largest fractions of heating and cooling for comment */
		atmdat.HCharHeatMax = MAX2(atmdat.HCharHeatMax,
			SumCTHeat_v/thermal.htot );

		atmdat.HCharCoolMax = MAX2(atmdat.HCharCoolMax,
			-SumCTHeat_v/thermal.htot);
	}

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

/*================================================================================*
 *================================================================================*/
static double HCTIon(
	/* ion is stage of ionization on C scale, 0 for atom */
	long int ion,		
	/* nelem is atomic number of element on C scale, 1 to 29 */
	/* HCTIon(1,5) is C+ + H+ => C++ + H */
	long int nelem)	
{
	double HCTIon_v, 
	  tused;

#	ifdef DEBUG_FUN
	fputs( "<+>HCTIon()\n", debug_fp );
#	endif
	/* H charge transfer ionization, using Jim Kingdon's ctdata.for */

	/* set up the rate coefficients if this is first call */
	if( !lgCTDataDefined )
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ," HCTIon doing 1-time init of charge transfer data\n");
		}
		lgCTDataDefined = TRUE;
		MakeHCTData();
	}

	/* check that data have been linked in,
	 * error here would mean that data below had not been loaded */
	ASSERT( CTIonData[2][0][0] > 0. );

	/* return zero for highly ionized species */
	if( ion >= 3 )
	{
		HCTIon_v = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->HCTIon()\n", debug_fp );
#		endif
		return( HCTIon_v );
	}

	/*begin sanity checks */
	/* check that ionization stage is ok for this element*/
	ASSERT( ion >= 0);
	ASSERT( ion <= nelem );

	/* now check the element is valid, this is routine HCTIon */
	ASSERT( nelem > 0 );
	ASSERT( nelem < LIMELM );

	/* may be no entry, first coefficient is zero in this case */
	if( CTIonData[nelem][ion][0] <= 0. )
	{
		HCTIon_v = 0.;
		
	}

	else
	{
		/* Make sure te is between temp. boundaries; set constant outside of range */
		tused = MAX2((double)phycon.te,CTIonData[nelem][ion][4]);
		tused = MIN2(tused,CTIonData[nelem][ion][5]);
		tused *= 1e-4;

		/* the interpolation equation */
		HCTIon_v = CTIonData[nelem][ion][0]*1e-9*(pow(tused,CTIonData[nelem][ion][1]))*
		  (1. + CTIonData[nelem][ion][2]*exp(CTIonData[nelem][ion][3]*tused))*
		  exp(-CTIonData[nelem][ion][6]/tused);
	}

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

/*================================================================================*
 *================================================================================*/
static double HCTRecom(
	/* ion is stage of ionization on C scale, 0 for rec to atom */
	long int ion,  	
	/* nelem is atomic number of element on C scale, 1 = he up to limelm */
	/* HCTRecom(1,5) would be C+2 + H => C+ + H+ */
	long int nelem)	
{
	double HCTRecom_v, 
	  tused;

#	ifdef DEBUG_FUN
	fputs( "<+>HCTRecom()\n", debug_fp );
#	endif
	/* 
	 * H charge transfer recombination using Jim Kingdon's block ctdata.for
	 */

	/* set up the rate coefficients if this is first call */
	if( !lgCTDataDefined )
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ," HCTIon doing 1-time init of charge transfer data\n");
		}
		lgCTDataDefined = TRUE;
		MakeHCTData();
	}

	/* this is check that data have been set up properly, will
	 * fail if arrays are not initialized properly */
	ASSERT( CTRecombData[1][0][0] > 0. );

	/* use Dalgarno estimate for highly ionized species, number reset with
	 * set charge transfer command */
	if( ion > 3 )
	{
		/* >>chng 96 nov 25, added this option, default is 1.92e-9
		 * Dalgarno's charge transfer */
		HCTRecom_v = atmdat.HCTAlex*((double)ion+1.);
		
#		ifdef DEBUG_FUN
		fputs( " <->HCTRecom()\n", debug_fp );
#		endif
		return( HCTRecom_v );
	}

	/* check that ion stage within bound for this atom */
	ASSERT( ion >= 0 && ion <= nelem );

	/* now check the element is valid, this is routine HCTIon */
	ASSERT( nelem > 0 && nelem < LIMELM );

	tused = MAX2((double)phycon.te,CTRecombData[nelem][ion][4]);
	tused = MIN2(tused,CTRecombData[nelem][ion][5]);
	tused *= 1e-4;

	if( tused == 0. )
	{
		HCTRecom_v = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->HCTRecom()\n", debug_fp );
#		endif
		return( HCTRecom_v );
	}

	/* the interpolation equation */
	HCTRecom_v = CTRecombData[nelem][ion][0]*1e-9*(pow(tused,CTRecombData[nelem][ion][1]))*
	  (1. + CTRecombData[nelem][ion][2]*sexp(-CTRecombData[nelem][ion][3]*tused));

	/* in sexp negative sign not typo - there are negative signs already
	 * in coefficient, and sexp has implicit negative sign */

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

/*================================================================================*
 *================================================================================*/
/*block data with Jim Kingdon's charge transfer data */
/* >>refer	H	CT	Kingdon, J, B., & Ferland, G.J. 1996, ApJS, 106, 205 */
/* 
 * first dimension is atomic number of atom, 0 for H 
 * second dimension is ionization stage,
 * 1=+0 for parent, etc
 * third dimension is atomic number of atom 
 * second dimension is ionization stage,
 * 1=+1 for parent, etc
 */

/* digital form of the fits to the charge transfer
 * ionization rate coefficients 
 *
 * Note: First parameter is in units of 1e-9!
 * Note: Seventh parameter is in units of 1e4 K */

/* digital form of the fits to the charge transfer
 * recombination rate coefficients (total)
 *
 * Note: First parameter is in units of 1e-9!
 * recombination 
 */

/* holds data for charge tranfer fits */
static void MakeHCTData(void)
{
	long int i, 
		j,
		nelem,
	  _r;

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

	/* >>chng 01 apr 24, zero out this block, as per PvH comments that
	 * translated block datas do not fully initialize arrays */
	/* first zero out entire arrays, since some may not have charge transfer data */
	for( nelem=0; nelem<LIMELM; ++nelem )
	{
		for( i=0; i<4; ++i )
		{
			for( j=0; j<7; ++j )
			{
				CTIonData[nelem][i][j] = 0.;
				CTRecombData[nelem][i][j] = 0.;
			}
			CTIonData[nelem][i][7] = 0.;
		}
	}

	/* 
	 * following are coefficients for charge transfer ionization,
	 * H+ + A => H + A+
	 */
	/* Lithium +0 */
	{ static double _itmp0[] = {2.84e-3 , 1.99 , 375.54 , -54.07 , 1e2 , 1e4 , 0.,
		-10.};

	for( i=1, _r = 0; i <= 8; i++ )
	{
		CTIonData[2][0][i-1] = _itmp0[_r++];
		}
	}

	/* C+0 ionization */
	{ static double _itmp1[] = {1.07e-6 , 3.15 , 176.43 , -4.29 , 1e3 , 1e5 , 0. ,2.34};
	for( i=1, _r = 0; i <= 8; i++ )
	{
		CTIonData[5][0][i-1] = _itmp1[_r++];
		}
	}
	{ static double _itmp2[] = {4.55e-3,-0.29,-0.92,-8.38,1e2,5e4,
	  1.086,-0.94};
	for( i=1, _r = 0; i <= 8; i++ )
	{
		CTIonData[6][0][i-1] = _itmp2[_r++];
		}
	}
	/* oxygen */
	{ static double _itmp3[] = {7.40e-2,0.47,24.37,-0.74,1e1,1e4,
	  0.023,-0.02};
	for( i=1, _r = 0; i <= 8; i++ )
	{
		CTIonData[7][0][i-1] = _itmp3[_r++];
		}
	}
	{ static double _itmp4[] = {3.34e-6,9.31,2632.31,-3.04,1e3,
	  2e4,0.0,-1.74};
	for( i=1, _r = 0; i <= 8; i++ )
	{
		CTIonData[10][0][i-1] = _itmp4[_r++];
		}
	}
	{ static double _itmp5[] = {9.76e-3,3.14,55.54,-1.12,5e3,3e4,
	  0.0,1.52};
	for( i=1, _r = 0; i <= 8; i++ )
	{
		CTIonData[11][0][i-1] = _itmp5[_r++];
		}
	}
	{ static double _itmp6[] = {7.60e-5,0.00,-1.97,-4.32,1e4,3e5,
	  1.670,-1.44};
	for( i=1, _r = 0; i <= 8; i++ )
	{
		CTIonData[11][1][i-1] = _itmp6[_r++];
		}
	}
	{ static double _itmp7[] = {0.92,1.15,0.80,-0.24,1e3,2e5,0.0,
	  0.12};
	for( i=1, _r = 0; i <= 8; i++ )
	{
		CTIonData[13][0][i-1] = _itmp7[_r++];
		}
	}
	/* Si+1 ionization */
	{ static double _itmp8[] = {2.26 , 7.36e-2 , -0.43 , -0.11 , 2e3 , 1e5 , 3.031
		,-2.72};
	for( i=1, _r = 0; i <= 8; i++ )
	{
		CTIonData[13][1][i-1] = _itmp8[_r++];
		}
	}
	{ static double _itmp9[] = {1.00e-5,0.00,0.00,0.00,1e3,1e4,
	  0.0,-3.24};
	for( i=1, _r = 0; i <= 8; i++ )
	{
		CTIonData[15][0][i-1] = _itmp9[_r++];
		}
	}
	{ static double _itmp10[] = {4.39,0.61,-0.89,-3.56,1e3,3e4,
	  3.349,-2.89};
	for( i=1, _r = 0; i <= 8; i++ )
	{
		CTIonData[23][1][i-1] = _itmp10[_r++];
		}
	}
	{ static double _itmp11[] = {2.83e-1,6.80e-3,6.44e-2,-9.70,
	  1e3,3e4,2.368,-2.04};
	for( i=1, _r = 0; i <= 8; i++ )
	{
		CTIonData[24][1][i-1] = _itmp11[_r++];
		}
	}
	{ static double _itmp12[] = {2.10,7.72e-2,-0.41,-7.31,1e4,1e5,
	  3.005,-2.56};
	for( i=1, _r = 0; i <= 8; i++ )
	{
		CTIonData[25][1][i-1] = _itmp12[_r++];
		}
	}
	{ static double _itmp13[] = {1.20e-2,3.49,24.41,-1.26,1e3,3e4,
	  4.044,-3.49};
	for( i=1, _r = 0; i <= 8; i++ )
	{
		CTIonData[26][1][i-1] = _itmp13[_r++];
		}
	}
	/* CT recombintion, A+n + H => A+n-1 + H+ */
	/* >>chng 01 may 03, first coef multiplied by 0.25, as per comment in
	 * >>refer	Li	CT	Stancil, P.C., & Zygelman, B., 1996, ApJ, 472, 102
	 * which corrected the error in 
	 * >>refer	He	CT	Zygelman, B., Dalgarno, A., Kimura, M., & Lane, N.F.,
	 * >>refercon	1989, Phys. Rev. A, 40, 2340
	 * this was used in the original Kingdon & Ferland paper so no correction required
	 * >>chng 04 apr 27, He was in error above as well, factor of 4, noted in 
	 * >>refer	He	CT	Stancil, P.C., Lepp, S., & Dalgarno, A. 1998, ApJ, 509, 1
	 */
	{ static double _itmp14[] = {/*7.47e-6*/1.87e-6,2.06,9.93,-3.89,6e3,1e5,
	  10.99};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[1][0][i-1] = _itmp14[_r++];
		}
	}
	{ static double _itmp15[] = {1.00e-5,0.,0.,0.,1e3,1e7,-40.81};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[1][1][i-1] = _itmp15[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[2][0][i-1] = 0.f;
		}
	{ static double _itmp16[] = {1.26,0.96,3.02,-0.65,1e3,3e4,3.02};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[2][1][i-1] = _itmp16[_r++];
		}
	}
	{ static double _itmp17[] = {1.00e-5,0.,0.,0.,2e3,5e4,-108.83};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[2][2][i-1] = _itmp17[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[3][0][i-1] = 0.f;
		}
	{ static double _itmp18[] = {1.00e-5,0.,0.,0.,2e3,5e4,-4.61};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[3][1][i-1] = _itmp18[_r++];
		}
	}
	{ static double _itmp19[] = {1.00e-5,0.,0.,0.,2e3,5e4,-140.26};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[3][2][i-1] = _itmp19[_r++];
		}
	}
	{ static double _itmp20[] = {5.17,0.82,-0.69,-1.12,2e3,5e4,
	  10.59};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[3][3][i-1] = _itmp20[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[4][0][i-1] = 0.f;
		}
	{ static double _itmp21[] = {2.00e-2,0.,0.,0.,1e3,1e9,2.46};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[4][1][i-1] = _itmp21[_r++];
		}
	}
	{ static double _itmp22[] = {1.00e-5,0.,0.,0.,2e3,5e4,-24.33};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[4][2][i-1] = _itmp22[_r++];
		}
	}
	/* B+4 recombinatino */
	{ static double _itmp23[] = {2.74 , 0.93 , -0.61 , -1.13 , 2e3 , 5e4 ,
	  11.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[4][3][i-1] = _itmp23[_r++];
		}
	}
	/* C+1 recombinatino */
	{ static double _itmp24[] = {4.88e-7 , 3.25 , -1.12 , -0.21 , 5.5e3 , 1e5 ,
		-2.34};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[5][0][i-1] = _itmp24[_r++];
		}
	}
	{ static double _itmp25[] = {1.67e-4,2.79,304.72,-4.07,5e3,
	  5e4,4.01};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[5][1][i-1] = _itmp25[_r++];
		}
	}
	{ static double _itmp26[] = {3.25,0.21,0.19,-3.29,1e3,1e5,5.73};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[5][2][i-1] = _itmp26[_r++];
		}
	}
	{ static double _itmp27[] = {332.46,-0.11,-9.95e-1,-1.58e-3,
	  1e1,1e5,11.30};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[5][3][i-1] = _itmp27[_r++];
		}
	}
	{ static double _itmp28[] = {1.01e-3,-0.29,-0.92,-8.38,1e2,
	  5e4,0.94};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[6][0][i-1] = _itmp28[_r++];
		}
	}
	{ static double _itmp29[] = {3.05e-1,0.60,2.65,-0.93,1e3,1e5,
	  4.56};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[6][1][i-1] = _itmp29[_r++];
		}
	}
	{ static double _itmp30[] = {4.54,0.57,-0.65,-0.89,1e1,1e5,
	  6.40};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[6][2][i-1] = _itmp30[_r++];
		}
	}
	/* N+4 recombination */
	{ static double _itmp31[] = { 2.95 , 0.55 , -0.39 , -1.07 , 1e3 , 1e6 ,
	  11.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[6][3][i-1] = _itmp31[_r++];
		}
	}
	{ static double _itmp32[] = {1.04,3.15e-2,-0.61,-9.73,1e1,1e4,
	  0.02};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[7][0][i-1] = _itmp32[_r++];
		}
	}
	{ static double _itmp33[] = {1.04,0.27,2.02,-5.92,1e2,1e5,6.65};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[7][1][i-1] = _itmp33[_r++];
		}
	}
	{ static double _itmp34[] = {3.98,0.26,0.56,-2.62,1e3,5e4,5.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[7][2][i-1] = _itmp34[_r++];
		}
	}
	{ static double _itmp35[] = {2.52e-1,0.63,2.08,-4.16,1e3,3e4,
	  8.47};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[7][3][i-1] = _itmp35[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[8][0][i-1] = 0.f;
		}
	{ static double _itmp36[] = {1.00e-5,0.,0.,0.,2e3,5e4,-21.37};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[8][1][i-1] = _itmp36[_r++];
		}
	}
	{ static double _itmp37[] = {9.86,0.29,-0.21,-1.15,2e3,5e4,
	  5.6};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[8][2][i-1] = _itmp37[_r++];
		}
	}
	{ static double _itmp38[] = {7.15e-1,1.21,-0.70,-0.85,2e3,5e4,
	  11.8};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[8][3][i-1] = _itmp38[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[9][0][i-1] = 0.f;
		}
	{ static double _itmp39[] = {1.00e-5,0.,0.,0.,5e3,5e4,-27.36};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[9][1][i-1] = _itmp39[_r++];
		}
	}
	{ static double _itmp40[] = {14.73,4.52e-2,-0.84,-0.31,5e3,
	  5e4,5.82};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[9][2][i-1] = _itmp40[_r++];
		}
	}
	{ static double _itmp41[] = {6.47,0.54,3.59,-5.22,1e3,3e4,8.60};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[9][3][i-1] = _itmp41[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[10][0][i-1] = 0.f;
		}
	{ static double _itmp42[] = {1.00e-5,0.,0.,0.,2e3,5e4,-33.68};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[10][1][i-1] = _itmp42[_r++];
		}
	}
	{ static double _itmp43[] = {1.33,1.15,1.20,-0.32,2e3,5e4,6.25};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[10][2][i-1] = _itmp43[_r++];
		}
	}
	{ static double _itmp44[] = {1.01e-1,1.34,10.05,-6.41,2e3,5e4,
	  11.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[10][3][i-1] = _itmp44[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[11][0][i-1] = 0.f;
		}
	{ static double _itmp45[] = {8.58e-5,2.49e-3,2.93e-2,-4.33,
	  1e3,3e4,1.44};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[11][1][i-1] = _itmp45[_r++];
		}
	}
	{ static double _itmp46[] = {6.49,0.53,2.82,-7.63,1e3,3e4,5.73};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[11][2][i-1] = _itmp46[_r++];
		}
	}
	{ static double _itmp47[] = {6.36,0.55,3.86,-5.19,1e3,3e4,8.60};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[11][3][i-1] = _itmp47[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[12][0][i-1] = 0.f;
		}
	{ static double _itmp48[] = {1.00e-5,0.,0.,0.,1e3,3e4,-5.23};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[12][1][i-1] = _itmp48[_r++];
		}
	}
	{ static double _itmp49[] = {7.11e-5,4.12,1.72e4,-22.24,1e3,
	  3e4,8.17};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[12][2][i-1] = _itmp49[_r++];
		}
	}
	{ static double _itmp50[] = {7.52e-1,0.77,6.24,-5.67,1e3,3e4,
	  8.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[12][3][i-1] = _itmp50[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[13][0][i-1] = 0.f;
		}
	/* Si+2 recombination */
	{ static double _itmp51[] = {6.77 , 7.36e-2 , -0.43 , -0.11 , 5e2 , 1e5 ,
	  2.72};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[13][1][i-1] = _itmp51[_r++];
		}
	}
	{ static double _itmp52[] = {4.90e-1,-8.74e-2,-0.36,-0.79,1e3,
	  3e4,4.23};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[13][2][i-1] = _itmp52[_r++];
		}
	}
	{ static double _itmp53[] = {7.58,0.37,1.06,-4.09,1e3,5e4,7.49};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[13][3][i-1] = _itmp53[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[14][0][i-1] = 0.f;
		}
	{ static double _itmp54[] = {1.74e-4,3.84,36.06,-0.97,1e3,3e4,
	  3.45};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[14][1][i-1] = _itmp54[_r++];
		}
	}
	{ static double _itmp55[] = {9.46e-2,-5.58e-2,0.77,-6.43,1e3,
	  3e4,7.29};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[14][2][i-1] = _itmp55[_r++];
		}
	}
	{ static double _itmp56[] = {5.37,0.47,2.21,-8.52,1e3,3e4,9.71};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[14][3][i-1] = _itmp56[_r++];
		}
	}
	{ static double _itmp57[] = {3.82e-7,11.10,2.57e4,-8.22,1e3,
	  1e4,-3.24};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[15][0][i-1] = _itmp57[_r++];
		}
	}
	{ static double _itmp58[] = {1.00e-5,0.,0.,0.,1e3,3e4,-9.73};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[15][1][i-1] = _itmp58[_r++];
		}
	}
	{ static double _itmp59[] = {2.29,4.02e-2,1.59,-6.06,1e3,3e4,
	  5.73};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[15][2][i-1] = _itmp59[_r++];
		}
	}
	{ static double _itmp60[] = {6.44,0.13,2.69,-5.69,1e3,3e4,8.60};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[15][3][i-1] = _itmp60[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[16][0][i-1] = 0.f;
		}
	{ static double _itmp61[] = {1.00e-5,0.,0.,0.,1e3,3e4,-10.21};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[16][1][i-1] = _itmp61[_r++];
		}
	}
	{ static double _itmp62[] = {1.88,0.32,1.77,-5.70,1e3,3e4,8.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[16][2][i-1] = _itmp62[_r++];
		}
	}
	{ static double _itmp63[] = {7.27,0.29,1.04,-10.14,1e3,3e4,
	  9.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[16][3][i-1] = _itmp63[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[17][0][i-1] = 0.f;
		}
	{ static double _itmp64[] = {1.00e-5,0.,0.,0.,1e3,3e4,-14.03};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[17][1][i-1] = _itmp64[_r++];
		}
	}
	{ static double _itmp65[] = {4.57,0.27,-0.18,-1.57,1e3,3e4,
	  5.73};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[17][2][i-1] = _itmp65[_r++];
		}
	}
	{ static double _itmp66[] = {6.37,0.85,10.21,-6.22,1e3,3e4,
	  8.60};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[17][3][i-1] = _itmp66[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[18][0][i-1] = 0.f;
		}
	{ static double _itmp67[] = {1.00e-5,0.,0.,0.,1e3,3e4,-18.02};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[18][1][i-1] = _itmp67[_r++];
		}
	}
	{ static double _itmp68[] = {4.76,0.44,-0.56,-0.88,1e3,3e4,
	  6.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[18][2][i-1] = _itmp68[_r++];
		}
	}
	{ static double _itmp69[] = {1.00e-5,0.,0.,0.,1e3,3e4,-47.3};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[18][3][i-1] = _itmp69[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[19][0][i-1] = 0.f;
		}
	{ static double _itmp70[] = {0.,0.,0.,0.,1e1,1e9,0.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[19][1][i-1] = _itmp70[_r++];
		}
	}
	{ static double _itmp71[] = {3.17e-2,2.12,12.06,-0.40,1e3,3e4,
	  6.6};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[19][2][i-1] = _itmp71[_r++];
		}
	}
	{ static double _itmp72[] = {2.68,0.69,-0.68,-4.47,1e3,3e4,
	  9.9};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[19][3][i-1] = _itmp72[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[20][0][i-1] = 0.f;
		}
	{ static double _itmp73[] = {0.,0.,0.,0.,1e1,1e9,0.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[20][1][i-1] = _itmp73[_r++];
		}
	}
	{ static double _itmp74[] = {7.22e-3,2.34,411.50,-13.24,1e3,
	  3e4,3.5};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[20][2][i-1] = _itmp74[_r++];
		}
	}
	{ static double _itmp75[] = {1.20e-1,1.48,4.00,-9.33,1e3,3e4,
	  10.61};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[20][3][i-1] = _itmp75[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[21][0][i-1] = 0.f;
		}
	{ static double _itmp76[] = {0.,0.,0.,0.,1e1,1e9,0.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[21][1][i-1] = _itmp76[_r++];
		}
	}
	{ static double _itmp77[] = {6.34e-1,6.87e-3,0.18,-8.04,1e3,
	  3e4,4.3};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[21][2][i-1] = _itmp77[_r++];
		}
	}
	{ static double _itmp78[] = {4.37e-3,1.25,40.02,-8.05,1e3,3e4,
	  5.3};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[21][3][i-1] = _itmp78[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[22][0][i-1] = 0.f;
		}
	{ static double _itmp79[] = {1.00e-5,0.,0.,0.,1e3,3e4,-1.05};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[22][1][i-1] = _itmp79[_r++];
		}
	}
	{ static double _itmp80[] = {5.12,-2.18e-2,-0.24,-0.83,1e3,
	  3e4,4.7};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[22][2][i-1] = _itmp80[_r++];
		}
	}
	{ static double _itmp81[] = {1.96e-1,-8.53e-3,0.28,-6.46,1e3,
	  3e4,6.2};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[22][3][i-1] = _itmp81[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[23][0][i-1] = 0.f;
		}
	{ static double _itmp82[] = {5.27e-1,0.61,-0.89,-3.56,1e3,3e4,
	  2.89};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[23][1][i-1] = _itmp82[_r++];
		}
	}
	{ static double _itmp83[] = {10.90,0.24,0.26,-11.94,1e3,3e4,
	  5.4};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[23][2][i-1] = _itmp83[_r++];
		}
	}
	{ static double _itmp84[] = {1.18,0.20,0.77,-7.09,1e3,3e4,6.6};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[23][3][i-1] = _itmp84[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[24][0][i-1] = 0.f;
		}
	{ static double _itmp85[] = {1.65e-1,6.80e-3,6.44e-2,-9.70,
	  1e3,3e4,2.04};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[24][1][i-1] = _itmp85[_r++];
		}
	}
	{ static double _itmp86[] = {14.20,0.34,-0.41,-1.19,1e3,3e4,
	  6.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[24][2][i-1] = _itmp86[_r++];
		}
	}
	{ static double _itmp87[] = {4.43e-1,0.91,10.76,-7.49,1e3,3e4,
	  7.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[24][3][i-1] = _itmp87[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[25][0][i-1] = 0.f;
		}
	{ static double _itmp88[] = {1.26,7.72e-2,-0.41,-7.31,1e3,1e5,
	  2.56};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[25][1][i-1] = _itmp88[_r++];
		}
	}
	{ static double _itmp89[] = {3.42,0.51,-2.06,-8.99,1e3,1e5,
	  6.3};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[25][2][i-1] = _itmp89[_r++];
		}
	}
	{ static double _itmp90[] = {14.60,3.57e-2,-0.92,-0.37,1e3,
	  3e4,10.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[25][3][i-1] = _itmp90[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[26][0][i-1] = 0.f;
		}
	{ static double _itmp91[] = {5.30,0.24,-0.91,-0.47,1e3,3e4,
	  2.9};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[26][1][i-1] = _itmp91[_r++];
		}
	}
	{ static double _itmp92[] = {3.26,0.87,2.85,-9.23,1e3,3e4,6.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[26][2][i-1] = _itmp92[_r++];
		}
	}
	{ static double _itmp93[] = {1.03,0.58,-0.89,-0.66,1e3,3e4,
	  10.51};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[26][3][i-1] = _itmp93[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[27][0][i-1] = 0.f;
		}
	{ static double _itmp94[] = {1.05,1.28,6.54,-1.81,1e3,1e5,3.0};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[27][1][i-1] = _itmp94[_r++];
		}
	}
	{ static double _itmp95[] = {9.73,0.35,0.90,-5.33,1e3,3e4,5.2};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[27][2][i-1] = _itmp95[_r++];
		}
	}
	{ static double _itmp96[] = {6.14,0.25,-0.91,-0.42,1e3,3e4,
	  10.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[27][3][i-1] = _itmp96[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[28][0][i-1] = 0.f;
		}
	{ static double _itmp97[] = {1.47e-3,3.51,23.91,-0.93,1e3,3e4,
	  3.44};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[28][1][i-1] = _itmp97[_r++];
		}
	}
	{ static double _itmp98[] = {9.26,0.37,0.40,-10.73,1e3,3e4,
	  5.6};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[28][2][i-1] = _itmp98[_r++];
		}
	}
	{ static double _itmp99[] = {11.59,0.20,0.80,-6.62,1e3,3e4,
	  9.};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[28][3][i-1] = _itmp99[_r++];
		}
	}
	for( i=1; i <= 7; i++ )
	{
		CTRecombData[29][0][i-1] = 0.f;
		}
	{ static double _itmp100[] = {1.00e-5,0.,0.,0.,1e3,3e4,-4.37};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[29][1][i-1] = _itmp100[_r++];
		}
	}
	{ static double _itmp101[] = {6.96e-4,4.24,26.06,-1.24,1e3,
	  3e4,7.8};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[29][2][i-1] = _itmp101[_r++];
		}
	}
	{ static double _itmp102[] = {1.33e-2,1.56,-0.92,-1.20,1e3,
	  3e4,11.73};
	for( i=1, _r = 0; i <= 7; i++ )
	{
		CTRecombData[29][3][i-1] = _itmp102[_r++];
	}
	}

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

/* ChargTranPun, punch charge transfer coef */
void ChargTranPun( FILE* ipPnunit , char* chPunch )
{
	long int j, jj;
	/* restore temp when done with this routine	*/
	float TempSave = phycon.te;

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

	/* this is the usual charge transfer option */
	if( strcmp( chPunch,"CHAR") == 0 )
	{
		/* charge exchange rate coefficients, entered with
		 * punch charge transfer command.  Queries Jim Kingdon's
		 * CT fits and routines to get H - He and higher CT rates
		 * rates are evaluated at the current temperature, which can
		 * be specified with the constant temperature commmand */
		/* first group is charge transfer recombination,
		 * the process A+x + H => A+x-1 + H+ */
		fprintf( ipPnunit, "#element\tion\n");
		for( j=1; j < LIMELM; j++ )
		{
			/* first number is atomic number, so add 1 to get onto physical scale */
			/* >>chng 00 may 09, caught by Jon Slavin */
			fprintf( ipPnunit, "%s\t", elementnames.chElementSym[j] );
			/*fprintf( ipPnunit, "%3ld\t", j+1 );*/

			for( jj=0; jj < j; jj++ )
			{
				fprintf( ipPnunit, "%.2e\t", 
				  HCTRecom(jj,j) );
			}
			fprintf( ipPnunit, "\n" );
		}

		/* second group is charge transfer ionization,
		 * the process A+x + H+ => A+x+1 + H */
		fprintf( ipPnunit, "\n#ionization rates, atomic number\n");
		for( j=1; j < LIMELM; j++ )
		{
			fprintf( ipPnunit, "%s\t", elementnames.chElementSym[j] );
			for( jj=0; jj < j; jj++ )
			{
				fprintf( ipPnunit, "%.2e\t", 
				  HCTIon(jj,j) );
			}
			fprintf( ipPnunit, "\n" );
		}
	}

	/* this is the charge transfer to produce output for agn,
	 * envoked with the punch charge transfer agn command */
	else if( strcmp( chPunch,"CHAG") == 0 )
	{
		/* this is boundary between two tables */
		double BreakEnergy = 100./13.0;
		double OHRate;
		float teinit = 5000.;
		float tefinal = 20000.;
		float te1;
		long int nelem, ion;

		te1 = teinit;
		fprintf(ipPnunit,"H ioniz\n X+i\\Te");
		while( te1 <= tefinal )
		{
			fprintf(ipPnunit,"\t%.0f K",te1);
			te1 *= 2.;
		}
		fprintf(ipPnunit,"\n");

		/* make sure rates already evaluated at least one time */
		ChargTranEval(&OHRate);

		/* loop over all elements, H charge transfer ionization */
		for( nelem=ipHELIUM; 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 )
				{
					/* skip high ionization CT */
					if( Heavy.Valence_IP_Ryd[nelem][ion] > BreakEnergy )
						break;
					/* most of these are actually zero */
					if( atmdat.HCharExcIonOf[nelem][ion] == 0 )
						continue;

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

					/* fully define the new temperature */
					phycon.te = teinit;
					tfidle(FALSE);
					while( phycon.te <= tefinal )
					{
						dense.IonLow[nelem] = 0;
						dense.IonHigh[nelem] = nelem+1;
						ChargTranEval(&OHRate);

						fprintf(ipPnunit,"\t%.2e",atmdat.HCharExcIonOf[nelem][ion]);
						phycon.te = phycon.te *2.f;
						tfidle(FALSE);
					}
					fprintf(ipPnunit,"\n");
				}
				fprintf(ipPnunit,"\n");
			}
		}

		te1 = teinit;
		fprintf(ipPnunit,"H recom\n X+i\\Te");
		while( te1 <= tefinal )
		{
			fprintf(ipPnunit,"\t%.0f K",te1);
			te1 *= 2.;
		}
		fprintf(ipPnunit,"\n");

		/* loop over all elements, H charge transfer recombination */
		for( nelem=ipHELIUM; 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 )
				{
					/* skip high ionization CT */
					if( Heavy.Valence_IP_Ryd[nelem][ion] > BreakEnergy )
						break;
					/* most of these are actually zero */
					if( atmdat.HCharExcRecTo[nelem][ion] == 0 )
						continue;

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

					/* fully define the new temperature */
					phycon.te = teinit;
					tfidle(FALSE);
					while( phycon.te <= tefinal )
					{
						dense.IonLow[nelem] = 0;
						dense.IonHigh[nelem] = nelem+1;
						ChargTranEval(&OHRate);

						fprintf(ipPnunit,"\t%.2e",atmdat.HCharExcRecTo[nelem][ion]);
						phycon.te = phycon.te *2.f;
						tfidle(FALSE);
					}
					fprintf(ipPnunit,"\n");
				}
				fprintf(ipPnunit,"\n");
			}
		}

#		if 0
		te1 = teinit;
		fprintf(ipPnunit,"He recom\n Elem\\Te");
		while( te1 <= tefinal )
		{
			fprintf(ipPnunit,"\t%.0f",te1);
			te1 *= 2.;
		}
		fprintf(ipPnunit,"\n");

		/* loop over all elements, H charge transfer recombination */
		for( nelem=ipHELIUM; 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 )
				{
					/* most of these are actually zero */
					if( atmdat.HeCharExcRecTo[nelem][ion] == 0 )
						continue;
					fprintf(ipPnunit,"%.2s%.2s", 
						elementnames.chElementSym[nelem],
						elementnames.chIonStage[ion]);

					/* fully define the new temperature */
					phycon.te = teinit;
					tfidle(FALSE);
					while( phycon.te <= tefinal )
					{
						dense.IonLow[nelem] = 0;
						dense.IonHigh[nelem] = nelem+1;
						ChargTranEval();

						fprintf(ipPnunit,"\t%.2e",atmdat.HeCharExcRecTo[nelem][ion]);
						phycon.te = phycon.te *2.;
						tfidle(FALSE);
					}
					fprintf(ipPnunit,"\n");
				}
				fprintf(ipPnunit,"\n");
			}
		}


		te1 = teinit;
		fprintf(ipPnunit,"He ioniz\n Elem\\Te");
		while( te1 <= tefinal )
		{
			fprintf(ipPnunit,"\t%.0f",te1);
			te1 *= 2.;
		}
		fprintf(ipPnunit,"\n");

		/* loop over all elements, H charge transfer recombination */
		for( nelem=ipHELIUM; 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 )
				{
					/* most of these are actually zero */
					if( atmdat.HeCharExcIonOf[nelem][ion] == 0 )
						continue;
					fprintf(ipPnunit,"%.2s%.2s", 
						elementnames.chElementSym[nelem],
						elementnames.chIonStage[ion]);

					/* fully define the new temperature */
					phycon.te = teinit;
					tfidle(FALSE);
					while( phycon.te <= tefinal )
					{
						dense.IonLow[nelem] = 0;
						dense.IonHigh[nelem] = nelem+1;
						ChargTranEval();

						fprintf(ipPnunit,"\t%.2e",atmdat.HeCharExcIonOf[nelem][ion]);
						phycon.te = phycon.te * 2.;
						tfidle(TRUE);
					}
					fprintf(ipPnunit,"\n");
				}
				fprintf(ipPnunit,"\n");
			}
		}
#		endif
	}
	else
	{
		fprintf( ioQQQ, " punch charge keyword insane\n" );
		puts( "[Stop in ChargTranPun]" );
		cdEXIT(EXIT_FAILURE);
	}

	phycon.te = TempSave;
	tfidle(FALSE);

#	ifdef DEBUG_FUN
	fputs( " <->ChargTranPun()\n", debug_fp );
#	endif
	return;
}
/*lint +e662 possible creation out of bounds pointer */
