/********************************************************************
 *                                                                  *
 *      Compute particle motions for electron blobs leaving an      *
 *  electron source shaped like a sphere at either cusp or coil     *
 *  face locations.  Every time step tracks more particles from     *
 *  source point and computes total electric and magnetic fields    *
 *  create from all previous sources as they interact with each     *
 *  other.                                                          *
 *                                                                  *
 *              Author = Mike Rosing                                *
 *              date = March 12, 2008                               *
 *                                                                  *
 *******************************************************************/

#include <stdio.h>
#include <math.h>

#define MAXSTEPS 25
#define NUMELMNTS ((MAXSTEPS+1)*(MAXSTEPS+2)*(MAXSTEPS+3)/6)

/*  MAXBLOBS is the storage limit for electron blobs to follow  */

#define MAXBLOBS 100000

/*  MAXWALL is the ratio of outer containment sphere radius to L (center
    of sphere to center of coil distance)  */

#define MAXWALL  (3.0)

/*  COILRADIUS is the ratio of coil radius to L  */

#define COILRADIUS (0.9)

/*  COILDIAMETER is the effective cross section of the coil, also
    in reference to L  */

#define COILDIAMETER  (0.01)

/*  ANGLEPS is the smallest allowed angle for a match.  I set it
    fairly large, but it could also be a function of MAXSTEPS.
*/

#define ANGLEPS  1e-4

/*  Cp is a constant that makes the fluid equation dimensionless.  It turns
    out it this combination of constants happens a lot, and arises any where
    the magnetic field does.  The formula for Cp is

         Cp^2 = 8*pi^2 * m * V / ( mu_0^2 * I_0^2 * e )

where m is electron mass, V is MaGrid potential, mu_0 is permiability of 
free space, I_0 is current in MaGrid coils and e is the charge on an
electron.  Thus the ratio of MaGrid voltage to square of current is an
important factor in the system.
*/

#define CONST_CP (0.02)

/*  SOURCEPOSITION is distance from center of polywell to center of
    electron source sphere.  It is in units of L.
*/

#define SOURCEPOSITION 2.6

/*  SOURCERADIUS is the radius of the electron source, in units of L  */

#define SOURCERADIUS 0.1

/*  SOURCEROWS are the number of rows to use on the electron source.  Total
    number of electron blobs will be SOURCEROWS*(SOURCEROWS+1)/2, but the
    number tracked will be 48 times that due to symmetry.
*/

#define SOURCEROWS  10
#define SOURCE_ELEMENTS (SOURCEROWS+1)*SOURCEROWS/2

/*  SOURCEVOLTAGE is the negative voltage applied to the electron source.
    It is relative to the MaGrid voltage so that the actual voltage on the
    electron source is SOURCEVOLTAGE * V.
*/

#define SOURCEVOLTAGE (-0.1)

/*  MAXTIMESTEPS is the number of times to do the calculation loop
    EMISSIONTIME is the number of time steps to allow electrons to leave 
    the electron source.
*/

#define MAXTIMESTEPS 10
#define EMISSIONTIME   5

/*  DELTANU is the dimensionless time scale.  It is

        DELTANU = e*mu0*I0/ ( 4*PI*m*L)  * delta_t

    For 200k Ampturns, L = 1m, a value of DELTANU = 1 corresponds to 1/3 of
    a nanosecond.
*/

#define DELTANU  (0.10)

/*  v0 is the velocity of the an electron accelerated to full MaGrid voltage
    assuming Newtonian energy transfer:

        v0 = sqrt(2*e*V/m)

    Since we need v0/(2*c) as a dimensionless parameter create

        BETA0 = (v0/2/c)^2

    for use in magnetic field calculations.
*/

#define BETA0  (0.05)

/*  to follow electron blobs, save position, velocity and fields that act on it.
    valid bit tells me if the particle has hit the wall - the only purpose is to
    avoid garbage collection activity so this could be used to save memory.
*/

typedef struct
{
    int valid;
    double r, thta, phi;
    double ux, uy, uz;
    double ex, ey, ez;
    double bx, by, bz;
}  Phase_Space;

/*  Create global buffers for data storage so integrals are easy */

double GridEfld[3*NUMELMNTS];
double SrcEfld[3*NUMELMNTS];
double MagFld[3*NUMELMNTS];
double ElcFld[3*NUMELMNTS];
int numangles[NUMELMNTS];
int GridCurrent[MAXTIMESTEPS];
int WallCurrent[MAXTIMESTEPS];
int SrcCurrent[MAXTIMESTEPS];
Phase_Space eblob[MAXBLOBS];

/*  subroutine to copy phase space elements from one place to another.
    Enter with pointers to source and destination.  Does not need to
    copy field values (yet).
*/

void copy_phase( Phase_Space *src, Phase_Space *dst)
{
    dst->valid = src->valid;
    dst->r = src->r;
    dst->thta = src->thta;
    dst->phi = src->phi;
    dst->ux = src->ux;
    dst->uy = src->uy;
    dst->uz = src->uz;
}

/*  Subroutine to find mirror points to phase space blobs.
    Enter with pointer to reflection plane normal vector, number of blobs
    to reflect and pointer to list of blobs.  Extends list in place and
    returns total # of blobs.
    Need room for 48 blobs total.
*/

int mirror_phase( double *normal, int nblb, Phase_Space *phspnts)
{
    int i, xtnd;
    double adotp, x, y, z, xprm, yprm, zprm, phick;

    xtnd = nblb;
    for(i=0; i<nblb; i++)
    {
	x = phspnts[i].r*cos(phspnts[i].thta);
	y = x*sin(phspnts[i].phi);
	x = x*cos(phspnts[i].phi);
	z = phspnts[i].r*sin(phspnts[i].thta);
	adotp = normal[0]*x + normal[1]*y + normal[2]*z;
	adotp *= 2.0;
	xprm = x - normal[0]*adotp;
	yprm = y - normal[1]*adotp;
	zprm = z - normal[2]*adotp;
	phspnts[xtnd].r = phspnts[i].r;
	phick = atan2(yprm, xprm);
	if (phick < 0.0) phick += 2.0*M_PI;
	phspnts[xtnd].phi = phick;
	phspnts[xtnd].thta = atan(zprm/sqrt(xprm*xprm + yprm*yprm));
	adotp = normal[0]*phspnts[i].ux;
	adotp += normal[1]*phspnts[i].uy;
	adotp += normal[2]*phspnts[i].uz;
	adotp *= 2.0;
	phspnts[xtnd].ux -= adotp*normal[0];
	phspnts[xtnd].uy -= adotp*normal[1];
	phspnts[xtnd].uz -= adotp*normal[2];
	phspnts[xtnd].valid = 1;
	xtnd++;
    }
    return xtnd;
}

/*  create 47 copies of one phase space point.
    Assumes storage space predefined.
*/

void sym_phase( double m[6][3], Phase_Space *src, Phase_Space *dst)
{
/*  do first two points independently  */

    copy_phase( src, dst);
    mirror_phase( m[0], 1, dst);
    copy_phase( &dst[1], &dst[2]);
    mirror_phase( m[1], 1, dst);

/*  copy three points into 3 more  */

    mirror_phase(m[2], 3, dst);

/*  next copy across x, y and z planes  */

    mirror_phase(m[3], 6, dst);
    mirror_phase(m[4], 12, dst);
    mirror_phase(m[5], 24, dst);
}

/*  Input a given i, j, k index (assumes r, theta, phi association),
    Returns 8 neighbor indecies.  Assumes maximum index is input.
*/

void neighborhood( int i, int j, int k, int *dexptr)
{
    int idex, kdex, jdex;
    int a, b, c;

    for(a=0; a<2; a++)
    {
	idex = (i-a)*(i-a+1)*(i-a+2)/6;
	for(b=0; b<2; b++)
	{
	    kdex = (k-b)*(k-b+1)/2;
	    for(c=0; c<2; c++)
	    {
		jdex = j-c;

/*  check to see if outside 1/48 grid.  Can only happen on "top side" cut  */

		if( jdex > k-b)
		    *dexptr = (idex + k*(k+1)/2 + j - 1)*3;
		else
		    *dexptr = (idex + kdex + jdex)*3;
		dexptr++;
	    }
	}
    }
}

/*  compute electric field generated by electron source on coil face.  */

main()
{
    FILE *readmag, *particles;
    int numread, dexptrs[8], numblobs, time;
    char filename[64], input[8];
    double magpot_amp;
    double across, r, phi, theta;
    int err, numval, i, j, k, idex, jdex, dex;
    double pi4, alpha, beta;
    double rprime, thetaprime, phiprime;
    double temp1, temp2, t, dist, dtr, dmr;
    int ip, jp, kp, npts, np, edex, kdex, passnum, abdex;
    double dsum, sum, dr, temp, acp;

    int srcnum, blobdex;
    double mirror[6][3], srcangle, coildistance, tparamtop, tparambot;
    double emag, cp, ctm, ctp, xprime, yprime, zprime, px, py, pz;
    double cpsi1, cpsi2, spsi1, spsi2, psicntr, psi[SOURCEROWS+2];
    double srcntr[(SOURCEROWS+1)*(SOURCEROWS+2)], rw, escale, bscale;
    double r0, r1, t00, t01, t10, t11, cp00, cp01, cp10, cp11;
    double ct00, ct01, ct10, ct11, st00, st01, st10, st11;
    double ex, ey, ez, bx, by, bz, w[8], r20, r21, rz;
    double dux, duy, duz, xb, yb, zb, xe, ye, ze;
    Phase_Space esrc[SOURCE_ELEMENTS], cpyblb[48];

/*  load MaGrid electric field from previous calculation  */

    sprintf(filename, "e_field_25.dat");
    readmag = fopen(filename, "rb");
    if( !readmag)
    {
	printf("can't read in potential data file %s.\n", filename);
	exit(0);
    }
    printf("opening potential data files and reading them in.\n");
    numread = fread(GridEfld, sizeof(double), 3*NUMELMNTS, readmag);
    if(numread < 3*NUMELMNTS)
	printf("\nWARNING: not all field file read in. \n\n");
    fclose(readmag);

/*  load MaGrid magnetic field from previous calculation  */

    sprintf(filename, "b_field_25.dat");
    readmag = fopen(filename, "rb");
    if( !readmag)
    {
	printf("can't read in magnetic potential data file %s.\n", filename);
	exit(0);
    }
    printf("opening magnetic potential data files and reading them in.\n");
    numread = fread(MagFld, sizeof(double), 3*NUMELMNTS, readmag);
    if(numread < 3*NUMELMNTS)
	printf("\nWARNING: not all field file read in. \n\n");
    fclose(readmag);

/*  Compute constants for main calculations  */

    pi4 = M_PI/4.0;          // ubiquitous constant
    srcangle = atan(SOURCERADIUS / SOURCEPOSITION);
    coildistance = sqrt(COILRADIUS*COILRADIUS + 1.0);
    tparamtop = MAXWALL - SOURCERADIUS;
    tparambot = SOURCEPOSITION*SOURCEPOSITION - tparamtop*tparamtop;
    emag = SOURCEVOLTAGE*SOURCERADIUS/tparamtop;
    emag *= MAXWALL;
    tparamtop *= SOURCERADIUS;

/*  compute electric field generated by electron source balls  */

    SrcEfld[0] = 0.0;
    SrcEfld[1] = 0.0;
    SrcEfld[2] = 0.0;
    dex = 3;
    for(i=1; i<MAXSTEPS; i++)
    {
	r = MAXWALL/MAXSTEPS*i;
	across = pi4/i;
	for(k=0; k<=i; k++)
	{
	    phi = across*k;
	    for(j=0; j<=k; j++)
	    {
		theta = across*j;

/* find cartesian equivelent position from 2D description */

		xb = r*cos(theta);
		yb = r*sin(theta) - SOURCEPOSITION;
		rz = sqrt(yb*yb + xb*xb);
		temp1 = tparamtop - SOURCEPOSITION*yb;
		temp2 = rz*rz - SOURCERADIUS*SOURCERADIUS;
		temp2 *= tparambot;
		t = (temp1 + sqrt(temp1*temp1 - temp2))/tparambot;
		rw = SOURCERADIUS + (MAXWALL - SOURCERADIUS)*t;

/*  find true 3D cartesian position for vector direction of E field  */

		xe = xb*cos(phi);
		ye = xb*sin(phi);
		ze = yb + SOURCEPOSITION;
		temp1 = emag/rz/rw/rw;
		SrcEfld[dex] = xe*temp1;
		dex++;
		SrcEfld[dex] = ye*temp1;
		dex++;
		SrcEfld[dex] = ze*temp1;
		dex++;
	    }
	}
    }

/*  combine MaGrid E field with electron source field  */

    dex = 0;
    for(i=0; i<MAXSTEPS; i++)
    {
	for(k=0; k<=i; k++)
	{
	    for(j=0; j<=k; j++)
	    {
		for(jdex=0; jdex<3; jdex++)
		{
		    ElcFld[dex] = GridEfld[dex] + SrcEfld[dex];
		    dex++;
		}
	    }
	}
    }

/*  Build table of centroids for source particles and their escape vectors  */

    srcnum = SOURCE_ELEMENTS*2;
    for(j=0; j<=SOURCEROWS; j++)
    {
	cpsi1 = 1.0 - (double)(1 + j)*(double)j/(double)srcnum;
	psi[j] = acos(cpsi1);
    }
    jdex = 0;
    for(j=1; j<=SOURCEROWS; j++)
    {
	cpsi1 = cos(psi[j-1]);
	cpsi2 = cos(psi[j]);
	spsi1 = sin(psi[j-1]);
	spsi2 = sin(psi[j]);
	psicntr = spsi2 - spsi1 + psi[j-1]*cpsi1 - psi[j]*cpsi2;
	psicntr /= cpsi1 - cpsi2;
	printf(" j = %d psi centroid is %lf \n", j, psicntr);
	for(k=0; k<j; k++)
	{
	    srcntr[jdex] = psicntr;
	    jdex++;
	    srcntr[jdex] = M_PI/8.0/j * (2.0*k + 1.0);
	    printf("%lf  ", srcntr[jdex]);
	    jdex++;
	}
	printf(" %d\n", jdex);
    }

/*  convert centroids to face coil positions  */

    rw = SOURCEPOSITION*SOURCEPOSITION + SOURCERADIUS*SOURCERADIUS;
    for(j=0; j<SOURCE_ELEMENTS; j++)
    {
	k = 2*j;
	cpsi1 = cos(srcntr[k+1]);  // cos(xi)
	cpsi2 = cos(srcntr[k]);    // cos(psi)
	spsi1 = sin(srcntr[k+1]);  // sin(xi)
	spsi2 = sin(srcntr[k]);    // sin(psi)
	temp1 = SOURCERADIUS*spsi2*cpsi1;
	temp2 = SOURCEPOSITION - SOURCERADIUS*cpsi2;
	esrc[j].phi = atan(temp1/temp2);  // phi
	temp1 = SOURCERADIUS*spsi1*spsi2;
	temp2 = rw - 2.0*SOURCERADIUS*SOURCEPOSITION*cpsi2;
	esrc[j].r = sqrt(temp2);      // radius
	temp2 = sqrt(temp2 - temp1*temp1);
	esrc[j].thta = atan(temp1/temp2);  // theta
	esrc[j].ux = - SOURCERADIUS*cpsi2;
	esrc[j].uy = SOURCERADIUS*spsi2*cpsi1;
	esrc[j].uz = SOURCERADIUS*spsi2*spsi1;
	esrc[j].valid = 1;
    printf("source %d data: r= %lf  theta = %lf  phi = %lf xi= %lf  psi = %lf\n", 
	   j, esrc[j].r, esrc[j].thta, esrc[j].phi, srcntr[k+1], srcntr[k]);
    printf("     ux = %lf  uy = %lf  uz = %lf\n", esrc[j].ux, esrc[j].uy, esrc[j].uz);
    }

    numblobs = 0;  // no source data yet

/*  initialize mirror reflection normal vectors */

    mirror[0][0] = 1.0/sqrt(2.0);
    mirror[0][1] = -mirror[0][0];
    mirror[0][2] = 0.0;
    mirror[1][0] = 0.0;
    mirror[1][1] = mirror[0][0];
    mirror[1][2] = mirror[0][1];
    mirror[2][0] = mirror[0][0];
    mirror[2][1] = 0.0;
    mirror[2][2] = mirror[0][1];
    mirror[3][0] = 0.0;
    mirror[3][1] = 0.0;
    mirror[3][2] = 1.0;
    mirror[4][0] = 1.0;
    mirror[4][1] = 0.0;
    mirror[4][2] = 0.0;
    mirror[5][0] = 0.0;
    mirror[5][1] = 1.0;
    mirror[5][2] = 0.0;

/*  initialize scale values based on input constants, for face source  */

    escale = SOURCEVOLTAGE*MAXWALL*SOURCERADIUS;
    escale /= 4*(MAXWALL - SOURCERADIUS);
    bscale = escale*CONST_CP*BETA0;
    escale /= SOURCE_ELEMENTS*2.0;

/*  begin loop over time stepping.
    First sweep across phase space and remove all particles striking surfaces.
*/

    for(time=0; time<MAXTIMESTEPS; time++)
    {
	printf("begin time step %d\n", time);
	rw = SOURCEPOSITION*SOURCEPOSITION;
	for(idex=0; idex<numblobs; idex++)
	{
	    if(eblob[idex].valid)
	    {
		if(eblob[idex].r > MAXWALL)
		{
		    WallCurrent[time]++;
		    eblob[idex].valid = 0;
		    printf("blob %d has hit wall at %lf\n", idex,eblob[idex].r); 
		    continue;
		}
		temp1 = eblob[idex].r - coildistance;
		if(temp1<0.0) temp1 = -temp1;
		if(temp1 <= COILDIAMETER)
		{
		    GridCurrent[time]++;
		    eblob[idex].valid = 0;
		    continue;
		}
		if((eblob[idex].thta <= srcangle) && (eblob[idex].phi <= srcangle))
		{
		    temp1 = 2.0*SOURCEPOSITION*eblob[idex].r;
		    temp1 *= cos(eblob[idex].thta)*cos(eblob[idex].phi);
		    temp1 = sqrt(rw + eblob[idex].r*eblob[idex].r - temp1);
		    if(temp1 <= SOURCERADIUS)
		    {
			SrcCurrent[time]++;
			eblob[idex].valid = 0;
		    }
		}
	    }
	}
	
/*  create new particles at sources if allowed  */

	if(time < EMISSIONTIME)
	{
	printf("begin emission\n");

	    for(j=0; j<srcnum/2; j++)
	    {
		copy_phase(&esrc[j], &eblob[numblobs]);
		numblobs++;
	    }
	}

	printf("numblobs = %d\n", numblobs);

/*  compute fields generated by blobs. Clear all fields for every blob. */

	for(blobdex=0; blobdex<numblobs; blobdex++)
	{
	    if(eblob[blobdex].valid)
	    {
		eblob[blobdex].ex = 0.0;
		eblob[blobdex].ey = 0.0;
		eblob[blobdex].ez = 0.0;
		eblob[blobdex].bx = 0.0;
		eblob[blobdex].by = 0.0;
		eblob[blobdex].bz = 0.0;
	    }
	}

/*    First generate 48 copies of each blob in storage.  */

	for(blobdex = 0; blobdex<numblobs; blobdex++)
	{
	    if(!eblob[blobdex].valid) continue;
	    sym_phase(mirror, &eblob[blobdex], cpyblb);
/*	    for(dex=0; dex<48; dex++)
	    {
		printf("r = %lf  phi = %lf  thta = %lf\n",
		       cpyblb[dex].r, cpyblb[dex].phi, cpyblb[dex].thta); 
		printf("ux = %lf  uy = %lf  uz = %lf\n",
		       cpyblb[dex].ux, cpyblb[dex].uy, cpyblb[dex].uz);
		printf("ex = %lf  ey = %lf  ez = %lf\n",
		       cpyblb[dex].ex, cpyblb[dex].ey, cpyblb[dex].ez); 
		printf("\n");
	    }
*/
	    for(dex=0; dex<numblobs; dex++)
	    {
		r = eblob[dex].r;
		r20 = r*r;
		phi = eblob[dex].phi;
		theta = eblob[dex].thta;
		xb = r*cos(theta);
		yb = xb*sin(phi);
		xb *= cos(phi);
		zb = r*sin(theta);
		for(j=0; j<48; j++)
		{

/*  compute distance between this point and all other blobs  */

		    if(!j && (dex == blobdex)) continue;  // don't do yourself!
		    cp = cos(phi - cpyblb[j].phi);
		    ctm = cos(theta - cpyblb[j].thta);
		    ctp = cos(theta + cpyblb[j].thta);
		    dtr = r*(ctm*(1.0 + cp) - ctp*(1.0 - cp));
		    dtr = cpyblb[j].r*(cpyblb[j].r - dtr) + r20;
		    dist = pow(dtr, 1.5);

/*  convert positions to cartesian form to get fields from distance  */

		    xprime = cpyblb[j].r*cos(cpyblb[j].thta);
		    yprime = xprime*sin(cpyblb[j].phi);
		    xprime *= cos(cpyblb[j].phi);
		    zprime = cpyblb[j].r*sin(cpyblb[j].thta);
		    px = (xb - xprime)/dist;
		    py = (yb - yprime)/dist;
		    pz = (zb - zprime)/dist;
		    eblob[dex].ex += px;
		    eblob[dex].ey += py;
		    eblob[dex].ez += pz;
		    eblob[dex].bx += cpyblb[j].uy*pz - cpyblb[j].uz*py;
		    eblob[dex].by += cpyblb[j].uz*px - cpyblb[j].ux*pz;
		    eblob[dex].bz += cpyblb[j].ux*py - cpyblb[j].uy*px;
/*		    printf("j=%d\n", j);
	    printf("r = %lf  phi = %lf  thta = %lf\n",
		   eblob[dex].r, eblob[dex].phi, eblob[dex].thta); 
	    printf("ux = %lf  uy = %lf  uz = %lf\n",
		   eblob[dex].ux, eblob[dex].uy, eblob[dex].uz);
	    printf("ex = %lf  ey = %lf  ez = %lf\n",
		   eblob[dex].ex, eblob[dex].ey, eblob[dex].ez); 
	    printf("\n");
	    if(j>2) exit(0);
*/
		}
	    }
	}

/*  scale all E and B fields according to dimensionless formula  */

	printf("begin scaling\n\n\n");

	for(dex=0; dex<numblobs; dex++)
	{
	    if(!eblob[dex].valid) continue;
	    eblob[dex].ex *= escale;
	    eblob[dex].ey *= escale;
	    eblob[dex].ez *= escale;
	    eblob[dex].bx *= bscale;
	    eblob[dex].by *= bscale;
	    eblob[dex].bz *= bscale;
	}

/*  compute change in velocity of each blob and move it.
    for each particle, estimate electric and magnetic field in its vicinity.
*/

	for(blobdex=0; blobdex<numblobs; blobdex++)
	{
	    if(!eblob[blobdex].valid) continue;
/*	    printf("r = %lf  phi = %lf  thta = %lf\n",
		   eblob[blobdex].r, eblob[blobdex].phi, eblob[blobdex].thta); 
	    printf("ux = %lf  uy = %lf  uz = %lf\n",
		   eblob[blobdex].ux, eblob[blobdex].uy, eblob[blobdex].uz);
	    printf("ex = %lf  ey = %lf  ez = %lf\n",
		   eblob[blobdex].ex, eblob[blobdex].ey, eblob[blobdex].ez); 
	    printf("bx = %lf  by = %lf  bz = %lf\n",
		   eblob[blobdex].bx, eblob[blobdex].by, eblob[blobdex].bz); 
	    printf("\n");
*/
	    xe = eblob[blobdex].r/MAXWALL*MAXSTEPS;
	    i = trunc(xe);
	    dex = trunc(xe + 0.5);
	    if( i==dex) dex++;
	    k = dex/pi4*eblob[blobdex].phi + 1;
	    if( k>dex) k = dex;
	    j = dex/pi4*eblob[blobdex].thta + 1;
	    if( j>k) j = k;

/*  compute indicies and weights to neighbors of this point  */

	    r0 = dex*MAXWALL/MAXSTEPS;
	    r1 = (dex-1)*MAXWALL/MAXSTEPS;
	    t00 = pi4/dex*j;
	    t01 = pi4/dex*(j-1);
	    t10 = pi4/(dex-1)*j;
	    t11 = pi4/(dex-1)*(j-1);
	    ct00 = cos(t00);
	    st00 = sin(t00);
	    ct01 = cos(t01);
	    st01 = sin(t01);
	    ct10 = cos(t10);
	    st10 = sin(t10);
	    ct11 = cos(t11);
	    st11 = sin(t11);
	    neighborhood(dex, j, k, dexptrs);
	    r = eblob[blobdex].r;
	    r20 = r*r;
	    r21 = r20 + r1*r1;
	    r20 += r0*r0;
	    temp1 = cos(eblob[blobdex].thta);
	    temp2 = sin(eblob[blobdex].thta);
	    t = eblob[blobdex].phi;
	    cp00 = cos(t - pi4/dex*k);
	    cp01 = cos(t - pi4/dex*(k-1));
	    cp10 = cos(t - pi4/(dex-1)*k);
	    cp11 = cos(t - pi4/(dex-1)*(k-1));
	    w[0] = r20 - r*r0*(temp1*ct00*cp00 + temp2*st00);
	    w[1] = r20 - r*r0*(temp1*ct00*cp01 + temp2*st00);
	    w[2] = r20 - r*r0*(temp1*ct01*cp00 + temp2*st01);
	    w[3] = r20 - r*r0*(temp1*ct01*cp01 + temp2*st01);
	    w[4] = r21 - r*r1*(temp1*ct10*cp10 + temp2*st10);
	    w[5] = r21 - r*r1*(temp1*ct10*cp11 + temp2*st10);
	    w[6] = r21 - r*r1*(temp1*ct11*cp10 + temp2*st11);
	    w[7] = r21 - r*r1*(temp1*ct11*cp11 + temp2*st11);
	    sum = 0.0;
	    for(i=0; i<8; i++) sum += w[i];
	    ex = 0.0;
	    ey = 0.0;
	    ez = 0.0;
	    bx = 0.0;
	    by = 0.0;
	    bz = 0.0;
	    for(i=0; i<8; i++)
	    {
		ex += ElcFld[dexptrs[i]]/w[i];
		ey += ElcFld[dexptrs[i]+1]/w[i];
		ez += ElcFld[dexptrs[i]+2]/w[i];
		bx += MagFld[dexptrs[i]]/w[i];
		by += MagFld[dexptrs[i]+1]/w[i];
		bz += MagFld[dexptrs[i]+2]/w[i];
	    }
	    ex *= sum;
	    ey *= sum;
	    ez *= sum;
	    bx *= sum;
	    by *= sum;
	    bz *= sum;

/*  compute change in velocity of this particle based on electric and magnetic
    field it is exposed to.
*/

	    dux = eblob[blobdex].uy*(bz + eblob[blobdex].bz);
	    dux -= eblob[blobdex].uz*(by + eblob[blobdex].by);
	    dux /= CONST_CP;
	    dux += ex + eblob[blobdex].ex;
	    dux *= escale*DELTANU;
	    duy = eblob[blobdex].uz*(bx + eblob[blobdex].bx);
	    duy -= eblob[blobdex].ux*(bz + eblob[blobdex].bz);
	    duy /= CONST_CP;
	    duy += ey + eblob[blobdex].ey;
	    duy *= escale*DELTANU;
	    duz = eblob[blobdex].ux*(by + eblob[blobdex].by);
	    duz -= eblob[blobdex].uy*(bx + eblob[blobdex].bx);
	    duz /= CONST_CP;
	    duz += ez + eblob[blobdex].ez;
	    duz *= escale*DELTANU;

/*  move particle to new position  */

	    eblob[blobdex].ux += dux;
	    eblob[blobdex].uy += duy;
	    eblob[blobdex].uz += duz;
	    px = eblob[blobdex].r*cos(eblob[blobdex].thta);
	    py = px*sin(eblob[blobdex].phi);
	    px *= cos(eblob[blobdex].phi);
	    pz = eblob[blobdex].r*sin(eblob[blobdex].thta);
	    px += eblob[blobdex].ux*DELTANU;
	    py += eblob[blobdex].uy*DELTANU;
	    pz += eblob[blobdex].uz*DELTANU;
	    r20 = px*px + py*py;
	    eblob[blobdex].r = sqrt(r20 + pz*pz);
	    phi = atan2(py, px);
	    if(phi < 0.0) eblob[blobdex].phi = 2.0*M_PI + phi;
	    else eblob[blobdex].phi = phi;
	    eblob[blobdex].thta = atan(pz/sqrt(r20));

/*  check to see if particle is outsize zone zero - swap it back if it is.
    Not exactly efficient, but goes through all possible combinations.  */

	    if(eblob[blobdex].thta < 0.0)
	    {
		eblob[blobdex].thta = -eblob[blobdex].thta;
		eblob[blobdex].uz = -eblob[blobdex].uz;
	    }
	    if(eblob[blobdex].phi > M_PI)
	    {
		eblob[blobdex].phi = 2.0*M_PI - eblob[blobdex].phi;
		eblob[blobdex].uy = -eblob[blobdex].uy;
	    }
	    if(eblob[blobdex].phi > M_PI/2.0)
	    {
		eblob[blobdex].phi = M_PI - eblob[blobdex].phi;
		eblob[blobdex].ux = -eblob[blobdex].ux;
	    }
	    if(eblob[blobdex].phi < eblob[blobdex].thta)
	    {
		temp1 =  eblob[blobdex].phi;
		eblob[blobdex].phi = eblob[blobdex].thta;
		eblob[blobdex].thta = temp1;
		temp1 = eblob[blobdex].uz;
		eblob[blobdex].uz = eblob[blobdex].uy;
		eblob[blobdex].uy = temp1;
	    }
	    if(eblob[blobdex].phi > M_PI/4.0)
	    {
		if(eblob[blobdex].thta > M_PI/2.0 - eblob[blobdex].phi)
		{
		    temp1 = eblob[blobdex].thta;
		    eblob[blobdex].thta = M_PI/2.0 - eblob[blobdex].phi;
		    eblob[blobdex].phi = M_PI/2.0 - temp1;
		    temp1 = eblob[blobdex].uz;
		    eblob[blobdex].uz = eblob[blobdex].ux;
		    eblob[blobdex].ux = temp1;
		}
		eblob[blobdex].phi = M_PI/2.0 - eblob[blobdex].phi;
		temp1 = eblob[blobdex].uy;
		eblob[blobdex].uy = eblob[blobdex].ux;
		eblob[blobdex].ux = temp1;
		
	    }
/*
	    printf("for particle %d we have:\n", blobdex);
	    printf("r = %lf  phi = %lf  thta = %lf\n",
		   eblob[blobdex].r, eblob[blobdex].phi, eblob[blobdex].thta); 
	    printf("ux = %lf  uy = %lf  uz = %lf\n",
		   eblob[blobdex].ux, eblob[blobdex].uy, eblob[blobdex].uz);
	    printf("ex = %lf  ey = %lf  ez = %lf\n",
		   eblob[blobdex].ex, eblob[blobdex].ey, eblob[blobdex].ez); 
	    printf("bx = %lf  by = %lf  bz = %lf\n",
		   eblob[blobdex].bx, eblob[blobdex].by, eblob[blobdex].bz); 
	    printf("\n");
*/
	}// blobdex

/*  save particle distibution to disk for later analysis  */
	printf("save particles to disk\n");

	sprintf(filename, "particles_%03d.dat", time);
	particles = fopen(filename, "w");
	if(!particles)
	{
	    printf("can't create  %s! \n", filename);
	    exit(0);
	}
	fwrite(&time, sizeof(int), 1, particles);
	fwrite(&numblobs, sizeof(int), 1, particles);
	fwrite(WallCurrent, sizeof(int), time+1, particles);
	fwrite(GridCurrent, sizeof(int), time+1, particles);
	fwrite(SrcCurrent, sizeof(int), time+1, particles);
	fwrite(eblob, sizeof(Phase_Space), numblobs, particles);
	fclose(particles);
    }// time
    printf("done!\n");
}