#include "peg_pic.h"

/****** Particle mover, case with b>0, Ey !=0. *****************/

void movebe() {
	int isp, i, j;
	float  vyold, atemp, s, Etemp, a[ng];

	for (isp=0; isp<nsp; isp++) {
		for (i=0; i< ng; i++) av_e[isp][i] +=e[i] ;

		if(dinj[isp]) {
			/* some particles are moving after nstep t.s.*/
			for (j=0; j< ng; j++) {
				a[j] = av_e[isp][j]*fr_E[isp];
			}

			for (i=np[isp]-1; i>=0; i--) {
				if (x[isp][i]<=xnc) {
					j = x[isp][i];
					s = x[isp][i] - j;
					atemp = a[j] + s*(a[j+1] - a[j]);
				} else {
					/* for i_t=2 */
					atemp =0.;
				}

				if (x[isp][i]<x_E) {
					/* for Ey field */
					Etemp= E1y[isp]*x[isp][i];
				} else if (x[isp][i]>(xnc-x_E)) {
					Etemp= E1y[isp]*(xnc-x[isp][i]);
				} else {
					Etemp= Ey[isp];
				}

				vx[isp][i] += atemp;
				vy[isp][i] += Etemp;
				vyold= vy[isp][i];
				vy[isp][i]= vyold*tkv[isp] -sz[isp]*vx[isp][i] +sx[isp]*vz[isp][i];
				vyold +=vy[isp][i];
				vx[isp][i] += tz[isp]*vyold + atemp;
				vy[isp][i] += Etemp;
				vz[isp][i] += -tx[isp]*vyold;
				x[isp][i] += nstep[isp]*vx[isp][i];
			}

			for (i=0; i< ng; i++) {
				av_e[isp][i]= 0.0;    /* init. of av. field */
			}

		} /* end of if(dinj[isp]>0) */
	}
}  /* END of movebe() */

/***************************************************************/
/****** Particle mover, case with b>0, Ey=0.   *****************/

void moveb() {
	int isp, i, j;
	float  vyold, atemp, s, a[ng];

	for (isp=0; isp<nsp; isp++) {
		for (i=0; i< ng; i++) av_e[isp][i] +=e[i] ;

		if(dinj[isp]) {
			/* some particles are moving after nstep t.s.*/
			for (j=0; j< ng; j++) {
				a[j] = av_e[isp][j]*fr_E[isp]; 
			}

			for (i=np[isp]-1; i>=0; i--) {
				if (x[isp][i]<=xnc)  {
					j = x[isp][i];
					s = x[isp][i] - j;
					atemp = a[j] + s*(a[j+1] - a[j]);
				} else {
					atemp =0.;                            /* for i_t=2 */
				}

				vx[isp][i] += atemp;
				vyold= vy[isp][i];
				vy[isp][i]= vyold*tkv[isp] -sz[isp]*vx[isp][i] +sx[isp]*vz[isp][i];
				vyold +=vy[isp][i];
				vx[isp][i] += tz[isp]*vyold + atemp;
				vz[isp][i] += -tx[isp]*vyold;
				x[isp][i] += nstep[isp]*vx[isp][i];
			}

			for (i=0; i< ng; i++) {
				av_e[isp][i]= 0.0;    /* init. of av. field */
			}

		} /* end of if(dinj[isp]>0) */
	}
}  /* END of moveb() */

/***************************************************************/
/****** Particle mover, case with b>0, Ex=Ey=0.   **************/

void moveb0() {
	int isp, i;
	float  vyold;

	for (isp=0; isp<nsp; isp++) {
		if(dinj[isp]) {
			/* some particles are moving after nstep t.s.*/
			for (i=np[isp]-1; i>=0; i--) {
				vyold= vy[isp][i];
				vy[isp][i]= vyold*tkv[isp] -sz[isp]*vx[isp][i] +sx[isp]*vz[isp][i];
				vyold +=vy[isp][i];
				vx[isp][i] += tz[isp]*vyold ;
				vz[isp][i] += -tx[isp]*vyold;
				x[isp][i] += nstep[isp]*vx[isp][i];
			}
		} /* end of if(dinj[isp]>0) */
	}
}  /* END of moveb0() */

/***************************************************************/
/****** Particle mover, case with b=Ey=0.   *****************/

void move1() {
	int isp, i, j;
	float atemp, s, a[ng];

	for (isp=0; isp<nsp; isp++)	{
		for (i=0; i< ng; i++) {
			av_e[isp][i] +=e[i] ;
		}

		if(dinj[isp]) {
			/* some particles are moving after nstep t.s.*/
			for (j=0; j< ng; j++) {
				a[j] = av_e[isp][j]*fr_E[isp];
			}

			for (i=np[isp]-1; i>=0; i--) {
				if (x[isp][i]<=xnc)	{
					j = x[isp][i];
					s = x[isp][i] - j;
					atemp = a[j] + s*(a[j+1] - a[j]);
				} else {
					atemp =0.;                            /* for i_t=2 */
				}

				vx[isp][i] += atemp;
				x[isp][i] += nstep[isp]*vx[isp][i];
			}

			for (i=0; i< ng; i++) {
				av_e[isp][i]= 0.0;    /* init. of av. field */
			}

		} /* end of if(dinj[isp]>0) */
	}
}  /* END of move1() */

/***************************************************************/
/****** Particle mover, case with b=Ex=Ey=0.   *****************/

void move0() {
	int isp, i;

	for (isp=0; isp<nsp; isp++) {
		if(dinj[isp]) {
//printf("jep, moving...np = %d\n", np[isp]);
			/* some particles are moving after nstep t.s.*/
			for (i=np[isp]-1; i>=0; i--) {
//printf("x=%g", x[isp][i]);
				x[isp][i] += nstep[isp]*vx[isp][i];
//printf("....vx=%f --> x=%f\n", vx[isp][i], x[isp][i]);

			}
		} 
	}
}  /* END of move0() */

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