/*HydColIon calculate hydrogenic ionization rates for all n, and Z*/
#include "cddefines.h"
#include "hydrogenic.h"
#include "phycon.h"
#include "tepowers.h"
#include "powi.h"
#include "sexp.h"
#include "ee1.h"

double HydColIon(
		long int n,	/* principal quantum number, > 1
						 * since only used for excited states */
		long int ipZ)/* charge, >=1 since only used for ions 
						 * ipZ = 1 is helium the least possible charge */
{
	long int i, 
	  _r;
	double H, 
	  HydColIon_v, 
	  Rnp, 
	  chim, 
	  e1, 
	  e2, 
	  e3, 
	  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( "<+>HydColIon()\n", debug_fp );
#	endif
	/*calculate hydrogenic ionization rates for all n, and Z
	 * Allen 1973, Astro. Quan. for low Te.
	 * Sampson and Zhang 1988, ApJ, 335, 516 for High Te.
	 * */

	/* data small/1e-37/ */

	/*begin  check */
	/* this routine only for ions, ipZ=0 is H, ipZ=1 he, etc */
	assert( ipZ > 0);
	/* this is quantum level, not for ground state (n=1), only
	 * n=2 or higher */
	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 = 8.617e-5*phycon.te;
	ryd = 13.60583;
	xn = (double)n;
	chim = POW2(ipZ+1.)*ryd/xn/xn;
	y = chim/tev;

	/*>>>chng 98 dec 17, following no longer the case,
	 * ee1 always returns e1 */
	/* E1(x) integral must be multiplied by e(-y) if y>1. */
	/*if( y < 1. )*/
	/*{*/
		/*e1 = ee1(y);*/
	/*}*/
	/*else*/
	/*{*/
		/*e1 = ee1(y)*hydro.hcbolt[ipZ][n];*/
	/*}*/
	e1 = ee1(y);
	e2 = hydro.hcbolt[ipZ][n] - y*e1;
	e3 = (hydro.hcbolt[ipZ][n] - y*e2)/2.;

	t1 = 1/xn*e1;
	t2 = 1./3./xn*(hydro.hcbolt[ipZ][n] - y*e3);
	t3 = 3.*H/xn/(3. - Rnp)*(y*e2 - 2.*y*e1 + hydro.hcbolt[ipZ][n]);
	t4 = 3.36*y*(e1 - e2);
	rate = 7.69415e-9*TePowers.sqrte*9.28278e-3*powi(xn/(ipZ+1),4)*g*y;
	rate *= t1 - t2 + t3 + t4;
	rate = MAX2(rate,small);

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

	rate2 = MAX2(rate2,small);

	/* Take the lowest of the two, they fit nicely together... */
	HydColIon_v = MIN2(rate,rate2);


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