/* 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 */
/*hydro_vs_deexcit compute collision strength for collisional deexcitation for hydrogen atom, from
 * >>refer	H1	collision	Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940
 *hydro_vs_ioniz generate hydrogenic collisional ionization rate coefficients 
 * for quantum number n */
/*Hion_coll_ioniz_ratecoef calculate hydrogenic ionization rates for ions 
 * with all n, and Z*/
/*Hion_colldeexc_cs collision strength for collisional de-excitation for any levels of hydrogenic Z*/
#include "cddefines.h"
#include "phycon.h"
#include "physconst.h"
#include "hydrooscilstr.h"
#include "iso.h"
#include "hydro_vs_rates.h"
#include "helike_cs.h"

static double hydro_vs_coll_str( double energy );
static double Therm_ave_coll_str_int_VS80( double EOverKT );

static long global_ipISO, global_nelem, global_ipHi, global_ipLo, global_Collider;
static double global_temp;
/* 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};

/*
 Neither of these rates can be modified to account for impact by non-electrons because they 
 are fits to thermally-averaged rates for electron impact...it can not be unravelled to 
 expose a projectile energy that can then be scaled to account for other projectiles.
 Instead, we have included the original cross section formula (eq 14) from 
 Vriens & Smeets (1980) below.
*/

/*hydro_vs_deexcit compute collision strength for collisional deexcitation for hydrogen atom, 
 * from Vriens and Smeets */
double hydro_vs_deexcit(
						long int ipISO,
						long int nelem,
						long int ipHi,
						long int ipLo
		/* long int ihigh, long int ilow */
		)
{
	double Anp, 
	  bn, 
	  Bnp,
	  delta, 
	  /*eni, 
	  enp,*/ 
	  Eni, 
	  Enp, 
	  gamma, 
	  hdexct_v, 
	  xlow, 
	  p, 
	  rate, 
	  ryd, 
	  s, 
	  tev;

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

	/* function written by Jason Ferguson to compute the Vriens and Smeets
	 * rate coeff. for collisional dexcitation between any two levels of H.
	 * valid for all ihigh, ilow
	 * at end converts this deexcitation rate to collision strength
	 * */
	/* p = (double)(ihigh); */
	p = (double)iso.quant_desig[ipISO][nelem][ipHi].n;
	/* xlow = (double)(ilow); */
	xlow = (double)iso.quant_desig[ipISO][nelem][ipLo].n;
	tev = EVRYD/TE1RYD*phycon.te;
	ryd = EVRYD;
	s = fabs(p-xlow);

	/* >>chng 04 dec 16, several things were incorrect.  There was no distinction between Enp and enp, 
	 * but the paper and the formulas have one.  There were also a factor ryd^2 missing in Bnp. */

	Enp = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]);
	Eni = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][ipLo];

	/*
	enp = Enp / tev;
	eni = Eni / tev;

	enp = -1./p/p + 1./xlow/xlow;
	eni = 1./xlow/xlow;
	*/
	/* This is an absorption oscillator strength. */
	Anp = 2./Enp*HydroOscilStr(xlow,p);

	bn = 1.4*log(xlow)/xlow - .7/xlow - .51/xlow/xlow + 1.16/xlow/xlow/xlow - 0.55/xlow/xlow/xlow/xlow;

	Bnp = 4.*ryd*ryd/p/p/p*(1./Enp/Enp + 4./3.*Eni/POW3(Enp) + bn*Eni*Eni/powi(Enp,4));

	delta = exp(-Bnp/Anp) + 0.06*s*s/p/xlow/xlow;

	gamma = ryd*log(1.+POW3(xlow)*tev/ryd)*(3. + 11.*POW2(s/xlow))/
		(6. + 1.6*p*s + .3/s/s + 0.8*(pow(p,1.5))/(pow(s,.5))*fabs(s-0.6));
	/* statistical weights taken care of in calling program:
	 * multiply this by xlow*xlow/p/p for true rate */
	rate = (1.6e-7)*sqrt(tev)/(tev + gamma)*(Anp*log(0.3*tev/ryd+delta) + Bnp);
	
	hdexct_v = rate;

	/* for collision strength
	 * stat weight taken care of in calling program */
	hdexct_v = rate*phycon.sqrte/COLL_CONST;

	ASSERT( hdexct_v >= 0. );

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

/*hydro_vs_deexcit compute collision rate for collisional deexcitation for hydrogen atom, 
 * from Vriens and Smeets */
double hydro_vs_excit(
						long int ipISO,
						long int nelem,
						long int ipHi,
						long int ipLo
		/* long int nhigh, long int nlow */
		)
{
	double Apn,
	  bp,
	  Bpn,
	  delta,
	  epn, 
	  Epi,
	  Epn,
	  gamma, 
	  /*hdexct_v, */
	  xlow, 
	  n, 
	  rate, 
	  ryd, 
	  s, 
	  tev;

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

	/* Computes the Vriens and Smeets
	 * rate coeff. for collisional excitation between any two levels of H.
	 * valid for all nhigh, nlow
	 * at end converts this excitation rate to collision strength	*/

	/* n = (double)(nhigh); */
	n = (double)iso.quant_desig[ipISO][nelem][ipHi].n;
	/* xlow = (double)(nlow); */
	xlow = (double)iso.quant_desig[ipISO][nelem][ipLo].n;

	tev = EVRYD/TE1RYD*phycon.te;
	ryd = EVRYD;
	s = fabs(n-xlow);

	/* >>chng 04 dec 16, several things were incorrect.  There was no distinction between Enp and enp, 
	 * but the paper and the formulas have one.  There were also a factor ryd^2 missing in Bnp. */

	Epn = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]);
	Epi = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][ipLo];

	epn = Epn / tev;

	/* This is an absorption oscillator strength. */
	Apn = 2./Epn*HydroOscilStr(xlow,n);

	bp = 1.4*log(xlow)/xlow - .7/xlow - .51/xlow/xlow + 1.16/xlow/xlow/xlow - 0.55/xlow/xlow/xlow/xlow;

	Bpn = 4.*ryd*ryd/n/n/n*(1./Epn/Epn + 4./3.*Epi/POW3(Epn) + bp*Epi*Epi/powi(Epn,4));

	delta = exp(-Bpn/Apn) + 0.06*s*s/n/xlow/xlow;

	gamma = ryd*log(1.+POW3(xlow)*tev/ryd)*(3. + 11.*POW2(s/xlow))/
		(6. + 1.6*n*s + .3/s/s + 0.8*(pow(n,1.5))/(pow(s,.5))*fabs(s-0.6));
	
	rate = (1.6e-7)*sqrt(tev)/(tev + gamma)*exp(-epn)*(Apn*log(0.3*tev/ryd+delta) + 
	  Bpn);
	
	/*hdexct_v = rate;

	hdexct_v = rate*phycon.sqrte/COLL_CONST; */

	ASSERT( rate >= 0. );

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


/* VS80 stands for Vriens and Smeets 1980 */
/* This routine calculates thermally-averaged collision strengths. */
double CS_VS80( long int ipISO, long int nelem, long int ipHi, long int ipLo, double temp, long int Collider )
{
	double coll_str;

	global_ipISO = ipISO;
	global_nelem = nelem;
	global_ipHi = ipHi;
	global_ipLo = ipLo;
	global_temp = temp;
	global_Collider = Collider;

#if	0
	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_VS80);
		}
		else
		{
			/* Do NOT average over Maxwellian */
			coll_str = hydro_vs_coll_str( EVRYD*global_temp/TE1RYD, FALSE );
		}
	}
	else
#endif
	{
		/* DO average over Maxwellian */
		coll_str =  qg32( 0.0, 1.0, Therm_ave_coll_str_int_VS80);
		coll_str += qg32( 1.0, 10.0, Therm_ave_coll_str_int_VS80);
	}

	return coll_str;
}

/* The integrand for calculating the thermal average of collision strengths */
static double Therm_ave_coll_str_int_VS80( double EOverKT )
{
	double integrand;
	
	integrand = exp( -1.*EOverKT ) * hydro_vs_coll_str( EOverKT * EVRYD*global_temp/TE1RYD );

	return integrand;
}

/*hydro_vs_deexcit compute collision strength for collisional deexcitation for hydrogen atom, 
 * from Vriens and Smeets */
static double hydro_vs_coll_str( double energy )
{
	double Apn, bp, Bpn, delta;
	double /* epi, epn,*/ Epi, Epn;
	double gamma, p, n;
	double ryd, s /* , tev, temp*/;
	double cross_section, coll_str, stat_weight;
	long ipISO, nelem, ipHi, ipLo, Collider;

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

	ipISO =	global_ipISO;
	nelem = global_nelem;
	ipHi = global_ipHi;
	ipLo = global_ipLo;
	/* temp = global_temp; */
	Collider = global_Collider;

	stat_weight = iso.stat[ipISO][nelem][ipLo];

	/* Computes the Vriens and Smeets
	 * rate coeff. for collisional dexcitation between any two levels of H.
	 * valid for all nhigh, nlow
	 * at end converts this excitation rate to collision strength	*/

	/*TODO	2	make these effective quantum numbers? */
	/* n = (double)(nhigh); */
	n = (double)iso.quant_desig[ipISO][nelem][ipHi].n;
	/* xlow = (double)(nlow); */
	p = (double)iso.quant_desig[ipISO][nelem][ipLo].n;

	/* tev = EVRYD/TE1RYD*temp; */
	ryd = EVRYD;
	s = fabs(n-p);

	Epn = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]);
	Epi = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][ipLo];

	/*
	epn = Epn / tev;
	epi = Epi / tev;
	*/

	/* This is an absorption oscillator strength. */
	Apn = 2./Epn*HydroOscilStr(p,n);

	bp = 1.4*log(p)/p - .7/p - .51/p/p + 1.16/p/p/p - 0.55/p/p/p/p;

	Bpn = 4.*ryd*ryd/n/n/n*(1./Epn/Epn + 4./3.*Epi/POW3(Epn) + bp*Epi*Epi/powi(Epn,4));

	delta = exp(-Bpn/Apn) - 0.4*Epn/ryd;

	gamma = ryd*(8. + 23.*POW2(s/p))/
		(8. + 1.1*n*s + 0.8/s/s + 0.4*(pow(n,1.5))/(pow(s,0.5))*fabs(s-1.0));
	
	/* Scale the energy to get an equivalent electron energy. */
	energy *= ColliderMass[ipELECTRON]/ColliderMass[Collider];

	/* cross section in units of PI*a_o^2 */
	if( energy/2/ryd+delta <= 0 )
	{
		cross_section = 0.;
	}
	else
	{
		cross_section = 2.*ryd/(energy + gamma)*(Apn*log(energy/2/ryd+delta) + Bpn);
		cross_section = MAX2( cross_section, 0. );
	}
	
	/* convert to collision strength */
	coll_str = cross_section * POW2( energy/EVRYD ) * stat_weight;
	/* coll_str /= (PI * Bohr_radius * Bohr_radius); */
	
#	ifdef DEBUG_FUN
	fputs( " <->hydro_vs_coll_str()\n", debug_fp );
#	endif

	return( coll_str );
}

/*hydro_vs_ioniz generate hydrogenic collisional ionization rate coefficients 
 * for quantum number n */
double hydro_vs_ioniz( long int ipISO, long int nelem, long int level )
{
	double coef, 
	  denom, 
	  eps, 
	  smeets_v, 
	  tryd;

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

	ASSERT( nelem <= 1 );

	/* a function written to calculate the rate coefficients 
	 * for hydrogen collisional ionizations from
	 * Jason Ferguson, summer 94
	 * valid for all n
	 * >>refer	H1	collision	Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940
	 * */
	tryd = EVRYD/TE1RYD*phycon.te;
	eps = TE1RYD/phycon.te*iso.xIsoLevNIonRyd[ipISO][nelem][level];
	denom = pow(eps,2.33) + 4.38*pow(eps,1.72) + 1.32*eps;
	
	coef = 9.56e-6/pow(tryd,1.5)/denom*iso.ConBoltz[ipISO][nelem][level];
	
	smeets_v = coef;

	ASSERT( smeets_v >= 0. );

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

/*Hion_coll_ioniz_ratecoef calculate hydrogenic ionization rates for all n, and Z*/
double Hion_coll_ioniz_ratecoef(
		/* the isoelectronic sequence */
		long int ipISO ,
		/* element, >=1 since only used for ions 
		 * nelem = 1 is helium the least possible charge */
		long int nelem,
		/* principal quantum number, > 1
		 * since only used for excited states */
		long int level)
{
	long int i,n, 
	  _r;
	double H, 
	  HydColIon_v, 
	  Rnp, 
	  charge,
	  chim, 
	  eone, 
	  etwo, 
	  ethree, 
	  g, 
	  rate, 
	  rate2, 
	  /* ryd, */
	  t1, 
	  t2, 
	  t3, 
	  t4, 
	  tev, 
	  xn, 
	  y;
	static double arrH[4], 
	  arrRnp[8], 
	  arrg[10];
	static double small = 0.;
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp0[] = {1.48,3.64,5.93,8.32};
		for( i=1, _r = 0; i <= 4; i++ )
		{
			arrH[i-1] = _itmp0[_r++];
			}
		}
		{ static double _itmp1[] = {2.20,1.90,1.73,1.65,1.60,1.56,1.54,
		  1.52};
		for( i=1, _r = 0; i <= 8; i++ )
		{
			arrRnp[i-1] = _itmp1[_r++];
			}
		}
		{ static double _itmp2[] = {0.8675,0.932,0.952,0.960,0.965,
		  0.969,0.972,0.975,0.978,0.981};
		for( i=1, _r = 0; i <= 10; i++ )
		{
			arrg[i-1] = _itmp2[_r++];
			}
		}
		_aini = 0;
	}

#	ifdef DEBUG_FUN
	fputs( "<+>Hion_coll_ioniz_ratecoef()\n", debug_fp );
#	endif
	/*calculate hydrogenic ionization rates for all n, and Z
	 * >>refer	HI	cs	Allen 1973, Astro. Quan. for low Te.
	 * >>refer	HI	cs	Sampson and Zhang 1988, ApJ, 335, 516 for High Te.
	 * */

	/*TODO	2	This routine, like the above one, needs to use proper energy.	*/
	charge = nelem - ipISO;
	/* this routine only for ions, nelem=0 is H, nelem=1 he, etc */
	ASSERT( charge > 0);
	/* this is quantum level, not for ground state (n=1), only
	 * n=2 or higher */
	n = iso.quant_desig[ipISO][nelem][level].n;
	ASSERT( n>1 );

	if( n > 4 )
	{
		H = 2.15*n;
	}
	else
	{
		H = arrH[n-1];
	}

	if( n > 8 )
	{
		Rnp = 1.52;
	}
	else
	{
		Rnp = arrRnp[n-1];
	}

	if( n > 10 )
	{
		g = arrg[9];
	}
	else
	{
		g = arrg[n-1];
	}

	tev = EVRYD/TE1RYD*phycon.te;
	/* ryd = EVRYD; */
	xn = (double)n;
	/* >>chng 19 dec 02, use proper energy!	*/
	/*chim = POW2(charge+1.)*ryd/xn/xn;*/
	chim = EVRYD * iso.xIsoLevNIonRyd[ipISO][nelem][level];
	y = chim/tev;

	eone = ee1(y);
	etwo = iso.ConBoltz[ipISO][nelem][level] - y*eone;
	ethree = (iso.ConBoltz[ipISO][nelem][level] - y*etwo)/2.;

	t1 = 1/xn*eone;
	t2 = 1./3./xn*(iso.ConBoltz[ipISO][nelem][level] - y*ethree);
	t3 = 3.*H/xn/(3. - Rnp)*(y*etwo - 2.*y*eone + iso.ConBoltz[ipISO][nelem][level]);
	t4 = 3.36*y*(eone - etwo);
	rate = 7.69415e-9*(double)phycon.sqrte*9.28278e-3*powi(xn/(charge+1),4)*g*y;
	rate *= t1 - t2 + t3 + t4;
	/* TODO	2	sometimes rate is negative...is that okay?	*/
	rate = MAX2(rate,small);

	rate2 = 2.1e-8*phycon.sqrte/chim/chim*dsexp(2.302585*5040.*
	  chim/(double)phycon.te);

	rate2 = MAX2(rate2,small);

	/* Take the lowest of the two, they fit nicely together... */
	/* >>chng 10 Sept 02, sometimes one of these is zero and the other is positive.
	 * in that case take the bigger one.	*/
	if( rate==0. || rate2==0. )
		HydColIon_v = MAX2(rate,rate2);
	else
		HydColIon_v = MIN2(rate,rate2);

	ASSERT( HydColIon_v >= 0. );

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


/*Hion_colldeexc_cs collision strength for collisional de-excitation for any levels of hydrogenic Z*/
double Hion_colldeexc_cs(long int ipHi, 
  long int ipLo, 
  /* element number, H is 0 */
  long int nelem,
  long int ipISO )
{
	long int i, 
	  nHi, 
	  nLo, 
	  _r;
	double A, 
	  An, 
	  D, 
	  H, 
	  HydColDwn_v, 
	  R, 
	  c, 
	  charge ,
	  e1, 
	  exp_int_2, 
	  gHi, 
	  gLo, 
	  rate, 
	  ryd, 
	  sexpy ,
	  tev, 
	  xn, 
	  xnp2, 
	  xp, 
	  y;
	static double arrAn[10], 
	  arrH[4], 
	  arrR[10];
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp0[] = {1.48,3.64,5.93,8.32};
		for( i=1, _r = 0; i <= 4; i++ )
		{
			arrH[i-1] = _itmp0[_r++];
			}
		}
		{ static double _itmp1[] = {1.30,0.59,0.38,0.286,0.229,0.192,
		  0.164,0.141,0.121,0.105};
		for( i=1, _r = 0; i <= 10; i++ )
		{
			arrAn[i-1] = _itmp1[_r++];
			}
		}
		{ static double _itmp2[] = {1.83,1.60,1.53,1.495,1.475,1.46,
		  1.45,1.45,1.46,1.47};
		for( i=1, _r = 0; i <= 10; i++ )
		{
			arrR[i-1] = _itmp2[_r++];
			}
		}
		_aini = 0;
	}

#	ifdef DEBUG_FUN
	fputs( "<+>Hion_colldeexc_cs()\n", debug_fp );
#	endif
	/* rate coefficient for collisional de-excitation for any levels of hydrogenic Z
	 * written by Jason Ferguson to compute the 
	 * >>refer	Hlike	cs	Sampson and Zhang 1988, ApJ, 335, 516
	 * rate coeff. for collisional excitation between any two levels of 
	 * hydrogenic atoms charge Z.  Valid for all ipLo, ipHi. */

	/* net charge, should be 2 for He+, where nelem==1 and ipISO==2 */
	charge = nelem + 1 - ipISO;

	/*begin sanity check */
	if( (charge < 0 || ipLo >= ipHi) || ipLo*ipHi < 0 )
	{
		fprintf( ioQQQ, "Hion_colldeexc_cs called with impossible parameters.\n" );
		fprintf( ioQQQ, "%ld %ld %f\n", ipHi, ipLo, charge );
		puts( "[Stop in Hion_colldeexc_cs]" );
		cdEXIT(EXIT_FAILURE);
	}
	/*end sanity check */

	/* we need to check to see what the gs and quantums are
	 * ipHi will always be 4 and greater. */
	nHi = ipHi;
	gHi = 2.*nHi*nHi;
	if( ipLo == 0 )
	{
		nLo = 1;
		gLo = 2.;
	}
	else if( ipLo == 1 )
	{
		nLo = 2;
		gLo = 2.;
	}
	else if( ipLo == 2 )
	{
		nLo = 2;
		gLo = 6.;
	}
	else
	{
		nLo = ipLo;
		gLo = 2.*nLo*nLo;
	}

	/* call hydrogen if charge = 1 */
	if( charge == 1 )
	{
		HydColDwn_v = hydro_vs_deexcit(ipISO, nelem, ipHi, ipLo)*gLo;
		
#		ifdef DEBUG_FUN
		fputs( " <->Hion_colldeexc_cs()\n", debug_fp );
#		endif
		return( HydColDwn_v );
	}

	if( nLo > 4 )
	{
		H = 2.15*nLo;
	}
	else
	{
		H = arrH[nLo-1];
	}
	if( nLo > 10 )
	{
		An = 1./(float)(nLo);
		R = 1.50;
	}
	else
	{
		An = arrAn[nLo-1];
		R = arrR[nLo-1];
	}

	xp = (double)(nHi);
	xn = (double)(nLo);
	xnp2 = POW2(xn/xp);
	tev = EVRYD/TE1RYD*phycon.te;
	ryd = EVRYD;

	y = POW2(charge)*ryd*(1./POW2(xn) - 1./POW2(xp))/tev;
	/*>>>chng 98 dec 17, redefine ee1 so always first exponential integral*/
	sexpy = dsexp(y);
	/* >>>chng 99 apr 15, as per Peter van Hoof comment, e1 underflowed due to 
	 * exp, then div by 0.  just set downward rate to zero when exp is this large */
	if( sexpy <= 0. )
	{
		/* do safe version */
		e1 = ee1_safe(y);
		
		A = 4.*powi(xn,4)/(1 - xnp2)*HydroOscilStr(xn,xp);
		D = A*H*(pow(1 - xnp2,R) - An*xnp2);
		c = 1.12*xn*A*(1 - xnp2);
		if( xp - xn == 1. )
		{
			c *= dsexp(0.006*powi(xn-1.,6)/charge);
		}

		rate = COLL_CONST/(POW2(xn*charge)*phycon.sqrte);
		rate *= (A + y*c - y*D)*e1 + D;

		/* Convert from excitation to de-excitation */
		/* >>>chng 99 apr 09, had not changed following, so low temp
		 * had too small a rate, and there was a discontinuity for he ii
		 * at T=30700 */
		HydColDwn_v = rate*gLo/gHi;
	}
	else
	{
		e1 = ee1(y);
		exp_int_2 = sexpy - y*e1;

		A = 4.*powi(xn,4)/(1 - xnp2)*HydroOscilStr(xn,xp);
		D = A*H*(pow(1 - xnp2,R) - An*xnp2);
		c = 1.12*xn*A*(1 - xnp2);
		if( xp - xn == 1. )
		{
			c *= dsexp(0.006*powi(xn-1.,6)/charge);
		}

		rate = COLL_CONST/(POW2(xn*charge)*phycon.sqrte);
		rate *= (A + y*c)*e1 + D*exp_int_2;

		/* Convert from excitation to de-excitation */
		/* >>>chng 99 apr 09, had not changed following, so low temp
		 * had too small a rate, and there was a discontinuity for he ii
		 * at T=30700 */
		HydColDwn_v = rate*gLo/gHi/sexpy;
	}

	HydColDwn_v = MAX2(HydColDwn_v,SMALLFLOAT);

	/* convert deexcitation rate to collision strength */
	HydColDwn_v = HydColDwn_v*gHi*phycon.sqrte/COLL_CONST;

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