/* 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 */
/*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"

/*hydro_vs_deexcit compute collision strength for collisional deexcitation for hydrogen atom, 
 * from Vriens and Smeets */
double hydro_vs_deexcit(long int ihigh, 
	  long int ilow)
{
	double anp, 
	  bn, 
	  bnp, 
	  delta, 
	  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);
	xlow = (double)(ilow);
	tev = EVRYD/TE1RYD*phycon.te;
	ryd = EVRYD;
	s = fabs(p-xlow);
	enp = -1./p/p + 1./xlow/xlow;
	eni = 1./xlow/xlow;
	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./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 nhigh, 
	  long int nlow)
{
	double apn,
	  bp,
	  bpn,
	  delta,
	  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 dexcitation between any two levels of H.
	 * valid for all nhigh, nlow
	 * at end converts this excitation rate to collision strength	*/

	n = (double)(nhigh);
	xlow = (double)(nlow);
	tev = EVRYD/TE1RYD*phycon.te;
	ryd = EVRYD;
	s = fabs(n-xlow);
	epn = -1./n/n + 1./xlow/xlow;
	epi = 1./xlow/xlow;
	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./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 );
}

/*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(nHi,nLo)*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 );
}
