/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*HeCollid evaluate collisional rates */
/*HeCSInterp interpolate on He1 collision strengths */
/*AtomCSInterp do the atom	*/
/*IonCSInterp do the ions	*/
/*CS_l_mixing_PS64 - find rate for l-mixing collisions by protons, for neutrals */
#include "cddefines.h" 
#include "iso.h"
#include "dense.h"
#include "helike.h"
#include "helike_cs.h"
#include "phycon.h"
#include "physconst.h"
#include "taulines.h"
#include "hydro_vs_rates.h"
#include "trace.h"
#include "ionbal.h"
#include "opacity.h"
#include "heavy.h"
#include "rfield.h"
#include "atmdat.h"
#include "math_complex.h"
#include "thirdparty.h"
#include "bevington.h"
#include "conv.h"
/*lint -e661 possible access of out-of-bounds pointer */


static double S62_Therm_ave_coll_str( double EProjectile_eV );
static double Therm_ave_coll_str_int( double EProjectileRyd);
static double collision_strength( double velOrEner, long int paramIsRedVel );
static double L_mix_integrand( double alpha );
static double StarkCollTransProb( long int n, long int l, long int lp, double alpha, double deltaPhi);
static double ellip_int_K( double k );

static long	int global_n, global_l, global_l_prime, global_s, global_z, global_Collider;
static double global_bmax, global_red_vel, global_an, global_collider_charge;
static double global_I_energy_eV, global_deltaE, global_temp, global_osc_str, global_stat_weight;
static double kTRyd;
/* These are masses relative to the proton mass of the electron, proton, and alpha particle. */
static double ColliderMass[3] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0};

/*===================================================================================*/
/* evaluate collisional rates */
void HeCollid( long int nelem)
{

	double factor1 , ConvLTEPOP, *ColIonizPerN, DimaRate, crate;
	long int ipLo , ipHi, n;
	/*long int ipFirstCollapsed;*/
	
	static double TeUsed[LIMELM]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
		0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0};

	static double TeUsedForCS[LIMELM]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
		0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0};

#	ifdef DEBUG_FUN
	fputs( "<+>HeCollid()\n", debug_fp );
#	endif
	
	if( (ColIonizPerN = (double*)MALLOC(sizeof(double)*
		(unsigned)(iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] + 1) ) )==NULL )
		BadMalloc();

	for( n=0; n <= ( iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] ); ++n)
	{
		ColIonizPerN[n] = 0.;
	}

#if	0
	/* >>chng 05 dec 30, reset TeUsedForCS if conv.nTotalIoniz is equal to zero. 
	 * this only affects things if we are doing a multiple model run */
	/* >>chng 06 jan 24, comment this out as per Ryan's instructions */
	if( !conv.nTotalIoniz )
		TeUsedForCS[nelem] = 0.;

	/* First calculate collision strengths or update if temp has changed much. */
	/*ipFirstCollapsed = iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem];*/
	if( (TeUsedForCS[nelem] == 0.) || 
		( TeUsedForCS[nelem]/phycon.te > 2.0 ) ||
		( TeUsedForCS[nelem]/phycon.te < 0.5 )  )
	{
		double cs_electron, cs_proton, cs_heplus;
		double cs_elec_old, cs_prot_old, cs_heplus_old;
		double old_Te;

		for( ipHi=ipHe2s3S; ipHi <iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
		{
			for( ipLo=ipHe1s1S; ipLo < ipHi; ipLo++ )
			{
				/* these are the current collision strengths, about to be the old ones */
  				cs_elec_old = EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs;
				cs_prot_old = helike.cs_proton[nelem][ipHi][ipLo];
				cs_heplus_old = helike.cs_heplus[nelem][ipHi][ipLo];

				/* get collision strength for electron, proton, and he+ impact */
  				cs_electron = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
				cs_proton = HeCSInterp( nelem , ipHi , ipLo, ipPROTON ); 
				cs_heplus = HeCSInterp( nelem , ipHi , ipLo, ipHE_PLUS );

				{
					/* option to print particulars of some line when called
					 * a prettier print statement is near where chSpin is defined below*/
					/*@-redef@*/
					enum {DEBUG_LOC=FALSE};
					/*@+redef@*/
					if( DEBUG_LOC && (nelem==ipHELIUM) && 
						(N_(ipHi)==N_(ipLo)) /* && (L_(ipHi)>=3) && (L_(ipLo)>=3) */
						&& (S_(ipHi)==S_(ipLo)) && (ipHi<43) )
					{
						fprintf(ioQQQ,"Lo %li n %li s %li l %li \t ", 
							ipLo, N_(ipLo), S_(ipLo), L_(ipLo) );
						fprintf(ioQQQ," Hi %li n %li s %li l %li \t", 
							ipHi, N_(ipHi), S_(ipHi), L_(ipHi) );
						fprintf(ioQQQ,"%.4e\t%.0f\t%.0f\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN,
							iso.stat[ipHE_LIKE][nelem][ipHi],
							iso.stat[ipHE_LIKE][nelem][ipLo],
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul,
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gf ,
							cs_electron,
							cs_proton,
							cs_heplus);
					}
				}

				if( TeUsedForCS[nelem] == 0. )
				{
					double saveTemp = phycon.te;
					
					old_Te = 2.0* phycon.te;
					/*TODO	2	rewrite HeCSInterp and all called routines to use a temp
					 * parameter, rather than phycon.te */
					phycon.te = old_Te;

  					cs_elec_old = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
					cs_prot_old = HeCSInterp( nelem , ipHi , ipLo, ipPROTON ); 
					cs_heplus_old = HeCSInterp( nelem , ipHi , ipLo, ipHE_PLUS );

					phycon.te = saveTemp;
				}
				else
				{
					old_Te = TeUsedForCS[nelem];
				}

				/* Find power laws for each collision strength, but if 
				 * the old or new collision strength is very small, assume no
				 * temperature dependence in the neighborhood.  */
				if( cs_electron <= 1.01E-10 || cs_elec_old <= 1.01E-10 )
				{
					helike.cs_elec_power[nelem][ipHi][ipLo] = 0.;
				}
				else
				{
					helike.cs_elec_power[nelem][ipHi][ipLo] = 
						log10( (double)cs_electron/(double)cs_elec_old )/
						log10( phycon.te/old_Te );
				}

				if( cs_proton <= 1.01E-10 || cs_prot_old <= 1.01E-10 )
				{
					helike.cs_prot_power[nelem][ipHi][ipLo] = 0.;
				}
				else
				{
					helike.cs_prot_power[nelem][ipHi][ipLo] = 
						log10( (double)cs_proton/(double)cs_prot_old )/
						log10( phycon.te/old_Te );
				}
				
				if( cs_heplus <= 1.01E-10 || cs_heplus_old <= 1.01E-10 )
				{
					helike.cs_heplus_power[nelem][ipHi][ipLo] = 0.;
				}
				else
				{
					helike.cs_heplus_power[nelem][ipHi][ipLo] = 
						log10( (double)cs_heplus/(double)cs_heplus_old )/
						log10( phycon.te/old_Te );
				}

				/* Put a very simple limit on these powers.  The range -15 to
				 * 15 should be very forgiving...We should expect any number outside
				 * that range to be bogus. */
				ASSERT( helike.cs_elec_power[nelem][ipHi][ipLo] < 15. );
				ASSERT( helike.cs_prot_power[nelem][ipHi][ipLo] < 15. );
				ASSERT( helike.cs_heplus_power[nelem][ipHi][ipLo] < 15. );
				ASSERT( helike.cs_elec_power[nelem][ipHi][ipLo] > -15. );
				ASSERT( helike.cs_prot_power[nelem][ipHi][ipLo] > -15. );
				ASSERT( helike.cs_heplus_power[nelem][ipHi][ipLo] > -15. );

				EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs = cs_electron;
				helike.cs_proton[nelem][ipHi][ipLo] = cs_proton;
				helike.cs_heplus[nelem][ipHi][ipLo] = cs_heplus;

			}
		}

		TeUsedForCS[nelem] = phycon.te;
	}
#endif

	/*lint -e777 */
	if( TeUsed[nelem] != phycon.te )
	{
	/*lint +e777 */
		TeUsed[nelem] = phycon.te;

		if( trace.lgTrace )
		{
			fprintf( ioQQQ, 
				"       HeCollid called nelem %li - will reeval Boltz fac, LTE dens\n",
				nelem );
		}

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

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

		iso.lgPopLTE_OK[ipHE_LIKE][nelem] = TRUE;
		/* fully define Boltzmann factors to continuum for model levels */
		for( ipLo=ipHe1s1S; ipLo < (iso.numLevels_max[ipHE_LIKE][nelem]); ipLo++ )
		{
			/* this Boltzmann factor is exp( +ioniz energy / Te ) */
			iso.ConBoltz[ipHE_LIKE][nelem][ipLo] = 
				dsexp(iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo]/phycon.te_ryd);

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

			if( iso.ConBoltz[ipHE_LIKE][nelem][ipLo] >= SMALLDOUBLE )
			{
				iso.PopLTE[ipHE_LIKE][nelem][ipLo] = 
					iso.stat[ipHE_LIKE][nelem][ipLo] / iso.ConBoltz[ipHE_LIKE][nelem][ipLo]* ConvLTEPOP;
				ASSERT( iso.PopLTE[ipHE_LIKE][nelem][ipLo] < BIGDOUBLE );
			}
			else
			{
				iso.PopLTE[ipHE_LIKE][nelem][ipLo] = 0.;
			}

			/* now check for any zeros - possible above underflowed to zero */
			if( iso.PopLTE[ipHE_LIKE][nelem][ipLo] <= 0. )
			{
				iso.lgPopLTE_OK[ipHE_LIKE][nelem] = FALSE;
			}
		}

		if( (TeUsedForCS[nelem] == 0.) || 
			( TeUsedForCS[nelem]/phycon.te > 1.15 ) ||
			( TeUsedForCS[nelem]/phycon.te < 0.85 ) )
		{
			for( ipHi=ipHe2s3S; ipHi <iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
			{
				for( ipLo=ipHe1s1S; ipLo < ipHi; ipLo++ )
				{

					EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs = 
						HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
					helike.cs_proton[nelem][ipHi][ipLo] = 
						HeCSInterp( nelem , ipHi , ipLo, ipPROTON );
					helike.cs_heplus[nelem][ipHi][ipLo] = 
						HeCSInterp( nelem , ipHi , ipLo, ipHE_PLUS );

					ASSERT( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs >= 0. );
					ASSERT( helike.cs_proton[nelem][ipHi][ipLo] >= 0. );
					ASSERT( helike.cs_heplus[nelem][ipHi][ipLo] >= 0. );
				}
			}
			TeUsedForCS[nelem] = phycon.te;
		}

		for( ipHi=ipHe2s3S; ipHi <iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
		{
			for( ipLo=ipHe1s1S; ipLo < ipHi; ipLo++ )
			{
				/********************************************************
				 ********************************************************
				 * NB - the collision strengths for proton and helium  *
				 * ion impact are multiplied by the ratio of collider  *
				 * to electron densities to precorrect for getting     *
				 * multiplied by the electron density when being put   *
				 * into rate matrix later                              *
				 ********************************************************
				 ********************************************************/

				EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ColUL = (float)(
					(
					/* due to electron impact */
					EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs+
					/* due to proton impact */
					helike.cs_proton[nelem][ipHi][ipLo]*
					(float)(dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr)+
					/* due to he+ impact */
					helike.cs_heplus[nelem][ipHi][ipLo]*
					(float)(dense.xIonDense[ipHELIUM][1]/dense.EdenHCorr)
					) / phycon.sqrte*COLL_CONST/(double)iso.stat[ipHE_LIKE][nelem][ipHi] );

#if	0
				/* collisional deexcitation rate, can be turned off with 
				 * atom he-like collisions off command */
				/*This is downward collision rate */
				EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ColUL = (float)(
					(
					/* due to electron impact */
					EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs*
					pow( phycon.te/TeUsedForCS[nelem], helike.cs_elec_power[nelem][ipHi][ipLo] ) +
					/* due to proton impact */
					helike.cs_proton[nelem][ipHi][ipLo]*
					(float)(dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr)*
					pow( phycon.te/TeUsedForCS[nelem], helike.cs_prot_power[nelem][ipHi][ipLo] ) +
					/* due to he+ impact */
					helike.cs_heplus[nelem][ipHi][ipLo]*
					(float)(dense.xIonDense[ipHELIUM][1]/dense.EdenHCorr)*
					pow( phycon.te/TeUsedForCS[nelem], helike.cs_heplus_power[nelem][ipHi][ipLo] )
					) / phycon.sqrte*COLL_CONST/(double)iso.stat[ipHE_LIKE][nelem][ipHi] );
#endif

				if( N_(ipHi) > iso.n_HighestResolved_max[ipHE_LIKE][nelem] &&
					N_(ipLo) <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
				{
					EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ColUL *=
						(float)( (2./3.)*(log((double)N_(ipHi))+2) );
				}

				/* define line Boltzmann factors 
				 * some energy-degenerate collisions have very small negative energy,
				 * when different l, same n, levels flip over */
				if( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyK > 0. )
				{
					iso.Boltzmann[ipHE_LIKE][nelem][ipHi][ipLo] = 
						sexp( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyK / phycon.te );
				}
				else
				{
					/* negative energy, assume h nu/kT is very small */
					iso.Boltzmann[ipHE_LIKE][nelem][ipHi][ipLo] = 1.;
				}
			}

			/* collisional ionization
			 * this is rate coefficient, not rate, for hydrogenic species */
			/*TODO	2	find ioniz rates for helium like species, these are hydrogenic rates */
			ASSERT( ipHi > 0 && ipHi < iso.numLevels_max[ipHE_LIKE][nelem]);
			if( nelem==ipHELIUM && helike.lgSetBenjamin )
			{
				/* atoms */
				
					double BenCollIonParamA[31] = {
						3.36E-09,9.36E-08,1.33E-07,1.58E-07,1.58E-07,1.58E-07,1.81E-07,5.36E-07,
						6.58E-07,7.23E-07,7.81E-07,7.81E-07,7.93E-07,1.62E-06,1.87E-06,2.00E-06,
						2.11E-06,2.11E-06,2.11E-06,2.11E-06,2.13E-06,3.61E-06,4.02E-06,4.21E-06,
						4.39E-06,4.40E-06,4.40E-06,4.40E-06,4.40E-06,4.40E-06,4.44E-06};

					double BenCollIonParamB[31] = {
						0.499,0.280,0.253,0.240,0.240,0.240,0.229,0.138,0.120,0.111,0.104,0.104,
						0.103,0.035,0.021,0.014,0.008,0.008,0.008,0.008,0.007,-0.047,-0.059,-0.064,
						-0.068,-0.068,-0.068,-0.068,-0.068,-0.068,-0.069};

					double BenCollIonParamC[31] = {
						28.625,5.655,4.733,4.331,4.331,4.331,4.037,2.304,2.072,1.972,1.895,
						1.895,1.880,1.298,1.207,1.167,1.135,1.135,1.134,1.134,1.128,0.866,
						0.821,0.802,0.785,0.785,0.784,0.784,0.784,0.784,0.781};

					ASSERT( ipHi < 31 );

					iso.ColIoniz[ipHE_LIKE][nelem][ipHi] = BenCollIonParamA[ipHi] * 
						pow( (double)phycon.te/10000., BenCollIonParamB[ipHi] ) *
						sexp( 10000.*BenCollIonParamC[ipHi]/phycon.te );
			}
			else 
			{
				if( nelem == ipHELIUM )
				{
					iso.ColIoniz[ipHE_LIKE][nelem][ipHi] = hydro_vs_ioniz( ipHE_LIKE, nelem, ipHi );
				}
				else
				{
					/* ions */
					iso.ColIoniz[ipHE_LIKE][nelem][ipHi] = Hion_coll_ioniz_ratecoef(ipHE_LIKE , nelem , ipHi);
				}

				if( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n > iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
				{
					iso.ColIoniz[ipHE_LIKE][nelem][ipHi] *= 2.*N_(ipHi);
				}

				/* zero collisional ionization if turned off.	*/
				iso.ColIoniz[ipHE_LIKE][nelem][ipHi] *= iso.lgColl_ionize[ipHE_LIKE];
			}

			ColIonizPerN[iso.quant_desig[ipHE_LIKE][nelem][ipHi].n] += iso.ColIoniz[ipHE_LIKE][nelem][ipHi];
		}

		helike.lgErrGenDone = TRUE;

		/* Here we arbitrary scale the highest level ionization to account for the fact
		 * that, if the atom is not full size, this level should be interacting with higher
		 * levels and not just the continuum.  We did add on collisional excitation terms instead
		 * but this caused a problem at low temperatures because the collisional ionization was 
		 * a sum of terms with different Boltzmann factors, while PopLTE had just one Boltzmann
		 * factor.  The result was a collisional recombination that had residual exponentials of
		 * the form exp(x/kT), which blew up at small T.	*/
		if( !iso.lgFullSize[ipHE_LIKE][nelem] && (nelem != ipHELIUM || !helike.lgSetBenjamin) )
			iso.ColIoniz[ipHE_LIKE][nelem][iso.numLevels_max[ipHE_LIKE][nelem]-1] *= 1.E5;

		/* Assert that total coll. ioniz. per n increases monotonically with increasing n.	*/
		if( iso.lgColl_ionize[ipHE_LIKE] && (nelem!=ipHELIUM || !helike.lgSetBenjamin ) )
		{
			for( n = 2; n < iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem]; n++  )
			{
				ASSERT( (ColIonizPerN[n] < ColIonizPerN[n+1]) || (ColIonizPerN[n]==0.) );
			}
		}

		/* this is to kill collisional excitations from n=1 and n=2 in the
		 * case b hummer and storey option */
		if( opac.lgCaseB_HummerStorey )
		{

			/* this kills collisional excitation and ionization from the ground
			 * 1s 1S level, and the two singlet n=2 levels, 2s 1S and 2p 1P */
			ipLo = ipHe1s1S;
			iso.ColIoniz[ipHE_LIKE][nelem][ipLo] = 0.;
			for( ipHi=ipLo+1; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
			{
				EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs = 0.;
				EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ColUL = 0.;
			}
			ipLo = ipHe2s1S;
			iso.ColIoniz[ipHE_LIKE][nelem][ipLo] = 0.;
			for( ipHi=ipLo+1; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
			{
				EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs = 0.;
				EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ColUL = 0.;
			}
			ipLo = ipHe2p1P;
			iso.ColIoniz[ipHE_LIKE][nelem][ipLo] = 0.;
			for( ipHi=ipLo+1; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
			{
				EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs = 0.;
				EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ColUL = 0.;
			}
		}
	}
	else
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, 
				"       HeCollid called nelem %li - no reeval Boltz fac, LTE dens\n",
				nelem );
		}
	}

	/* atmdat_coll_ion( Atom num (H=1), Num elec, >=0 */
	DimaRate = atmdat_coll_ion( nelem+1, 2 , phycon.te );

	crate = DimaRate*dense.EdenHCorr*iso.lgColl_ionize[ipHE_LIKE];
	/* NB ColIoniz does not include secondaries since not total, but CollidRate does,
	 * since is total */
	iso.ColIoniz[ipHE_LIKE][nelem][ipHe1s1S] = DimaRate*iso.lgColl_ionize[ipHE_LIKE];
	
	/* collisional ionization rate from ground with only thermal */
	/* this part of routine is always reevaluated, even if c is constant,
	 * this is important since density is in crate */
	ionbal.CollIonRate_Ground[nelem][nelem-1][0] = crate;

	/* cooling due to collisional ionization, which only includes thermal */
	ionbal.CollIonRate_Ground[nelem][nelem-1][1] = (crate*
		rfield.anu[Heavy.ipHeavy[nelem][nelem-1]-1]* EN1RYD);

	/* collisional ionizaiton of atomic helium */
	if( nelem == ipHELIUM && !helike.lgSetBenjamin )
	{
		double c2sion , c2pion;
		/* >>refer	he1	cs	Seaton, M.S. 1964, Plan Sp Sci 12, 55. */
		c2sion = 8.8e-10*(double)phycon.sqrte*iso.ConBoltz[ipHE_LIKE][nelem][ipHe2s3S];
		c2pion = 3.*8.8e-10*(double)phycon.sqrte*iso.ConBoltz[ipHE_LIKE][nelem][ipHe2p3P0];
		iso.ColIoniz[ipHE_LIKE][nelem][ipHe2s3S] = c2sion*iso.lgColl_ionize[ipHE_LIKE];
		/* Zhang and Sampson collision excitation data has that excitations out of the separate
		 * terms in 2^3P are equal.  Since coll. ionization is really just collisional excitation,
		 * the same thing should happen here.  So, this treatment should be correct.	*/
		iso.ColIoniz[ipHE_LIKE][nelem][ipHe2p3P0] = c2pion*iso.lgColl_ionize[ipHE_LIKE]/3.;
		iso.ColIoniz[ipHE_LIKE][nelem][ipHe2p3P1] = c2pion*iso.lgColl_ionize[ipHE_LIKE]/3.;
		iso.ColIoniz[ipHE_LIKE][nelem][ipHe2p3P2] = c2pion*iso.lgColl_ionize[ipHE_LIKE]/3.;
	}
	{
		/* prints table of wavelengths of lines to ground */
		/*@-redef@*/
		enum {PRINTIT=FALSE};
		/*@+redef@*/
		if( PRINTIT )
		{

			nelem = ipHELIUM;
			for( ipHi=1; ipHi<iso.numLevels_max[ipHE_LIKE][ipHELIUM]; ++ipHi )
			{
				for( ipLo=0; ipLo<ipHi; ++ipLo )
				{
					fprintf(ioQQQ,"Z %li Lo %li n %li s %li l %li \t ", 
						nelem+1, ipLo,
						N_(ipLo),
						S_(ipLo),
						L_(ipLo) );
					fprintf(ioQQQ," Hi %li n %li s %li l %li \t", 
						ipHi, 
						N_(ipHi),
						S_(ipHi),
						L_(ipHi) );
					fprintf(ioQQQ,"%.4e\t%.4e\tcs\t%.4e\n",
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN,
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul ,
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs);
				}
			}
			exit(1);
		}
	}


	free( ColIonizPerN );

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


/* Choose either AtomCSInterp or IonCSInterp */
float HeCSInterp(long int nelem,
				 long int ipHi,
				 long int ipLo,
				 long int Collider )
{
	float cs, factor1;

	/* This variable is for diagnostic purposes:
	 * a string used in the output to designate where each cs comes from.	*/	
	const char *where = "      ";

	if( iso.lgColl_excite[ipHE_LIKE] == FALSE )
		return (float)1E-10;
	
	if( nelem == ipHELIUM )
		cs = AtomCSInterp( nelem, ipHi , ipLo, &factor1, &where, Collider );
	else
		cs = IonCSInterp( nelem, ipHi , ipLo, &factor1, &where, Collider );

	putError(nelem, ipHi, ipLo, IPCOLLIS, 0.15f );

#if	0
	if( N_(ipLo)<=2 && N_(ipHi)<=5 )
		putError(nelem, ipHi, ipLo, IPCOLLIS, (float)BRAY_STD_DEV );
	else
	{
		putError(nelem, ipHi, ipLo, IPCOLLIS, 0.f );
	}
	else if( N_(ipHi)==N_(ipLo) )
		putError(nelem, ipHi, ipLo, IPCOLLIS, (float)LMIX_STD_DEV );
	else 
		putError(nelem, ipHi, ipLo, IPCOLLIS, (float)VRIENS_STD_DEV );
#endif

	ASSERT( cs >= 0.f );
	
	/* in many cases the correction factor for split states has already been made,
	 * if not then factor is still negative */
	/* Remove the second test here when IonCSInterp is up to par with AtomCSInterp */
	ASSERT( factor1 >= 0.f || nelem!=ipHELIUM );
	if( factor1 < 0.f )
	{
		ASSERT( helike.lgCS_Vriens );

		factor1 = 1.f;
	}

	/* take factor into account, usually 1, ratio of stat weights if within 2 3P 
	 * and with collisions from collapsed to resolved levels */
	cs *= factor1;

	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		
		if( DEBUG_LOC && ( nelem==ipOXYGEN ) && (cs > 0.f) && (iso.quant_desig[ipHE_LIKE][nelem][ipHi].n == 2) 
			&& ( iso.quant_desig[ipHE_LIKE][nelem][ipLo].n <= 2 ) )
			fprintf(ioQQQ,"%li\t%li\t%li\t-\t%li\t%li\t%li\t%.2e\t%s\n", 
				iso.quant_desig[ipHE_LIKE][nelem][ipLo].n , iso.quant_desig[ipHE_LIKE][nelem][ipLo].s ,
				iso.quant_desig[ipHE_LIKE][nelem][ipLo].l ,	iso.quant_desig[ipHE_LIKE][nelem][ipHi].n ,
				iso.quant_desig[ipHE_LIKE][nelem][ipHi].s , iso.quant_desig[ipHE_LIKE][nelem][ipHi].l , cs,where);
	}

	return MAX(cs, (float)1E-10);
}

float AtomCSInterp(long int nelem,
				   long int ipHi,
				   long int ipLo,
				   float *factor1,
				   const char **where,
				   long int Collider )
{
	long ipArray;
	float cs;

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

	ASSERT( nelem == ipHELIUM );

	/* init values, better be positive when we exit */
	cs = -1.f; 

	/* this may be used for splitting up the collision strength within 2 3P when
	 * the lower level is withint 2 3P, and for collisions between resolved and collapsed levels.
	 * It may be set somewhere in routine, so set to negative value here as flag saying not set */
	*factor1 = -1.f;

	/* for most of the helium iso sequence, the order of the J levels within 2 3P 
	 * in increasing energy, is 0 1 2 - the exception is atomic helium itself,
	 * which is swapped, 2 1 0 */

	/* this branch is for upper and lower levels within 2p3P */
	if( ipLo >= ipHe2p3P0 && ipHi <= ipHe2p3P2 && Collider==ipELECTRON )
	{
		*factor1 = 1.f;
		/* atomic helium, use Berrington private comm cs */
		
		/* >>refer	he1	cs	Berrington, Keith, 2001, private communication - email follows
		> Dear Gary,
		> I could not find any literature on the He fine-structure
		> problem (but I didn't look very hard, so there may be 
		> something somewhere). However, I did a quick R-matrix run to 
		> see what the magnitude of the collision strengths are... At 
		> 1000K, I get the effective collision strength for 2^3P J=0-1, 
		>  0-2, 1-2 as 0.8,0.7,2.7; for 10000K I get 1.2, 2.1, 6.0
		*/
		/* indices are the same and correct, only thing special is that energies are in inverted order...was right first time.	*/
		if( ipLo == ipHe2p3P0 && ipHi == ipHe2p3P1 )
		{
			cs = 1.2f;
		}
		else if( ipLo == ipHe2p3P0 && ipHi == ipHe2p3P2 )
		{
			cs = 2.1f;
		}
		else if( ipLo == ipHe2p3P1 && ipHi == ipHe2p3P2 )
		{
			cs = 6.0f;
		}
		else
		{
			cs = 1.0f;
			TotalInsanity();
		}
		
		*where = "Berr  ";
		/* statistical weights included */
	}
	/* >>chng 02 feb 25, Bray data should come first since it is the best we have.	*/
	/* this branch is the Bray et al data, for n <= 5, where quantal calcs exist 
	 * must exclude ipLo >= ipHe2p1P because they give no numbers for those	*/
	else if( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n <= 5 && 
		( ipHi < iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ) &&
		helike.HeCS[nelem][ipHi][ipLo][0] >= 0.f && Collider== ipELECTRON )
	{
		ASSERT( *factor1 == -1.f );
		ASSERT( ipLo < ipHi );
		ASSERT( ipHe2p3P0 == 3 );

		/* ipLo is within 2^3P	*/
		if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
		{
			/* *factor1 is ratio of statistical weights of level to term */
			
			/* ipHe2p3P0, ipHe2p3P1, ipHe2p3P2 have indices 3,4,and 5, but j=0,1,and 2.	*/
			*factor1 = (2.f*((float)ipLo-3.f)+1.f) / 9.f;
			
			/* ipHi must be above ipHe2p3P2 since 2p3Pj->2p3Pk taken care of above	*/
			ASSERT( ipHi > ipHe2p3P2 );
		}
		/* ipHi is within 2^3P	*/
		else if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
		{
			ASSERT( ipLo < ipHe2p3P0 );

			*factor1 = (2.f*((float)ipHi-3.f)+1.f) / 9.f;
		}
		/* neither are within 2^3P...no splitting necessary	*/
		else 
		{
			*factor1 = 1.f;
		}

		/* SOME OF THESE ARE NOT N-CHANGING!	*/
		/* Must be careful about turning each one on or off.	*/
		
		/* this is the case where we have quantal calculations */
		/* >>refer	He1	cs	Bray, I., Burgess, A., Fursa, D.V., & Tully, J.A., 2000, A&AS, 146, 481-498 */
		/* check whether we are outside temperature array bounds,
		 * and use extreme value if we are */
		if( phycon.alogte <= helike.CSTemp[0] )
		{
			cs = helike.HeCS[nelem][ipHi][ipLo][0];
		}
		else if( phycon.alogte >= helike.CSTemp[helike.nCS-1] )
		{
			cs = helike.HeCS[nelem][ipHi][ipLo][helike.nCS-1];
		}
		else
		{
			float flow; 

			/* find which array element within the cs vs temp array */
			ipArray = (long)((phycon.alogte - helike.CSTemp[0])/(helike.CSTemp[1]-helike.CSTemp[0]));
			ASSERT( ipArray < helike.nCS );
			ASSERT( ipArray >= 0 );
			/* when taking the average, this is the fraction from the lower temperature value */
			flow = (float)( (phycon.alogte - helike.CSTemp[ipArray])/
				(helike.CSTemp[ipArray+1]-helike.CSTemp[ipArray])) ;
			ASSERT( (flow >= 0.f) && (flow <= 1.f) );

			cs = helike.HeCS[nelem][ipHi][ipLo][ipArray] * (1.f-flow) +
				helike.HeCS[nelem][ipHi][ipLo][ipArray+1] * flow;
		}

		*where = "Bray ";

		/* options to kill collisional excitation and/or l-mixing	*/
		if ( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n == iso.quant_desig[ipHE_LIKE][nelem][ipLo].n )
			/* iso.lgColl_l_mixing turned off with atom he-like l-mixing collisions off command */
			cs *= (float)iso.lgColl_l_mixing[ipHE_LIKE];
		else
		{
			/* iso.lgColl_excite turned off with atom he-like collisional excitation off command */
			cs *= (float)iso.lgColl_excite[ipHE_LIKE];
			/* This kills collisions to n=5, for comparison to Benjamin et al 99.	*/
			if( ( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n >= 5 ) && helike.lgSetBenjamin )
				cs = 0.f;
		}

		ASSERT( cs >= 0.f );
		/* statistical weights included */
	}
	/* this branch, n-same, l-changing collision, and not case of He with quantal data */
	else if( (iso.quant_desig[ipHE_LIKE][nelem][ipHi].n == iso.quant_desig[ipHE_LIKE][nelem][ipLo].n ) &&
		(iso.quant_desig[ipHE_LIKE][nelem][ipHi].s == iso.quant_desig[ipHE_LIKE][nelem][ipLo].s ) )
	{
		ASSERT( *factor1 == -1.f );
		*factor1 = 1.f;
		
		/* ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n >= 3 ); */
		ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] );

		if( (iso.quant_desig[ipHE_LIKE][nelem][ipLo].l <=2) &&
			abs(iso.quant_desig[ipHE_LIKE][nelem][ipHi].l - iso.quant_desig[ipHE_LIKE][nelem][ipLo].l)== 1 )
		{
			/* Use the method given in 
			 * >>refer He	CS	Seaton, M. J. 1962, Proc. Phys. Soc. 79, 1105 
			 * statistical weights included */
			cs = (float)CS_l_mixing_S62(ipHE_LIKE, nelem, ipLo, ipHi, (double)phycon.te, Collider); 

			/*fprintf( ioQQQ,"ipHi\t%li\tipLo\t%li\tcs\t%e\n", ipHi, ipLo, cs );*/
		}
		else if( helike.lgCS_Vrinceanu == TRUE )
		{
			if( iso.quant_desig[ipHE_LIKE][nelem][ipLo].l >=3 &&
				iso.quant_desig[ipHE_LIKE][nelem][ipHi].l >=3 )
			{
				/* Use the method given in 
				 * >>refer He	CS	Vrinceanu, D. \& Flannery, M. R. 2001, PhysRevA 63, 032701 
				 * statistical weights included */
				cs = (float)CS_l_mixing_VF01( nelem,
					iso.quant_desig[ipHE_LIKE][nelem][ipLo].n,
					iso.quant_desig[ipHE_LIKE][nelem][ipLo].l,
					iso.quant_desig[ipHE_LIKE][nelem][ipHi].l,
					iso.quant_desig[ipHE_LIKE][nelem][ipLo].s,
					(double)phycon.te,
					Collider );
			}
			else
			{
				cs = 0.f;
			}
		}
		/* this branch, l changing by one */
		else if( abs(iso.quant_desig[ipHE_LIKE][nelem][ipHi].l - iso.quant_desig[ipHE_LIKE][nelem][ipLo].l)== 1)
		{
			/* >>refer	He	cs	Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */
			/* statistical weights included */
			cs = (float)CS_l_mixing_PS64( nelem, ipLo, ipHi, Collider);
		}
		else
		{
			/* l changes by more than 1, but same-n collision */
			cs = 0.f;
		}

		/* ipLo is within 2^3P	*/
		if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
		{
			*factor1 = (2.f*((float)ipLo-3.f)+1.f) / 9.f;
		}

		/* ipHi is within 2^3P	*/
		if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
		{
			*factor1 = (2.f*((float)ipHi-3.f)+1.f) / 9.f;
		}

		*where = "lmix  ";
		cs *= (float)iso.lgColl_l_mixing[ipHE_LIKE];
	}

	/* this is an atomic n-changing collision with no quantal calculations */
	else if ( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n != iso.quant_desig[ipHE_LIKE][nelem][ipLo].n )
	{
		ASSERT( *factor1 == -1.f );
		/* this is an atomic n-changing collision with no quantal calculations */
		/* gbar g-bar goes here */

		/* >>chng 02 jul 25, add option for fits to quantal cs data */
		if( helike.lgCS_Vriens )
		{
			/* this is Vriens & Smeets collision strength for neutral H */
			/* There is no possibility for modifying this routine for an arbitrary
			 * collider, so the below routine thermally averages the original Vriens
			 * & Smeets cross-section, modified for other colliders, and converted to
			 * a collision strength. */
			/*
			cs = (float)hydro_vs_deexcit(
							ipHE_LIKE, 
							nelem,
							ipHi,
							ipLo ) * iso.stat[ipHE_LIKE][nelem][ipLo];
			*/
			
			/* >>refer He CS	Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940
			 * statistical weight is NOT included in the routine */
			cs = (float)CS_VS80(
							ipHE_LIKE,
							nelem,
							ipHi,
							ipLo,
							phycon.te,
							Collider ) * iso.stat[ipHE_LIKE][nelem][ipLo];

			*factor1 = 1.f;
			*where = "Vriens";
		}
		else if( helike.lgCS_None )
		{
			cs = 0.f;
			*factor1 = 1.f;
			*where = "no gb";
		}
		else if( helike.lgCS_new )
		{
			*factor1 = 1.f;
			/* Don't know if stat weights are included in this, but they're probably
			 * wrong anyway since they are based in part on the former (incorrect)
			 * implementation of Vriens and Smeets rates */
			
			/* two different fits, allowed and forbidden */
			if( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul > 1. )
			{
				/* permitted lines - large A */
				double x = 
					log10(MAX2(34.7,EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN));

				if( helike.lgCS_new == 1 )
				{
					/* this is the broken power law fit, passing through both quantal
					 * calcs at high energy and asymptotically goes to VS at low energies */
					if( x < 4.5 )
					{
						/* low energy fit for permitted transitions */
						cs = (float)pow( 10. , -1.45*x + 6.75);
					}
					else
					{
						/* higher energy fit for permitted transitions */
						cs = (float)pow( 10. , -3.33*x+15.15);
					}
				}
				else if( helike.lgCS_new == 2 )
				{
					/* single parallel fit for permitted transitions, runs parallel to VS */
					cs = (float)pow( 10. , -2.3*x+10.3);
				}
				else
					TotalInsanity();
			}
			else
			{
				/* fit for forbidden transitions */
				if( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN < 25119.f )
				{
					cs = 0.631f; 
				}
				else
				{
					cs = (float)pow(10., 
						-3.*log10(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN)+12.8);
				}
			}

			*where = "newgb";

			/* option to print the resulting collision strengths
			if( nelem==1 )
				fprintf(ioQQQ,"%.3e\t%.3e\t%.3e\n", 
				EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul ,
				EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN,
				cs ); */
		}
		else
			TotalInsanity();

		/* ipLo is within 2^3P	*/
		if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
		{
			*factor1 = (2.f*((float)ipLo-3.f)+1.f) / 9.f;
		}

		/* ipHi is within 2^3P	*/
		if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
		{
			*factor1 = (2.f*((float)ipHi-3.f)+1.f) / 9.f;
		}

		/* options to turn off collisions */
		cs *= (float)iso.lgColl_excite[ipHE_LIKE];

		/* This kills collisions to n=5, for comparison to Benjamin et al 99.	*/
		if( ( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n >= 5 ) && helike.lgSetBenjamin )
			cs = 0.f;
	}
	else
	{
		/* If spin changing collisions are prohibited in the l-mixing routine,
		 * they will fall here, and will have been assigned no collision strength.	
		 * Assign zero for now.	*/
		ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].s != iso.quant_desig[ipHE_LIKE][nelem][ipLo].s );
		cs = 0.f;
		*factor1 = 1.f;
	}

	ASSERT( cs >= 0.f );

	/* Change things to be more like the models by Benjamin, Skillman, and Smits (1999). */
	if( helike.lgSetBenjamin )
	{
		if( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n >=5 && 
			iso.quant_desig[ipHE_LIKE][nelem][ipLo].n != iso.quant_desig[ipHE_LIKE][nelem][ipHi].n )
		{
			ASSERT( cs == 0.f );
		}
		else if( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n >=5 &&
			iso.quant_desig[ipHE_LIKE][nelem][ipLo].s != iso.quant_desig[ipHE_LIKE][nelem][ipHi].s )
		{
			cs = 0.f;
		}
	}

	/*putError(nelem,ipHi,ipLo,IPCOLLIS,-1);*/
	
#	ifdef DEBUG_FUN
	fputs( " <->AtomCSInterp()\n", debug_fp );
#	endif

	return(cs);

}

/* IonCSInterp interpolate on collision strengths for element nelem */
float IonCSInterp( long nelem , long ipHi , long ipLo, float *factor1, const char **where, long Collider  )
{
	float cs;

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

	ASSERT( nelem > ipHELIUM );
	ASSERT( nelem < LIMELM );

	/* init values, better be positive when we exit */
	cs = -1.f; 

	/* this may be used for splitting up the collision strength for collisions between
	 * resolved and collapsed levels.  It may be set somewhere in routine, so set to 
	 * negative value here as flag saying not set */
	*factor1 = -1.f;


	/* >>chng 02 mar 04,  the approximation here for transitions within 2p3P was not needed,
	 * because the Zhang data give these transitions.  They are of the same order, but are 
	 * specific to the three transitions	*/

	/* this branch is ground to n=2 or from n=2 to n=2, for ions only	*/
	/*>>refer Helike	CS	Zhang, Honglin, & Sampson, Douglas H. 1987, ApJS 63, 487	*/
	if( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n==2 
		&& iso.quant_desig[ipHE_LIKE][nelem][ipLo].n<=2 )
	{
		*where = "Zhang";
		*factor1 = 1.;
		
		/* Collisions from gound	*/
		if( ipLo == ipHe1s1S )
		{
			switch( ipHi )
			{
			case 1:	/* to 2tripS	*/
				cs = 0.25f/(float)POW2(nelem+1.);
				break;
			case 2: /* to 2singS	*/
				cs = 0.4f/(float)POW2(nelem+1.);
				break;
			case 3: /* to 2tripP0	*/
				cs = 0.15f/(float)POW2(nelem+1.);
				break;
			case 4: /* to 2tripP1	*/
				cs = 0.45f/(float)POW2(nelem+1.);
				break;
			case 5: /* to 2tripP2	*/
				cs = 0.75f/(float)POW2(nelem+1.);
				break;
			case 6: /* to 2singP	*/
				cs = 1.3f/(float)POW2(nelem+1.);
				break;
			default:
				TotalInsanity();
				break;
			}
			cs *= (float)iso.lgColl_excite[ipHE_LIKE];
		}
		/* collisions from 2tripS to n=2	*/
		else if( ipLo == ipHe2s3S )
		{
			switch( ipHi )
			{
			case 2: /* to 2singS	*/
				cs = 2.75f/(float)POW2(nelem+1.);
				break;
			case 3: /* to 2tripP0	*/
				cs = 60.f/(float)POW2(nelem+1.);
				break;
			case 4: /* to 2tripP1	*/
				cs = 180.f/(float)POW2(nelem+1.);
				break;
			case 5: /* to 2tripP2	*/
				cs = 300.f/(float)POW2(nelem+1.);
				break;
			case 6: /* to 2singP	*/
				cs = 5.8f/(float)POW2(nelem+1.);
				break;
			default:
				TotalInsanity();
				break;
			}
			cs *= (float)iso.lgColl_l_mixing[ipHE_LIKE];
		}
		/* collisions from 2singS to n=2	*/
		else if( ipLo == ipHe2s1S )
		{
			switch( ipHi )
			{
			case 3: /* to 2tripP0	*/
				cs = 0.56f/(float)POW2(nelem+1.);
				break;
			case 4: /* to 2tripP1	*/
				cs = 1.74f/(float)POW2(nelem+1.);
				break;
			case 5: /* to 2tripP2	*/
				cs = 2.81f/(float)POW2(nelem+1.);
				break;
			case 6: /* to 2singP	*/
				cs = 190.f/(float)POW2(nelem+1.);
				break;
			default:
				TotalInsanity();
				break;
			}
			cs *= (float)iso.lgColl_l_mixing[ipHE_LIKE];
		}
		/* collisions from 2tripP0 to n=2	*/
		else if( ipLo == ipHe2p3P0 )
		{
			switch( ipHi )
			{
			case 4: /* to 2tripP1	*/
				cs = 8.1f/(float)POW2(nelem+1.);
				break;
			case 5: /* to 2tripP2	*/
				cs = 8.2f/(float)POW2(nelem+1.);
				break;
			case 6: /* to 2singP	*/
				cs = 3.9f/(float)POW2(nelem+1.);
				break;
			default:
				TotalInsanity();
				break;
			}
			cs *= (float)iso.lgColl_l_mixing[ipHE_LIKE];
		}
		/* collisions from 2tripP1 to n=2	*/
		else if( ipLo == ipHe2p3P1 )
		{
			switch( ipHi )
			{
			case 5: /* to 2tripP2	*/
				cs = 30.f/(float)POW2(nelem+1.);
				break;
			case 6: /* to 2singP	*/
				cs = 11.7f/(float)POW2(nelem+1.);
				break;
			default:
				TotalInsanity();
				break;
			}
			cs *= (float)iso.lgColl_l_mixing[ipHE_LIKE];
		}
		/* collisions from 2tripP2 to n=2	*/
		else if( ipLo == ipHe2p3P2 )
		{
			/* to 2singP	*/
			cs = 19.5f/(float)POW2(nelem+1.) * (float)iso.lgColl_l_mixing[ipHE_LIKE];
		}
		else
			TotalInsanity();

		/* statistical weights included */
	}

	/* this branch, n-same, l-changing collisions */
	else if( (iso.quant_desig[ipHE_LIKE][nelem][ipHi].n == iso.quant_desig[ipHE_LIKE][nelem][ipLo].n ) &&
		(iso.quant_desig[ipHE_LIKE][nelem][ipHi].s == iso.quant_desig[ipHE_LIKE][nelem][ipLo].s ) )
	{
		ASSERT( *factor1 == -1.f );
		*factor1 = 1.f;
		
		/* ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n >= 3 ); */
		ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] );

		if( (iso.quant_desig[ipHE_LIKE][nelem][ipLo].l <=2) &&
			abs(iso.quant_desig[ipHE_LIKE][nelem][ipHi].l - iso.quant_desig[ipHE_LIKE][nelem][ipLo].l)== 1 )
		{
			/* Use the method given in 
			 * >>refer He	CS	Seaton, M. J. 1962, Proc. Phys. Soc. 79, 1105 
			 * statistical weights included */
			cs = (float)CS_l_mixing_S62(ipHE_LIKE, nelem, ipLo, ipHi, (double)phycon.te, Collider); 

			/*fprintf( ioQQQ,"ipHi\t%li\tipLo\t%li\tcs\t%e\n", ipHi, ipLo, cs );*/
		}
		else if( helike.lgCS_Vrinceanu == TRUE )
		{
			if( iso.quant_desig[ipHE_LIKE][nelem][ipLo].l >=3 &&
				iso.quant_desig[ipHE_LIKE][nelem][ipHi].l >=3 )
			{
				/* Use the method given in 
				 * >>refer He	CS	Vrinceanu, D. \& Flannery, M. R. 2001, PhysRevA 63, 032701 
				 * statistical weights included */
				cs = (float)CS_l_mixing_VF01( nelem,
					iso.quant_desig[ipHE_LIKE][nelem][ipLo].n,
					iso.quant_desig[ipHE_LIKE][nelem][ipLo].l,
					iso.quant_desig[ipHE_LIKE][nelem][ipHi].l,
					iso.quant_desig[ipHE_LIKE][nelem][ipLo].s,
					(double)phycon.te,
					Collider );
			}
			else
			{
				cs = 0.f;
			}
		}
		/* this branch, l changing by one */
		else if( abs(iso.quant_desig[ipHE_LIKE][nelem][ipHi].l - iso.quant_desig[ipHE_LIKE][nelem][ipLo].l)== 1)
		{
			/* >>refer	He	cs	Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */
			/* statistical weights included */
			cs = (float)CS_l_mixing_PS64( nelem, ipLo, ipHi, Collider);
		}
		else
		{
			/* l changes by more than 1, but same-n collision */
			cs = 0.f;
		}

		/* ipLo is within 2^3P	*/
		if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
		{
			*factor1 = (2.f*((float)ipLo-3.f)+1.f) / 9.f;
		}

		/* ipHi is within 2^3P	*/
		if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
		{
			*factor1 = (2.f*((float)ipHi-3.f)+1.f) / 9.f;
		}

		*where = "lmix  ";
		cs *= (float)iso.lgColl_l_mixing[ipHE_LIKE];
	}

	/* this branch, n changing collisions for ions */
	else if( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n != iso.quant_desig[ipHE_LIKE][nelem][ipLo].n )
	{
		long int nlo = iso.quant_desig[ipHE_LIKE][nelem][ipLo].n;
		
		/* Even spin-changing collisions are included here.  Is this okay?	*/

		/* for this routine 2s is 1, grnd is 0 */
		if( nlo == 1)
			nlo = 0;

		/* ionic n-changing collision */
		cs = (float)Hion_colldeexc_cs(iso.quant_desig[ipHE_LIKE][nelem][ipHi].n, nlo, nelem, ipHE_LIKE );
		cs *= iso.lgColl_excite[ipHE_LIKE];

		*where = "hydro";
		/* statistical weights included */
	}
	else
	{
		/* what's left are deltaN=0, spin changing collisions.
		 * These have not been accounted for.	*/
		/* Make sure what comes here is what we think it is.	*/
		ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n == iso.quant_desig[ipHE_LIKE][nelem][ipLo].n );
		ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].s != iso.quant_desig[ipHE_LIKE][nelem][ipLo].s );
		cs = 0.f;
		*where = "spin ";
	}

	ASSERT( cs >= 0.f );

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

	/*putError(nelem,ipHi,ipLo,IPCOLLIS,-1);*/

	return(cs);
}


/*CS_l_mixing_S62 - find rate for l-mixing collisions by protons, for neutrals */
/* The S62 stands for Seaton 1962 */
double CS_l_mixing_S62(
	long ipISO,
	long nelem /* the chemical element, 1 for He */,
	long ipLo /* lower level, 0 for ground */,
	long ipHi /* upper level, 0 for ground */,
	double temp,
	long Collider)
{
	/* >>refer	He	l-mixing	Seaton, M.J., 1962, Proc. Phys. Soc. */
	double coll_str, deltaE;

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

	global_temp = temp;
	global_stat_weight = iso.stat[ipISO][nelem][ipLo];
	/* >>chng 05 sep 06, RP  - update energies of excited states */
	/* global_deltaE = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] -
		iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]); */
	global_deltaE = EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg/EN1EV;
	deltaE = global_deltaE;
	global_I_energy_eV = iso.xIsoLevNIonRyd[ipISO][nelem][0];
	global_Collider = Collider;

	ASSERT( TRANS_PROB_CONST*POW2(deltaE/WAVNRYD/EVRYD) > 0. );

	global_osc_str = EmisLines[ipISO][nelem][ipHi][ipLo].Aul/
		(TRANS_PROB_CONST*POW2(deltaE/WAVNRYD/EVRYD));

	/* This returns a thermally averaged collision strength */
	/* The very wide range of energies is needed to get the full
	 * effective velocity range for each impactor.  */
	coll_str =  qg32(    0.0,     1.0, S62_Therm_ave_coll_str);
	coll_str += qg32(    1.0,    10.0, S62_Therm_ave_coll_str);
	coll_str += qg32(   10.0,   100.0, S62_Therm_ave_coll_str);
	coll_str += qg32(  100.0,  1000.0, S62_Therm_ave_coll_str);
	coll_str += qg32( 1000.0, 10000.0, S62_Therm_ave_coll_str);
	coll_str += qg32(10000.0,100000.0, S62_Therm_ave_coll_str);
		
	ASSERT( coll_str > 0. );

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

	/*putError(nelem,ipHi,ipLo,ipCollis,-1);*/

	return coll_str;
}

/* The integrand for calculating the thermal average of collision strengths */
static double S62_Therm_ave_coll_str( double proj_energy )
{

	double integrand, cross_section, /*Rnot,*/ osc_strength, coll_str, zOverB2;
	double deltaE, /* betanot, */ betaone, zeta_of_betaone, /* cs1, */ cs2, temp, stat_weight;
	double I_energy_eV, Dubya;
	/*double Bohr_radius = 5.29177249e-9;*/
	long i, Collider;

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

	deltaE = global_deltaE;
	osc_strength = global_osc_str;
	temp = (double)global_temp;
	stat_weight = global_stat_weight;
	/* Which one to use? */
	I_energy_eV = global_I_energy_eV;
	I_energy_eV = EVRYD;
	Collider = global_Collider;

	/* Rnot = 1.1229*H_BAR/sqrt(ELECTRON_MASS*deltaE*EN1EV)/Bohr_radius; in units of Bohr_radius */

	proj_energy *= ColliderMass[ipELECTRON]/ColliderMass[Collider];
	/* The deltaE here is to make sure that the collider has no less
	 * than the energy difference between the initial and final levels. */
	proj_energy += deltaE;
	Dubya = 0.5*(2.*proj_energy-deltaE);
	ASSERT( Dubya > 0. );

	/* betanot = sqrt(proj_energy/I_energy_eV)*(deltaE/2./Dubya)*Rnot; */

	ASSERT( I_energy_eV > 0. );
	ASSERT( osc_strength > 0. );

	zOverB2 = 0.5*POW2(Dubya/deltaE)*deltaE/I_energy_eV/osc_strength;
	/* This correctly takes account of the mass in equation 31. */
	zOverB2 *= POW2(ColliderMass[ipELECTRON]/ColliderMass[Collider]);

	ASSERT( zOverB2 > 0. );
		
	if( zOverB2 > 100. )
	{
		betaone = sqrt( 1./zOverB2 );
	}
	else if( zOverB2 < 0.54 )
	{
		/* Low betaone approximation */
		betaone = (1./3.)*(log(PI)-log(zOverB2)+1.28);
		if( betaone > 2.38 )
		{
			/* average with this over approximation */
			betaone += 0.5*(log(PI)-log(zOverB2));
			betaone *= 0.5;
		}
	}
	else
	{
		long ip_zOverB2 = 0;
		double zetaOVERbeta2[101] = {
			1.030E+02,9.840E+01,9.402E+01,8.983E+01,8.583E+01,8.200E+01,7.835E+01,7.485E+01,
			7.151E+01,6.832E+01,6.527E+01,6.236E+01,5.957E+01,5.691E+01,5.436E+01,5.193E+01,
			4.961E+01,4.738E+01,4.526E+01,4.323E+01,4.129E+01,3.943E+01,3.766E+01,3.596E+01,
			3.434E+01,3.279E+01,3.131E+01,2.989E+01,2.854E+01,2.724E+01,2.601E+01,2.482E+01,
			2.369E+01,2.261E+01,2.158E+01,2.059E+01,1.964E+01,1.874E+01,1.787E+01,1.705E+01,
			1.626E+01,1.550E+01,1.478E+01,1.409E+01,1.343E+01,1.280E+01,1.219E+01,1.162E+01,
			1.107E+01,1.054E+01,1.004E+01,9.554E+00,9.094E+00,8.655E+00,8.234E+00,7.833E+00,
			7.449E+00,7.082E+00,6.732E+00,6.397E+00,6.078E+00,5.772E+00,5.481E+00,5.202E+00,
			4.937E+00,4.683E+00,4.441E+00,4.210E+00,3.989E+00,3.779E+00,3.578E+00,3.387E+00,
			3.204E+00,3.031E+00,2.865E+00,2.707E+00,2.557E+00,2.414E+00,2.277E+00,2.148E+00,
			2.024E+00,1.907E+00,1.795E+00,1.689E+00,1.589E+00,1.493E+00,1.402E+00,1.316E+00,
			1.235E+00,1.157E+00,1.084E+00,1.015E+00,9.491E-01,8.870E-01,8.283E-01,7.729E-01,
			7.206E-01,6.712E-01,6.247E-01,5.808E-01,5.396E-01};

		ASSERT( zOverB2 >= zetaOVERbeta2[100] );

		/* find beta in the table */
		for( i=0; i< 100; ++i )
		{
			if( ( zOverB2 < zetaOVERbeta2[i] ) && ( zOverB2 >= zetaOVERbeta2[i+1] ) )
			{
				/* found the temperature - use it */
				ip_zOverB2 = i;
				break;
			}
		}

		ASSERT( (ip_zOverB2 >=0) && (ip_zOverB2 < 100) );

		betaone = (zOverB2 - zetaOVERbeta2[ip_zOverB2]) / 
			(zetaOVERbeta2[ip_zOverB2+1] - zetaOVERbeta2[ip_zOverB2]) *
			(pow(10., ((double)ip_zOverB2+1.)/100. - 1.) - pow(10., ((double)ip_zOverB2)/100. - 1.)) +
			pow(10., (double)ip_zOverB2/100. - 1.);

	}

	zeta_of_betaone = zOverB2 * POW2(betaone);
	
	/* cs1 = betanot * gsl_sf_bessel_K0_e(betanot) * gsl_sf_bessel_K1_e(betanot); */
	cs2 = 0.5*zeta_of_betaone + betaone *
		gsl_sf_bessel_K0_e(betaone) * gsl_sf_bessel_K1_e(betaone);
		
	/* cross_section = MIN(cs1, cs2); */
	cross_section = cs2;
	
	/* cross section in units of PI/a_o^s */
	cross_section *= 8. * (I_energy_eV/deltaE) * osc_strength * (I_energy_eV/proj_energy);
	
	/* convert to collision strength */
	coll_str = cross_section * stat_weight * POW2( proj_energy/EVRYD );
	
	integrand = exp( -1.*(proj_energy-deltaE)*EVDEGK/temp ) * coll_str;

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

	return integrand;
}

/*CS_l_mixing_PS64 - find rate for l-mixing collisions by protons, for neutrals */
double CS_l_mixing_PS64(
	long nelem /* the chemical element, 1 for He */,
	long ipLo /* lower level, 0 for ground */,
	long ipHi /* upper level, 0 for ground */,
	long Collider)
{
	/* >>refer	He	l-mixing	Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */
	double cs, Dul,
		TwoLogDebye, TwoLogRc1, 
		factor1, factor2, 
		/* factor3, radint, TwoLogRc2, */
		bestfactor,	factorpart,
		reduced_mass, reduced_mass_2_emass,
		rate,  tau /* ,RMS , deltaE */;
	/* only do proton collisions for now */
	const double ChargIncoming = 1.;

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

	ASSERT( ipHi > ipLo );
	/* In this routine, three different cutoff radii are calculated, and as per PS64,
	 * the least of these is selected.  Must be careful to take least positive result!	*/
	
	/* This reduced mass is in grams.	*/
	reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/
		(dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT;
	/* reduced mass in gm, only do proton collisions for now
	 * this mass always appears relative to the electron mass, so define it that way */
	reduced_mass_2_emass = reduced_mass / ELECTRON_MASS;

	/* This is the lifetime of ipLo.	*/
	tau = helike.Lifetime[nelem][ipLo];

	/* equation 46 of PS64 */
	/* min on density added to prevent this from becoming large and negative
	 * at very high n_e - Pengelly & Seaton did not apply this above
	 * 1e12 cm^-3 */
	/* This is actually 2 times the log of the Debye radius.	*/
	TwoLogDebye = 1.68 + log10( phycon.te / MIN2(1e12 , dense.eden ) );

	/* This is cutoff = 0.72v(tau), where tau is the lifetime of the level nl.	*/
	TwoLogRc1 = 10.95 + log10( phycon.te * tau * tau / reduced_mass_2_emass );
	
#if	0
	/* this is root mean squared velocity */
	/* >>chng 02 mar 17 replaced numbers with constant.	*/
	/* >>chng 03 nov 11 this was wrong...RMSv was two order of magnitude too small!	*/
	RMSv = sqrt( 2. * BOLTZMANN * phycon.te / reduced_mass );
	/* This is the difference between energy levels, in eV.	*/
	deltaE = EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN * WAVNRYD * EVRYD;

	/* This is cutoff = 1.12hbar*v/deltaE.	*/
	/* If deltaE is zero or negative, set this cutoff to zero, so that it will be excluded below.	*/
	if ( deltaE > 0. )
		TwoLogRc2 = 2. * log10( 1.12 * HBAReV * RMSv / deltaE );
	else
		TwoLogRc2 = 0.;

	/* Calculate the hydrogenic radial integral */
	if (iso.quant_desig[ipHE_LIKE][nelem][ipHi].l > iso.quant_desig[ipHE_LIKE][nelem][ipLo].l)
		radint = hri(iso.quant_desig[ipHE_LIKE][nelem][ipHi].n,iso.quant_desig[ipHE_LIKE][nelem][ipHi].l,
			iso.quant_desig[ipHE_LIKE][nelem][ipLo].n, iso.quant_desig[ipHE_LIKE][nelem][ipLo].l, nelem);
	else if (iso.quant_desig[ipHE_LIKE][nelem][ipHi].l < iso.quant_desig[ipHE_LIKE][nelem][ipLo].l)
		radint = hri(iso.quant_desig[ipHE_LIKE][nelem][ipHi].n,iso.quant_desig[ipHE_LIKE][nelem][ipLo].l,
			iso.quant_desig[ipHE_LIKE][nelem][ipLo].n, iso.quant_desig[ipHE_LIKE][nelem][ipHi].l, nelem);
	else 
		radint = 0.;

	ASSERT( radint > 0. );

	Dul = 8./3. * radint * radint * ChargIncoming * ChargIncoming;
#endif

	/* NB - this only works for helium-like ions, the nelem in pow2 is the
	 * residual nuclear charge equation 44 of PS64 */
	Dul = POW2( ChargIncoming / (double)nelem ) * 6. * POW2( (double)iso.quant_desig[ipHE_LIKE][nelem][ipLo].n ) *
		( POW2((double)iso.quant_desig[ipHE_LIKE][nelem][ipLo].n) - 
		POW2((double)iso.quant_desig[ipHE_LIKE][nelem][ipLo].l) - iso.quant_desig[ipHE_LIKE][nelem][ipLo].l - 1);

	ASSERT( Dul > 0. );
	ASSERT( phycon.te  / Dul / reduced_mass_2_emass > 0. );

	factorpart = (11.54 + log10( phycon.te  / Dul / reduced_mass_2_emass ) );

	if( (factor1 = factorpart + TwoLogDebye) <= 0.)
		factor1 = BIGDOUBLE;
	if( (factor2 = factorpart + TwoLogRc1) <= 0.)
		factor2 = BIGDOUBLE;

#if	0
	if( (TwoLogRc2 == 0.) || ( (factor3 = factorpart + TwoLogRc2) <= 0.) )
		factor3 = BIGDOUBLE;
	
	factor3 = BIGDOUBLE;
#endif

	/* Now we find the least positive result.	*/
	bestfactor = MIN2(factor1,factor2);
	/* It would be suspect if bestfactor gets above the double digit range,
	 * so let's make sure less than 100.	*/
	ASSERT( (bestfactor>0.) && (bestfactor<100.) );

	/* This is the rate coefficient.   Units: cm^3 s-1	*/
	rate = 9.93e-6 * sqrt( reduced_mass_2_emass  ) * Dul / phycon.sqrte * bestfactor;
	/* this is the TOTAL rate from nl to nl+/-1. So assume we can
	 * divide by two to get the rate nl to either nl+1 or nl-1. */
	rate /= 2.;

	/* convert rate to collision strength */
	cs = rate / COLL_CONST * phycon.sqrte * iso.stat[ipHE_LIKE][nelem][ipLo] ;

#if	0
	/* convert from collsion strength to cross_section */
	cs *= PI/(2.*iso.quant_desig[ipHE_LIKE][nelem][ipLo].l+1.);
	cs /= POW2( reduced_mass * RMSv / H_BAR );
#endif

	ASSERT( cs > 0. );

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

	/*putError(nelem,ipHi,ipLo,ipCollis,-1);*/

	return cs ;
}

/*CS_l_mixing - find collision strength for l-mixing collisions for neutrals */
/* The VF stands for Vrinceanu & Flannery 2001 */
/* >>refer	He	l-mixing	Vrinceanu, D. & Flannery, M. R. 2001, PhysRevA 63, 032701	*/
/* >>refer	He	l-mixing	Hezel, T. P., Burkhardt, C. E., Ciocca, M., He, L-W., */
/* >>refercon	Leventhal, J. J. 1992, Am. J. Phys. 60, 329 */
double CS_l_mixing_VF01(long int nelem,
						long int n,
						long int l,
						long int lp,
						long int s,
						double temp,
						long int Collider )
{

	double coll_str;

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

	global_z = nelem;
	global_n = n;
	global_l = l;
	global_l_prime = lp;
	global_s = s;
	global_temp = temp;
	global_collider_charge = 1.;
	global_Collider = Collider;

	ASSERT( l != 0 );
	ASSERT( lp != 0 );

	kTRyd = temp / TE1RYD;
	if( helike.lgCS_therm_ave == FALSE )
	{
		/* Must do some thermal averaging for densities greater
		 * than about 10000 and less than about 1e10,
		 * because kT gives significantly different results.
		 * Still, do sparser integration than is done below */
		if( (dense.eden > 10000.) && (dense.eden < 1E10 ) )
		{
			coll_str =  qg32( 0.0, 6.0, Therm_ave_coll_str_int);
		}
		else
		{
			/* Do NOT average over Maxwellian */
			coll_str = collision_strength( kTRyd, FALSE );
		}
	}
	else
	{
		/* DO average over Maxwellian */
		coll_str =  qg32( 0.0, 1.0, Therm_ave_coll_str_int);
		coll_str += qg32( 1.0, 10.0, Therm_ave_coll_str_int);
	}

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

	return coll_str;
}

/* The integrand for calculating the thermal average of collision strengths */
static double Therm_ave_coll_str_int( double EOverKT )
{
	double integrand;
	
#	ifdef DEBUG_FUN
	fputs( " <+>Therm_ave_coll_str_int()\n", debug_fp );
#	endif
	
	integrand = exp( -1.*EOverKT ) * collision_strength( EOverKT * kTRyd, FALSE );

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

	return integrand;
}

static double collision_strength( double velOrEner, long int paramIsRedVel )
{

	double cross_section, coll_str, RMSv, aveRadius, reduced_vel, E_Proj_Ryd;
	double ConstantFactors, reduced_mass, CSIntegral, stat_weight;
	double ColliderCharge = global_collider_charge;
	double quantum_defect1, quantum_defect2, omega_qd1, omega_qd2, omega_qd;
	double reduced_b_max, reduced_b_min, alphamax, alphamin, step, alpha1 /*, alpha2*/;
	double Bohr_rad = 5.29177249e-9;

	long nelem, n, l, lp, s, ipLo, ipHi, Collider = global_Collider;

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

	nelem = global_z;
	n = global_n;
	l = global_l;
	lp = global_l_prime;
	s = global_s;
	stat_weight = (2.*l + 1.) * (2.*s + 1.);

	/*ASSERT( n != 0 );
	ASSERT( l != 0 );
	ASSERT( lp != 0 );*/
	/* these shut up the lint, already done above */
	assert( n != 0 );
	assert( l != 0 );
	assert( lp != 0 );
	
	/* This reduced mass is in grams.	*/
	reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/
		(dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT;
	ASSERT( reduced_mass > 0. );

	/* this is root mean squared velocity */
	/* use this as projectile velocity for thermally-averaged cross-section? */
	/*RMSv = sqrt( 2. * BOLTZMANN * global_temp / reduced_mass );*/
	aveRadius = (Bohr_rad/nelem)*POW2((double)n);
	ASSERT( aveRadius < 1.e-4 );
	/* >>chng 05 jul 14, as per exchange with Ryan Porter & Peter van Hoof, avoid
	 * roundoff error and give ability to go beyond zinc */
	/*ASSERT( aveRadius >=  Bohr_rad );*/
	ASSERT( aveRadius > 3.9/LIMELM * Bohr_rad );
	global_an = aveRadius;
	
	RMSv = nelem*POW2(ELEM_CHARGE_ESU)/n/H_BAR;
	ASSERT( RMSv > 0. );

	ASSERT( ColliderMass[Collider] > 0. );

	if( paramIsRedVel == TRUE )
	{
		/* velOrEner is a reduced velocity */
		reduced_vel = velOrEner;
		/* The proton mass is needed here because the ColliderMass array is
		 * expressed in units of the proton mass, but here we need absolute mass. */
		E_Proj_Ryd = 0.5 * POW2( reduced_vel * RMSv ) * ColliderMass[Collider] *
			PROTON_MASS / EN1RYD;
	}
	else
	{	
		/* velOrEner is a projectile energy in Rydbergs */
		E_Proj_Ryd = velOrEner;
		reduced_vel = sqrt( 2.*E_Proj_Ryd*EN1RYD/ColliderMass[Collider]/PROTON_MASS )/RMSv;
	}
	
	/* put limits on the reduced velocity.   These limitis should be more than fair. */
	ASSERT( reduced_vel > 1.e-10 );
	ASSERT( reduced_vel < 1.e10 );

	global_red_vel = reduced_vel;

	/* Factors outside integral	*/
	ConstantFactors = 4.5*PI*POW2(ColliderCharge*aveRadius/reduced_vel);

	/* Reduced here means in units of aveRadius: */
	reduced_b_min = 1.5 * ColliderCharge / reduced_vel;
	ASSERT( reduced_b_min > 1.e-10 );
	ASSERT( reduced_b_min < 1.e10 );
	
	/* Debye radius: appears to be too large, results in 1/v^2 variation. */
	/* reduced_b_max = sqrt( BOLTZMANN*global_temp/ColliderCharge/dense.eden ) 
		/ (PI2*ELEM_CHARGE_ESU); */

	ipLo = QuantumNumbers2Index[nelem][n][l][s];
	ipHi = QuantumNumbers2Index[nelem][n][l][s];

	{
		quantum_defect1  = (double)n- (double)nelem/sqrt(iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo]);
		quantum_defect2  = (double)n- (double)nelem/sqrt(iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHi]);
		
		/* The magnitude of each quantum defect must be between zero and one. */
		ASSERT( fabs(quantum_defect1)  < 1.0 );
		ASSERT( fabs(quantum_defect1)  > 0.0 );
		ASSERT( fabs(quantum_defect2)  < 1.0 );
		ASSERT( fabs(quantum_defect2)  > 0.0 );
	}

	/* The quantum defect precession frequencies */
	omega_qd1 = fabs( 5.* quantum_defect1 * (1.-0.6*POW2((double)l/(double)n)) / POW3( (double)n ) / (double)l );
	omega_qd2 = fabs( 5.* quantum_defect2 * (1.-0.6*POW2((double)l/(double)n)) / POW3( (double)n ) / (double)l );
	/* Take the average for the two levels, for reciprocity. */
	omega_qd = 0.5*( omega_qd1 + omega_qd2 );
	
	ASSERT( omega_qd > 0. );

	reduced_b_max = sqrt( 1.5 * ColliderCharge * n / omega_qd )/aveRadius;
	reduced_b_max = MAX2( reduced_b_max, reduced_b_min );
	ASSERT( reduced_b_max > 0. );
	global_bmax = reduced_b_max*aveRadius;
	alphamin = 1.5*ColliderCharge/(reduced_vel * reduced_b_max);
	alphamax = 1.5*ColliderCharge/(reduced_vel * reduced_b_min);

	ASSERT( alphamin > 0. );
	ASSERT( alphamax > 0. );

	alphamin = MAX2( alphamin, 1.e-30 );
	alphamax = MAX2( alphamax, 1.e-20 );

	CSIntegral = 0.;

	if( alphamax > alphamin )
	{

		step = (alphamax - alphamin)/5.;
		alpha1 = alphamin;
		CSIntegral += qg32(  alpha1, alpha1+step, L_mix_integrand);
		CSIntegral += qg32(  alpha1+step, alpha1+4.*step, L_mix_integrand);
	}

	/* Calculate cross section */
	cross_section = ConstantFactors * CSIntegral;

	/* convert to collision strength */
	coll_str = cross_section * stat_weight / PI / Bohr_rad / Bohr_rad;
	coll_str *= 0.25 * POW2( E_Proj_Ryd );

	coll_str = MAX2( SMALLFLOAT, coll_str);

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

	return coll_str;

}

static double L_mix_integrand( double alpha )
{
	double integrand, deltaPhi, b;

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

	ASSERT( alpha >= 1.e-30 );
	ASSERT( global_bmax > 0. );
	ASSERT( global_red_vel > 0. );

	/* >>refer He l-mixing Kazansky, A. K. & Ostrovsky, V. N. 1996, JPhysB: At. Mol. Opt. Phys. 29, 3651 */
	b = 1.5*global_collider_charge*global_an/(global_red_vel * alpha);
	if( b < global_bmax )
	{
		deltaPhi = -1.*PI + 2.*asin(b/global_bmax);
	}
	else
	{
		deltaPhi = 0.;
	}
	integrand = 1./alpha/alpha/alpha;
	integrand *= StarkCollTransProb( global_n, global_l, global_l_prime, alpha, deltaPhi );
	
#	ifdef DEBUG_FUN
	fputs( " <->L_mix_integrand()\n", debug_fp );
#	endif
	
	return integrand;
}

static double StarkCollTransProb( long n, long l, long lp, double alpha, double deltaPhi)
{
	double probability;
	double cosU1, cosU2, sinU1, sinU2, cosChiOver2, sinChiOver2, cosChi, A, B;
	
#	ifdef DEBUG_FUN
	fputs( " <+>StarkCollTransProb()\n", debug_fp );
#	endif

	ASSERT( alpha > 0. );

	/* These are defined on page 11 of VF01 */ 
	cosU1 = 2.*POW2((double)l/(double)n) - 1.;
	cosU2 = 2.*POW2((double)lp/(double)n) - 1.;

	sinU1 = sqrt( 1. - cosU1*cosU1 );
	sinU2 = sqrt( 1. - cosU2*cosU2 );


	cosChiOver2 = (1. + alpha*alpha*cos( sqrt(1.+alpha*alpha) * deltaPhi ) )/(1.+alpha*alpha);
	sinChiOver2 = sqrt( 1. - cosChiOver2*cosChiOver2 );
	cosChi = 2. * POW2( cosChiOver2 ) - 1.;

	if( l == 0 )
	{
		if( -1.*cosU2 - cosChi < 0. )
		{
			probability = 0.;
		}
		else
		{
			/* Here is the initial state S case */
			ASSERT( sinChiOver2 > 0. );
			ASSERT( sinChiOver2*sinChiOver2 > POW2((double)lp/(double)n) );
			/* This is equation 35 of VF01.  There is a factor of hbar missing in the denominator of the
			 * paper, but it's okay if you use atomic units (hbar = 1). */
			probability = (double)lp/(POW2((double)n)*sinChiOver2*sqrt( POW2(sinChiOver2) - POW2((double)lp/(double)n) ) );
		}
	}
	else
	{
		double OneMinusCosChi = 1. - cosChi;

		if( OneMinusCosChi == 0. )
		{
			double hold = sin( deltaPhi / 2. );
			/* From approximation at bottom of page 10 of VF01. */
			OneMinusCosChi = 8. * alpha * alpha * POW2( hold );
		}

		if( OneMinusCosChi == 0. )
		{
			probability = 0.;
		}
		else
		{
			/* Here is the general case */
			A = (cosU1*cosU2 - sinU1*sinU2 - cosChi)/OneMinusCosChi;
			B = (cosU1*cosU2 + sinU1*sinU2 - cosChi)/OneMinusCosChi;
	
			ASSERT( B > A );

			/* These are the three cases of equation 34. */
			if( B <= 0. )
			{
				probability = 0.;
			}
			else
			{
				ASSERT( POW2( sinChiOver2 ) > 0. );

				probability = 2.*lp/(PI* /*H_BAR* */ n*n*POW2( sinChiOver2 ));

				if( A < 0. )
				{
					probability *= ellip_int_K( sqrt(B/(B-A)) );
					probability /= sqrt( B - A );
				}
				else
				{
					probability *= ellip_int_K( sqrt((B-A)/B) );
					probability /= sqrt( B );
				}
			}
		}

	}

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

	return probability;

}


/* This is a modified version of two routines
 * in the GNU science library, version 1.4:
 *		gsl_sf_ellint_Kcomp_e
 *		gsl_sf_ellint_RF_e
 * The routines include references:
 *		[Carlson, Numer. Math. 33 (1979) 1, (4.5)]
 *		[Abramowitz+Stegun, 17.3.33]
 */
static double ellip_int_K( double k )
{
	double result;

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

	ASSERT( k*k < 1.0 );
  
	if(k*k >= 1.0 - sqrt(DBL_EPSILON))
	{
		/* [Abramowitz+Stegun, 17.3.33] */
		const double y = 1.0 - k*k;
		const double a[] = { 1.38629436112, 0.09666344259, 0.03590092383 };
		const double b[] = { 0.5, 0.12498593597, 0.06880248576 };
		const double ta = a[0] + y*(a[1] + y*a[2]);
		const double tb = -log(y) * (b[0] * y*(b[1] + y*b[2]));
		result = ta + tb;
	}
	else
	{
		double y = 1.0 - k*k;
		
		/*const double lolim = 5.0 * DBL_MIN;
		const double uplim = 0.2 * DBL_MAX;*/

		const double c1 = 1.0 / 24.0;
		const double c2 = 3.0 / 44.0;
		const double c3 = 1.0 / 14.0;
		double xn = 0.0;
		double yn_ = y;
		double zn = 1.0;
		double mu =0., xndev=0., yndev=0., zndev=0., e2_, e3, s;
		/* yn and e2 hide VS.net internal symbol yn and cddefines symbol e2, respectively. 
		 * An underscore at the end prevents this collision. */
		
		ASSERT( y >= 5.0 * DBL_MIN );
		ASSERT( y < 0.2 * DBL_MAX );

		while(1)
		{
			double epslon, lamda;
			double xnroot, ynroot, znroot;
			mu = (xn + yn_ + zn) / 3.0;
			ASSERT( mu > 0. );
			xndev = 2.0 - (mu + xn) / mu;
			yndev = 2.0 - (mu + yn_) / mu;
			zndev = 2.0 - (mu + zn) / mu;
			epslon = MAX3(fabs(xndev), fabs(yndev), fabs(zndev));
			if (epslon < 0.001) break;
			xnroot = sqrt(xn);
			ynroot = sqrt(yn_);
			znroot = sqrt(zn);
			lamda = xnroot * (ynroot + znroot) + ynroot * znroot;
			xn = (xn + lamda) * 0.25;
			yn_ = (yn_ + lamda) * 0.25;
			zn = (zn + lamda) * 0.25;
			ASSERT( xn>0. && yn_>0. && zn>0. );
		}
		e2_ = xndev * yndev - zndev * zndev;
		e3 = xndev * yndev * zndev;
		s = 1.0 + (c1 * e2_ - 0.1 - c2 * e3) * e2_ + c3 * e3;
		result = s / sqrt(mu);
	}

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

	return result;

}
/*lint +e661 possible access of out-of-bounds pointer */

