/*hdexct compute collision strength for collisional deexcitation for hydrogen atom, 
 * from Vriens and Smeets */
#include "powi.h"
#include "cddefines.h"
#include "tepowers.h"
#include "phycon.h"
#include "hydrooscilstr.h"
#include "hdexct.h"

double hdexct(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( "<+>hdexct()\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 = 8.61716e-5*phycon.te;
	ryd = 13.60583;
	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*TePowers.sqrte/8.629e-6;

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