#include "pdp1.h"
#include "xgrafix.h"


// Function void imp_move(int isp) {{{
void imp_move(int isp)
{
	register int i, j;
	float k, s, atemp, ke_temp;

	k = sp_k[isp];
	for (j=0; j<ng; j++) {
		a[j] = e[j]*a_scale[isp];
		if (theRunWithXFlag) sp_j[isp][j] = sp_ke[isp][j] = 0.0;
	}

	for (i= 0; i<np[isp]; i++) {
		//***** Post-push to time n ****
		j = x[isp][i];
		s = x[isp][i] - j;
		atemp= a[j] + s*(a[j+1] - a[j]);

		x[isp][i] += 0.5*k*atemp;
		vx[isp][i]+= atemp;

		//***** Get diagnostics at time n ***
		j = x[isp][i];
		s = x[isp][i] - j;

		if (theRunWithXFlag) {
			ke_temp= vx[isp][i]*vx[isp][i] +vy[isp][i]*vy[isp][i] +vz[isp][i]*vz[isp][i];
			sp_j[isp][j]   += (1-s)*vx[isp][i];
			sp_j[isp][j+1] += s*vx[isp][i];
			sp_ke[isp][j]  += (1-s)*ke_temp;
			sp_ke[isp][j+1]+= s*ke_temp;
		}

		//***** Pre-push to time n+1  ******
		atemp= a[j] + s*(a[j+1] - a[j]);
		x[isp][i] += k*(vx[isp][i] +0.5*atemp);
		vx[isp][i]+= atemp;
	}

	//***** Fix diagnostics and end points and rescale ***
	if (theRunWithXFlag) {
		sp_j[isp][0]  *= 2.;
		sp_j[isp][nc] *= 2.;
		sp_ke[isp][0] *= 2.;
		sp_ke[isp][nc]*= 2.;
		for (j=0; j<ng; j++) jdote[isp][j] += sp_j[isp][j]*e[j];
	}
} // }}}

// Function void exp_move(int isp) {{{
void exp_move(int isp)
{
	register int i, j;
	float k, s, vtemp, atemp, ke_temp, vxold, vyold, vzold;

	vtemp = 0.0;

	k = sp_k[isp];
	for (j=0; j<ng; j++) {
		a[j] = e[j]*a_scale[isp];
		if (theRunWithXFlag) sp_j[isp][j] = sp_ke[isp][j] = 0.0;
	}
	//********** If external B field ********* {{{
	if (b>0.0) {
		for (i=np[isp]-1; i>=0; i--) {
			j = x[isp][i];
			s = x[isp][i] - j;

			if (theRunWithXFlag)
				vtemp = vx[isp][i];

			atemp = a[j] + s*(a[j+1] - a[j]);
			vx[isp][i] += atemp;

			//***** Boris rotation *******
			vxold= vx[isp][i] +tz[isp]*vy[isp][i];
			vyold= vy[isp][i] -tz[isp]*vx[isp][i] +tx[isp]*vz[isp][i];
			vzold= vz[isp][i] -tx[isp]*vy[isp][i];

			vx[isp][i] += sz[isp]*vyold;
			vy[isp][i] += -sz[isp]*vxold +sx[isp]*vzold;
			vz[isp][i] += -sx[isp]*vyold;

			//***** Advence velocity and position *****
			vx[isp][i] += atemp;
			x[isp][i]  += k*vx[isp][i];

			//***** Get diagnostics at time n ***
			if (theRunWithXFlag) {
				vtemp += vx[isp][i];
				ke_temp= vx[isp][i]*vx[isp][i] +vy[isp][i]*vy[isp][i] +vz[isp][i]*vz[isp][i];
				sp_j[isp][j]   += (1-s)*vtemp;
				sp_j[isp][j+1] += s*vtemp;
				sp_ke[isp][j]  += (1-s)*ke_temp;
				sp_ke[isp][j+1]+= s*ke_temp;
			}
		}
	} // }}}
	//********** If NO external B field ********* {{{
	else {
		for (i=np[isp]-1; i>=0; i--) {
			j = x[isp][i];
			s = x[isp][i] - j;

			if (theRunWithXFlag)
				vtemp = vx[isp][i];

			//***** Advence velocity and position *****
			vx[isp][i]+= a[j] + s*(a[j+1] - a[j]);
			x[isp][i] += k*vx[isp][i];

			//***** Get diagnostics at time n ***
			if (theRunWithXFlag) {
				vtemp += vx[isp][i];
				ke_temp= vx[isp][i]*vx[isp][i] +vy[isp][i]*vy[isp][i] +vz[isp][i]*vz[isp][i];
				sp_j[isp][j]   += (1-s)*vtemp;
				sp_j[isp][j+1] += s*vtemp;
				sp_ke[isp][j]  += (1-s)*ke_temp;
				sp_ke[isp][j+1]+= s*ke_temp;
			}
		}
	} // }}}

	//***** Fix diagnostics and end points and rescale ***
	if (theRunWithXFlag) {
		sp_j[isp][0]  *= 2.;
		sp_j[isp][nc] *= 2.;
		sp_ke[isp][0] *= 2.;
		sp_ke[isp][nc]*= 2.;
		for (j=0; j<ng; j++) jdote[isp][j] += sp_j[isp][j]*e[j];
	}
} // }}}


// Modeline for ViM {{{1
// vim:set ts=4:
// vim600:fdm=marker fdl=0 fdc=3:
// 1}}}