/*HydroRecCool hydrogen recombination cooling, called by HydroCool */
#include "cddefines.h"
#include "physconst.h"
#include "hydrogenic.h"
#include "phycon.h"
#include "powi.h"
#include "tepowers.h"

double HydroRecCool(
	/* n is the prin quantum number on the physical scale */
	long int n , 
	/* ipZ is the charge on the C scale, 0 is hydrogen */
	long int ipZ)
{
	long int nm1; /* save n - 1*/

	double fac, 
	  hclf_v, 
	  x;
	static double  a[15]={-26.6446988,-26.9066506,-27.0619832,-27.1826903,
	  -27.2783527,-27.3595949,-27.569406,-27.611159,-27.654748,-27.70479,
	  -27.745398,-27.776126,-27.807261,-27.833093,-27.866134};
	static double  b[15]={-0.40511045,-0.41644707,-0.45834359,-0.49137946,
	  -0.51931762,-0.54971231,-0.18555807,-0.29204736,-0.36741085,
	  -0.45843009,-0.49753713,-0.51418739,-0.52287028,-0.52445456,
	  -0.52292803};
	static double  c[15]={11.29232731,11.71035693,12.89309608,13.85569087,
	  14.67354775,15.56090323,6.147461,9.0304953,11.094731,13.630431,
	  14.721959,15.172335,15.413946,15.458123,15.428761};
	static double  d[15]={.067257375,.07638384,.089925637,.102252192,
	  .111016071,.119518918,0.0093832482,0.028119606,0.039357697,0.050378417,
	  0.051644049,0.051367182,0.04938724,0.050139066,0.043085968};
	static double  e[15]={-1.99108378,-2.26898352,-2.65163846,-3.02333001,
	  -3.29462338,-3.56633674,-1.0019228,-1.5128672,-1.8327058,-2.1866371,
	  -2.2286257,-2.1932699,-2.1205891,-2.1317169,-1.9175186};
	static double  f[15]={-0.0050802618,-0.005849291,-0.0074854405,-0.0085677543,
	  -0.0093067267,-0.0098455637,0.040903604,0.037491802,0.035618861,
	  0.034132954,0.032418252,0.02947883,0.027393564,0.027607009,0.02433868};
	static double  g[15]={.166267838,.196780541,.242675042,.282237211,
	  .310204623,.335160025,-0.81087376,-0.73435108,-0.69164333,-0.64907209,
	  -0.61216299,-0.55239109,-0.51048669,-0.51963194,-0.4504203};
	static double  h[15]={.00020528663,.00027588492,.00033980652,.00041445515,
	  .00046423276,.0005121808,-0.0011986559,-0.0011333973,-0.0010992935,
	  -0.0010878727,-0.0010412678,-0.00095539899,-0.00089141547,-0.00089294364,
	  -0.00079179756};
	static double  i[15]={-0.0071357493,-0.0099630604,-0.01178647,-0.014696455,
	  -0.01670318,-0.01887373,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.};

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

	/* confirm that n is 1 or greater. */
	assert( n > 0 );

	/* this is log of (temperature divided by charge squared) since sqlogz is log10 Z^2 */
	x = TePowers.telogn[0] - TePowers.sqlogz[ipZ];

	/* this will be called for very silly low temperatures, due to high charge.
	 * evaluate at lowest fitted temp, and then extrapolate using known form.
	 * this is ok since these very high charge species do not contribute much
	 * recombination cooling */

	if( n > 15 || x < 0.2 )
	{
		/* use scale factor from actual recombination rate for this level, this ion,
		 * fac accounts for decreasing efficiency high up
		 * fac is 0.38 at n=15, and 0.32 at 25 */
		fac = 0.219 + (2.58 - 2.586/POW2((double)n));
		fac = 0.35*0.69* pow( (15./(double)n ) , 0.67 );
		hclf_v = hydro.hrec[ipZ][n][ipRecRad]*phycon.te*BOLTZMANN* fac;
		
#		ifdef DEBUG_FUN
		fputs( " <->HydroRecCool()\n", debug_fp );
#		endif
		return( hclf_v );
	}

	/* low temp (smal x) caught above, but still here for historical reasons */
	if( x < 0.2 || x > 10. )
	{
		fprintf( ioQQQ, " HydroRecCool called with invalid temperature=%e ipZ=%li\n", 
		  phycon.te , ipZ );
		puts( "[Stop in HydroRecCool]" );
		exit(1);
	}

	/* convert onto c array for n*/
	nm1 = n - 1;

	if( ipZ == 0 )
	{
		/* simple case, most often called, for hydrogen itself */
		fac = (a[nm1] + 
		  c[nm1]*TePowers.telogn[0] + 
		  e[nm1]*TePowers.telogn[1] + 
		  g[nm1]*TePowers.telogn[2] + 
		  i[nm1]*TePowers.telogn[3])/
		  (1. + b[nm1]*TePowers.telogn[0] + 
		  d[nm1]*TePowers.telogn[1] + 
		  f[nm1]*TePowers.telogn[2] + 
		  h[nm1]*TePowers.telogn[3]);
	}
	else
	{
		/* hydrogenic ions, expand as powers in t-z2 */
		fac = (a[nm1] + 
		  c[nm1]*x + 
		  e[nm1]*POW2( x) + 
		  g[nm1]*POW3( x) + 
		  i[nm1]*powi( x,4) ) /
		  (1. + b[nm1]*x  + 
		  d[nm1]*POW2( x ) + 
		  f[nm1]*POW3( x ) + 
		  h[nm1]*powi( x ,4) );
	}

	hclf_v = pow(10.,fac)*POW3(ipZ+1.);

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