#include "peg_pic.h"

/** Calculation of the primary (en1) and secondary (en2) electron     **/
/** energies after the ionization collision. "i" is the neutral type. **/
/** If "Hydion=0, then uses the original expression from XPDP1 code;  **/
/** Else uses the exact formula.                                      **/ 

void Escat(float *en1, float *en2, int i) {
	float e_p, engy;

	if (Hydion[i]) {
		if (Hydion[i]<4) {
			/* e+H (D, T) collision */
			e_p = *en1/ionengy0[i];
			*en2 = (feh_ecr(e_p)-1.)*ionengy0[i];
			*en1 -= (*en2 + ionengy0[i]);
		} else if (Hydion[i]==4) {
			/* e+O2 collision */
			engy = *en1 - ionengy0[i];
			*en2 = 10.0*tan(frand()*atan(engy/20.0));
			*en1 = engy - *en2;
		} else {
			/* other ionizations have to be added here **/
			printf("Data for ionization is missing\n");
			exit(1);
		}
	} else {
		engy = *en1 - ionengy0[i];
		*en2 = 10.0*tan(frand()*atan(engy/20.0));
		*en1 = engy - *en2;
	} 
} /* END of Escat() */

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

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

float sigma3(int jj, float energy) {
	int 		i, j;
	float		deng;
	static float 	**sion, **sion1;
	static int 	init_flag=1;
	char  	tchar[80];

	/********  Initialization   **********/
	if(init_flag) {
		sion= (float **)malloc(Nnum*sizeof(float *));
		sion1= (float **)malloc(Nnum*sizeof(float *));
		
		for (i=0; i<Nnum; i++) {
			sion[i]= (float *)malloc(200*sizeof(float));
			sion1[i]= (float *)malloc(NEMAX*sizeof(float));
		}

		for (j=0; j<Nnum; j++) {
			if (Hydion[j]) {
				dec1(tchar, &i, 0, Hydion[j], 2, j, 0);


				if (!i) {
					/* reading from the table */
					tcrossec(sion[j], tchar, 200, 0.1);
					tcrossec(sion1[j], tchar, NEMAX, 1);
				} else {
					/* reading from the file */
					for (i=0; i<200; i++) {
						deng= 0.1*i;
						sion[j][i] = (*fcrs)(deng);
					}
					
					for (i=0; i<NEMAX; i++) {
						sion1[j][i] = (*fcrs)((float)i);
					}
				}
			} else {
				/* approximated cross-sections */
				for (i=0; i<200; i++) {
					deng = 0.1*i;
					sion[j][i] = forig_i(deng, j);
				}
				
				for (i=0; i<NEMAX; i++) {
					sion1[j][i] = forig_i((float)i, j);
				}
			}

			if (dshort>0.) {
				for (i=0; i<200; i++) {
					sion[j][i] *=dshort;
				}
				
				for (i=0; i<NEMAX; i++) {
					sion1[j][i] *=dshort;
				}
			}
		}
		init_flag =0;
	}
	/**************************************/

	if (energy<20) {
		i= 10*energy;
		return(sion[jj][i]);
	} else {
		i= energy + 0.5;
		
		if(NEMAX <= i) {
			i= NEMAX-1;
		}
		
		return(sion1[jj][i]);
	}

} /* END of sigma3 */
/*******************************************/
