#include <stdlib.h>

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

//**************************************************************
// Routine to adjust (initialize, re-pack and inject) particles
//  to the desired boundary conditions

// PG for volume source

void sterm(int count);
static float Svel[2048][NSMAX] ;
static int npold, ionsp;

// Function void adjust(int isp) {{{
void adjust(int isp)
{
	static int secsp, init_flag=1;
	static float extra[NSMAX][2];

	register int i, ii, j, k;
	int nnp, secountl, secountr;
	int nreflux[NSMAX], ninjected, s;
	float dum, del_t, vmag, theta;
	// extern float vnext();
	extern float vnext(int, int);

	secountl = secountr = 0;

	// INITIALIZE array for computing positions of injected particles
	if (init_flag) {
		secsp = secondary -1;

		// "enter" is now the no. of particles injected each "dt"
		for (i=0; i<nsp; i++) {
			if (fabs(enter[i][0]) > 0.)
				extra[i][0] = 0.5123123;
			else
				extra[i][0] = 0.;
			if (fabs(enter[i][1]) > 0.)
				extra[i][1] = 0.5123123;
			else
				extra[i][1] = 0.;
		}
		init_flag = 0;
	}

	if (psource && (isp == ionsp)) {
		npold = np[ionsp]; // Save the number of ions before adjust
	}


	jwall[isp] = 0.0;
	if (secondary)
		secountl = secountr = 0;
	nreflux[isp]= 0;    // nreflux[isp] equals zero

	if (np[isp] > 0) {
		nnp = np[isp] -1;
		i = 0;

		// eliminate "outsiders", allow for secondary electron emission,
		// and if it left thru LH plate, ADD charge there to sigma
		// (plate surface density).
		do {
			if (x[isp][i] >= xnc) {
				x[isp][i] = x[isp][nnp];
				vx[isp][i] = vx[isp][nnp];
				vy[isp][i] = vy[isp][nnp];
				vz[isp][i] = vz[isp][nnp];
				nnp--;
				if (reflux)  nreflux[isp]++;
				else if (secondary) {
					if (frand() < seec[isp]) {
						secountr++;
					}
				}
			}
			else if (x[isp][i] <= 0) {
				//********* LHS wall diagnostics   *************** {{{
				if (theRunWithXFlag) {
					dum = (vx[isp][i]*vx[isp][i] + vy[isp][i]*vy[isp][i]
							+vz[isp][i]*vz[isp][i] - emin[isp])/de[isp];
					s = dum;
					if (s<nbin[isp]-1 && dum>=0) {
						dum -= s;
						fe[isp][s]  += (!s) ? 2*(1-dum) : 1-dum;
						fe[isp][s+1]+= (s==nbin[isp]-2) ? 2*dum : dum;
					}

					dum = -atan(sqrt(vy[isp][i]*vy[isp][i]+vz[isp][i]*vz[isp][i])
							/vx[isp][i])/dtheta[isp];
					s = dum;
					dum -= s;
					ftheta[isp][s]  += (!s) ? 2*(1-dum) : 1-dum;
					ftheta[isp][s+1]+= (s==nbin[isp]-2) ? 2*dum : dum;
				} // }}}
				x[isp][i]  = x[isp][nnp];
				vx[isp][i] = vx[isp][nnp];
				vy[isp][i] = vy[isp][nnp];
				vz[isp][i] = vz[isp][nnp];
				nnp--;
				jwall[isp] += q[isp];
				if (secondary) {
					if (frand() < seec[isp]) {
						secountl++;
						jwall[secsp] -= q[secsp];
					}
				}
			}
			else {
				//*****  MID system diagnostics  *********** {{{
				if (theRunWithXFlag) {
					if (xs_mid[isp] <= x[isp][i] && x[isp][i] <= xf_mid[isp]) {
						s = (vx[isp][i]*vx[isp][i]
								+ vy[isp][i] * vy[isp][i]
								+ vz[isp][i] * vz[isp][i]
								- emin_mid[isp]) / de_mid[isp];

						if (0 <= s && s< nbin_mid[isp])
							fe_mid[isp][s] += 1;
					}
				} // }}}
				i++;
			}
		} while (i <= nnp);
		np[isp] = nnp + 1;
	}

	// INJECT new particles at walls, one species at a time {{{

	if (inject[isp] || nreflux[isp]) {
		for (k=0; k<2; k++) {
			extra[isp][k] += enter[isp][k];
			if (reflux && k == 1) extra[isp][k] += nreflux[isp];
			// "extra" = the fraction of "dt" when emitted
			ninjected = extra[isp][k];

			while (extra[isp][k] >= 1.0) {
				del_t= ((int)extra[isp][k])/(ninjected +1.);
				ii = np[isp];
				np[isp]++;

				if (ii > maxnp[isp]) {    // Move array boundaries here
					printf("ADJUST: too many particles, species %d",isp);
					exit(1);
				}

				// Choose V's
				vx[isp][ii] = vnext(k,isp);
				vmag= vperp0[isp] +vperpt[isp]*sqrt(2*fabs(log(frand()+DBL_MIN)));
				theta= 2*M_PI*frand();
				vy[isp][ii] = vmag*sin(theta);
				vz[isp][ii] = vmag*cos(theta);

				// Adjust Vx for effect of electric field
				vx[isp][ii] += (del_t-0.5)*e[k*nc]*a_scale[isp];
				x[isp][ii] = k*nc + del_t*vx[isp][ii];
				if (!k) jwall[isp] -= q[isp];
				extra[isp][k] -= 1.0;
			}
		}
	} // }}}

	if (secondary) {
		i = np[secsp];
		np[secsp] += secountl +secountr;

		if (np[secsp] > maxnp[isp])
			printf("ADJUST(Secondaries): too many particles. MUST EXIT!");

		for (j=i; j< i+secountl; j++) {
			x[secsp][j]= frand();
			vx[secsp][j]= vnext(0,secsp);
			vmag= vperp0[secsp]
				+ vperpt[secsp] * sqrt(2*fabs(log(frand()+DBL_MIN)));
			theta= 2*M_PI*frand();
			vy[secsp][j] = vmag*sin(theta);
			vz[secsp][j] = vmag*cos(theta);
		}
		i += secountl;
		for (j=i; j< i+secountr; j++) {
			x[secsp][j]= xnc - frand();
			vx[secsp][j]= vnext(1,secsp);
			vmag= vperp0[secsp]
				+ vperpt[secsp] * sqrt(2*fabs(log(frand()+DBL_MIN)));
			theta= 2*M_PI*frand();
			vy[secsp][j] = vmag*sin(theta);
			vz[secsp][j] = vmag*cos(theta);
		}
	}
	jwall[isp] /= sp_k[isp]*dt;

	// try to add a volume source of ionization
	if (psource && (isp == ionsp)) {
		if (t>tstrt)
			sterm(np[ionsp]-npold); // Create ion/elec pair for each ion lost ?
	}

} // }}}

//***********************************************************
// Function to give bit-reversed velocities to the injection
// scheme of ADJUST, one at a time

// Function float vnext(int k, int isp) {{{
float vnext(int k, int isp)
{
	static int iv[NSMAX][2], qinit=1;
	static double svf[NSMAX][2];
	static float vel[1024][NSMAX][2];
	int iiv, ksign, ii, i;
	double v1, v2, ddv, vv, vf, df, sumvf, ss, ds;
	// extern double fdist(), revers();
	extern double fdist(double, int, int);
	extern double revers(unsigned int);

	// At t=0., set up iv arrays
	if (qinit) {
		for (i=0; i<nsp; i++) iv[i][0]= iv[i][1]= 0;
		qinit= 0;
	}

	// For a COLD beam, injection is simple:
	if (vt[isp][k]==0.) return (v0[isp][k]);

	if (iv[isp][k]==0) {
		// Integrate velocity-distribution function
		// using Simpsons's rule (to normalize it)
		if (k==0) {
			v1 = max(vc[isp][0], v0[isp][0]-6.*vt[isp][0]);
			v2 = max(v0[isp][0]+6.*vt[isp][0], v1+3.*vt[isp][0]);
		}
		else {
			v1 = min(vc[isp][1], v0[isp][1]+6.*vt[isp][1]);
			v2 = min(v0[isp][1]-6.*vt[isp][1], v1-3.*vt[isp][1]);
		}

		ddv = (v2-v1)/2048.;
		vv = v1;
		svf[isp][k] = vv*fdist(vv,k,isp)/2.;

		for (i=0; i<1023; i++) {
			vv += ddv;
			df = fdist(vv,k,isp);
			svf[isp][k] += vv*(df+df);
			vv += ddv;
			df = fdist(vv,k,isp);
			svf[isp][k] += vv*df;
		}   // end "for (i...)" loop

		df = fdist(v2,k,isp)/2.;
		svf[isp][k] = (svf[isp][k] - v2*df)*2.*ddv/3.;
	}   // end "if (iv[...)" loop

	iiv = iv[isp][k]%1024;
	if (iiv==0) {
		ksign =  1 - 2*k;
		// ksign = +1(k=0)  or -1(k=1)
		sumvf = 0.;
		ii = 0;
		if (k==0) vv = max(vc[isp][0], v0[isp][0]-6.*vt[isp][0]);
		else vv = min(vc[isp][1], v0[isp][1]+6.*vt[isp][1]);
		ds = svf[isp][k]/1024.;
		iiv = iv[isp][k]/1024.;
		ss = revers(iiv)*ds;

L11: ddv = ksign*vt[isp][k]/128.;
	 vf = vv*fdist(vv,k,isp);

	 ddv=ksign*min(fabs(ddv),.45*ds/fabs(vf+DBL_MIN));
	 // Vary the integration step to aim for the next
	 //	desired value for the integral

	 do {
		 df = ddv*(vf + 4.*(vv+ddv)*fdist(vv+ddv,k,isp) +
				 (vv+ddv+ddv)*fdist(vv+ddv+ddv,k,isp))/3.;
		 ddv *= 0.5;
	 } while (df >= ds);
	 // If integral increases too much, go back w. dv/2 step
	 ddv *= 2.;

	 vv += 2.*ddv;
	 sumvf += df;
	 if (sumvf < ss) goto L11;

	 vel[ii][isp][k] = vv - (sumvf-ss)*ddv/df;
	 // Interpolate to find the desired velocity

	 ii += 1;
	 ss += ds;
	 if (ii < 1024) goto L11;

	 if (iv[isp][k]==0) iv[isp][k] = 1;
	 // Skip the very first velocity
	}  // end "if (iiv..)"  loop

	ii = iv[isp][k]%1024;
	iiv = Index[ii];

	iv[isp][k] += 1;
	if (iv[isp][k] < 0) iv[isp][k] = 1;
	// Keep track of the total number of particles injected

	return (vel[iiv][isp][k]);
} // }}}

// Function void sterm(int count) {{{
void sterm(int count)
{
	extern double fdist();
	int i, ix, j, k, in, isp, initnp[NSMAX][2];
	float vmag, theta, r;

	//**********************************************
	// Load one species and one direction at a time



	for (isp=0; isp<nsp; isp++) {
		if ((np[isp] - count) > maxnp[isp]) {
			for (k=0; k<2; k++) {
				printf("np[isp] initnp[isp][%d] maxnp[isp] %d %d %d \n  ",
						k, np[isp], initnp[isp][k], maxnp[isp]);
			}
			puts("source term: too many particles, species ");
			putchar(isp+49);
			exit(1);
		}
	}

	if (count < 0) { // add only when there is a loss
		for (j=0; j<abs(count); j++) {
			r=frand();

			// Loading the ions (isp=1)  randomly in the plasma
			isp=1;
			ix = np[isp];
			x[isp][ix]= (xnc-10e-6)*r;
			in=2048.0*frand();
			vx[isp][ix] = Svel[in][isp];
			vmag= vperp0[isp]+vperpt[isp]*sqrt(2*fabs(log(frand()+DBL_MIN)));
			theta= 2*M_PI*frand();
			vy[isp][ix] = vmag*sin(theta);
			vz[isp][ix] = vmag*cos(theta);

			if (dde) {
				for (i=0; i<1; i++) {
					x[isp][ix] += dde*sin((2*M_PI*x[isp][ix]/xnc));
				}
			}
			np[isp]++;


			//	 printf ("added ion %f %e\n", x[isp][ix], vx[isp][ix]);

			// Loading the electrons  (isp=0)  randomly in the plasma
			isp=0;
			ix = np[isp];
			x[isp][ix]= (xnc-10e-6)*r;
			in=2048.0*frand();
			vx[isp][ix] = Svel[in][isp];
			vmag= vperp0[isp]+vperpt[isp]*sqrt(2*fabs(log(frand()+DBL_MIN)));
			theta= 2*M_PI*frand();
			vy[isp][ix] = vmag*sin(theta);
			vz[isp][ix] = vmag*cos(theta);

			if (dde) {
				for (i=0; i<1; i++) {
					x[isp][ix] += dde*sin((2*M_PI*x[isp][ix]/xnc));
				}
			}
			np[isp]++;
			//printf ("added electron %f %e\n", x[isp][ix], vx[isp][ix]);

		}   // end of THERMAL loader  "if ( VT>0.."

		//	printf ("\n\n");
	}

} // }}}

// Function void svel(void) {{{
void svel(void)
{
	extern double fdist();
	double vv;
	int i, k, ksign, ii, isp;
	float v1, v2, ddv, sf, ddf, sumf, ss, ds, ff;


	// old version had ionsp declared global, do it here now

	ionsp = ionspecies-1;

	for (isp=0; isp<nsp; isp++) {
		for (k=0; k<2; k++) {
			ksign = 1 -2*k;

			// Integrate the [dist fcn] from    v1 to v2,
			// using Simpson's rule, to normalize them

			if (k==0) {
				v1 = max(vc[isp][0], v0[isp][0]-6.*vt[isp][0] );
				v2 = max(v0[isp][0]+6.*vt[isp][0], v1+3.*vt[isp][0]);
			}
			else {
				v1 = min(vc[isp][1], v0[isp][1]+6.*vt[isp][1] );
				v2 = min(v0[isp][1]-6.*vt[isp][1], v1-3.*vt[isp][1]);
			}

			ddv = (v2-v1)/2048.;
			vv = v1;
			sf = fdist(vv,k,isp)/2.;

			for (i=0; i<1024; i++) {
				vv += ddv;
				ddf = fdist(vv, k, isp);
				sf += ddf + ddf;
				vv += ddv;
				ddf = fdist(vv, k, isp);
				sf += ddf;
			}
			ddf = fdist(v2, k, isp)/2;
			sf -= ddf;
			sf = (sf + sf)*fabs(ddv)/3.;
			// Integrate again, varying the stepsize to aim for the values
			// of the integral at which particle velocities are desired

			sumf = 0.;
			ii = 0;
			vv = v1;
			ds = sf/1024.;
			ss = ds/2.;
			do {
				do {
					ff = fdist(vv, k, isp);
					ddv = ksign*vt[isp][k]/128.;
					if (ff > 0.)
						ddv = ksign*min(fabs(ddv), .45*ds/ff);
					do {
						ddf = fdist(vv, k, isp) + 4.*fdist(vv + ddv, k, isp)
							+ fdist(vv+ddv+ddv, k, isp);
						ddf *= fabs(ddv)/3.;
						if (ddf >= ds) ddv *= 0.5;
					} while (ddf >= ds);

					vv += ddv+ddv;
					sumf += ddf;
				} while (sumf < ss);

				// When a desired velocity has been passed, interpolate
				//   to find it, and start looking for the next velocity
				Svel[ii][isp] = vv - (sumf - ss)*ddv/ddf;
				Svel[2047-ii][isp]=-Svel[ii][isp];
				ii++;
				ss += ds;
			} while (ii < 1024);

			// Compute the number of particles to be loaded, and load the
			// 1024 velocities uniformly in space, in bit-rev order

		}   // end  "for ( k=0,1...)" loop
	}  // end for (isp=0 thru nsp-1..)" loop

} // }}}

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