#include "peg_pic.h"

/** Cross section of electron-neutral particle excitation collisions **/
/** "jj" is the neutral type, "energy" is the electron energy and    **/
/** "exctp" is the type of the excitation collision.		     **/ 

float sigma2(int jj, float energy, int exctp) {
	int 		i, j, k;
	float		*tmpcs, *tmpcs1, deng;
	static float 	***sext, ***sext1;
	static int 	init_flag=1;
	char  	tchar[80];

	/**********  Initialization   **********/
	if(init_flag) {
		tmpcs = (float *)malloc(1000*sizeof(float)); 
		tmpcs1= (float *)malloc(NEMAX*sizeof(float));
		sext= (float ***)malloc(Nnum*sizeof(float **));
		sext1= (float ***)malloc(Nnum*sizeof(float **));
		
		for (i=0; i<Nnum; i++) {
			sext[i]= (float **)malloc((Nex+1)*sizeof(float *));
			sext1[i]= (float **)malloc((Nex+1)*sizeof(float *));
			for (j=0; j<=Nex; j++) {
				sext[i][j]= (float *)malloc(1000*sizeof(float));
				sext1[i][j]= (float *)malloc(NEMAX*sizeof(float));
			}
		}

		for (j=0; j<Nnum; j++) {
			if (Hydex[j]) {
				/* exact cross-sections */
				for (i=0; i<1000; i++) {
					tmpcs[i] = 0.;
				}
				
				for (i=0; i<NEMAX; i++) {
					tmpcs1[i] = 0.;
				}
				
				for (k=0; k<=Nexel[j]; k++)	{
					dec1(tchar, &i, 0, Hydex[j], 1, j, k);

					if (!i) {
						/* reading from a table */
						tcrossec(sext[j][k], tchar, 1000, 0.02);
						tcrossec(sext1[j][k], tchar, NEMAX, 1);
					} else {
						/* reading from a file */
						for (i=0; i<1000; i++) {
							deng= 0.02*i;
							sext[j][k][i] = (*fcrs)(deng);
						}
						
						for (i=0; i<NEMAX; i++) {
							sext1[j][k][i] = (*fcrs)((float)i);
						}
					}

					if (Nexel[j]) {
						for (i=0; i<1000; i++) {
							tmpcs[i] += sext[j][k][i];
						}
						
						for (i=0; i<NEMAX; i++) {
							tmpcs1[i] += sext1[j][k][i];
						}
					}
				}

				for (k=1; k<Nexel[j]; k++) {
					for (i=0; i<1000; i++) {
						sext[j][k+1][i] += sext[j][k][i];
					}
					
					for (i=0; i<NEMAX; i++) {
						sext1[j][k+1][i] += sext1[j][k][i];
					}
				}


				if (Nexel[j]) {
					for (i=0; i<1000; i++) {
						sext[j][0][i] = tmpcs[i];
						sext[j][1][i] = sext[j][0][i]-tmpcs[i];
					}
					
					for (i=0; i<NEMAX; i++) {
						sext1[j][0][i] = tmpcs1[i];
						sext1[j][1][i] = sext1[j][0][i]-tmpcs1[i];
					}
				}
			} else {
				/* approximated cross-sections */
				for (i=0; i<1000; i++) {
					deng= 0.02*i;
					sext[j][0][i] = forig_ex(deng, j);
				}
				
				for (i=0; i<NEMAX; i++) {
					sext1[j][0][i] = forig_ex((float)i, j);
				}
			}

			if (dshort>0.) {
				for (k=0; k<=Nexel[j]; k++) {
					for (i=0; i<1000; i++) {
						sext[j][k][i] *= dshort; 
					}
					
					for (i=0; i<NEMAX; i++) {
						sext1[j][k][i] *= dshort;
					}
				}
			}
		}
		
		free(tmpcs);
		init_flag =0;
	}
	/* end of if (init_flag) */

	if (energy<20) {
		i= 50*energy;
		return(sext[jj][exctp][i]);
	} else {
		i= energy +.5;
		
		if(NEMAX <= i) {
			i= NEMAX-1;
		}
		
		return(sext1[jj][exctp][i]);
	}

}  /* END of sigma2 */

/*****************************************/
