/* 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 */
/*HydroCollid evaluate collision rate for model hydrogen atom */
/*HCSAR_interp interpolate on collision strengths */
/*C6cs123 line collision rates for lower levels of hydrogenic carbon, n=1,2,3 */
/*Ca20cs123 line collision rates for lower levels of hydrogenic calcium, n=1,2,3 */
/*Hydcs123 Hydrogenic de-excitation collision rates n=1,2,3 */
/*H1cs123 hydrogen collision data levels involving 1s,2s,2p,3. */
/*Ne10cs123 line collision rates for lower levels of hydrogenic neon, n=1,2,3 */
/*He2cs123 line collision strengths for lower levels of helium ion, n=1,2,3, by K Korista */
/*Fe26cs123 line collision rates for lower levels of hydrogenic iron, n=1,2,3 */
#include "cddefines.h"
/* the dec/compac ccc compiler becomes confused with pointer arith - disable warning */
#ifdef __DECC
#    pragma message disable (BADBOUNDCHK)
#endif
#include "physconst.h"
#include "taulines.h"
#include "iso.h"
#include "ionbal.h"
#include "trace.h"
#include "dense.h"
#include "phycon.h"
#include "rfield.h"
#include "tfidle.h"
#include "heavy.h"
#include "secondaries.h"
#include "opacity.h"
#include "converge.h"
#include "hydro_vs_rates.h"
#include "atmdat.h"
#include "hydrogenic.h"

/*Fe26cs123 line collision rates for lower levels of hydrogenic iron, n=1,2,3 */
static double Fe26cs123(long int i, 
  long int j);

/*He2cs123 line collision strengths for lower levels of helium ion, n=1,2,3, by K Korista */
static double He2cs123(long int i, 
  long int j);

/*H1cs123 hydrogen collision data levels involving 1s,2s,2p,3. */
static double H1cs123(long int i, 
  long int j, 
  long int chType);

/*Hydcs123 Hydrogenic de-excitation collision rates n=1,2,3 */
static double Hydcs123(long int ilow, 
  long int ihigh, 
  long int iz, 
  long int chType);

/*C6cs123 line collision rates for lower levels of hydrogenic carbon, n=1,2,3 */
static double C6cs123(long int i, 
  long int j);

/*Ca20cs123 line collision rates for lower levels of hydrogenic calcium, n=1,2,3 */
static double Ca20cs123(long int i, 
  long int j);

/*HCSAR_interp interpolate on collision strengths */
static float HCSAR_interp( int ipLo , int ipHi );

/*Ne10cs123 line collision rates for lower levels of hydrogenic neon, n=1,2,3 */
static double Ne10cs123(long int i, 
  long int j);

 /* >>refer	H1	cs	Anderson, H., Ballance, C.P., Badnell, N.R., & Summers, H.P., 
 * >>refercon	2000, J Phys B, 33, 1255; erratum, 2002 */
#define NHCSTE 8
#define NHCS 6
static float HCSTE[NHCSTE] = {5802.,11604.,34812.,58020.,116040.,174060.,232080.,290100.};
/* these were the original data, replaced with set below on 
>>chng 02 may 06, repair erratum on Anderson et al. H cs data */
#if 0
static float HCS[NHCSTE*NHCS*(NHCS-1)] =
{0., 0., 0., 0., 0., 0., 0., 0., /* 1s - 1s */
2.600E-01f, 2.960E-01f, 3.250E-01f, 3.370E-01f, 3.560E-01f, 3.680E-01f, 3.750E-01f, 3.800E-01f,/* 1s-2s*/
4.270E-01f, 5.360E-01f, 8.570E-01f, 1.150E+00f, 1.750E+00f, 2.130E+00f, 2.350E+00f, 2.460E+00f,/* 1s-2p*/
2.372E-01f, 2.605E-01f, 3.413E-01f, 4.104E-01f, 5.063E-01f, 5.351E-01f, 5.335E-01f, 5.236E-01f,/* 1s-3*/
9.990E-02f, 1.144E-01f, 1.564E-01f, 1.857E-01f, 2.263E-01f, 2.422E-01f, 2.445E-01f, 2.428E-01f,/* 1s-4*/
7.325E-02f, 8.263E-02f, 9.641E-02f, 1.064E-01f, 1.202E-01f, 1.233E-01f, 1.219E-01f, 1.191E-01f,/* 1s-5*/
0., 0., 0., 0., 0., 0., 0., 0., /* 2s - 1s */
0., 0., 0., 0., 0., 0., 0., 0., /* 2s - 2s */
0., 0., 0., 0., 0., 0., 0., 0., /* 2s - 2p */
5.860E+00f, 7.640E+00f, 1.419E+01f, 2.015E+01f, 3.156E+01f, 3.969E+01f, 4.598E+01f, 5.114E+01f,/* 2s - 3 */
2.276E+00f, 2.747E+00f, 4.365E+00f, 5.609E+00f, 7.803E+00f, 9.267E+00f, 1.038E+01f, 1.126E+01f,/* 2s - 4 */
1.610E+00f, 1.897E+00f, 2.421E+00f, 2.799E+00f, 3.512E+00f, 4.023E+00f, 4.450E+00f, 4.790E+00f,/* 2s - 5 */
0., 0., 0., 0., 0., 0., 0., 0., /* 2p - 1s */
0., 0., 0., 0., 0., 0., 0., 0., /* 2p - 2s */
0., 0., 0., 0., 0., 0., 0., 0., /* 2p - 2p */
2.228E+01f, 2.818E+01f, 5.026E+01f, 7.193E+01f, 1.178E+02f, 1.542E+02f, 1.828E+02f, 2.078E+02f,/* 2p - 3 */
8.772E+00f, 1.042E+01f, 1.600E+01f, 2.058E+01f, 2.898E+01f, 3.482E+01f, 3.923E+01f, 4.280E+01f,/* 2p - 4 */
6.161E+00f, 7.209E+00f, 9.091E+00f, 1.049E+01f, 1.308E+01f, 1.488E+01f, 1.626E+01f, 1.736E+01f,/* 2p - 5 */
0., 0., 0., 0., 0., 0., 0., 0., /* 3 - 1s */
0., 0., 0., 0., 0., 0., 0., 0., /* 3 - 2s */
0., 0., 0., 0., 0., 0., 0., 0., /* 3 - 2p */
0., 0., 0., 0., 0., 0., 0., 0., /* 3 - 3 */
8.648E+01f, 1.366E+02f, 3.252E+02f, 4.834E+02f, 7.715E+02f, 9.692E+02f, 1.123E+03f, 1.247E+03f,/* 3 - 4 */
6.229E+01f, 8.263E+01f, 1.255E+02f, 1.527E+02f, 1.951E+02f, 2.211E+02f, 2.400E+02f, 2.544E+02f,/* 3 - 5 */
0., 0., 0., 0., 0., 0., 0., 0., /* 4 - 1s */
0., 0., 0., 0., 0., 0., 0., 0., /* 4 - 2s */
0., 0., 0., 0., 0., 0., 0., 0., /* 4 - 2p */
0., 0., 0., 0., 0., 0., 0., 0., /* 4 - 3 */
0., 0., 0., 0., 0., 0., 0., 0., /* 4 - 4 */
4.625E+02f, 8.070E+02f, 1.707E+03f, 2.269E+03f, 3.122E+03f, 3.675E+03f, 4.083E+03f, 4.414E+03f /* 4 - 5 */
};
#endif

/*>>chng 02 may 06, repair erratum on Anderson et al. H cs data - only first few were updated
 * once it became clear that cs had not changed at any possible temperatures - changes where
 * at T > 1e5 */
static float HCS[NHCSTE*NHCS*(NHCS-1)] =
{0., 0., 0., 0., 0., 0., 0., 0., /* 1s - 1s */
2.600E-01f, 2.960E-01f, 3.260E-01f, 3.390E-01f, 3.730E-01f, 4.060E-01f, 4.360E-01f, 4.610E-01f,/* 1s-2s*/
4.290E-01f, 5.290E-01f, 8.530E-01f, 1.150E+00f, 1.810E+00f, 2.350E+00f, 2.810E+00f, 3.200E+00f,/* 1s-2p*/
2.372E-01f, 2.605E-01f, 3.413E-01f, 4.140E-01f, 5.500E-01f, 6.510E-01f, 7.300E-01f, 7.976E-01f,/* 1s-3*/
9.990E-02f, 1.144E-01f, 1.564E-01f, 1.857E-01f, 2.263E-01f, 2.422E-01f, 2.445E-01f, 2.428E-01f,/* 1s-4*/
7.325E-02f, 8.263E-02f, 9.641E-02f, 1.064E-01f, 1.202E-01f, 1.233E-01f, 1.219E-01f, 1.191E-01f,/* 1s-5*/
0., 0., 0., 0., 0., 0., 0., 0., /* 2s - 1s */
0., 0., 0., 0., 0., 0., 0., 0., /* 2s - 2s */
0., 0., 0., 0., 0., 0., 0., 0., /* 2s - 2p */
5.860E+00f, 7.640E+00f, 1.419E+01f, 2.015E+01f, 3.156E+01f, 3.969E+01f, 4.598E+01f, 5.114E+01f,/* 2s - 3 */
2.276E+00f, 2.747E+00f, 4.365E+00f, 5.609E+00f, 7.803E+00f, 9.267E+00f, 1.038E+01f, 1.126E+01f,/* 2s - 4 */
1.610E+00f, 1.897E+00f, 2.421E+00f, 2.799E+00f, 3.512E+00f, 4.023E+00f, 4.450E+00f, 4.790E+00f,/* 2s - 5 */
0., 0., 0., 0., 0., 0., 0., 0., /* 2p - 1s */
0., 0., 0., 0., 0., 0., 0., 0., /* 2p - 2s */
0., 0., 0., 0., 0., 0., 0., 0., /* 2p - 2p */
2.228E+01f, 2.818E+01f, 5.026E+01f, 7.193E+01f, 1.178E+02f, 1.542E+02f, 1.828E+02f, 2.078E+02f,/* 2p - 3 */
8.772E+00f, 1.042E+01f, 1.600E+01f, 2.058E+01f, 2.898E+01f, 3.482E+01f, 3.923E+01f, 4.280E+01f,/* 2p - 4 */
6.161E+00f, 7.209E+00f, 9.091E+00f, 1.049E+01f, 1.308E+01f, 1.488E+01f, 1.626E+01f, 1.736E+01f,/* 2p - 5 */
0., 0., 0., 0., 0., 0., 0., 0., /* 3 - 1s */
0., 0., 0., 0., 0., 0., 0., 0., /* 3 - 2s */
0., 0., 0., 0., 0., 0., 0., 0., /* 3 - 2p */
0., 0., 0., 0., 0., 0., 0., 0., /* 3 - 3 */
8.648E+01f, 1.366E+02f, 3.252E+02f, 4.834E+02f, 7.715E+02f, 9.692E+02f, 1.123E+03f, 1.247E+03f,/* 3 - 4 */
6.229E+01f, 8.263E+01f, 1.255E+02f, 1.527E+02f, 1.951E+02f, 2.211E+02f, 2.400E+02f, 2.544E+02f,/* 3 - 5 */
0., 0., 0., 0., 0., 0., 0., 0., /* 4 - 1s */
0., 0., 0., 0., 0., 0., 0., 0., /* 4 - 2s */
0., 0., 0., 0., 0., 0., 0., 0., /* 4 - 2p */
0., 0., 0., 0., 0., 0., 0., 0., /* 4 - 3 */
0., 0., 0., 0., 0., 0., 0., 0., /* 4 - 4 */
4.625E+02f, 8.070E+02f, 1.707E+03f, 2.269E+03f, 3.122E+03f, 3.675E+03f, 4.083E+03f, 4.414E+03f /* 4 - 5 */
};

void HydroCollid(
	/* element atomic number on c scale */
	long int nelem)
{
	long int i, 
	  ipHi, 
	  ipLo,
	  n;
	double CStemp, 
	  factor ,
	  ConvLTEPOP;
	/*lint -e785 too few initializers for aggregate, but it defaults to all zeros */
	static float told[LIMELM] = {0.};
	/*lint +e785 */

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

	/* check that we were called with valid charge */
	ASSERT( nelem >= 0);
	ASSERT( nelem < LIMELM );

	/* evaluate collision rate for model hydrogen atom */
	/* nelem is atomic number minus 1, so H is 0*/

	{
		/*@-redef@*/
		enum {AGN=FALSE};
		/*@+redef@*/
		if( AGN  )
		{
			float tear[3]={10000.f,15000.f,20000.f};
			for( i=0;i<3;++i)
			{
				phycon.te = tear[i];
				tfidle(FALSE);
				for( ipHi=1;ipHi<4;++ipHi )
				{
					for( ipLo=0; ipLo<ipHi; ++ipLo )
					{
						if( !(ipLo==1 && ipHi==2 ) )
						{
							fprintf(ioQQQ,"%.3f\t", HCSAR_interp(ipLo,ipHi) ) ;
						}
					}
				}
				fprintf(ioQQQ," \n");
			}
			cdEXIT(EXIT_SUCCESS);
		}
	}

	/***********************************************************************
	 *                                                                     *
	 * get secondary excitation by suprathermal electrons                  *
	 *                                                                     *
	 ***********************************************************************/

	if( nelem == 0 )
	{
		/* following are for hydrogen ATOMS, ions are in next else */
		/* secondary excitation of Lyman lines */
		secondaries.Hx12[0][ipH2p] = (float)(secondaries.x12tot*2./3./1.261);
		secondaries.Hx12[0][ipH2s] = (float)(secondaries.x12tot*1./3./1.261);
		secondaries.Hx12[0][3] = (float)(secondaries.x12tot*0.165/1.261);
		secondaries.Hx12[0][4] = (float)(secondaries.x12tot*0.059/1.261);

		/* rest are just a wild guess */
		for( i=5; i < iso.numLevels[ipH_LIKE][nelem]; i++ )
		{
			secondaries.Hx12[0][i] = (float)(secondaries.Hx12[0][i-1]/2.);
		}
	}
	else if( nelem==1 )
	{
		/* following for ions, guess from Janev et al. 
		 * only evaluated for He, and there will be a problem if
		 * hydrogenic helium not present, but hydrogenic species
		 * of heavier elements are present.  I can see how this
		 * would ever happen */
		secondaries.Hx12[1][ipH2p] = (float)(secondaries.x12tot/3.*2./3.);
		secondaries.Hx12[1][ipH2s] = (float)(secondaries.x12tot/3./3.);
		/* this is just a wild guess */
		for( i=3; i < iso.numLevels[ipH_LIKE][nelem]; i++ )
		{
			secondaries.Hx12[1][i] = (float)(secondaries.Hx12[1][i-1]/2.);
		}
	}

	/* skip most of this sub if temperature has not changed by much,
	 * the check on conv.nTotalIoniz is to make sure that we redo this
	 * on the very first call in a grid calc - it is 0 on the first call */
	/* loc1 bug 99 apr 08, tried commenting out first block so always reeval, no joy */
	if( fabs(phycon.te-told[nelem])/phycon.te < 0.00002 && 
		EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].ColUL > 0. && conv.nTotalIoniz )
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, 
				"       HydroCollid called nelem %li - no reeval Boltz fac, LTE dens\n",
				nelem );
		}
		/* (but still eval col below, search for 
		 * following code executed even when temperature did not change" ) */
	}
	else
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, 
				"       HydroCollid called nelem %li - will reeval Boltz fac, LTE dens\n",
				nelem );
		}
		told[nelem] = phycon.te;

		/***********************************************************************
		 *                                                                     *
		 * Boltzmann factors for all levels, ionization and                    *
		 * exciation between levels                                            *
		 *                                                                     *
		 ***********************************************************************/

		/* define boltz facs between adjacent levels, will go from ipHi to one below */
		for( ipHi=ipH2s; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
		{
			/* HdeltaTe is array of line temps, upper to lower
			 * their boltzmann factors for this temp,
			 * the array addressing for hlbolt is correct, and backwards from
			 * other similar variables - [up][lo] */
			iso.Boltzmann[ipH_LIKE][nelem][ipHi][ipHi-1] = 
				exp(-(double)EmisLines[ipH_LIKE][nelem][ipHi][ipHi-1].EnergyK/(double)phycon.te);
		}

		/* the highest level, up to the continuum */
		/* >>chng 01 sep 23, got rid of HCionTe */
		iso.ConBoltz[ipH_LIKE][nelem][iso.numLevels[ipH_LIKE][nelem]-1] =
			exp(-iso.xIsoLevNIonRyd[ipH_LIKE][nelem][iso.numLevels[ipH_LIKE][nelem]-1]/(double)phycon.te_ryd);

		/* define remaining line Boltzmann factors in terms of existing factors between levels */
		for( ipLo=ipH1s; ipLo < (iso.numLevels[ipH_LIKE][nelem] - 2); ipLo++ )
		{
			for( ipHi=ipLo + 2; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
			{
				iso.Boltzmann[ipH_LIKE][nelem][ipHi][ipLo] = 
					iso.Boltzmann[ipH_LIKE][nelem][ipHi-1][ipLo]*
				  iso.Boltzmann[ipH_LIKE][nelem][ipHi][ipHi-1];
			}
		}

		for( i=ipH1s; i < (iso.numLevels[ipH_LIKE][nelem] - 1); i++ )
		{
			iso.ConBoltz[ipH_LIKE][nelem][i] = 
				iso.Boltzmann[ipH_LIKE][nelem][iso.numLevels[ipH_LIKE][nelem]-1][i]*
			  iso.ConBoltz[ipH_LIKE][nelem][iso.numLevels[ipH_LIKE][nelem]-1];
			/* this is to prevent close-to-overflow problems in model atom */
			if( iso.ConBoltz[ipH_LIKE][nelem][i] < SMALLDOUBLE )
			{
				iso.ConBoltz[ipH_LIKE][nelem][i] = 0.;
			}
		}

		/***********************************************************************
		 *                                                                     *
		 * LTE abundances for all levels, ionization and                       *
		 * exciation between levels                                            *
		 *                                                                     *
		 ***********************************************************************/

		/* following factor is actually 4.1412957e-16 (old=4.14158E-16), 
		 * but e- stat weight is in */
		/* ConvLTEPOP = 2.0708e-16/phycon.te32;*/
		/* >>chng 99 jun 02, use codata and infinite mass nuc,
		 * actually correct reduced mass is below */
		/* HION_LTE_POP	is planck^2 / (2 pi m_e k ), must raise this to 3/2 when used */
		factor = HION_LTE_POP*dense.AtomicWeight[nelem]/
			(dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT) ;

		/* term in () is stat weight of electron * ion */
		ConvLTEPOP = pow(factor,1.5)/(2.*iso.stat_ion[ipH_LIKE])/phycon.te32;

		for( n=ipH1s; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
		{
			/* >>>chng 99 mar 16, had gone to gt 0 previous week, but problems
			 * with very very small hcbolt, change back to 1e-100 */
			if( iso.ConBoltz[ipH_LIKE][nelem][n] > 1e-100 )
			{
				iso.PopLTE[ipH_LIKE][nelem][n] = 
					iso.stat[ipH_LIKE][nelem][n] / iso.ConBoltz[ipH_LIKE][nelem][n]* ConvLTEPOP;
			}
			else
			{
				iso.PopLTE[ipH_LIKE][nelem][n] = 0.;
			}
		}

		/* now check for any zeros - if present then matrix cannot be used */
		iso.lgPopLTE_OK[ipH_LIKE][nelem] = TRUE;
		for( n=ipH1s; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
		{
			if( iso.PopLTE[ipH_LIKE][nelem][n] <= 0. )
			{
				iso.lgPopLTE_OK[ipH_LIKE][nelem] = FALSE;
				break;
			}
		}

		/***********************************************************************
		 *                                                                     *
		 * collisional ionization, and three body recombination                *
		 * Vriens and Smeets collisional ionization data set                   *
		 * use for hydrogen itselt                                             *
		 *                                                                     *
		 ***********************************************************************/
		if( nelem == ipHYDROGEN )
		{
			/* atomic hydrogen
			 * this is total to 2s and 2p, will be split below
			 * HColOn is set to 0 with hydrogen collisions off command, usually 1 */
			iso.ColIoniz[ipH_LIKE][nelem][ipH2p] = iso.lgColl_ionize[ipH_LIKE]* 
				hydro_vs_ioniz(ipH_LIKE, nelem, ipH2p);

			/* now do rest of levels, up to second highest */
			for( n=3; n < iso.numLevels[ipH_LIKE][nelem]-1; n++ )
			{
				/* HColOn is option to turn off collisions, "hydrogen collis off" comnd */
				iso.ColIoniz[ipH_LIKE][nelem][n] = iso.lgColl_ionize[ipH_LIKE] * 
					hydro_vs_ioniz(ipH_LIKE, nelem, n);
			}
			/* always leave highest level coupled to continuum - HColIonOn not there*/
			iso.ColIoniz[ipH_LIKE][nelem][iso.numLevels[ipH_LIKE][nelem]-1] = 
				hydro_vs_ioniz( ipH_LIKE, nelem, iso.numLevels[ipH_LIKE][nelem]-1 );
		}
		else
		{
			/* hydrogenic ions */
			/* this is total to 2s and 2p */
			iso.ColIoniz[ipH_LIKE][nelem][ipH2p] = Hion_coll_ioniz_ratecoef(ipH_LIKE,nelem,ipH2p)*iso.lgColl_ionize[ipH_LIKE];
			for( n=3; n < iso.numLevels[ipH_LIKE][nelem]-1; n++ )
			{
				/* iso.lgColl_ionize is option to turn off collisions, "atom h-like collis off" comnd */
				iso.ColIoniz[ipH_LIKE][nelem][n] = Hion_coll_ioniz_ratecoef(ipH_LIKE,nelem,n)*iso.lgColl_ionize[ipH_LIKE];
			}
			/* always leave highest level coupled to continuum - lgColl_ionize not there*/
			iso.ColIoniz[ipH_LIKE][nelem][iso.numLevels[ipH_LIKE][nelem]-1] = Hion_coll_ioniz_ratecoef(ipHYDROGEN,nelem,iso.numLevels[ipH_LIKE][nelem]-1);
		}

		/* total 2s and 2p was set in previous branch, now split into 2s and 2p */
		iso.ColIoniz[ipH_LIKE][nelem][ipH2p] *= 0.75;
		iso.ColIoniz[ipH_LIKE][nelem][ipH2s] = iso.ColIoniz[ipH_LIKE][nelem][ipH2p]/3.;

		/* collisional ionization from ground */
		iso.ColIoniz[ipH_LIKE][nelem][ipH1s] = atmdat_coll_ion(nelem+1,1,phycon.te);

		/* option to kill collisional ionization with Hydro collis off */
		iso.ColIoniz[ipH_LIKE][nelem][ipH1s] *= iso.lgColl_ionize[ipH_LIKE];
		iso.ColIoniz[ipH_LIKE][nelem][ipH2s] *= iso.lgColl_ionize[ipH_LIKE];
		iso.ColIoniz[ipH_LIKE][nelem][ipH2p] *= iso.lgColl_ionize[ipH_LIKE];

		/* hummer and storey option for case b kills coll ioniz from 1 and 2 */
		if( opac.lgCaseB_HummerStorey )
		{
			iso.ColIoniz[ipH_LIKE][nelem][ipH1s] = 0.;
			iso.ColIoniz[ipH_LIKE][nelem][ipH2p] = 0.;
			iso.ColIoniz[ipH_LIKE][nelem][ipH2s] = 0.;
		}

		/***********************************************************************
		 * 
		 * collisional strengths for all hydrogen lines
		 * 
		 ************************************************************************/

		/* evaluate collisional deexcitation rates,
		 * Hydcs123 evaluates downward collision strength, which is stored in
		 * hydrogen lines vector
		 * Hydcs123, upon call 1=1s, 2=2s, 3=2p, 4=3 */
		for( ipLo=ipH1s; ipLo < (iso.numLevels[ipH_LIKE][nelem] - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
			{
				if( ipHi == 1 )
				{
					/* 2s - 1s */
					if( nelem==0 )
					{
						CStemp = HCSAR_interp(ipLo,ipHi);
					}
					else
					{
						CStemp = Hydcs123(1,2,nelem,'e');
					}
				}
				else if( ipHi == 2 )
				{
					if( ipLo == 0 )
					{
						/* 2p - 1s */
						if( nelem==0 )
						{
							CStemp = HCSAR_interp(ipLo,ipHi);
						}
						else
						{
							CStemp = Hydcs123(1,3,nelem,'e');
						}
					}
					/* this is 2s - 2p */
					else if( ipLo == 1 )
					{
						/* option to kill 2s2p collisions */
						if( iso.lgColl_l_mixing[ipH_LIKE] )
						{
							/* want total collision strength for 2s - 2p - including
							 * both electrons and protons
							 * first is 2s - 2p for electrons ONLY 
							 * final term will be multiplied by edenHcor so fix here
							 * so that final is correct */
							CStemp = (Hydcs123(2,3,nelem,'e')*dense.eden + 
								Hydcs123(2,3,nelem,'p')*dense.xIonDense[ipHYDROGEN][1]) / dense.EdenHCorr;
#							if 0							 
							/* 2s - 2p for electrons ONLY */
							CStemp = Hydcs123(2,3,nelem,'e');
							/* 2s - 2p for protons ONLY */
							HCol2s2pp = Hydcs123(2,3,nelem,'p');
							/* trick, since will later be multiplied by electron,
							 * not proton density */
							/*CStemp += HCol2s2pp*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/
							/* NB - all collision rates will be multiplied by this number */
							CStemp += HCol2s2pp*dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr;
#							endif
						}
						else
						{
							CStemp = 0.;
						}
					}
					else
					{
						/* this is not possible */
						CStemp = -DBL_MAX;
					}
				}

				else if( ipHi == 3 )
				{
					/* 3 - 1s */
					if( nelem==0 )
					{
						CStemp = HCSAR_interp(ipLo,ipHi);
					}
					else
					{
						if( ipLo == 0 )
						{
							/* 3 - 1 */
							CStemp = Hydcs123(1,4,nelem,'e');
						}
						else if( ipLo == 1 )
						{
							/* 3 - 2s */
							CStemp = Hydcs123(2,4,nelem,'e');
						}
						else if( ipLo == 2 )
						{
							/* 3 - 2p */
							CStemp = Hydcs123(3,4,nelem,'e');
						}
						else
						{
							/* this is not possible */
							CStemp = -DBL_MAX;
						}
					}
				}
				else if( nelem==ipHYDROGEN && ipHi <= 5 )
				{
					CStemp = HCSAR_interp(ipLo,ipHi);
				}
				else
				{
					/* highly excited levels */
					CStemp = Hion_colldeexc_cs(ipHi,ipLo,nelem , ipH_LIKE );
				}
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].cs = (float)CStemp;
			}
		}

		/* save 2s-2p collisions since we always want to include, even when other
		 * collisions are turned off, since 2s-2p act as one level.  there is a
		 * separate no 2s2p command, which is implemented below */
		CStemp = EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].cs;

		/* kill all collisional excitation if this set with 
		 * atom h-like collisional excitation */
		if( iso.lgColl_excite[ipH_LIKE] == 0 )
		{
			for( ipLo=ipH1s; ipLo < (iso.numLevels[ipH_LIKE][nelem]-1); ipLo++ )
			{
				for( ipHi=ipLo + 1; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
				{
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].cs = 0.;
				}
			}
		}

		/* case b hummer&storey does not have any collisions from n=1 or 2 
		 * just kill those collisions here */
		if( opac.lgCaseB_HummerStorey )
		{
			for( ipLo=ipH1s; ipLo <= ipH2p; ipLo++ )
			{
				for( ipHi=3; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
				{
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].cs = 0.;
				}
			}
		}

		/* reset the 2s-2p collsions unless turned off with atom h-like collisions l-mixing
		 * command, since we usually want to include these even when line excit
		 * is turned off */
		if( iso.lgColl_l_mixing[ipH_LIKE] )
		{
			EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].cs = (float)CStemp;
		}

		/***********************************************************************
		 *                                                                     *
		 * collisional deexcitation for all lines                              *
		 *                                                                     *
		 ***********************************************************************/

		for( ipLo=ipH1s; ipLo < (iso.numLevels[ipH_LIKE][nelem] - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
			{
				/* collisional deexcitation rate */
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].ColUL = 
					(float)(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].cs/
				  phycon.sqrte*COLL_CONST/EmisLines[ipH_LIKE][nelem][ipHi][ipLo].gHi );
				/*This is downward collision rate */
			}
		}

		if( (trace.lgTrace && trace.lgIsoTraceFull[ipH_LIKE]) && (nelem == trace.ipIsoTrace[ipH_LIKE]) )
		{
			fprintf( ioQQQ, "       HydroCollid rate n->n-1:" );
			for( n=ipH1s; n < (iso.numLevels[ipH_LIKE][nelem] - 1); n++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n+1][n].ColUL ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid col ion cof:" );
			for( n=ipH1s; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e",  iso.ColIoniz[ipH_LIKE][nelem][n] ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dn2 1s:" );
			for( n=ipH2s; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][ipH1s].ColUL ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dn2 2s:" );
			for( n=ipH2p; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][ipH2s].ColUL ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dn2 2p:" );
			for( n=3; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][ipH2p].ColUL ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dwn2 3:" );
			for( n=4; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][3].ColUL ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dwn2 4:" );
			for( n=5; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][4].ColUL ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dwn2 5:" );
			for( n=6; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][5].ColUL) );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dwn2 6:" );
			for( n=7; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][6].ColUL ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid rate dwn2 7:" );
			for( n=8; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][7].ColUL) );
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid boltz fac c:" );
			for( n=ipH1s; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", iso.ConBoltz[ipH_LIKE][nelem][n] ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "       HydroCollid HLTE(n,Z):  " );
			for( n=ipH1s; n < iso.numLevels[ipH_LIKE][nelem]; n++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e", iso.PopLTE[ipH_LIKE][nelem][n] ));
			}
			fprintf( ioQQQ, "\n" );
		}
	}

	/* this must be reevaluated every time since eden can change when Te does not */
	/* save into main array - collisional ionization by thermal electrons */
	ionbal.CollIonRate_Ground[nelem][nelem][0] = 
		(iso.ColIoniz[ipH_LIKE][nelem][ipH1s]*dense.EdenHCorr) ;

	/* cooling due to collisional ionization, which only includes thermal electrons */
	ionbal.CollIonRate_Ground[nelem][nelem][1] = (iso.ColIoniz[ipH_LIKE][nelem][ipH1s]*dense.EdenHCorr*
		rfield.anu[Heavy.ipHeavy[nelem][nelem]-1]* EN1RYD);

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


#define HCSAR(ilo , ihi , nte ) (*(HCS+(nte)+(ilo)*48+(ihi)*8))

/*HCSAR_interp interpolate on collision strengths */
static float HCSAR_interp( int ipLo , int ipHi )
{

	static int ip=1;
	float cs;

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

	if( ipLo==1 && ipHi==2 )
	{
		fprintf(ioQQQ,"HCSAR_interp was called for the 2s-2p transition, which it cannot do\n");
		cdEXIT(EXIT_FAILURE);
	}
	if( phycon.te <= HCSTE[0] )
	{
		cs = HCSAR( ipLo , ipHi , 0 );
	}
	else if( phycon.te >= HCSTE[NHCSTE-1] )
	{
		cs = HCSAR( ipLo , ipHi , NHCSTE-1 );
	}
	else
	{
		/* the ip index is most likely correct since it points to the last temperature */
		if( !(HCSTE[ip-1] < phycon.te ) && ( phycon.te <= HCSTE[ip]) )
		{
			/* we must find the temperature in the array */
			for( ip=1; ip<NHCSTE ; ++ip )
			{
				if( (HCSTE[ip-1] < phycon.te ) && ( phycon.te <= HCSTE[ip]) )
					break;
			}
		}
		/* we now have the index */
		cs = HCSAR( ipLo , ipHi , ip-1 ) + 
			(HCSAR( ipLo , ipHi , ip ) - HCSAR( ipLo , ipHi , ip-1 ) ) / (HCSTE[ip]-HCSTE[ip-1] ) *
			(phycon.te - HCSTE[ip-1] );
		if( cs <= 0.)
		{
			fprintf(ioQQQ," insane cs returned by HCSAR_interp, values are\n");
			fprintf(ioQQQ,"%.3f %.3f \n", HCSAR( ipLo , ipHi , ip-1 ),HCSAR( ipLo , ipHi , ip ) );
		}
	}

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

/*Hydcs123 Hydrogenic de-excitation collision strengths bewteen levels n=1,2,3,
 * for any charge.  routine only called by hydrocollid to fill in hydroline arrays
 * with collision strengths */
static double Hydcs123(
	 /* lower principal quantum number, 1, 2, or 3, in this routine 
	  * 1 is 1s, 2 is 2s, 3 is 2p, and 4 is 3 */
	 long int ipLow, 
	 /* upper principal quantum nubmer, 2, 3, or 4 */
	 long int ipHi, 
	 /* charge, 0 for hydrogen, 1 for helium, etc */
	 long int nelem, 
	 /* = 'e' for electron collisions, ='p' for proton */
	 long int chType)
{
	long int i, 
	  j, 
	  k, 
	  _r;
	double C, 
	  D, 
	  EE, 
	  expq ,
	  Hydcs123_v, 
	  Ratehigh, 
	  Ratelow, 
	  TeUse, 
	  gLo, 
	  gHi, 
	  q, 
	  rate, 
	  slope, 
	  temp, 
	  temphigh, 
	  templow, 
	  tev, 
	  x, 
	  QuanNLo, 
	  QuanNUp, 
	  Charge, 
	  ChargeSquared, 
	  zhigh, 
	  zlow;
	static double ae[5], 
	  ap[5], 
	  be[5], 
	  bp[5], 
	  ce[5], 
	  cp[5], 
	  de[5], 
	  dp[5], 
	  ee[5], 
	  ep[5];
	static double A[2]={4.4394,0.0};
	static double B[2]={0.8949,0.8879};
	static double C0[2]={-0.6012,-0.2474};
	static double C1[2]={-3.9710,-3.7562};
	static double C2[2]={-4.2176,2.0491};
	static double D0[2]={2.930,0.0539};
	static double D1[2]={1.7990,3.4009};
	static double D2[2]={4.9347,-1.7770};
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp0[] = {-2113.113,729.0084,1055.397,854.632,
		  938.9912};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			ap[i-1] = _itmp0[_r++];
			}
		}
		{ static double _itmp1[] = {-6783.515,-377.7190,724.1936,493.1107,
		  735.7466};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			bp[i-1] = _itmp1[_r++];
			}
		}
		{ static double _itmp2[] = {-3049.719,226.2320,637.8630,388.5465,
		  554.6369};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			cp[i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {3514.5153,88.60169,-470.4055,-329.4914,
		  -450.8459};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			dp[i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {0.005251557,0.009059154,0.008725781,
		  0.009952418,0.01098687};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			ep[i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {-767.5859,-643.1189,-461.6836,-429.0543,
		  -406.5285};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			ae[i-1] = _itmp5[_r++];
			}
		}
		{ static double _itmp6[] = {-1731.9178,-1442.548,-1055.364,
		  -980.3079,-930.9266};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			be[i-1] = _itmp6[_r++];
			}
		}
		{ static double _itmp7[] = {-939.1834,-789.9569,-569.1451,-530.1974,
		  -502.0939};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			ce[i-1] = _itmp7[_r++];
			}
		}
		{ static double _itmp8[] = {927.4773,773.2008,564.3272,524.2944,
		  497.7763};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			de[i-1] = _itmp8[_r++];
			}
		}
		{ static double _itmp9[] = {-0.002528027,-0.003793665,-0.002122103,
		  -0.002234207,-0.002317720};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			ee[i-1] = _itmp9[_r++];
			}
		}
		_aini = 0;
	}

#	ifdef DEBUG_FUN
	fputs( "<+>Hydcs123()\n", debug_fp );
#	endif
	/* Hydrogenic de-excitation collision rates n=1,2,3 
	 * >>refer	h1	cs	Callaway, J. 1983, Phys Let A, 96, 83
	 * >>refer	h1	cs	Zygelman, B., & Dalgarno, A. 1987, Phys Rev A, 35, 4085 
	 * for 2p-2s only.
	 * The fit from Callaway is in nuclear charge for 1s - 2s,2p only.
	 * For transtions involving level 3, interpolation in Z involving
	 * the functions He2cs123,C6cs123,Ne10cs123,Ca20cs123, Fe26cs123.
	 *
	 * The fits from ZD are for 2p-2s for Z=2,6,12,16,18 only other charges are
	 * interpolated, both electron and proton rates are included,
	 * the variable chType is either 'e' or 'p'.
	 *
	 * ipLow is the lower level and runs from 1 to 3 (1s, 2s, 2p)
	 * ipHi is the upper level and runs from 2 to 4 (2s, 2p, 3) */

	/* for Callaway fit: */
	/* for Zygelman and Dalgarno: */

	/* first entry is 2p, then 2s */

	/* fit in nuclear charge Z for 2p-2s collisions in hydrogenic species
	 * equation is a+bx+cx^2ln(x)+dexp(x)+eln(x)/x^2, where x=te/Z^2 in a.u.
	 * first are the proton rates: */
	/* these are electron rates: */

	/* charge is on scale with H=1, He=2, etc */
	ASSERT( nelem >= 0 );
	ASSERT( nelem < LIMELM );

	/* these are the pointers to upper and lower levels.  1=1s, 2=2s, 3=2p, 4=3 */
	ASSERT( ipLow > 0);
	ASSERT( ipLow <= 3);
	ASSERT( ipHi > 1 );
	ASSERT( ipHi <=4 );

	/* set quantum numbers and stat. weights of the transitions: */
	if( ipHi == 4 )
	{
		/* upper is n=3 then set level, stat. weight */
		QuanNUp = 3.;
		gHi = 18.;
		/* following will be set here even though it is not used in this case,
		 * to prevent good compilers from falsing on i not set,
		 * there is assert when used to make sure it is ok */
		i = -1;
	}
	else if( ipHi == 3 )
	{
		/* upper is nl=2p then set level, stat. weight */
		QuanNUp = 2.;
		gHi = 6.;
		/* used to point within vectors defined above */
		i = 0;
	}
	else if( ipHi == 2 )
	{
		/* upper is nl=2s then set level, stat. weight */
		QuanNUp = 2.;
		gHi = 2.;
		/* used to point within vectors defined above */
		i = 1;
	}
	else
	{
		fprintf( ioQQQ, " Insane levels in Hydcs123\n" );
		puts( "[Stop in hydcs123]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* which lower level? */
	if( ipLow == 1 )
	{
		/* lower is n=1 then set level, stat. weight */
		QuanNLo = 1.;
		gLo = 2.;
	}
	else if( ipLow == 2 )
	{
		/* lower is nl=2s then set level, stat. weight */
		QuanNLo = 2.;
		gLo = 2.;
	}
	else if( ipLow == 3 )
	{
		QuanNLo = 2.;
		gLo = 6.;
	}
	else
	{
		fprintf( ioQQQ, " Insane levels in Hydcs123\n" );
		puts( "[Stop in hydcs123]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( nelem == 0 )
	{
		/* atomic hydrogen, use special routine */
		Hydcs123_v = H1cs123(ipLow,ipHi,chType);
		
#		ifdef DEBUG_FUN
		fputs( " <->Hydcs123()\n", debug_fp );
#		endif
		return( Hydcs123_v );
	}

	/* following is for charged species, since H itself returned above */

	/* following is the physical charge */
	Charge = (double)(nelem + 1);
	/* square of charge */
	ChargeSquared = Charge*Charge;

	if( ipLow == 2 && ipHi == 3 )
	{
		/*************************************** this is 2p-2s:
		 * series of if statements determines which entries Charge is between. */
		if( nelem < 1 )
		{
			/* this can't happen since routine returned above when ip=1 and 
			 * special atomic hydrogen routine called */
			fprintf( ioQQQ, " insane charge given to Hydcs123\n" );
			puts( "[Stop in hydcs123]" );
			cdEXIT(EXIT_FAILURE);
		}
		else if( nelem == 1 )
		{
			zlow = 2.;
			j = 1;
			zhigh = 2.;
			k = 1;
		}
		/* Li through C */
		else if( nelem <= 5 )
		{
			zlow = 2.;
			j = 1;
			zhigh = 6.;
			k = 2;
		}
		else if( nelem <= 11 )
		{
			zlow = 6.;
			j = 2;
			zhigh = 12.;
			k = 3;
		}
		else if( nelem <= 15 )
		{
			zlow = 12.;
			j = 3;
			zhigh = 16.;
			k = 4;
		}
		else if( nelem <= 17 )
		{
			zlow = 16.;
			j = 4;
			zhigh = 18.;
			k = 5;
		}
		/* following changed to else from else if, 
		 * to prevent false comment in good compilers */
		/*else if( nelem > 18 )*/
		else 
		{
			zlow = 18.;
			j = 5;
			zhigh = 18.;
			k = 5;
		}

		/* convert Te to a.u./Z^2
		 * determine rate at the low Charge */
		x = EVRYD/TE1RYD*phycon.te/27.211396/zlow/zlow;
		TeUse = MIN2(x,0.80);
		x = MAX2(0.025,TeUse);

		/* what type of collision are we dealing with? */
		if( chType == 'e' )
		{
			/* electron collisions */
			Ratelow = ae[j-1] + be[j-1]*x + ce[j-1]*x*x*log(x) + de[j-1]*
			  exp(x) + ee[j-1]*log(x)/x/x;
		}
		else if( chType == 'p' )
		{
			Ratelow = ap[j-1] + bp[j-1]*x + cp[j-1]*x*x*log(x) + dp[j-1]*
			  exp(x) + ep[j-1]*log(x)/x/x;
		}
		else 
		{
			/* this can't happen */
			fprintf( ioQQQ, " insane collision species given to Hydcs123\n" );
			puts( "[Stop in hydcs123]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* determine rate at the high Charge */
		x = EVRYD/TE1RYD*phycon.te/27.211396/zhigh/zhigh;
		TeUse = MIN2(x,0.80);
		x = MAX2(0.025,TeUse);
		if( chType == 'e' )
		{
			Ratehigh = ae[k-1] + be[k-1]*x + ce[k-1]*x*x*log(x) + 
			  de[k-1]*exp(x) + ee[k-1]*log(x)/x/x;
		}
		else
		{
			Ratehigh = ap[k-1] + bp[k-1]*x + cp[k-1]*x*x*log(x) + 
			  dp[k-1]*exp(x) + ep[k-1]*log(x)/x/x;
		}
		/* linearly interpolate in charge */
		/*lint -e777 test floats for equality */
		if( zlow == zhigh )
		/*lint +e777 test floats for equality */
		{
			rate = Ratelow;
		}
		else
		{
			slope = (Ratehigh - Ratelow)/(zhigh - zlow);
			rate = slope*(Charge - zlow) + Ratelow;
		}
		rate = rate/ChargeSquared/Charge*1.0e-7;
		/* must convert to cs and need to know the valid temp range */
		templow = 0.025*27.211396*TE1RYD/EVRYD*ChargeSquared;
		temphigh = 0.80*27.211396*TE1RYD/EVRYD*ChargeSquared;
		TeUse = MIN2((double)phycon.te,temphigh);
		temp = MAX2(TeUse,templow);
		Hydcs123_v = rate*gHi*sqrt(temp)/COLL_CONST;

	}
	else if( ipHi == 4 )
	{
		/* n = 3
		 * for the rates involving n = 3, must do something different. */
		if( nelem < 1 )
		{
			fprintf( ioQQQ, " insane charge given to Hydcs123\n" );
			puts( "[Stop in hydcs123]" );
			cdEXIT(EXIT_FAILURE);
		}
		else if( nelem == 1 )
		{
			zlow = 2.;
			Ratelow = He2cs123(ipLow,ipHi);
			zhigh = 2.;
			Ratehigh = Ratelow;
		}
		else if( nelem > 1 && nelem <= 5 )
		{
			zlow = 2.;
			Ratelow = He2cs123(ipLow,ipHi);
			zhigh = 6.;
			Ratehigh = C6cs123(ipLow,ipHi);
		}
		else if( nelem > 5 && nelem <= 9 )
		{
			zlow = 6.;
			Ratelow = C6cs123(ipLow,ipHi);
			zhigh = 10.;
			Ratehigh = Ne10cs123(ipLow,ipHi);
		}
		else if( nelem > 9 && nelem <= 19 )
		{
			zlow = 10.;
			Ratelow = Ne10cs123(ipLow,ipHi);
			zhigh = 20.;
			Ratehigh = Ca20cs123(ipLow,ipHi);
		}
		else if( nelem > 19 && nelem <= 25 )
		{
			zlow = 20.;
			Ratelow = Ca20cs123(ipLow,ipHi);
			zhigh = 26.;
			Ratehigh = Fe26cs123(ipLow,ipHi);
		}
		/*>>>chng 98 dec 17, to else to stop comment from good compilers*/
		/*else if( nelem > 26 )*/
		else
		{
			Charge = 26.;
			zlow = 26.;
			Ratelow = Fe26cs123(ipLow,ipHi);
			zhigh = 26.;
			Ratehigh = Ratelow;
		}

		/* linearly interpolate */
		/*lint -e777 test floats for equality */
		if( zlow == zhigh )
		/*lint +e777 test floats for equality */
		{
			rate = Ratelow;
		}
		else
		{
			slope = (Ratehigh - Ratelow)/(zhigh - zlow);
			rate = slope*(Charge - zlow) + Ratelow;
		}
		Hydcs123_v = rate;

	}
	else
	{
		/* this branch 1-2s, 1-2p */
		if( nelem == 1 )
		{
			/* this brance for helium, then return */
			Hydcs123_v = He2cs123(ipLow,ipHi);
			
#			ifdef DEBUG_FUN
			fputs( " <->Hydcs123()\n", debug_fp );
#			endif
			return( Hydcs123_v );
		}

		/* electron temperature in eV */
		tev = phycon.te / EVDEGK;
		/* energy in eV for hydrogenic species and these quantum numbers */
		EE = ChargeSquared*EVRYD*(1./QuanNLo/QuanNLo - 1./QuanNUp/QuanNUp);
		/* EE/kT for this transion */
		q = EE/tev;
		TeUse = MIN2(q,10.);
		/* q is now EE/kT but between 1 and 10 */
		q = MAX2(1.,TeUse);
		expq = exp(q);

		/* i must be 0 or 1 */
		ASSERT( i==0 || i==1 );
		C = C0[i] + C1[i]/Charge + C2[i]/ChargeSquared;
		D = D0[i] + D1[i]/Charge + D2[i]/ChargeSquared;

		/* following code changed so that ee1 always returns e1,
		 * orifinal version only returned e1 for x < 1 */
		/* use disabled e1: */
		/*if( q < 1. )*/
		/*{*/
			/*rate = (B[i-1] + D*q)*exp(-q) + (A[i-1] + C*q - D*q*q)**/
			  /*ee1(q);*/
		/*}*/
		/*else*/
		/*{*/
			/*rate = (B[i-1] + D*q) + (A[i-1] + C*q - D*q*q)*ee1(q);*/
		/*}*/
		/*rate *= 8.010e-8/2./ChargeSquared/tev*sqrt(tev);*/
		/* convert to de-excitation */
		/*if( q < 1. )*/
		/*{*/
			/*rate = rate*exp(q)*gLo/gHi;*/
		/*}*/
		/*else*/
		/*{*/
			/*rate = rate*gLo/gHi;*/
		/*}*/

		/*>>>chng 98 dec 17, ee1 always returns e1 */
		rate = (B[i] + D*q)/expq + (A[i] + C*q - D*q*q)*
			  ee1(q);
		rate *= 8.010e-8/2./ChargeSquared/tev*sqrt(tev);
		/* convert to de-excitation */
		rate *= expq*gLo/gHi;

		/* convert to cs */
		Hydcs123_v = rate*gHi*phycon.sqrte/COLL_CONST;
	}

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

/*C6cs123 line collision rates for lower levels of hydrogenic carbon, n=1,2,3 */
static double C6cs123(long int i, 
  long int j)
{
	double C6cs123_v, 
	  TeUse, 
	  t, 
	  x;
	static double a[3]={-92.23774,-1631.3878,-6326.4947};
	static double b[3]={-11.93818,-218.3341,-849.8927};
	static double c[3]={0.07762914,1.50127,5.847452};
	static double d[3]={78.401154,1404.8475,5457.9291};
	static double e[3]={332.9531,5887.4263,22815.211};

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

	/* These are fits to Table 5 of
	 * >>refer	c6	cs	Aggarwal, K.M., & Kingston, A.E. 1991, J Phys B, 24, 4583
	 * C VI collision rates for 1s-3, 2s-3, and 2p-3, 
	 * principal quantum numbers n and l.
	 *
	 * i is the lower level and runs from 1 to 3 (1s, 2s, 2p)
	 * j is the upper level and runs from 2 to 4 (2s, 2p, 3)
	 * 1s-2s,2p is not done here.
	 * check temperature: fits only good between 3.8 < log Te < 6.2
	 */
	/* arrays for fits of 3 transitions see the code below for key: */

	TeUse = MAX2(phycon.te,6310.);
	t = MIN2(TeUse,1.6e6);
	x = log10(t);

	if( i == 1 && j == 2 )
	{
		/* 1s - 2s (first entry) */
		fprintf( ioQQQ, " Carbon VI 2s-1s not done in C6cs123\n" );
		puts( "[Stop in c6cs123]" );
		cdEXIT(EXIT_FAILURE);
	}

	else if( i == 1 && j == 3 )
	{
		/* 1s - 2p (second entry) */
		fprintf( ioQQQ, " Carbon VI 2p-1s not done in C6cs123\n" );
		puts( "[Stop in c6cs123]" );
		cdEXIT(EXIT_FAILURE);
	}

	else if( i == 1 && j == 4 )
	{
		/* 1s - 3 (first entry) */
		C6cs123_v = a[0] + b[0]*x + c[0]*x*x*sqrt(x) + d[0]*log(x) + 
		  e[0]*log(x)/x/x;
		/* C6Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else if( i == 2 && j == 4 )
	{
		/* 2s - 3 (second entry)         */
		C6cs123_v = a[1] + b[1]*x + c[1]*x*x*sqrt(x) + d[1]*log(x) + 
		  e[1]*log(x)/x/x;
		/* C6Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else if( i == 3 && j == 4 )
	{
		/* 2p - 3s (third entry) */
		C6cs123_v = a[2] + b[2]*x + c[2]*x*x*sqrt(x) + d[2]*log(x) + 
		  e[2]*log(x)/x/x;
		/* C6Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else
	{

		fprintf( ioQQQ, "  insane levels for C VI n=1,2,3 !!!\n" );
		puts( "[Stop in c6cs123]" );
		cdEXIT(EXIT_FAILURE);
	}

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

/*Ca20cs123 line collision rates for lower levels of hydrogenic calcium, n=1,2,3 */
static double Ca20cs123(long int i, 
  long int j)
{
	double Ca20cs123_v, 
	  TeUse, 
	  t, 
	  x;
	static double a[3]={-12.5007,-187.2303,-880.18896};
	static double b[3]={-1.438749,-22.17467,-103.1259};
	static double c[3]={0.008219688,0.1318711,0.6043752};
	static double d[3]={10.116516,153.2650,717.4036};
	static double e[3]={45.905343,685.7049,3227.2836};

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

	/* 
	 * These are fits to Table 5 of
	 * >>refer	ca20	cs	Aggarwal, K.M., & Kingston, A.E. 1992, J Phys B, 25, 751
	 * Ca XX collision rates for 1s-3, 2s-3, and 2p-3, 
	 * principal quantum numbers n and l.
	 *
	 * i is the lower level and runs from 1 to 3 (1s, 2s, 2p)
	 * j is the upper level and runs from 2 to 4 (2s, 2p, 3)
	 * 1s-2s,2p is not done here.
	 * check temperature: fits only good between 5.0 < log Te < 7.2
	 */

	/* arrays for fits of 3 transitions see the code below for key: */

	TeUse = MAX2(phycon.te,1.0e5);
	t = MIN2(TeUse,1.585e7);
	x = log10(t);

	if( i == 1 && j == 2 )
	{
		/* 1s - 2s (first entry) */
		fprintf( ioQQQ, " Ca XX 2s-1s not done in Ca20cs123\n" );
		puts( "[Stop in ca20cs123]" );
		cdEXIT(EXIT_FAILURE);
	}

	else if( i == 1 && j == 3 )
	{
		/* 1s - 2p (second entry) */
		fprintf( ioQQQ, " Ca XX 2p-1s not done in Ca20cs123\n" );
		puts( "[Stop in ca20cs123]" );
		cdEXIT(EXIT_FAILURE);
	}

	else if( i == 1 && j == 4 )
	{
		/* 1s - 3 (first entry) */
		Ca20cs123_v = a[0] + b[0]*x + c[0]*x*x*sqrt(x) + d[0]*log(x) + 
		  e[0]*log(x)/x/x;
		/* Ca20Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else if( i == 2 && j == 4 )
	{
		/* 2s - 3 (second entry)         */
		Ca20cs123_v = a[1] + b[1]*x + c[1]*x*x*sqrt(x) + d[1]*log(x) + 
		  e[1]*log(x)/x/x;
		/* Ca20Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else if( i == 3 && j == 4 )
	{
		/* 2p - 3s (third entry) */
		Ca20cs123_v = a[2] + b[2]*x + c[2]*x*x*sqrt(x) + d[2]*log(x) + 
		  e[2]*log(x)/x/x;
		/* Ca20Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else
	{

		fprintf( ioQQQ, "  insane levels for Ca XX n=1,2,3 !!!\n" );
		puts( "[Stop in ca20cs123]" );
		cdEXIT(EXIT_FAILURE);
	}

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

/*Ne10cs123 line collision rates for lower levels of hydrogenic neon, n=1,2,3 */
static double Ne10cs123(long int i, 
  long int j)
{
	double Ne10cs123_v, 
	  TeUse, 
	  t, 
	  x;
	static double a[3]={3.346644,151.2435,71.7095};
	static double b[3]={0.5176036,20.05133,13.1543};
	static double c[3]={-0.00408072,-0.1311591,-0.1099238};
	static double d[3]={-3.064742,-129.8303,-71.0617};
	static double e[3]={-11.87587,-541.8599,-241.2520};

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

	/*  These are fits to Table 5 of
	 *  >>refer	ne10	cs	Aggarwal, K.M., & Kingston, A.E. 1991, PhyS, 44, 517
	 *  Ne X collision rates for 1s-3, 2s-3, and 2p-3, 
	 *  principal quantum numbers n and l.
	 *
	 *  i is the lower level and runs from 1 to 3 (1s, 2s, 2p)
	 *  j is the upper level and runs from 2 to 4 (2s, 2p, 3)
	 *  1s-2s,2p is not done here.
	 *  check temperature: fits only good between 3.8 < log Te < 6.2
	 * */
	/* arrays for fits of 3 transitions see the code below for key: */

	TeUse = MAX2(phycon.te,6310.);
	t = MIN2(TeUse,1.6e6);
	x = log10(t);

	if( i == 1 && j == 2 )
	{
		/* 1s - 2s (first entry) */
		fprintf( ioQQQ, " Neon X 2s-1s not done in Ne10cs123\n" );
		puts( "[Stop in ne10cs123]" );
		cdEXIT(EXIT_FAILURE);
	}

	else if( i == 1 && j == 3 )
	{
		/* 1s - 2p (second entry) */
		fprintf( ioQQQ, " Neon X 2p-1s not done in Ne10cs123\n" );
		puts( "[Stop in ne10cs123]" );
		cdEXIT(EXIT_FAILURE);
	}

	else if( i == 1 && j == 4 )
	{
		/* 1s - 3 (first entry) */
		Ne10cs123_v = a[0] + b[0]*x + c[0]*x*x*sqrt(x) + d[0]*log(x) + 
		  e[0]*log(x)/x/x;
		/*  Ne10Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else if( i == 2 && j == 4 )
	{
		/* 2s - 3 (second entry)         */
		Ne10cs123_v = a[1] + b[1]*x + c[1]*x*x*sqrt(x) + d[1]*log(x) + 
		  e[1]*log(x)/x/x;
		/*  Ne10Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else if( i == 3 && j == 4 )
	{
		/* 2p - 3s (third entry) */
		Ne10cs123_v = a[2] + b[2]*x + c[2]*x*x*sqrt(x) + d[2]*log(x) + 
		  e[2]*log(x)/x/x;
		/*  Ne10Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else
	{

		fprintf( ioQQQ, "  insane levels for Ne X n=1,2,3 !!!\n" );
		puts( "[Stop in ne10cs123]" );
		cdEXIT(EXIT_FAILURE);
	}

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

/*H1cs123 hydrogen collision data levels involving 1s,2s,2p,3. */
static double H1cs123(long int i, 
  long int j, 
  long int chType)
{
	long int k, 
	  l, 
	  _r;
	double H1cs123_v, 
	  T, 
	  TeUse, 
	  Tmax1, 
	  Tmax2, 
	  Tmax3, 
	  Tmin, 
	  bk, 
	  cs, 
	  cs3d, 
	  cs3p, 
	  cs3s, 
	  ev, 
	  rate, 
	  ryd, 
	  te03, 
	  te05, 
	  te10, 
	  te20;
	static double b[10][11];
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp0[] = {2.5501041e-01,9.973588e-02,-5.404749e-02,
		  2.520311e-02,-4.05306e-03,3.508995e-02,28.056,7.2945,0.0,
		  0.0,3.6105477e-01,1.73103797,3.7297254e-01,-4.3240114e-01,
		  1.4296095e-01,-1.670572e-02,0.0,0.0,0.0,0.0,6.461453e-02,
		  -1.1160826e-01,7.6957269e-01,-2.18918875,3.14216481,-2.23498961,
		  6.2651998e-01,1.12515e-02,28.056,7.2945,9.895325e-02,2.3240113e-01,
		  2.0112555e-01,8.7813331e-01,-3.34239518,3.7075526,-1.37320827,
		  0.0,0.0,0.0,5.742134e-02,-6.09531e-03,9.3848281e-01,-3.40393911,
		  5.50742903,-4.25196238,1.2687104,7.24136e-03,28.056,7.2945,
		  1.2168,4.7349,6.5942,-14.9419,10.7255,-3.3407,0.3823,0.0,
		  0.0,0.0,1.9985,13.2632,11.9120,-13.0364,5.6720,-1.2155,0.1034,
		  0.0,0.0,0.0,0.9382,32.8948,-17.9286,-0.8772,5.8667,-2.4439,
		  0.3147,0.0,0.0,0.0,2.0097,3.7515,-7.9734,10.9359,-7.1074,
		  2.1359,-0.2405,0.0,0.0,0.0,6.5316,26.8258,-18.2834,2.22,
		  4.1958,-2.0916,0.289,0.0,0.0,0.0,8.9833,126.0513,-21.9017,
		  20.7445,-20.3788,7.3643,-0.9082,0.0,0.0,0.0};

		for( k=1, _r = 0; k <= 11; k++ )
		{
			for( l=1; l <= 10; l++ )
			{
				b[l-1][k-1] = _itmp0[_r++];
				}
			}
		}
		_aini = 0;
	}

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

	/* These are fits given by 
	 *  >>refer	h1	cs	Callaway, J. 1994, At. Data Nucl. Data Tables, 57, 9
	 * H I collision strengths for 1s-2s, 1s-2p, 1s-3, 2s-3, and 2p-3, 
	 * principle quantum numbers n and l.
	 * CLOUDY sums angular momentum (s,p,d) components of level n=3.
	 *
	 * SEE IMPORTANT NOTE BELOW ABOUT EST. VALID ENERGY RANGES OF FITS.
	 ************************************************************************
	 *
	 * i is the lower level and runs from 1 to 3 (1s, 2s, 2p)
	 * j is the upper level and runs from 2 to 4 (2s, 2p, 3)
	 * */
	/* Callaway's fits use 10 variables and there are 11 transitions: */

	/* across are coef. to each transtion, down are the transitions */

	/* some constants
	 * */
	ev = 1.60184e-12;
	bk = 1.380622e-16;
	ryd = EVRYD /* 13.60583 */;

	/* convert temperature T (K) to energy in Rydbergs
	 * */
	T = bk*phycon.te/(ev*ryd);

	/**********************************************************
	 *************-----------IMPORTANT------------*************
	 **********************************************************
	 *
	 * The Callaway fits to effective cross-sections are high-order 
	 * polynomials which become insane not too far from the  
	 * unstated upper energy limits, given below. 
	 * Fits are more stable at low energy and here the validity 
	 * of cross-sections limits the validity of the fits.
	 *
	 * Energy limits of validity:
	 * minimum energy = 0.01 Rydbergs (about 1580K)
	 * */
	Tmin = 0.01;

	/* 1s - 2s,2p maximum energy = 4 Rydbergs (~ 632,000K)
	 * */
	Tmax1 = 4.0;

	/* 1s - 3 max energy is 1.2667 Rydberg (200,000K)
	 * */
	Tmax2 = 1.2667;

	/* 2s,p - 3 max energy is 1.0 Rydberg (158,888K)
	 * */
	Tmax3 = 1.0;

	/**********************************************************
	 * initialize cs */
	cs = 0.;

	/**************fits begin here**************
	 * */
	if( i == 1 && j == 2 )
	{

		/* 1s - 2s (first entry)
		 *
		 * check temperature:  */
		TeUse = MIN2(T,Tmax1);
		T = MAX2(TeUse,Tmin);

		cs = 0.;
		for( k=0; k < 5; k++ )
		{
			cs += b[k][0]*powi(T,k);
		}

		H1cs123_v = cs + b[5][0]*log(b[6][0]*T)*exp(-b[7][0]*T);
	}

	else if( i == 1 && j == 3 )
	{

		/* 1s - 2p (second entry)
		 *
		 * check temperature:  */
		TeUse = MIN2(T,Tmax1);
		T = MAX2(TeUse,Tmin);

		cs = 0.;
		for( k=0; k < 6; k++ )
		{
			cs += b[k][1]*powi(T,k);
		}
		H1cs123_v = cs;
	}

	else if( i == 1 && j == 4 )
	{

		/* 1s - 3s (third entry)
		 *
		 * check temperature:  */
		TeUse = MIN2(T,Tmax2);
		T = MAX2(TeUse,Tmin);

		cs3s = 0.;
		for( k=0; k < 7; k++ )
		{
			cs3s += b[k][2]*powi(T,k);
		}

		cs3s += b[7][2]*log(b[8][2]*T)*exp(-b[9][2]*T);

		/* 1s - 3p (fourth entry) */
		cs3p = 0.;
		for( k=0; k < 7; k++ )
		{
			cs3p += b[k][3]*powi(T,k);
		}

		/* 1s - 3d (fifth entry) */
		cs3d = 0.;
		for( k=0; k < 7; k++ )
		{
			cs3d += b[k][4]*powi(T,k);
		}

		cs3d += b[7][4]*log(b[8][4]*T)*exp(-b[9][4]*T);

		/* Summed cs: 1s - 3  */
		H1cs123_v = cs3s + cs3p + cs3d;
	}

	else if( i == 2 && j == 4 )
	{

		/* 2s - 3s (sixth entry)
		 *
		 * check temperature:  */
		TeUse = MIN2(T,Tmax3);
		T = MAX2(TeUse,Tmin);

		cs3s = 0.;
		for( k=0; k < 7; k++ )
		{
			cs3s += b[k][5]*powi(T,k);
		}

		/* 2s - 3p (seventh entry) */
		cs3p = 0.;
		for( k=0; k < 7; k++ )
		{
			cs3p += b[k][6]*powi(T,k);
		}

		/* 2s - 3d (eighth entry) */
		cs3d = 0.;
		for( k=0; k < 7; k++ )
		{
			cs3d += b[k][7]*powi(T,k);
		}

		/* Summed cs: 2s - 3  */
		H1cs123_v = cs3s + cs3p + cs3d;
	}

	else if( i == 3 && j == 4 )
	{
		/* 2p - 3s (ninth entry)
		 *
		 * check temperature:  */
		TeUse = MIN2(T,Tmax3);
		T = MAX2(TeUse,Tmin);

		cs3s = 0.;
		for( k=0; k < 7; k++ )
		{
			cs3s += b[k][8]*powi(T,k);
		}

		/* 2p - 3p (tenth entry) */
		cs3p = 0.;
		for( k=0; k < 7; k++ )
		{
			cs3p += b[k][9]*powi(T,k);
		}

		/* 2p - 3d (eleventh entry) */
		cs3d = 0.;
		for( k=0; k < 7; k++ )
		{
			cs3d += b[k][10]*powi(T,k);
		}

		/* Summed cs: 2p - 3 */
		H1cs123_v = cs3s + cs3p + cs3d;
	}

	else if( i == 2 && j == 3 )
	{
		/* 2p - 2s 
		 * */
		te20 = pow(phycon.te,.2f);
		te03 = pow(phycon.te,.03f);
		te10 = pow(phycon.te,.10f);
		te05 = pow(phycon.te,.05f);
		if( chType == 'e' )
		{
			rate = 5.738e-4/(te20*te20/te03);
		}
		else
		{
			rate = 6.290e-4/(te10*te05);
		}
		/*  convert to a cs */
		H1cs123_v = rate*sqrt(phycon.te)*6./COLL_CONST;

	}
	else
	{
		fprintf( ioQQQ, "  insane levels for H I n=1,2,3 !!!\n" );
		fprintf( ioQQQ, "%4ld%4ld\n", i, j );
		puts( "[Stop in h1cs123]" );
		H1cs123_v = 0.;
		cdEXIT(EXIT_FAILURE);
	}

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

/*He2cs123 line collision strengths for lower levels of helium ion, n=1,2,3, by K Korista */
static double He2cs123(long int i, 
  long int j)
{
	double He2cs123_v, 
	  cs3d, 
	  cs3p, 
	  cs3s, 
	  t;
	static double a[11]={0.12176209,0.32916723,0.46546497,0.044501688,
	  0.040523277,0.5234889,1.4903214,1.4215094,1.0295881,4.769306,
	  9.7226127};
	static double b[11]={0.039936166,2.9711166e-05,-0.020835863,3.0508137e-04,
	  -2.004485e-15,4.41475e-06,1.0622666e-05,2.0538877e-06,0.80638448,
	  2.0967075e-06,7.6089851e-05};
	static double c[11]={143284.77,0.73158545,-2.159172,0.43254802,
	  2.1338557,8.9899702e-06,-2.9001451e-12,1.762076e-05,52741.735,
	  -2153.1219,-3.3996921e-11};

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

	/* These are fits to Table 2
	 * >>refer	he2	cs	Aggarwal, K.M., Callaway, J., Kingston, A.E., Unnikrishnan, K.
	 * >>refercon	1992, ApJS, 80, 473
	 * He II collision rates for 1s-2s, 1s-2p, 1s-3, 2s-3, and 2p-3, 
	 * principal quantum numbers n and l.
	 * CLOUDY sums angular momentum (s,p,d) components of level n=3.
	 *
	 * i is the lower level and runs from 1 to 3 (1s, 2s, 2p)
	 * j is the upper level and runs from 2 to 4 (2s, 2p, 3)
	 * check temperature: fits only good between 5,000K and 500,000K
	 * */
	/* array for fits of 11 transitions see the code below for key: */

	t = phycon.te;
	if( t < 5000. )
	{
		t = 5000.;
	}
	else if( t > 5.0e05 )
	{
		t = 5.0e05;
	}

	/**************fits begin here**************
	 * */
	if( i == 1 && j == 2 )
	{
		/* 1s - 2s (first entry) */
		He2cs123_v = a[0] + b[0]*exp(-t/c[0]);
		/*  He2Rate123 = cs*COLL_CONST/SQRT(te)/2.
		 * */
	}
	else if( i == 1 && j == 3 )
	{
		/* 1s - 2p (second entry) */
		He2cs123_v = a[1] + b[1]*pow(t,c[1]);
		/*  He2Rate123 = cs*COLL_CONST/SQRT(te)/6.
		 * */
	}
	else if( i == 1 && j == 4 )
	{
		/* 1s - 3s (third entry) */
		cs3s = a[2] + b[2]*log(t) + c[2]/log(t);

		/* 1s - 3p (fourth entry) */
		cs3p = a[3] + b[3]*pow(t,c[3]);

		/* 1s - 3d (fifth entry) */
		cs3d = a[4] + b[4]*pow(t,c[4]);

		/* 1s - 3 sum */
		He2cs123_v = cs3s + cs3p + cs3d;
		/*  He2Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else if( i == 2 && j == 4 )
	{
		/* 2s - 3s (sixth entry)         */
		cs3s = (a[5] + c[5]*t)/(1 + b[5]*t);

		/* 2s - 3p (seventh entry) */
		cs3p = a[6] + b[6]*t + c[6]*t*t;

		/* 2s - 3d (eighth entry) */
		cs3d = (a[7] + c[7]*t)/(1 + b[7]*t);

		/* 2s - 3 sum */
		He2cs123_v = cs3s + cs3p + cs3d;
		/*  He2Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else if( i == 3 && j == 4 )
	{
		/* 2p - 3s (ninth entry) */
		cs3s = a[8] + b[8]*exp(-t/c[8]);

		/* 2p - 3p (tenth entry) */
		cs3p = a[9] + b[9]*t + c[9]/t;

		/* 2p - 3d (eleventh entry) */
		cs3d = a[10] + b[10]*t + c[10]*t*t;

		/* 2p - 3 sum */
		He2cs123_v = cs3s + cs3p + cs3d;
		/*  He2Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else
	{

		/**************fits end here************** */
		fprintf( ioQQQ, "  insane levels for He II n=1,2,3 !!!\n" );
		puts( "[Stop in he2cs123]" );
		cdEXIT(EXIT_FAILURE);
	}

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

/*Fe26cs123 line collision rates for lower levels of hydrogenic iron, n=1,2,3 */
static double Fe26cs123(long int i, 
  long int j)
{
	double Fe26cs123_v, 
	  TeUse, 
	  t, 
	  x;
	static double a[3]={-4.238398,-238.2599,-1211.5237};
	static double b[3]={-0.4448177,-27.06869,-136.7659};
	static double c[3]={0.0022861,0.153273,0.7677703};
	static double d[3]={3.303775,191.7165,972.3731};
	static double e[3]={15.82689,878.1333,4468.696};

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

	/*  These are fits to Table 5 of
	 *  >>refer	fe26	cs	Aggarwal, K.M., & Kingston, A.E. 1993, ApJS, 85, 187
	 *  Fe XXVI collision rates for 1s-3, 2s-3, and 2p-3, 
	 *  principal quantum numbers n and l.
	 *
	 *  i is the lower level and runs from 1 to 3 (1s, 2s, 2p)
	 *  j is the upper level and runs from 2 to 4 (2s, 2p, 3)
	 *  1s-2s,2p is not done here.
	 *  check temperature: fits only good between 5.2 < log Te < 7.2
	 *  */
	/*  arrays for fits of 3 transitions see the code below for key: */

	TeUse = MAX2(phycon.te,1.585e5);
	t = MIN2(TeUse,1.585e7);
	x = log10(t);

	if( i == 1 && j == 2 )
	{
		/* 1s - 2s (first entry) */
		fprintf( ioQQQ, " Fe XXVI 2s-1s not done in Fe26cs123\n" );
		puts( "[Stop in fe26cs123]" );
		cdEXIT(EXIT_FAILURE);
	}

	else if( i == 1 && j == 3 )
	{
		/* 1s - 2p (second entry) */
		fprintf( ioQQQ, " Fe XXVI 2p-1s not done in Fe26cs123\n" );
		puts( "[Stop in fe26cs123]" );
		cdEXIT(EXIT_FAILURE);
	}

	else if( i == 1 && j == 4 )
	{
		/* 1s - 3 (first entry) */
		Fe26cs123_v = a[0] + b[0]*x + c[0]*x*x*sqrt(x) + d[0]*log(x) + 
		  e[0]*log(x)/x/x;
		/*  Fe26Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else if( i == 2 && j == 4 )
	{
		/* 2s - 3 (second entry)         */
		Fe26cs123_v = a[1] + b[1]*x + c[1]*x*x*sqrt(x) + d[1]*log(x) + 
		  e[1]*log(x)/x/x;
		/*  Fe26Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else if( i == 3 && j == 4 )
	{
		/* 2p - 3s (third entry) */
		Fe26cs123_v = a[2] + b[2]*x + c[2]*x*x*sqrt(x) + d[2]*log(x) + 
		  e[2]*log(x)/x/x;
		/*  Fe26Rate123 = cs*COLL_CONST/SQRT(te)/18.
		 * */
	}
	else
	{

		fprintf( ioQQQ, "  insane levels for Ca XX n=1,2,3 !!!\n" );
		puts( "[Stop in fe26cs123]" );
		cdEXIT(EXIT_FAILURE);
	}

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


