/*Hydcs123 Hydrogenic de-excitation collision strengths bewteen levels n=1,2,3,
 * for any charge.  routine only called by hydrocollid to fill in hydroline arrays
 * with collision strengths */
#include "cddefines.h"
#include "hydrogenic.h"
#include "physconst.h"
#include "phycon.h"
#include "tepowers.h"
#include "ee1.h"
/*C6cs123 line collision rates for lower levels of hydrogenic carbon, n=1,2,3 */
double C6cs123(long int i, 
  long int j);

/*Ca20cs123 line collision rates for lower levels of hydrogenic calcium, n=1,2,3 */
double Ca20cs123(long int i, 
  long int j);


double Hydcs123(
	 /* lower principal quantum number, 1, 2, or 3, in this routine 
	  * 1 is 1s, 2 is 2s, 3 is 2p, and 4 is 3 */
	 long int ipLow, 
	 /* upper principal quantum nubmer, 2, 3, or 4 */
	 long int ipHi, 
	 /* charge, 0 for hydrogen, 1 for helium, etc */
	 long int ipZ, 
	 /* = 'e' for electron collisions, ='p' for proton */
	 long int chType)
{
	long int i, 
	  j, 
	  k, 
	  _r;
	double C, 
	  D, 
	  E, 
	  expq ,
	  Hydcs123_v, 
	  Ratehigh, 
	  Ratelow, 
	  TeUse, 
	  gLo, 
	  gHi, 
	  q, 
	  rate, 
	  slope, 
	  temp, 
	  temphigh, 
	  templow, 
	  tev, 
	  x, 
	  QuanNLo, 
	  QuanNUp, 
	  Charge, 
	  ChargeSquared, 
	  zhigh, 
	  zlow;
	static double ae[5], 
	  ap[5], 
	  be[5], 
	  bp[5], 
	  ce[5], 
	  cp[5], 
	  de[5], 
	  dp[5], 
	  ee[5], 
	  ep[5];
	static double A[2]={4.4394,0.0};
	static double B[2]={0.8949,0.8879};
	static double C0[2]={-0.6012,-0.2474};
	static double C1[2]={-3.9710,-3.7562};
	static double C2[2]={-4.2176,2.0491};
	static double D0[2]={2.930,0.0539};
	static double D1[2]={1.7990,3.4009};
	static double D2[2]={4.9347,-1.7770};
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp0[] = {-2113.113,729.0084,1055.397,854.632,
		  938.9912};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			ap[i-1] = _itmp0[_r++];
			}
		}
		{ static double _itmp1[] = {-6783.515,-377.7190,724.1936,493.1107,
		  735.7466};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			bp[i-1] = _itmp1[_r++];
			}
		}
		{ static double _itmp2[] = {-3049.719,226.2320,637.8630,388.5465,
		  554.6369};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			cp[i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {3514.5153,88.60169,-470.4055,-329.4914,
		  -450.8459};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			dp[i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {0.005251557,0.009059154,0.008725781,
		  0.009952418,0.01098687};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			ep[i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {-767.5859,-643.1189,-461.6836,-429.0543,
		  -406.5285};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			ae[i-1] = _itmp5[_r++];
			}
		}
		{ static double _itmp6[] = {-1731.9178,-1442.548,-1055.364,
		  -980.3079,-930.9266};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			be[i-1] = _itmp6[_r++];
			}
		}
		{ static double _itmp7[] = {-939.1834,-789.9569,-569.1451,-530.1974,
		  -502.0939};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			ce[i-1] = _itmp7[_r++];
			}
		}
		{ static double _itmp8[] = {927.4773,773.2008,564.3272,524.2944,
		  497.7763};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			de[i-1] = _itmp8[_r++];
			}
		}
		{ static double _itmp9[] = {-0.002528027,-0.003793665,-0.002122103,
		  -0.002234207,-0.002317720};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			ee[i-1] = _itmp9[_r++];
			}
		}
		_aini = 0;
	}

#	ifdef DEBUG_FUN
	fputs( "<+>Hydcs123()\n", debug_fp );
#	endif
	/* Hydrogenic de-excitation collision rates n=1,2,3 
	 * >>refer Callaway, J. 1983, Phys Let A, 96, 83
	 * >>refer Zygelman, B., & Dalgarno, A. 1987, Phys Rev A, 35, 4085 
	 * for 2p-2s only.
	 * The fit from Callaway is in nuclear charge for 1s - 2s,2p only.
	 * For transtions involving level 3, interpolation in Z involving
	 * the functions He2cs123,C6cs123,Ne10cs123,Ca20cs123, Fe26cs123.
	 *
	 * The fits from ZD are for 2p-2s for Z=2,6,12,16,18 only other charges are
	 * interpolated, both electron and proton rates are included,
	 * the variable chType is either 'e' or 'p'.
	 *
	 * ipLow is the lower level and runs from 1 to 3 (1s, 2s, 2p)
	 * ipHi is the upper level and runs from 2 to 4 (2s, 2p, 3) */

	/* for Callaway fit: */
	/* for Zygelman and Dalgarno: */

	/* first entry is 2p, then 2s */

	/* fit in nuclear charge Z for 2p-2s collisions in hydrogenic species
	 * equation is a+bx+cx^2ln(x)+dexp(x)+eln(x)/x^2, where x=te/Z^2 in a.u.
	 * first are the proton rates: */
	/* these are electron rates: */

	/* charge is on scale with H=1, He=2, etc */
	assert( ipZ >= 0 );
	assert( ipZ < LIMELM );

	/* these are the pointers to upper and lower levels.  1=1s, 2=2s, 3=2p, 4=3 */
	assert( ipLow > 0);
	assert( ipLow <= 3);
	assert( ipHi > 1 );
	assert( ipHi <=4 );

	/* set quantum numbers and stat. weights of the transitions: */
	if( ipHi == 4 )
	{
		/* upper is n=3 then set level, stat. weight */
		QuanNUp = 3.;
		gHi = 18.;
		/* following will be set here even though it is not used in this case,
		 * to prevent good compilers from falsing on i not set,
		 * there is assert when used to make sure it is ok */
		i = -1;
	}
	else if( ipHi == 3 )
	{
		/* upper is nl=2p then set level, stat. weight */
		QuanNUp = 2.;
		gHi = 6.;
		/* used to point within vectors defined above */
		i = 0;
	}
	else if( ipHi == 2 )
	{
		/* upper is nl=2s then set level, stat. weight */
		QuanNUp = 2.;
		gHi = 2.;
		/* used to point within vectors defined above */
		i = 1;
	}
	else
	{
		fprintf( ioQQQ, " Insane levels in Hydcs123\n" );
		puts( "[Stop in hydcs123]" );
		exit(1);
	}

	/* which lower level? */
	if( ipLow == 1 )
	{
		/* lower is n=1 then set level, stat. weight */
		QuanNLo = 1.;
		gLo = 2.;
	}
	else if( ipLow == 2 )
	{
		/* lower is nl=2s then set level, stat. weight */
		QuanNLo = 2.;
		gLo = 2.;
	}
	else if( ipLow == 3 )
	{
		QuanNLo = 2.;
		gLo = 6.;
	}
	else
	{
		fprintf( ioQQQ, " Insane levels in Hydcs123\n" );
		puts( "[Stop in hydcs123]" );
		exit(1);
	}

	if( ipZ == 0 )
	{
		/* atomic hydrogen, use special routine */
		Hydcs123_v = H1cs123(ipLow,ipHi,chType);
		
#		ifdef DEBUG_FUN
		fputs( " <->Hydcs123()\n", debug_fp );
#		endif
		return( Hydcs123_v );
	}

	/* following is for charged species, since H itself returned above */

	/* following is the physical charge */
	Charge = ipZ + 1;
	/* square of charge */
	ChargeSquared = Charge*Charge;

	if( ipLow == 2 && ipHi == 3 )
	{
		/*************************************** this is 2p-2s:
		 * series of if statements determines which entries Charge is between. */
		if( ipZ < 1 )
		{
			/* this can't happen since routine returned above when ip=1 and 
			 * special atomic hydrogen routine called */
			fprintf( ioQQQ, " insane charge given to Hydcs123\n" );
			puts( "[Stop in hydcs123]" );
			exit(1);
		}
		else if( ipZ == 1 )
		{
			zlow = 2.;
			j = 1;
			zhigh = 2.;
			k = 1;
		}
		/* Li through C */
		else if( ipZ <= 5 )
		{
			zlow = 2.;
			j = 1;
			zhigh = 6.;
			k = 2;
		}
		else if( ipZ <= 11 )
		{
			zlow = 6.;
			j = 2;
			zhigh = 12.;
			k = 3;
		}
		else if( ipZ <= 15 )
		{
			zlow = 12.;
			j = 3;
			zhigh = 16.;
			k = 4;
		}
		else if( ipZ <= 17 )
		{
			zlow = 16.;
			j = 4;
			zhigh = 18.;
			k = 5;
		}
		/*>>> following changed to else from else if, 
		 * to prevent false comment in good compilers */
		/*else if( ipZ > 18 )*/
		else 
		{
			zlow = 18.;
			j = 5;
			zhigh = 18.;
			k = 5;
		}

		/* convert Te to a.u./Z^2
		 * determine rate at the low Charge */
		x = 8.6171e-5*phycon.te/27.211396/zlow/zlow;
		TeUse = MIN2(x,0.80);
		x = MAX2(0.025,TeUse);

		/* what type of collision are we dealing with? */
		if( chType == 'e' )
		{
			/* electron collisions */
			Ratelow = ae[j-1] + be[j-1]*x + ce[j-1]*x*x*log(x) + de[j-1]*
			  exp(x) + ee[j-1]*log(x)/x/x;
		}
		else if( chType == 'p' )
		{
			Ratelow = ap[j-1] + bp[j-1]*x + cp[j-1]*x*x*log(x) + dp[j-1]*
			  exp(x) + ep[j-1]*log(x)/x/x;
		}
		else 
		{
			/* this can't happen */
			fprintf( ioQQQ, " insane collision species given to Hydcs123\n" );
			puts( "[Stop in hydcs123]" );
			exit(1);
		}

		/* determine rate at the high Charge */
		x = 8.6171e-5*phycon.te/27.211396/zhigh/zhigh;
		TeUse = MIN2(x,0.80);
		x = MAX2(0.025,TeUse);
		if( chType == 'e' )
		{
			Ratehigh = ae[k-1] + be[k-1]*x + ce[k-1]*x*x*log(x) + 
			  de[k-1]*exp(x) + ee[k-1]*log(x)/x/x;
		}
		else
		{
			Ratehigh = ap[k-1] + bp[k-1]*x + cp[k-1]*x*x*log(x) + 
			  dp[k-1]*exp(x) + ep[k-1]*log(x)/x/x;
		}
		/* linearly interpolate in charge */
		if( zlow == zhigh )
		{
			rate = Ratelow;
		}
		else
		{
			slope = (Ratehigh - Ratelow)/(zhigh - zlow);
			rate = slope*(Charge - zlow) + Ratelow;
		}
		rate = rate/ChargeSquared/Charge*1.0e-7;
		/* must convert to cs and need to know the valid temp range */
		templow = 0.025*27.211396/8.6171e-5*ChargeSquared;
		temphigh = 0.80*27.211396/8.6171e-5*ChargeSquared;
		TeUse = MIN2(phycon.te,temphigh);
		temp = MAX2(TeUse,templow);
		Hydcs123_v = rate*gHi*sqrt(temp)/8.629e-6;

	}
	else if( ipHi == 4 )
	{
		/* n = 3
		 * for the rates involving n = 3, must do something different. */
		if( ipZ < 1 )
		{
			fprintf( ioQQQ, " insane charge given to Hydcs123\n" );
			puts( "[Stop in hydcs123]" );
			exit(1);
		}
		else if( ipZ == 1 )
		{
			zlow = 2.;
			Ratelow = He2cs123(ipLow,ipHi);
			zhigh = 2.;
			Ratehigh = Ratelow;
		}
		else if( ipZ > 1 && ipZ <= 5 )
		{
			zlow = 2.;
			Ratelow = He2cs123(ipLow,ipHi);
			zhigh = 6.;
			Ratehigh = C6cs123(ipLow,ipHi);
		}
		else if( ipZ > 5 && ipZ <= 9 )
		{
			zlow = 6.;
			Ratelow = C6cs123(ipLow,ipHi);
			zhigh = 10.;
			Ratehigh = Ne10cs123(ipLow,ipHi);
		}
		else if( ipZ > 9 && ipZ <= 19 )
		{
			zlow = 10.;
			Ratelow = Ne10cs123(ipLow,ipHi);
			zhigh = 20.;
			Ratehigh = Ca20cs123(ipLow,ipHi);
		}
		else if( ipZ > 19 && ipZ <= 25 )
		{
			zlow = 20.;
			Ratelow = Ca20cs123(ipLow,ipHi);
			zhigh = 26.;
			Ratehigh = Fe26cs123(ipLow,ipHi);
		}
		/*>>>chng 98 dec 17, to else to stop comment from good compilers*/
		/*else if( ipZ > 26 )*/
		else
		{
			Charge = 26.;
			zlow = 26.;
			Ratelow = Fe26cs123(ipLow,ipHi);
			zhigh = 26.;
			Ratehigh = Ratelow;
		}

		/* linearly interpolate */
		if( zlow == zhigh )
		{
			rate = Ratelow;
		}
		else
		{
			slope = (Ratehigh - Ratelow)/(zhigh - zlow);
			rate = slope*(Charge - zlow) + Ratelow;
		}
		Hydcs123_v = rate;

	}
	else
	{
		/* this branch 1-2s, 1-2p */
		if( ipZ == 1 )
		{
			/* this brance for helium, then return */
			Hydcs123_v = He2cs123(ipLow,ipHi);
			
#			ifdef DEBUG_FUN
			fputs( " <->Hydcs123()\n", debug_fp );
#			endif
			return( Hydcs123_v );
		}

		/* electron temperature in eV */
		tev = phycon.te / EVDEGK;
		/* energy in eV for hydrogenic species and these quantum numbers */
		E = ChargeSquared*EVRYD*(1./QuanNLo/QuanNLo - 1./QuanNUp/QuanNUp);
		/* E/kT for this transion */
		q = E/tev;
		TeUse = MIN2(q,10.);
		/* q is now E/kT but between 1 and 10 */
		q = MAX2(1.,TeUse);
		expq = exp(q);

		/* i must be 0 or 1 */
		assert( i==0 || i==1 );
		C = C0[i] + C1[i]/Charge + C2[i]/ChargeSquared;
		D = D0[i] + D1[i]/Charge + D2[i]/ChargeSquared;

		/* following code changed so that ee1 always returns e1,
		 * orifinal version only returned e1 for x < 1 */
		/* use disabled e1: */
		/*if( q < 1. )*/
		/*{*/
			/*rate = (B[i-1] + D*q)*exp(-q) + (A[i-1] + C*q - D*q*q)**/
			  /*ee1(q);*/
		/*}*/
		/*else*/
		/*{*/
			/*rate = (B[i-1] + D*q) + (A[i-1] + C*q - D*q*q)*ee1(q);*/
		/*}*/
		/*rate *= 8.010e-8/2./ChargeSquared/tev*sqrt(tev);*/
		/* convert to de-excitation */
		/*if( q < 1. )*/
		/*{*/
			/*rate = rate*exp(q)*gLo/gHi;*/
		/*}*/
		/*else*/
		/*{*/
			/*rate = rate*gLo/gHi;*/
		/*}*/

		/*>>>chng 98 dec 17, ee1 always returns e1 */
		rate = (B[i] + D*q)/expq + (A[i] + C*q - D*q*q)*
			  ee1(q);
		rate *= 8.010e-8/2./ChargeSquared/tev*sqrt(tev);
		/* convert to de-excitation */
		rate *= expq*gLo/gHi;

		/* convert to cs */
		Hydcs123_v = rate*gHi*TePowers.sqrte/8.629e-6;
	}


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

/*C6cs123 line collision rates for lower levels of hydrogenic carbon, n=1,2,3 */
double C6cs123(long int i, 
  long int j)
{
	double C6cs123_v, 
	  TeUse, 
	  t, 
	  x;
	static double a[3]={-92.23774,-1631.3878,-6326.4947};
	static double b[3]={-11.93818,-218.3341,-849.8927};
	static double c[3]={0.07762914,1.50127,5.847452};
	static double d[3]={78.401154,1404.8475,5457.9291};
	static double e[3]={332.9531,5887.4263,22815.211};

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

	/* These are fits to Table 5 of
	 * >>refer Aggarwal, K.M., & Kingston, A.E. 1991, J Phys B, 24, 4583
	 * C VI collision rates for 1s-3, 2s-3, and 2p-3, 
	 * principal quantum numbers n and l.
	 *
	 * i is the lower level and runs from 1 to 3 (1s, 2s, 2p)
	 * j is the upper level and runs from 2 to 4 (2s, 2p, 3)
	 * 1s-2s,2p is not done here.
	 * check temperature: fits only good between 3.8 < log Te < 6.2
	 */
	/* arrays for fits of 3 transitions see the code below for key: */

	TeUse = MAX2(phycon.te,6310.);
	t = MIN2(TeUse,1.6e6);
	x = log10(t);

	if( i == 1 && j == 2 )
	{
		/* 1s - 2s (first entry) */
		fprintf( ioQQQ, " Carbon VI 2s-1s not done in C6cs123\n" );
		puts( "[Stop in c6cs123]" );
		exit(1);

	}
	else if( i == 1 && j == 3 )
	{
		/* 1s - 2p (second entry) */
		fprintf( ioQQQ, " Carbon VI 2p-1s not done in C6cs123\n" );
		puts( "[Stop in c6cs123]" );
		exit(1);

	}
	else if( i == 1 && j == 4 )
	{
		/* 1s - 3 (first entry) */
		C6cs123_v = a[0] + b[0]*x + c[0]*x*x*sqrt(x) + d[0]*log(x) + 
		  e[0]*log(x)/x/x;
		/* C6Rate123 = cs*8.629e-6/SQRT(te)/18.
		 * */
	}
	else if( i == 2 && j == 4 )
	{
		/* 2s - 3 (second entry)         */
		C6cs123_v = a[1] + b[1]*x + c[1]*x*x*sqrt(x) + d[1]*log(x) + 
		  e[1]*log(x)/x/x;
		/* C6Rate123 = cs*8.629e-6/SQRT(te)/18.
		 * */
	}
	else if( i == 3 && j == 4 )
	{
		/* 2p - 3s (third entry) */
		C6cs123_v = a[2] + b[2]*x + c[2]*x*x*sqrt(x) + d[2]*log(x) + 
		  e[2]*log(x)/x/x;
		/* C6Rate123 = cs*8.629e-6/SQRT(te)/18.
		 * */
	}
	else
	{

		fprintf( ioQQQ, "  insane levels for C VI n=1,2,3 !!!\n" );
		puts( "[Stop in c6cs123]" );
		exit(1);

	}

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

/*Ca20cs123 line collision rates for lower levels of hydrogenic calcium, n=1,2,3 */

double Ca20cs123(long int i, 
  long int j)
{
	double Ca20cs123_v, 
	  TeUse, 
	  t, 
	  x;
	static double a[3]={-12.5007,-187.2303,-880.18896};
	static double b[3]={-1.438749,-22.17467,-103.1259};
	static double c[3]={0.008219688,0.1318711,0.6043752};
	static double d[3]={10.116516,153.2650,717.4036};
	static double e[3]={45.905343,685.7049,3227.2836};

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

	/* 
	 * These are fits to Table 5 of
	 * >>refer Aggarwal, K.M., & Kingston, A.E. 1992, J Phys B, 25, 751
	 * Ca XX collision rates for 1s-3, 2s-3, and 2p-3, 
	 * principal quantum numbers n and l.
	 *
	 * i is the lower level and runs from 1 to 3 (1s, 2s, 2p)
	 * j is the upper level and runs from 2 to 4 (2s, 2p, 3)
	 * 1s-2s,2p is not done here.
	 * check temperature: fits only good between 5.0 < log Te < 7.2
	 */

	/* arrays for fits of 3 transitions see the code below for key: */

	TeUse = MAX2(phycon.te,1.0e5);
	t = MIN2(TeUse,1.585e7);
	x = log10(t);

	if( i == 1 && j == 2 )
	{
		/* 1s - 2s (first entry) */
		fprintf( ioQQQ, " Ca XX 2s-1s not done in Ca20cs123\n" );
		puts( "[Stop in ca20cs123]" );
		exit(1);

	}
	else if( i == 1 && j == 3 )
	{
		/* 1s - 2p (second entry) */
		fprintf( ioQQQ, " Ca XX 2p-1s not done in Ca20cs123\n" );
		puts( "[Stop in ca20cs123]" );
		exit(1);

	}
	else if( i == 1 && j == 4 )
	{
		/* 1s - 3 (first entry) */
		Ca20cs123_v = a[0] + b[0]*x + c[0]*x*x*sqrt(x) + d[0]*log(x) + 
		  e[0]*log(x)/x/x;
		/* Ca20Rate123 = cs*8.629e-6/SQRT(te)/18.
		 * */
	}
	else if( i == 2 && j == 4 )
	{
		/* 2s - 3 (second entry)         */
		Ca20cs123_v = a[1] + b[1]*x + c[1]*x*x*sqrt(x) + d[1]*log(x) + 
		  e[1]*log(x)/x/x;
		/* Ca20Rate123 = cs*8.629e-6/SQRT(te)/18.
		 * */
	}
	else if( i == 3 && j == 4 )
	{
		/* 2p - 3s (third entry) */
		Ca20cs123_v = a[2] + b[2]*x + c[2]*x*x*sqrt(x) + d[2]*log(x) + 
		  e[2]*log(x)/x/x;
		/* Ca20Rate123 = cs*8.629e-6/SQRT(te)/18.
		 * */
	}
	else
	{

		fprintf( ioQQQ, "  insane levels for Ca XX n=1,2,3 !!!\n" );
		puts( "[Stop in ca20cs123]" );
		exit(1);

	}

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