#include "peg_pic.h"


/***************************************************************/
/** Scattering of particles; uses primary energy "energy" and **/
/** velocity "vel". "e-flag" is the flag for electron elastic **/
/** collisions. If =1, then the energy gained by neutral *******/
/** particle during the elastic collison is substracted from  **/
/** the electron energy. "itype" is the type of colliding ******/
/** neutral particle, and "ang" is the coefficient for the *****/
/** scattaring angle (see the manual)  *************************/

void newvel(float energy, float vel, float *vx, float *vy, float *vz, int itype, int e_flag, float ang) {
	int      i;
	float    mag;
	double   ux, uy, uz, up, c_th, s_th, sc, ss, dux, duy, duz;
	static float *mm;
	static int   init_flag=1;

	if(energy < 1e-30) {
		c_th = 1.;
	} else {
		if (ang>1.e-6) {
			// anizotropic scattering
			energy *= ang;
			c_th = (energy +2 -2*pow(energy +1,frand()))/energy;
		} else {
			// izotropic scattering
			c_th = 1. - 2.0*frand();
		}
	}

	s_th = sqrt(fabs(1. - c_th*c_th));
	mag=frand()*TWOPI;

	if (e_flag) {
		/* for electron elastic collisions */
		if (init_flag) {
			mm = (float *)malloc(Nnum*sizeof(float));
			init_flag = 0;

			for (i=0; i<Nnum; i++) {
				mm[i] = 2*m[ecolsp]/m[icolsp[i]];
			}
		}
		
		vel *= sqrt(1 - mm[itype]*(1-c_th));
	}

	ux = *vx;
	uy = *vy;
	uz = *vz;
	up = sqrt(ux*ux + uy*uy);

	if (up) {
		s_th /= up;
		ss = s_th*sin(mag);
		sc = s_th*cos(mag);


		dux = ux*uz*sc - uy*ss + ux*c_th;
		duy = uy*uz*sc + ux*ss + uy*c_th;
		duz = -up*up*sc + uz*c_th;
		/*
  printf("2 vx= %e, vy= %e, vz= %e, up= %e, sin= %e, cos= %e\n", dux , duy, duz, up, ss, sc);
		 */
	} else {
		dux = s_th*cos(mag);
		duy = s_th*sin(mag);
		duz = c_th;
	}

	*vx = vel*dux;
	*vy = vel*duy;
	*vz = vel*duz;

} /* END of newvel() */

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