#include "peg_pic.h"

/* Functions used for charged-neutral particle collision cross-sections */

/****************** e+neutral collisions ****************/

/*******************************************************/
/* Integral elastic cross-sections (original input) */

float forig_el(float en, int j) {
	float  mag;

	if (en <= elsengy0[j]) {
		mag = 0;

	} else if (en > elsengy0[j] && en <= elsengy1[j]) {
		mag = selsmax[j]*(en - elsengy0[j])/(elsengy1[j] - elsengy0[j]);
		
	} else if (en > elsengy1[j] && en <= elsengy2[j]) {
		mag = selsmax[j];

	} else {
		mag = selsmax[j]*elsengy2[j]/log(elsengy2[j])*log(en)/en;
	}

	return(mag);
}

/*******************************************************/
/* Integral excitation cross-sections (original input) */


float forig_ex(float en, int j) {
	float  mag;

	if (en <= ext0[j]) {
		mag = 0;
	} else if (en > ext0[j] && en <= extengy1[j]) {
		mag = sextmax[j]*(en - ext0[j])/(extengy1[j] - ext0[j]);
	} else if (en > extengy1[j] && en <= extengy2[j]) {
		mag = sextmax[j]; 
	} else {
		mag = sextmax[j]*extengy2[j]/log(extengy2[j])*log(en)/en;
	}

	return(mag);
}

/*******************************************************/
/* Integral ionization cross-sections (original input) */

float forig_i(float en, int j) {
	float mag;

	if (en <= ionengy0[j]) {
		mag = 0.;
	} else if (en > ionengy0[j] && en <= ionengy1[j]) {
		mag = sionmax[j]*(en - ionengy0[j])/(ionengy1[j] - ionengy0[j]);
	} else if (en > ionengy1[j] && en <= ionengy2[j]) {
		mag = sionmax[j];
	} else {
		mag = sionmax[j]*ionengy2[j]/log(ionengy2[j])*log(en)/en;
	}

	return(mag);
}

/*******************************************************/
/** Integral ionization cross-sections for e+H, D, T ***/

float feh_i(float en) {
	float ionengy = 13.6, sionmax  = 5.511e-20, eps_p, mag;

	if (en <= ionengy) {
		return(0.);
	}

	eps_p = 1. + en/ionengy;
	mag = sionmax*(1. - 1/(eps_p-1.) + log(eps_p-1.)*(0.18102 -
			1./eps_p + 0.028706/(eps_p*eps_p) - 2.0055/(eps_p*eps_p*eps_p) +
			1.1819/(eps_p*eps_p*eps_p*eps_p) - 0.36406/(eps_p*eps_p*eps_p*eps_p*eps_p)))
			/(1.+eps_p);

	return(mag);
}

/************************************************************/
/** Energy of the new electron after e+H, D, T ionization **/
/** "e_p" is the normalised energy of the primary electron **/

float feh_ecr(float e_p) {
	float rengy;

	do rengy = 0.5*(e_p - 1.)*frand() + 1.;
	while ((frand()*(1.- 1./e_p + 1./e_p/e_p +
			0.49911*log(e_p))) > ((1./rengy/rengy+1./(e_p-rengy+1.)/(e_p-rengy+1.)-
					1./rengy/(e_p-rengy+1.) + log(e_p)*(-0.014353+0.75206/rengy-
							0.29548/(rengy*rengy) + 0.056884/(rengy*rengy*rengy))/
							(rengy*rengy*rengy))));

	return(rengy);
} 

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