/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*Badnell_rec_init This code is written by Terry Yun, 2005 *
 * It reads rate coefficient fits into 3D arrays and output array.out for testing *
 * The testing can be commented out */
/*Badnell_DR_rate_eval This code is written by Terry Yun, 2005 *
 * It interpolates the rate coefficients in a given temperature.*
   It recieves ATOMIC_NUM_BIG, NELECTRONS values, temperature and returns the rate coefficient*
   It returns
        '-2': initial <= final
              init < 0 or init >302 or final < 0 or final > 302
        '-1': the transition is not defined
        '99': unknown invalid entries                         */ 
#include "cddefines.h"
#include "phycon.h"
#include "path.h"
#include "ionbal.h"

#define MAX_FIT_PAR_DR 9
/* one larger than largest atomic number in Nigel's data file */
#define ATOMIC_NUM_BIG 55
/* number of electrons */
#define NELECTRONS 12
static double C[ATOMIC_NUM_BIG][NELECTRONS][MAX_FIT_PAR_DR]; 
static double E[ATOMIC_NUM_BIG][NELECTRONS][MAX_FIT_PAR_DR];    
static int defn_C[ATOMIC_NUM_BIG][NELECTRONS];
static int defn_E[ATOMIC_NUM_BIG][NELECTRONS];
static int nn_counter[ATOMIC_NUM_BIG][NELECTRONS];

#define MAX_FIT_PAR_RR 6
static double par[ATOMIC_NUM_BIG][NELECTRONS][MAX_FIT_PAR_RR]; 
static int defn[ATOMIC_NUM_BIG][NELECTRONS];

/* #define PRINT_DR */
/* #define PRINT_RR */

#ifdef	PRINT_DR
#define FILE_NAME_OUT "array.out"
#endif

#ifdef	PRINT_RR
#define FILE_NAME_OUT "array.out"
#endif

/*Badnell_rec_init This code is written by Terry Yun, 2005 *
 * It reads rate coefficient fits into 3D arrays and output array.out for testing *
 * The testing can be commented out */
void Badnell_rec_init( void )
{

	double par_C[MAX_FIT_PAR_DR];
	double par_E[MAX_FIT_PAR_DR];
	char chLine[INPUT_LINE_LENGTH];
	int z_value, n_value;
	int i, j, k;
	int count, number;
	double temp_par[MAX_FIT_PAR_RR];
	int M_state, W_state;

	/*char chFilename[FILENAME_PATH_LENGTH_2] ;*/
#	define NBLOCK	2
	int data_begin_line[NBLOCK];/*it tells you where the data set begins(begins with 'Z')*/
	int length_of_line;	/*this variable checks for a blank line*/
	FILE *ioDATA;
	char chFilename[FILENAME_PATH_LENGTH_2] ;
	int yr, mo, dy;
	char *chs;

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

	#ifdef PRINT_DR
	FILE *ofp;
	ofp= fopen(FILE_NAME_OUT, "w");
	if( ofp == NULL )
	{
		printf("Can't open %s\n", FILE_NAME_OUT);
		exit(1);
	}
	#endif

	#ifdef PRINT_RR
	FILE *ofp;
	ofp= fopen(FILE_NAME_OUT, "w");
	if( ofp == NULL )
	{
		printf("Can't open %s\n", FILE_NAME_OUT);
		exit(1);
	}
	#endif

	for( i=0; i<NBLOCK; ++i )
	{
		/* set to really large negative number - crash if used before being redefined */
		data_begin_line[i] = INT_MIN;
	}

	/* check on path is file not here and path set */
	/* path was parsed in getset */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "badnell_dr.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "badnell_dr.dat" );
	}

	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		printf("Badnell_rec_init cannot open %s\n", chFilename);
		exit(1);
	}

	count = 0;
	number = 0;

	/*Find out the number line where the data starts 
	 * there are two main blocks of data and each starts with a Z in column 2 */
	while( fgets(chLine, (int)sizeof(chLine), ioDATA) != NULL )
	{
		count++;

		if( chLine[2]=='Z' )
		{
			/* number has to be 0 or 1, and indicates the first or second block of data
			 * count is the line number for the start of that block */
			data_begin_line[number] = count;
			ASSERT( number < NBLOCK );
			number++;
		}
	}

	/*set a flag for a undefined data*/
	for( i=0; i<ATOMIC_NUM_BIG; i++ )
	{
		for( j=0; j<NELECTRONS; j++ )
		{
			defn_C[i][j] = 0;
			defn_E[i][j] = 0;   
		}
	}

	/*set fitting coeffs to zero initially*/
	for( i=0; i<ATOMIC_NUM_BIG; i++ )
	{
		for( j=0; j<NELECTRONS; j++ )
		{
			for( k=0; k<MAX_FIT_PAR_DR; k++ )
			{
				C[i][j][k] = 0;
				E[i][j][k] = 0;  
			}
		}
	}

	count = 0;
	/*Start from beginning to read in again*/  
	fseek(ioDATA, 0, SEEK_SET);
	/* read magic number for DR data */
	fgets(chLine, (int)sizeof(chLine), ioDATA);
	count++;
	{
		int dr_yr = 2006, dr_mo = 1, dr_dy = 27;
		chs = strchr(chLine, ')');
		++chs;
		sscanf(chs, "%4i%2i%2i",&yr, &mo, &dy);
		if((yr != dr_yr) || (mo != dr_mo) || (dy != dr_dy))
		{
			printf("The version of %s I found (%i %i %i) is not the current version (%i %i %i).\n", chFilename, yr, mo, dy, dr_yr, dr_mo, dr_dy);
			cdEXIT(EXIT_FAILURE);
		}
	}

	while( fgets(chLine, (int)sizeof(chLine), ioDATA) != NULL )
	{
		count++;
		length_of_line = (int)strlen(chLine);

		/*reading in coeffs C */
		if( count > data_begin_line[0] && count < data_begin_line[1] && length_of_line >3 )
		{
			/*set array par_C to zero */
			for( i=0; i<MAX_FIT_PAR_DR; i++ )
			{
				par_C[i] = 0;
			}
			sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf%lf%lf%lf", &z_value, &n_value, &M_state, &W_state, &par_C[0], &par_C[1], &par_C[2], &par_C[3], &par_C[4], &par_C[5], &par_C[6], &par_C[7], &par_C[8]);

			if(M_state == 1)
			{
				/*Set a flag to '1' when the indices are defined */
				defn_C[z_value][n_value] = 1;

				/*counting the number of coeffs */
				nn_counter[z_value][n_value] = 8;   
				while( par_C[nn_counter[z_value][n_value]] == 0 )
				{
					--nn_counter[z_value][n_value];
				}
				++nn_counter[z_value][n_value];

			/*assign the values into array */
				for( i=0; i<nn_counter[z_value][n_value]; i++ )
					C[z_value][n_value][i] = par_C[i];
			}
		}
	}

	/*starting again to read in E values */
	fseek(ioDATA, 0, SEEK_SET);
	count = 0; 
	while( fgets(chLine, (int)sizeof(chLine), ioDATA) != NULL )
	{
		count++;
		length_of_line = (int)strlen(chLine);
		if( count > data_begin_line[1] && length_of_line > 3 )
		{

		  /*set array par_E to zero*/
			for( i=0; i<MAX_FIT_PAR_DR; i++ )
			{
				par_E[i] = 0;
			}   
			sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf%lf%lf%lf", &z_value, &n_value, &M_state, &W_state, &par_E[0], &par_E[1], &par_E[2], &par_E[3], &par_E[4], &par_E[5], &par_E[6], &par_E[7], &par_E[8]);

			if(M_state == 1)
			{
				/*counting the number of coeffs*/
				nn_counter[z_value][n_value] = 8;
				while( par_E[nn_counter[z_value][n_value]] == 0 )
				{
					--nn_counter[z_value][n_value];
				}
				++nn_counter[z_value][n_value];

				/*Set a flag to '1' when the indices are defined*/
				defn_E[z_value][n_value] = 1;

				/*assign the values into array*/
				for( i=0; i<nn_counter[z_value][n_value]; i++ )
					E[z_value][n_value][i] = par_E[i];
			}
		}
	}

	/*output coeffs for defined values for testing */
	#ifdef PRINT_DR
	for( i=0; i<ATOMIC_NUM_BIG; i++ )
	{
		for( j=0; j<NELECTRONS; j++ )
		{
			if( defn_C[i][j] == 1 )
			{
				fprintf(ofp, "%i %i %e %e %e %e %e %e %e %e %e\n", i, j, C[i][j][0], C[i][j][1], C[i][j][2], C[i][j][3], C[i][j][4], C[i][j][5], C[i][j][6], C[i][j][7], C[i][j][8]);
			}
		}
	}
	for( i=0; i<ATOMIC_NUM_BIG; i++ )
	{
		for( j=0; j<NELECTRONS; j++ )
		{
			if( defn_E[i][j] == 1 )
			{
				fprintf(ofp, "%i %i %e %e %e %e %e %e %e %e %e\n", i, j, E[i][j][0], E[i][j][1], E[i][j][2], E[i][j][3], E[i][j][4], E[i][j][5], E[i][j][6], E[i][j][7], E[i][j][8]);
			}
		}
	}
	fclose(ofp);
   #endif
   
	/*checking for the match of defn_C and defn_E - Both have to be defined*/
	for( i=0; i<ATOMIC_NUM_BIG; i++ )
	{
		for( j=0; j<NELECTRONS; j++ )
		{
			if( defn_C[i][j] != defn_E[i][j] )
			{
				printf("C%i, E%i: %i %i\n", i, j, defn_C[i][j], defn_E[i][j]);
				printf("Critical error in defn\n");
			}
		}
	}

	/* fclose(ioDATA); */

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

	/***************************/
	/* addition: code for RR(Radiative recombination) */

	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "badnell_rr.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "badnell_rr.dat" );
	}

	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		printf("Badnell_rec_init cannot open %s\n", chFilename);
		exit(1);
	}

	/*set a flag for a undefined data*/
	for( i=0; i<ATOMIC_NUM_BIG; i++ )
	{
		for( j=0; j<NELECTRONS; j++ )
		{
			defn[i][j] = 0;
		}
	}

	/*set fitting coeffs to zero initially*/
	for( i=0; i<ATOMIC_NUM_BIG; i++ )
	{
		for( j=0; j<NELECTRONS; j++ )
		{
			for( k=0; k<MAX_FIT_PAR_RR; k++ )
			{
				par[i][j][k] = 0;
			}
		}
	}

	/* read magic number for RR data */
	{
		int rr_yr = 2005, rr_mo = 12, rr_dy = 27;
		fgets(chLine, (int)sizeof(chLine), ioDATA);
		chs = strchr(chLine, ')');
		++chs;
		sscanf(chs, "%4i%2i%2i", &yr, &mo, &dy);
		if((yr != rr_yr)||(mo != rr_mo)||(dy != rr_dy))
		{
			printf("The version of %s I found (%i %i %i) is not the current version (%i %i %i).\n", chFilename, yr, mo, dy, rr_yr, rr_mo, rr_dy);
			cdEXIT(EXIT_FAILURE);
		}
	}

	while( fgets(chLine, (int)sizeof(chLine), ioDATA) != NULL )
	{
		/*reading in coeffs */
		/*set array par to zero */
		for( i=0; i<MAX_FIT_PAR_RR; i++ )
		{
			temp_par[i] = 0;
		}
		if(chLine[0] != '#')
		{
			sscanf(chLine, "%i%i%i%i%lf%lf%lf%lf%lf%lf", &z_value, &n_value, &M_state, &W_state, &temp_par[0], &temp_par[1], &temp_par[2], &temp_par[3], &temp_par[4], &temp_par[5]);

			if(M_state == 1)
			{
				/*Set a flag to '1' when the indices are defined */  
				defn[z_value][n_value] = 1;
				/*assign the values into array */
				for( i=0; i<MAX_FIT_PAR_RR; i++ )
					par[z_value][n_value][i] = temp_par[i];
			}
		}
	}

	/*output coeffs for defined values for testing */
	#ifdef PRINT_RR
	count = 0;
	for( i=0; i<ATOMIC_NUM_BIG; i++ )
	{
		for( j=0; j<NELECTRONS; j++ )
		{
			if( defn[i][j] == 1 )
			{
				fprintf(ofp, "%i %i %e %e %e %e %e %e\n", i, j, par[i][j][0], par[i][j][1], par[i][j][2], par[i][j][3], par[i][j][4], par[i][j][5]);
				count++;
			}
		}
	}
	fprintf(ofp, "total lines are %i ", count);

	fclose(ofp);
	#endif

	fclose(ioDATA);

	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			long int nelem;
			for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
			{
				int ion;
				fprintf(ioQQQ,"\nDEBUG rr rec\t%li",nelem);
				for( ion=0; ion<=nelem; ++ion )
				{
					fprintf(ioQQQ,"\t%.2e", Badnell_RR_rate_eval(nelem+1 , nelem-ion ) );
				}
				fprintf(ioQQQ,"\n");
				fprintf(ioQQQ,"DEBUG dr rec\t%li",nelem);
				for( ion=0; ion<=nelem; ++ion )
				{
					fprintf(ioQQQ,"\t%.2e", Badnell_DR_rate_eval(nelem+1 , nelem-ion ) );
				}
				fprintf(ioQQQ,"\n");
			}
			cdEXIT(EXIT_FAILURE);
		}
	}

	return;
}

/*Badnell_DR_rate_eval This code is written by Terry Yun, 2005 *
 * It interpolates the rate coefficients in a given temperature.*
   It recieves ATOMIC_NUM_BIG, NELECTRONS values, temperature and returns the rate coefficient*
   It returns
        '-2': initial <= final
              init < 0 or init >302 or final < 0 or final > 302
        '-1': the transition is not defined
        '99': unknown invalid entries                         */ 
double Badnell_DR_rate_eval(
	/* atomic number on physics scale - He - 2 */
	int z_val, 
	/* number of core electrons before capture of free electron */
	int n_val )
{

	double k, sum;
	int i;

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

	/*Invalid entries returns '-2':if the z_values are smaller than equal to the n_values */
	if( z_val <= n_val )
	{
		k = -2;
	}
	/*Invalid intries returns '-2' if z_val and n_val are out of the range*/
	else if( z_val <2 || z_val > 30 || n_val < 1 || n_val > 11 )
	{
		k = -2;
	}
	/*undefined z and n returns '-1'*/
	else if( defn_C[z_val][n_val] == 0 )
	{
		k = -1;
	}
	else if( defn_C[z_val][n_val] == 1 )
	{
		sum = 0;
		i = 0;
		/* loop over all non-zero terms */
		for(i=0; i<nn_counter[z_val][n_val]; ++i  )
		{
			sum += (C[z_val][n_val][i] * sexp( E[z_val][n_val][i]/phycon.te));
		}

		/*k = pow(phycon.te, -1.5) * sum;*/
		k = sum / phycon.te32;
	}
	/*unkown invalid entries returns '-99'*/
	else
	{
		k = -99;
	}

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

	return k;

}

double Badnell_RR_rate_eval(
	/* atomic number on physics scale - He - 2 */
	int z_val, 
	/* number of core electrons before capture of free electron */
	int n_val )
{
	double k;
	double B, D, F;

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

	/*Invalid entries returns '-2':if the z_values are smaller than equal to the n_values */
	if( z_val <= n_val )
	{
		k = -2;
	}
	/*Invalid intries returns '-2' if z_val and n_val are out of the range*/
	else if( z_val <1 || z_val > 30 || n_val < 0 || n_val > 11 )
	{
		k = -2;
	}
	/*undefined z and n returns '-1'*/
	else if( defn[z_val][n_val] == 0 )
	{
		k = -1;
	}
	/* coeffs:A=par[0], B=par[1], T0=par[2], T1=par[3], C=par[4], T2=par[5] */
	else if( defn[z_val][n_val] == 1 )
	{

		/* k=A*[(T/T0)^1/2*(1+(T/T0)^1/2)^1-B*(1+(T/T1)^1/2)^1+B]^-1 
		where B = B + C*exp(-T2/T) */
		double temp;

		temp = -par[z_val][n_val][5]/phycon.te; /* temp = (-T2/T) */
		B = par[z_val][n_val][1] + par[z_val][n_val][4]*exp(temp);
		D = sqrt(phycon.te/par[z_val][n_val][2]); /* D = (T/T0)^1/2 */
		F = sqrt(phycon.te/par[z_val][n_val][3]); /* F = (T/T1)^1/2 */
		k = par[z_val][n_val][0]/(D*pow((1.+D),(1.-B))*pow((1.+F),(1.+B)));

	}

	/*unkown invalid entries returns '-99'*/
	else
	{
		k = -99;
	}

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

	return k;
}
