#include "peg_pic.h"

/** Cross section of electron-neutral particle elastic collisions **/

float sigma1(int jj, float energy) {
	int 		i, j;
	float 	deng;
	static float 	**sels, **sels1;
	static int 	init_flag=1;
	char  	tchar[80];

	/********  Initialization   *************/
	if(init_flag) {
		sels= (float **)malloc(Nnum*sizeof(float *));
		sels1= (float **)malloc(Nnum*sizeof(float *));
		
		for (i=0; i<Nnum; i++) {
			sels[i]= (float *)malloc(1000*sizeof(float));
			sels1[i]= (float *)malloc(NEMAX*sizeof(float));
		}

		for (j=0; j<Nnum; j++) {
			if (Hyd[j]) {                   /* exact cross-sections */
				dec1(tchar, &i, 0, Hyd[j], 0, 0, 0);

				if (!i) {
					/* reading from the table */
					tcrossec(sels[j], tchar, 1000, 0.01);
					tcrossec(sels1[j], tchar, NEMAX, 1);
				} else  {
					/* reading from the file */
					for (i=0; i<1000; i++) {
						deng= 0.01*i;
						sels[j][i] = (*fcrs)(deng);
					}
					
					for (i=0; i<NEMAX; i++) {
						sels1[j][i] = (*fcrs)((float)i);
					}
				}
			} else {
				/* approximated cross-sections */
				for (i=0; i<1000; i++) {
					deng= 0.01*i;
					sels[j][i] = forig_el(deng, j);
				}    
				
				for (i=0; i<NEMAX; i++) {
					sels1[j][i] = forig_el((float)i, j);
				}
			}

			if (dshort>0.) {
				for (i=0; i<1000; i++) sels[j][i] *=dshort;
				for (i=0; i<NEMAX; i++) sels1[j][i] *=dshort;
			}
		}
		init_flag =0;
	} /* end of if (init_flag) */

	if (energy<10) {
		i= 100*energy;
		return(sels[jj][i]);
	} else {
		i= energy + 0.5;
		
		if(NEMAX <= i) {
			i= NEMAX-1;
		}
		
		return(sels1[jj][i]);
	}

} /* END of sigma1 */

/*******************************************************************/
/**** Cross section of ion-neutral particle elastic collisions *****/

float sigma4(int ii, int k, float energy) {
	int i, j, p;     
	float deng;
	static float ***siel1, ***siel2, ***siel3;
	static int init_flag=1;
	char  tchar[80];                                                                                      

	/**********  Initialization  **************/
	if(init_flag) {
		p=0;
		
		for (i=0; i<Nnum; i++) {
			for (j=0; j<Nnum; j++) {
				if (HH[i][j]) {
					p +=1;
				}
			}
		}

		if (p) {
			/* allocate space for exact cross-sections */
			siel1= (float ***)malloc(Nnum*sizeof(float **));
			siel2= (float ***)malloc(Nnum*sizeof(float **));
			siel3= (float ***)malloc(Nnum*sizeof(float **));

			for (i=0; i<Nnum; i++) {
				siel1[i]= (float **)malloc(Nnum*sizeof(float *));
				siel2[i]= (float **)malloc(Nnum*sizeof(float *)); 
				siel3[i]= (float **)malloc(Nnum*sizeof(float *));

				for (j=0; j<Nnum; j++) {
					siel1[i][j]= (float *)malloc(1000*sizeof(float));
					siel2[i][j]= (float *)malloc(200*sizeof(float));
					siel3[i][j]= (float *)malloc(NEMAX*sizeof(float));
				}
			}                                                                 
		}   /* end of if (p) loop */

		for (i=0; i<Nnum; i++) {
			for (j=0; j<Nnum; j++) {
				if (HH[i][j]) {
					/* exact cross-sections */
					dec1(tchar, &p, 1, HH[i][j], 0, 0, 0);

					if (!p) {
						/* reading from the table */
						tcrossec(siel1[i][j], tchar, 1000, 1.e-4);
						tcrossec(siel2[i][j], tchar, 200, 0.1);
						tcrossec(siel3[i][j], tchar, NEMAX, 1);
					} else {
						/* reading from the file */
						for (p=0; p<1000; p++) {
							deng= 1.e-4*p;
							siel1[i][j][p] = (*fcrs)(deng);
						}

						for (p=0; p<200; p++) {
							deng= 0.1*p;
							siel2[i][j][p] = (*fcrs)(deng);
						}

						for (p=0; p<NEMAX; p++) {
							siel3[i][j][p] = (*fcrs)((float)p);
						}
					}

					if (dshort>0.) {
						for (i=0; i<Nnum; i++) {
							for (j=0; j<Nnum; j++) {
								for (p=0; p<1000; p++) {
									siel1[i][j][p] *= dshort;
								}
								
								for (p=0; p<200; p++) {
									siel2[i][j][p] *= dshort;
								}
								
								for (p=0; p<NEMAX; p++) {
									siel3[i][j][p] *= dshort;
								}
							}
						}
					}
				} else if (dshort>0.) {
					/* approximated cross-sections */
					ascat[i][j] *=dshort;
					bscat[i][j] *=dshort;
				}
			}
		}      /* end of for (i=0; i<Nnum; i++) loop */

		init_flag =0;

	}  /* end of initialization */

	if (HH[ii][k]) {
		if (energy<0.1) {
			i= 10000*energy;
			//printf("ii= %d, k= %d, i= %d, t= %g\n", ii,k,i,t);
			//printf("i= %f\n", siel2[ii][k][i]);
			return(siel1[ii][k][i]);
		}

		if (energy<20.0) {
			i= 10*energy;
			//printf("ii= %d, k= %d, i= %d, t= %g\n", ii,k,i,t);
			//printf("i= %f\n", siel2[ii][k][i]);
			return(siel2[ii][k][i]);
		} else if (energy<NEMAX) {
			//printf("ii= %d, k= %d, i= %d, t= %g\n", ii,k,i,t);
			//printf("i= %f\n", siel2[ii][k][i]);
			i= energy;
			return(siel3[ii][k][i]);
		} else {
			return(0.);
		}
	} else {
		return(ascat[ii][k] + bscat[ii][k]/(sqrt(energy) +1e-30));
	}
}
/* END of sigma4 */                                                                                     

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

