/*fe2par evaluate FeII partition function */
#include "cddefines.h"
#include "fe2par.h"
#define	NFE2PR	61
void locate(double[],long,double,long*);

double fe2par(double te)
{
	long int i;
	double du, 
	  fe2par_v, 
	  slope;
	static double fe2pt[NFE2PR]={1.000e00,5.000e00,2.000e01,5.000e01,
	  1.000e02,2.000e02,3.000e02,5.000e02,7.000e02,1.000e03,1.200e03,
	  1.500e03,2.000e03,2.500e03,3.000e03,3.500e03,4.000e03,4.500e03,
	  5.000e03,5.500e03,6.000e03,6.500e03,7.000e03,7.500e03,8.000e03,
	  8.500e03,9.000e03,9.500e03,1.000e04,1.100e04,1.200e04,1.300e04,
	  1.400e04,1.500e04,1.600e04,1.700e04,1.800e04,1.900e04,2.000e04,
	  2.200e04,2.400e04,2.600e04,2.800e04,3.000e04,3.500e04,4.000e04,
	  4.500e04,5.000e04,5.500e04,6.000e04,6.500e04,7.000e04,7.500e04,
	  8.000e04,8.500e04,9.000e04,9.500e04,1.000e05,1.050e05,1.100e05,
	  1.150e05};
	static double fe2pf[NFE2PR]={1.000e01,1.000e01,1.000e01,1.000e01,
	  1.003e01,1.056e01,1.159e01,1.403e01,1.639e01,1.960e01,2.157e01,
	  2.426e01,2.816e01,3.146e01,3.431e01,3.684e01,3.915e01,4.132e01,
	  4.341e01,4.547e01,4.755e01,4.967e01,5.185e01,5.411e01,5.646e01,
	  5.891e01,6.146e01,6.412e01,6.688e01,7.273e01,7.903e01,8.576e01,
	  9.295e01,1.006e02,1.087e02,1.174e02,1.265e02,1.362e02,1.465e02,
	  1.689e02,1.938e02,2.214e02,2.517e02,2.846e02,3.781e02,4.855e02,
	  6.038e02,7.301e02,8.616e02,9.959e02,1.131e03,1.266e03,1.399e03,
	  1.530e03,1.658e03,1.782e03,1.903e03,2.020e03,2.134e03,2.243e03,
	  2.349e03};

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

	/* function to evaluate partition function for FeII
	 *
	 * Temperature (K) */

	/* Partition Fct. */

	if( te < 1. )
	{
		fe2par_v = 10.;
	}
	else if( te > fe2pt[NFE2PR-1] )
	{
		fe2par_v = fe2pf[NFE2PR-1];
	}
	else
	{
		locate(fe2pt,NFE2PR,te,&i);
		slope = (fe2pf[i] - fe2pf[i-1])/(fe2pt[i] - fe2pt[i-1]);
		du = slope*(te - fe2pt[i-1]);
		fe2par_v = fe2pf[i-1] + du;
	}

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

/*locate locate value within array  */
void locate(double xx[], 
	  long int n, 
	  double x, 
	  long int *j)
{
	long int jl, 
	  jm, 
	  ju;

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

	/* This gives the index J of the array XX(J) such that XX(J) < X < XX(J+1), 
	 * where X (in our case the temperature) is provided. 
	 * when x exactly equals xx(j), returns xx(j+1) */
	jl = 0;
	ju = n + 1;
	/*10   if(ju-jl.gt.1)then */
	while( ju - jl > 1 )
	{
		jm = (ju + jl)/2;
		if( (xx[n-1] > xx[0]) == (x > xx[jm-1]) )
		{
			jl = jm;
		}
		else
		{
			ju = jm;
		}
		/* go to 10 */
	}
	*j = jl;

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