#include "peg_pic.h"

/************************************************************************/
/********* Some functions used in different routines ********************/


/**************** Particle energy ****************************/

float Ep(int ii, int jj) {
	return(vx[ii][jj]*vx[ii][jj] + vy[ii][jj]*vy[ii][jj] + vz[ii][jj]*vz[ii][jj]);
}

/*******************************************************************/
/** Decide which data table (tchr) or which function (datp()) ******/
/** has to be used for charged-neutral particle collision **********/
/** cross-sections. "ctp" is the collision tyle, elastic, ...; *****/
/** "dattp" is the flag for data type: if=0, cross-sections are ****/
/** given as a table, else they are given as a function. "chtp" ****/
/** is the type of the collided charged particle and ntp is either */
/** the type of the collided neutral particles if charged particles */
/** are the electrons, or the collision type for the ion-neutral   */
/** collisions (H+H, H+D, and so on). "ctp" is the collision number*/
/** corresponding to this process. "extp" is used for the excitation*/
/** collision type. */ 

void dec1(char *tchr, int *dattp, int chtp, int ntp, int ctp, int cpt, int extp) {

	if (chtp==0) {
		/* the charged particles are electrons */
		if (ntp<=3) {
			/* neutrals are H, D, or T */
			if (!ctp) {
				/* elastic collisions */
				*dattp = 0;
				strcpy(tchr, "eHel.txt");
			} else if (ctp==1) {
				/* excitation collisions */
				*dattp = 0;
				strcpy(tchr, "eHex.txt");
			} else {
				/* ionization collision */
				*dattp= 1;
				fcrs = feh_i;
			}
		} else if (ntp==4) {
			/* neutrals are O2 mol. */
			if (!ctp) {
				/* elastic collisions */
				*dattp = 0;
				strcpy(tchr, "eO2_mt_PIC.txt");
			} else if (ctp==1)  {
				/* excitation collisions */
				*dattp = 0;

				if (extp==0) {
					strcpy(tchr, "eO2_v1_PIC.txt");   // v=1, excit.
				} else if (extp==1) {
					strcpy(tchr, "eO2_dis_PD_PIC.txt"); // 3P 1D dissociat.
				} else if (extp==2) {
					strcpy(tchr, "eO2_v2_PIC.txt"); // v=2, excit.
				} else if (extp==3) {
					strcpy(tchr, "eO2_dis_PP_PIC.txt"); // 3P 3P dissociat.
				} else if (extp==4) {
					strcpy(tchr, "eO2_dis_DD_PIC.txt"); // 1D 1D dissociat.
				} else if (extp==5) {
					strcpy(tchr, "eO2_csigm_PIC.txt");// Metastable excitation (c1, A3) 
				} else if (extp==6) {
					strcpy(tchr, "eO2_sing_PIC.txt"); // Metastable excitation (a1)
				} else if (extp==7) {
					strcpy(tchr, "eO2_dis_3p3P_PIC.txt"); // ? 3p^3P dissociative excit.
				} else if (extp==8) {
					strcpy(tchr, "eO2_bsigm_PIC.txt"); // Metastable excitation (b1)
				} else if (extp==9) {
					strcpy(tchr, "eO2_rot_PIC.txt"); // Rotational excitation
				} else if (extp==10) {
					strcpy(tchr, "eO2_dis_3p5P_PIC.txt"); // 3P 3p^5P dissociative excit.
				}
			} else {
				/* ionization collision */
				*dattp= 0;
				strcpy(tchr, "eO2_i_PIC.txt");
			}
		} else {
			/* new e+N collisions have to be added here */
			printf("Data for electron-neutral (= %d) particle collisions is missing\n", ntp);
			exit(1);
		}
		/* end of if (chtp==0) loop */
	} else {
		/* the charged particles are ions */
		if (ntp==1) {
			/* HH collisions */
			*dattp = 0;

			if (!ctp) {
				/* elastic collisions */
				strcpy(tchr, "HHel.txt");
			} else {
				/* charge exchange collisions */
				strcpy(tchr, "HHcx.txt");
			}
		} else if (ntp==2) {
			/* DD collisions */
			*dattp = 0;

			if (!ctp) {
				/* elastic collisions */
				strcpy(tchr, "DDel.txt");
			} else  {
				/* charge exchange collisions */
				strcpy(tchr, "DDcx.txt");
			}
		} else if (ntp==3) {
			/* TT collisions */
			*dattp = 0;

			if (!ctp) {
				/* elastic collisions */
				strcpy(tchr, "TTel.txt");
			} else {
				/* charge exchange collisions */
				strcpy(tchr, "TTcx.txt");
			}
		} else if (ntp==4) {
			/* HD collisions */
			*dattp = 0;

			if (!ctp) {
				/* elastic collisions */
				strcpy(tchr, "HDel.txt");
			} else {
				/* charge exchange collisions */
				strcpy(tchr, "HDcx.txt");
			}
		} else if (ntp==5) {
			/* HT collisions */
			*dattp = 0;

			if (!ctp)  {
				/* elastic collisions */
				strcpy(tchr, "HTel.txt");
			} else {
				/* charge exchange collisions */
				strcpy(tchr, "HTcx.txt");
			}
		} else if (ntp==6) {
			/* DT collisions */
			*dattp = 0;

			if (!ctp) {
				/* elastic collisions */
				strcpy(tchr, "DTel.txt");
			} else {
				/* charge exchange collisions */
				strcpy(tchr, "DTcx.txt");
			}
		} else if (ntp==7) {
			/* O2+O2 collisions */
			*dattp = 0;

			if (!ctp) {
				/* elastic collisions */
				strcpy(tchr, "zero.txt");
			} else {
				/* charge exchange collisions */
				strcpy(tchr, "O2_O2+_CX_new.txt");
			}
		} else {
			/* new i+N collisions have to be added here */
			printf("Data for ion-neutral particle collisionis is missing\n");
			exit(1);
		}
	}   /* end of if (chtp!=0) loop */
}
/* End of dec1() */


/*******************************************************************/
/** Calculation of matrix for cross-sections using tabulated data **/
/** Function reads the file "atdata/crsarr" and criates ************/
/** cross-section table tarr[i] (i=0,..., imax). "de" is the *******/
/** energy grid **/

void tcrossec(float *tarr, char *crsarr, int imax, float de) {
	int    i, j, jmax;
	float  en, *enarr, *crsarr1;
	char   nm[80];
	FILE   *fp;

	strcpy(nm, dargv3);
	strcat(nm, crsarr);

	fp=fopen(nm,"r");

	if (!fp) {
		puts("\nCan't find input file with cross-sections");
		exit(1);
	}
	
	fscanf(fp,"%d\n", &jmax);
	enarr= (float*)malloc((jmax+1)*sizeof(float));
	crsarr1= (float*)malloc((jmax+1)*sizeof(float));

	for (i=0;i<=jmax;i++) {
		fscanf(fp,"%g %g\n", &enarr[i], &crsarr1[i]);
	}

	fclose(fp);

	for (i=0; i<imax; i++) {
		en= de*i;
		if (en <= enarr[0]) {
			tarr[i] = crsarr1[0];
		} else if (en>enarr[jmax]) {
			tarr[i] = crsarr1[jmax];
		} else {
			j = 0;
			
			while(en>enarr[j]) {
				j++;
			}
			
			tarr[i] = (crsarr1[j-1]*(enarr[j]-en) + crsarr1[j]*(en-enarr[j-1]))/
			(enarr[j] - enarr[j-1]);
		}
	}
}
/* End of tcrossec()  */

/********************************************************************/
/** Function for plotting data "vel" [0,nmax] to the ASCII file fx **/

void prtest(float *vel, int nmax) {
	int   ii;
	float mag;
	FILE  *mydmp;
	char  *Haha="fx";

	mydmp = fopen(Haha, "w");

	for(ii=0; ii<VVMAX; ii++) {
		mag = vel[ii];
		fprintf(mydmp,"%f\n", mag);
	}
	
	fclose(mydmp);

} /* End of prtest() */

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