#include <stdlib.h>

#include "pdp1.h"


// Function void load(int initnp[][2]) {{{
void load(int initnp[][2])
{
	// extern double fdist();
	extern double fdist(double, int, int);
	double vv;
	int i, ix, k, isp, ksign, ii;
	float ddx, vmag, theta, v1, v2, ddv, sf, ddf, sumf,
	ss, ds, ff, vel[1024];

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

	for (isp=0; isp<nsp; isp++) {
		np[isp]= 0;

		for (k=0; k<2; k++) {
			if (!initnp[isp][k]) continue;

			ddx= xnc/initnp[isp][k];
			if (np[isp] + initnp[isp][k] > maxnp[isp]) {
				puts("LOAD: too many particles, species ");
				exit(1);
			}
			if (vt[isp][k]==0. ) {    /* Loader for a COLD beam   */
				for (i=0; i<initnp[isp][k]; i++) {
					ix = i +np[isp];
					if (!i)  x[isp][ix]= .5*ddx;
					else    x[isp][ix] = x[isp][ix-1] +ddx;
					vx[isp][ix] = v0[isp][k];
					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<initnp[isp][k]; i++) {
						ix = np[isp] +i;
						x[isp][ix] += dde*sin(2*M_PI*x[isp][ix]/xnc);
					}
				}
			}
			else {    /* Loader for THERMAL distribution */
				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.;

L12:    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;
		if (sumf < ss) goto L12;

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

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

		for (i=0; i<initnp[isp][k]; i++) {
			ix = i +np[isp];
			if (!i)  x[isp][ix]= .5*ddx;
			else    x[isp][ix] = x[isp][ix-1] + ddx;
			vx[isp][ix] = vel[Index[(i+k)%1024]];
			/* vx[isp][ix] = vel[Index[ix%1024]]; */
			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<initnp[isp][k]; i++) {
				ix = np[isp] +i;
				x[isp][ix] += dde*sin(2*M_PI*x[isp][ix]/xnc);
			}
		}
			}   /* end of THERMAL loader  "if ( VT>0.." */
			np[isp] += initnp[isp][k];

		}   /* 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}}}