/*******************************************************************
 *                                                                 *
 *      The purpose of this program is to compute single particle  *
 *  orbits in the polywell structure.  There are several ideas     *
 *  incorporated:  The equations are all unitless so scaling can   *
 *  can be applied later.  The symmetry of the polywell was shown  *
 *  by Indrek on talk-polywell.org to be 1/48, not just 1/8.  This *
 *  allows a much finer grid as well as the ability to store data  *
 *  in ram even on a 1GB machine.  While the equations of motion   *
 *  are computed in cartesian coordinates, the E and B fields are  *
 *  stored in spherical coordinates. They are unusual coordinates  *
 *  with "theta" being in the x,y plane and "phi" going up from    *
 *  the x,y plane to the point at z.  The zeroth quadrant of this  *
 *  spherical system holds all the E and B data, and Indrek's      *
 *  vector form is used to find the lookup point within that block.*
 *                                                                 *
 *     Very simple particle tracking is done using dimensionless   *
 *  steps in time and space and velocity.  Two 3D equations are    *
 *  used to compute the change in velocity and then the change in  *
 *  position based on that velocity.  Use some kind of Monte-Carlo *
 *  method to explore the full grid.                               *
 *                                                                 *
 *                        Author = Mike Rosing                     *
 *                         date = Jan. 1, 2008                     *
 *                                                                 *
 ******************************************************************/

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

#define MAXDIM 400

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

#define MAXWALL  (3.0)

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

/*  globals set in main only once  */

double dr;
double bfield[NUMELMNTS], efield[NUMELMNTS];

/*  subroutine to return offset into E or B field buffers.
    Enter with x, y, z position vector.
    Returns with offset to start of Fx, Fy, Fz data or -1 if
      input point is out of bounds.
*/

int bufroffset(double x, double y, double z, int *mapping)
{
    unsigned long mask1, mask2, mask3, mask;
    int offset, i, j, k;
    double r, theta, phi, xm, ym, zm;
    double dtheta, dphi;

/*  First map all octants back to zeroth octant  */

    if (x < 0) 
    {
	x = -x;
	mapping[3] = -1;
    }
    else mapping[3] = 1;
    if (y < 0) 
    {
	y = -y;
	mapping[4] = -1;
    }
    else mapping[4] = 1;
    if (z < 0) 
    {
	z = -z;
	mapping[5] = -1;
    }
    else mapping[5] = 1;

/*  find sextant within base octant.  Each bit position is 
    for a particular sextant  */

    if (x < y) mask1 = 0xE;
    else mask1 = 0x31;
    if (x < z) mask2 = 0x1C;
    else mask2 = 0x23;
    if( y < z) mask3 = 0x38;
    else mask3 = 0x7;
    mask = mask1 & mask2 & mask3;
    if( !mask)
    {
	printf("Error in bufroffset: x = %lf  y = %lf  z = %lf\n",
	       x, y, z);
	exit(0);
    }
    if( mask & 1)
    {
	xm = x;
	ym = y;
	zm = z;
	mapping[0] = 0;
	mapping[1] = 1;
	mapping[2] = 2;
    }
    else if (mask & 2)
    {
	xm = y;
	ym = x;
	zm = z;
	mapping[0] = 1;
	mapping[1] = 0;
	mapping[2] = 2;
    }
    else if (mask & 4)
    {
	xm = y;
	ym = z;
	zm = x;
	mapping[0] = 1;
	mapping[1] = 2;
	mapping[2] = 0;
    }
    else if (mask & 8)
    {
	xm = z;
	ym = y;
	zm = x;
	mapping[0] = 2;
	mapping[1] = 1;
	mapping[2] = 0;
    }
    else if (mask & 0x10)
    {
	xm = z;
	ym = x;
	zm = y;
	mapping[0] = 2;
	mapping[1] = 0;
	mapping[2] = 1;
    }
    else if (mask & 0x20)
    {
	xm = x;
	ym = z;
	zm = y;
	mapping[0] = 0;
	mapping[1] = 2;
	mapping[2] = 1;
    }
    r = sqrt(xm*xm + ym*ym + zm*zm);
    theta = atan2(ym, xm);
    phi = atan2(zm,  sqrt(xm*xm + ym*ym));
    i = lrint(r/dr);
    if (i > MAXDIM) return -1;
    dtheta = 4.0*i/M_PI;
    j = lrint(theta*dtheta);
    if (j > MAXDIM) return -1;
    dphi = dtheta;
    k = lrint(phi*dphi);
    if (k > MAXDIM) return -1;
    offset = (i*(i+1)*(i+2)/6 + j*(j+1)/2 + k)*3;
    return offset;
}

main()
{
    FILE  *bfile, *efile, *sqgle, *vggle;
    double x, y, z, ux, uy, uz;
    int numread, index, i, map[6];
    double deltanu, param, volts, amps;
    double uxp, uyp, uzp, me, ce;
    double bx, by, bz;
    double magridfreq, nu, gridamp;

    me = 9.10938189e-31;    // mass electron
    ce = 1.602176462e-19;   // charge of electron

    printf("opening data files and reading them in.\n");
    efile = fopen("e_field.dat", "rb");
    if(!efile)
    {
	printf("can't open e_field.dat file.\n");
	exit(0);
    }
    bfile = fopen("b_field.dat", "rb");
    if(!bfile)
    {
	printf("can't open b_field.dat file.\n");
	exit(0);
    }
    numread = fread(bfield, sizeof(double), NUMELMNTS, bfile);
    if(numread < NUMELMNTS)
	printf("\nWARNING: not all B field file read in. %d  \n\n", 
	       numread);
    printf("\n");
    numread = fread(efield, sizeof(double), NUMELMNTS, efile);
    if(numread < NUMELMNTS)
	printf("\nWARNING: not all E field file read in. \n\n");
    fclose(efile);
    fclose(bfile);
    printf("Voltage on grid: ");
    scanf ("%lf", &volts);
    printf("Amp turns in grid: ");
    scanf("%lf", &amps);
    param = sqrt(2.0*me*volts/ce)/amps/(2.0e-7);
    printf("containment parameter = %lf\n", param);
    dr = MAXWALL/MAXDIM;

/*  pick starting conditions  */

    deltanu = 0.005;
    x = 0.0;
    y = 0.0;
    z = 0.0;
    ux = 0.001;
    uy = 0.001;
    uz = 0.001;
    magridfreq = 5.0;

    sqgle = fopen("squiggle.dat", "w");
    if(!sqgle)
    {
	printf("can't create squiggle data file\n");
	exit(0);
    }
    vggle = fopen("vwiggle.dat", "w");
    if(!vggle)
    {
	printf("can't create vwiggle data file\n");
	exit(0);
    }
    nu = 0.0;
    for(i=0; i<MAXSTEPS; i++)
    {
	index = bufroffset( x, y, z, map);
	if( index < 0)
	{
	    printf("Hit wall.\n");
	    printf("%lg  %lg %lg\n", x, y, z);
	    printf("%lg  %lg %lg\n", ux, uy, uz);
	    break;
	}
	x += deltanu * ux;
	y += deltanu * uy;
	z += deltanu * uz;
	bx = map[3]*bfield[index+map[0]];
	by = map[4]*bfield[index+map[1]];
	bz = map[5]*bfield[index+map[2]];
	uxp = ux;
	uyp = uy;
	uzp = uz;
	gridamp = param*cos(magridfreq*nu);
	ux -= deltanu*(gridamp*map[3]*efield[index+map[0]] + uyp*bz - uzp*by);
	uy -= deltanu*(gridamp*map[4]*efield[index+map[1]] + uzp*bx - uxp*bz);
	uz -= deltanu*(gridamp*map[5]*efield[index+map[2]] + uxp*by - uyp*bx);
	nu += deltanu;
	fprintf(sqgle, "%lg %lg %lg\n", x, y, z);
	fprintf(vggle, "%lg %lg %lg\n", ux, uy, uz);
    }
    printf("went for %d steps\n", i);
    fclose(sqgle);
    fclose(vggle);
}