/*HydColDwn collision strength for collisional de-excitation for any levels of hydrogenic Z*/
#include "cddefines.h"
#include "phycon.h"
#include "powi.h"
#include "sexp.h"
#include "tepowers.h"
#include "hydrooscilstr.h"
#include "hdexct.h"
#include "ee1.h"
#include "hydrogenic.h"

double HydColDwn(long int ipHi, 
  long int ipLo, 
  long int iz)
{
	long int i, 
	  nHi, 
	  nLo, 
	  _r;
	double A, 
	  An, 
	  D, 
	  H, 
	  HydColDwn_v, 
	  R, 
	  c, 
	  e1, 
	  e2, 
	  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( "<+>HydColDwn()\n", debug_fp );
#	endif
	/*rate coefficient for collisional de-excitation for any levels of hydrogenic Z
	 *written by Jason Ferguson to compute the 
	 * >>refer 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.
	 * */

	/*begin sanity check */
	if( (iz < 1 || ipLo >= ipHi) || ipLo*ipHi < 0 )
	{
		fprintf( ioQQQ, "HydColDwn called with impossible parameters.\n" );
		fprintf( ioQQQ, "%5ld%5ld%5ld\n", ipHi, ipLo, iz );
		puts( "[Stop in hydrocoldwn]" );
		exit(1);
	}
	/*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 iz = 1 */
	if( iz == 1 )
	{
		HydColDwn_v = hdexct(nHi,nLo)*gLo;
		
#	ifdef DEBUG_FUN
		fputs( " <->HydColDwn()\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 = 8.61716e-5*phycon.te;
	ryd = 13.60583;

	y = POW2((double)iz)*ryd*(1./POW2(xn) - 1./POW2(xp))/tev;
	/*>>>chng 98 dec 17, redefine ee1 so always first exponential integral*/
	/*if( y < 1. )*/
	/*{*/
		/*e1 = ee1(y);*/
		/*e2 = sexp(y) - y*e1;*/
	/*}*/
	/*else*/
	/*{*/
		/*e1 = ee1(y);*/
		/*e2 = 1 - y*e1;*/
	/*}*/
	sexpy = sexp(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. )
	{
		HydColDwn_v = 0.;
	}
	else
	{
		e1 = ee1(y);
		e2 = 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 *= sexp(0.006*powi(xn-1.,6)/iz);
		}

		rate = 8.63e-6/(POW2(xn*iz)*TePowers.sqrte);
		rate *= (A + y*c)*e1 + D*e2;

		/* 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 */
		/*if( y < 1 )*/
		/*{*/
			HydColDwn_v = rate*gLo/gHi/sexpy;
		/*}*/
		/*else*/
		/*{*/
			/*HydColDwn_v = rate*gLo/gHi;*/
		/*}*/
	}

	HydColDwn_v = MAX2(HydColDwn_v,1.e-38);

	/* convert deexcitation rate to collision strength */
	HydColDwn_v = HydColDwn_v*gHi*TePowers.sqrte/8.629e-6;

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