#include <string.h>
#include "cddefines.h"  
#include "hsrate.h"  

double HSRate( 
	  /* upper and lower quantum numbers*/
	  long int iHi, 
	  long int iLo, 
	  /* charge, 1 or 2 for now */
	  long int ipZ, 
	  /* temperature to interpolate, for H between 500-30,000K*/
	  double TempIn, 
	  /* density to interpolate */
	  double DenIn)

/* general utility to interpolate line emissivities from the Storey & Hummer tables
 of case B emissivities.  
 iHi, iLo, two long integers are the principal quantum
 numbers, and are upper and lower levels in any order.  
 ipZ nuclear charge, 1 or 2 have data now
 TempIn = temperature, and must lie within the range of the table, which depends on
 the ion charge, and is 500 - 30,000K for hydrogen.  
 DenIn is the density and must not exceed the high density limit to the table.  

 routine returns -1 if conditions outside temperature range, or
 above density range of tabulated case b results
 If desired density is low limit, lower limit is used instead 
*/

{
	long int 
		ipTemp, /*pointer to temperature*/
		ipDens, /*pointer to density*/
		ipDensHi ,
		ipTempHi ;
	int ipUp , ipLo , ip ;
	double x1 , x2 , yy1 , yy2 , z1 , z2 , za , zb ,z;

	/*make sure ipZ is 1 or 2*/
	if( ipZ<1 || ipZ> 2 )
	{
		printf("HSRate called with improper ipZ, was%li and must be 1 or 2",ipZ);
		exit(-67);
	}
	/* now back ipZ back one, to be on c scale */
	--ipZ;

	/*===========================================================*/
	/* following is option to have n in either order, result is that ipUp and ipLo
	 * will be the upper and lower levels */
	if( iHi > iLo )
	{
		ipUp = iHi;  ipLo = iLo ;
	}
	else if( iHi < iLo )
	{
		ipUp = iLo ; ipLo = iHi ;
	}
	else
	{ 
		printf("HSRate called with indices equal, %ld  %ld  \n",iHi,iLo);
		exit(1);
	}

	/* now check that they are in range of the model atom*/
	if( ipLo <2 ) 
	{
		printf(" HSRate called with lower quantum number less than 2, = %i\n",
			ipLo);
		exit(1);
	}

	if( ipUp >25 ) 
	{
		printf(" HSRate called with upper quantum number greater than 25, = %i\n",
			ipUp);
		exit(1);
	}

	/*===========================================================*/
	/*bail if above high density limit */
	if( DenIn > CaseBHS.Density[ipZ][CaseBHS.nDensity[ipZ]-1] ) 
	{
		/* this is flag saying bogus results */
		return -1;
	}

	if( DenIn <= CaseBHS.Density[ipZ][0] )
	{
		/* this case, desired density is below lower limit to table,
		 * just use the lower limit */
		ipDens = 0;
	}
	else
	{
		/* this case find where within table density lies */
		for( ipDens=0; ipDens < CaseBHS.nDensity[ipZ]-1; ipDens++ )
		{
			if( DenIn >= CaseBHS.Density[ipZ][ipDens] && 
				DenIn < CaseBHS.Density[ipZ][ipDens+1] ) break ;
		}
	}


	/*===========================================================*/
	/* confirm within temperature range*/
	if( TempIn < CaseBHS.ElecTemp[ipZ][0] || 
		TempIn > CaseBHS.ElecTemp[ipZ][CaseBHS.ntemp[ipZ]-1] ) 
	{
		/* this is flag saying bogus results */
		return -1;
	}

	/* find where within grid this temperature lies */ 
	for( ipTemp=0; ipTemp < CaseBHS.ntemp[ipZ]-1; ipTemp++ )
	{
		if( TempIn >= CaseBHS.ElecTemp[ipZ][ipTemp] && 
			TempIn < CaseBHS.ElecTemp[ipZ][ipTemp+1] ) break ;
	}

	/*===========================================================*/
	/*we now have the pointers within the temperature array*/

	if( ipDens+1 > CaseBHS.nDensity[ipZ]-1 )
	{ 
		/* special case, when cell is highest density point */
		ipDensHi = CaseBHS.nDensity[ipZ]-1;
	}
	else if( DenIn < CaseBHS.Density[ipZ][0])
	{
		 /* or density below lower limit to table, set both bounds to 0 */
		ipDensHi = 0;
	}
	else
	{ 
		ipDensHi = ipDens+1;
	}

	/*special case, if cell is highest temperature point*/
	if( ipTemp+1 > CaseBHS.ntemp[ipZ]-1 )
	{ 
		ipTempHi = CaseBHS.ntemp[ipZ]-1;
	}
	else
	{ 
		ipTempHi = ipTemp+1;
	}

	x1 = log10( CaseBHS.Density[ipZ][ipDens] );
	x2 = log10( CaseBHS.Density[ipZ][ipDensHi] );

	yy1 = log10( CaseBHS.ElecTemp[ipZ][ipTemp] );
	yy2 = log10( CaseBHS.ElecTemp[ipZ][ipTempHi] );

	/*now generate the pointer to the array, expression from storey code -1 for c*/
	ip = (((CaseBHS.ncut[ipZ]-ipUp)*(CaseBHS.ncut[ipZ]+ipUp-1))/2)+ipLo - 1;

	/*pointer must lie within line array*/
	assert( ip<NLINEHS ); 
	assert( ip >= 0 );

	/* interpolate on emission rate*/
	z1 = log10( CaseBHS.Emiss[ipZ][ipTemp][ipDens][ip]) ;
	z2 = log10( CaseBHS.Emiss[ipZ][ipTemp][ipDensHi][ip]) ;

	if(x2 == x1 ) 
	{
		za = z2;
	}
	else 
	{
		za = z1 + log10( DenIn / CaseBHS.Density[ipZ][ipDens] ) * (z2-z1)/(x2-x1);
	}

	z1 = log10( CaseBHS.Emiss[ipZ][ipTempHi][ipDens][ip]) ;
	z2 = log10( CaseBHS.Emiss[ipZ][ipTempHi][ipDensHi][ip]) ;

	if(x2 == x1 ) 
	{
		zb = z2;
	}
	else 
	{
		zb = z1 + log10( DenIn / CaseBHS.Density[ipZ][ipDens] ) * (z2-z1)/(x2-x1);
	}

	if( yy2 == yy1 ) 
	{
		z = zb ;
	}
	else 
	{
		z = za + log10( TempIn / CaseBHS.ElecTemp[ipZ][ipTemp] ) * (zb-za)/(yy2-yy1);
	}

	return ( pow(10.,z) );

}