#include "peg_pic.h"

/* calculation of the effective ionization probability */

float efsigma(float energy,float density) {
	int          i, j, n, k;
	static int   rows, cols, init_flag=1;
	static float *dens, **seff, emin, emin1, *sv_excit;
	float	     coef, dum, dum1;
	char	     nm[80];
	FILE         *fp; 

	if (init_flag) {
		/* Effective ionization rate coefficiets are in efi_rc.txt */
		strcpy(nm, dargv3);
		strcat(nm, "efi_rc.txt"); 
		fp=fopen(nm,"r");

		if (!fp) {
			puts("\nCan't find input file for ef. Ionization (tableEfI)");
			exit(1);
		}
		
		fscanf(fp,"%d %d %f\n", &rows, &cols, &emin);
		dens=(float*)malloc(cols*sizeof(float));
		
		for(i=0;i<cols;i++) {
			fscanf(fp,"%g",dens+i);
		}

		seff=(float**)malloc(cols*sizeof(float*));

		for(i=0;i<cols;i++) {
			seff[i]=(float*)malloc(rows*sizeof(float));
		}

		sv_excit=(float*)malloc(rows*sizeof(float));

		for(i=0;i<rows;i++) {
			for(j=0;j<cols;j++) {
				fscanf(fp,"%g",&seff[j][i]);
				if (dshort1) seff[j][i] *=dshort;
			}
		}
		
		fclose(fp);

		for(i=0;i<cols;i++) {
			dens[i] /= (nc2p/area/dx);
		}

		emin1 = 9*emin+0.5; /* N0=10*emin + 0.5 (-emin), 10=1/delta_E */

		/**********************************************************************/
		/**        Calculation of <sigma_excit.xV>, (averaged over           **/
		/**                        the maxwell-distribution)                 **/

		coef= 6.6924e4;    /* 0.1*sqrt(1.6022e-19*2/pi/m_e), 0.1=del_E*/ 
		j= emin*10 + 1;
		
		for(i=0;i<rows;i++) { 
			sv_excit[i] =0.;
		}

		for(i=1;i<rows;i++) {
			if (i<j) {
				dum= 0.1*i;                 /* delta_E = 0.1 eV */
			} else {
				dum= emin + (i-j) + 1 ;     /* delta_E = 1 eV */
			}

			k= 10*NEMAX;
			
			for(n=1;n<k;n++) {                    /* Maxwell averaging */
				dum1= 0.1*n;
				sv_excit[i] +=(sigma2(0, dum1, 0)*dum1*exp(-dum1/dum));
			}

			sv_excit[i] *= (coef/dum/sqrt(dum));
		}

		cols--;
		rows--;
		init_flag=0;
	}

	/*  printf("s= %e\n", seff[0][200]); */

	if (density-dens[0]<=0) {
		return (0.);
	} else if (density-dens[cols]>0) {
		printf("Too high density! Density=%g",density);
		exit(1);
	} else {
		if (energy<emin) {
			i= energy*10+0.5;
		} else { 
			i= energy+emin1;
			i=max(i,rows);
		} 

		j=1;
		
		while(density-dens[j]>0) {
			j++;
		}

		return ( ((seff[j][i]*(density-dens[j-1])+
				seff[j-1][i]*(dens[j]-density) ) /
				(dens[j]-dens[j-1]) - seff[0][i]) /sv_excit[i] );
	}

	return (0.); /* just to avoid "warning" */

} /* End of EFSIGMA */
