/********************************************************************
 *                                                                  *
 *      Compute potential from a Maxwellian electron density        *
 *  function.                                                       *
 *  Experimentally based approximate electron distribution in       *
 *  polywell.  See notes on "Electron Fluid in a Polywell take 2"   *
 *                                                                  *
 *       density ~ 1/(exp(-a*r*r) + b*r*r)                          *
 *                                                                  *
 *              Author = Mike Rosing                                *
 *                date = Feb 22, 2008                               *
 *                                                                  *
 *******************************************************************/

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

#define MAXSTEPS 25

/*  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)
#define NUMELMNTS ((MAXSTEPS+1)*(MAXSTEPS+2)*(MAXSTEPS+3)/6)

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

double GridPot[NUMELMNTS];
double ElecPot0[NUMELMNTS], ElecPot1[NUMELMNTS];
double MagPot[3*NUMELMNTS];
double energy[NUMELMNTS];
double angles[2*48*NUMELMNTS];
int numangles[NUMELMNTS];

/* ELECTRON_TEMP is the temperature to MaGrid voltage ratio.  It is 
   really an energy ratio defined by

              ELECTRON_TEMP*e*V = k*T  

where V is the MaGrid voltage, k is Boltzmann's constant and T is the
electron temperature.
*/

#define ELECTRON_TEMP  0.3

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

/*  The density distribution is arbitrarily defined to start out with

         density(r) = 1/exp(-alpha*r^2) + beta*r^2

Call the first coefficient DENSITY_A and the second DENSITY_B.  These 
values are chosen to give a flat center and close to zero near the walls.
*/

#define DENSITY_A  9.0
#define DENSITY_B  10.0

/*  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.05)

/*  Find all symmetry points in every sextant and octant to a given point.
    The sextant split line is along theta=phi = PI/4, and the octants are
    split along the (x,y), (x,z) and (y,z) planes.  All points are at same
    radius so only two angles are needed for each point.
    Input is phi, theta in 0,0 (octant, sextant)
    Output is vector list with return value of number of points in the list.
    First entry includes input value followed by N-1 points in phi, theta
    order.
*/

int symmetry( double *source, double *vector)
{
    double phi, theta, pi2, phicheck, thcheck, vec[96];
    int i, j, k, n, check;

    pi2 = M_PI/2.0;
    phi = source[0];
    theta = source[1];
    for(i=0; i<6; i++)  // outer loop over sextants
    {
	k = 16*i;
	for(j=0; j<16; j+=2)  // inner loop over octants
	{
	    vec[k+j] = phi;
	    vec[k+j+1] = theta;
	    if(j==6 || j==14) theta = - theta;
	    phi += pi2;
	    if(phi > 2.0*M_PI) phi -= 2.0*M_PI;
	    phicheck = phi - 2.0*M_PI;
	    if(phicheck < 0.0) phicheck = -phicheck;
	    if(phicheck < ANGLEPS) phi = 0.0;
	}
	switch (i)
	{
	    case 0:
		phi = pi2 - source[0];
		break;
	    case 1:
		phi = pi2 - source[1];
		theta = source[0];
		break;
	    case 2:
		phi = pi2 - source[1];
		theta = pi2 - source[0];
		break;
	    case 3:
		phi = source[1];
		theta = pi2 - source[0];
		break;
	    case 4:
		phi = source[1];
		theta = source[0];
		break;
	    default:
		break;
	}
    }

/*  see if any of these points are on special case boundary  */

    check = 0;
    phi = source[0];
    theta = source[1];
    phicheck = phi - pi2/2.0;
    if (phicheck < 0.0) phicheck = -phicheck;
    if (phicheck < ANGLEPS)  check = 1;
    if (theta < ANGLEPS) check = 1;
    thcheck = phi - theta;
    if (thcheck < 0.0) thcheck = -thcheck;
    if( thcheck < ANGLEPS) check = 1;
    if (check)
    {

/*  brute force simple check to remove duplicate points  */

	vector[0] = phi;
	vector[1] = theta;
	n = 1;
	for( i=2; i<96; i+=2)
	{
	    check = 1;
	    j = 0;
	    while( (check > 0) && (j < n))
	    {
		k = 2*j;
		phicheck = vec[i] - vector[k];
		if (phicheck < 0.0) phicheck = -phicheck;
		thcheck = vec[i+1] - vector[k+1];
		if (thcheck < 0.0) thcheck = -thcheck;
		if ((phicheck < ANGLEPS) && (thcheck < ANGLEPS))
		    check = 0;
		j++;
	    }
	    if(check)
	    {
		vector[2*n] = vec[i];
		vector[2*n+1] = vec[i+1];
		n++;
	    }
	}
	return n;
    }
    else
    {
	for(i=0; i<96; i+=2)
	{
	    vector[i] = vec[i];
	    vector[i+1] = vec[i+1];
	}
	return 48;
    }
}

/*  General integral over 1/48 volume segment.
    Input is pointer to table of NUMELMNTS which is
    summed over the volume segment.  Note that 
    r^2 dr  cos(theta) dtheta  dphi is computed here, not
    in input table.
    Returns integral value by accounting for duplicate faces.
*/

double integrate_sphere( double *func)
{
    int i, j, k, idex, kdex;
    double  dr, theta, acp, sum, dsum;

    sum = func[0];
    dr = pow(MAXWALL/MAXSTEPS, 3.0);
    dr *= M_PI*M_PI/16.0;

/*  sum over lines which have very limited numbers  */
/*  x, y, z axies first  */

    dsum = 0.0;
    for(i=1; i<MAXSTEPS; i++)
    {
	idex = i*(i+1)*(i+2)/6;
	dsum += func[idex];
    }
    dsum *= 6.0;
    sum += dsum;

/*  center of each octant where sextants meet  */

    dsum = 0.0;
    for(i=1; i<MAXSTEPS; i++)
    {
	idex = i*(i+1)*(i+2)/6 + i*(i+1)/2 + i;
	dsum += func[idex]*cos(M_PI/4.0);
    }
    dsum *= 8.0;
    sum += dsum;

/*  center lines on each face between octants  */

    dsum = 0.0;
    for(i=1; i<MAXSTEPS; i++)
    {
	idex = i*(i+1)*(i+2)/6 + i*(i+1)/2;
	dsum += func[idex];
    }
    dsum *= 12.0;
    sum += dsum;

/*  face planes  */

    dsum = 0.0;
    for(i=2; i<MAXSTEPS; i++)
    {
	idex = i*(i+1)*(i+2)/6;
	acp = M_PI/4.0/i;
	for(k=1; k<i; k++)
	{
	    kdex = k*(k+1)/2 + idex;
	    dsum += func[kdex];
	    theta = acp*k;
	    dsum += func[kdex + i]*cos(theta);
	}
	for(j=1; j<i; j++)
	{
	    theta = acp*j;
	    kdex = j + i*(i+1)/2 + idex;
	    dsum += func[kdex]*cos(theta);
	}
    }
//    dsum *= 24.0;
    dsum *= 16.0;  //  Gives right answer - why??
    sum += dsum;

/*  core integral, no special cases  */

    dsum = 0.0;
    for(i=2; i<MAXSTEPS; i++)
    {
	idex = i*(i+1)*(i+2)/6;
	acp = M_PI/4.0/i;
	for(k=1; k<i; k++)
	{
	    kdex = k*(k+1)/2 + idex;
	    for(j=1; j<k; j++)
	    {
		theta = acp*j;
		dsum += func[kdex+j]*cos(theta);
	    }
	}
    }
    dsum *= 48.0;
    sum += dsum;
    return sum*dr;
}

main()
{
    FILE *readmag, *fi_pot, *denscheck;
    int numread;
    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, scale;
    double rprime, thetaprime, phiprime;
    double anglprime[2], sympts[48], t1, t2, t3, dist, dtr, dmr;
    int ip, jp, kp, npts, np, edex, kdex, passnum, abdex;
    double delta_volume, dsum, sum, dr, temp, acp;
    double abmax, abval;

/*  load grid potential from previous calculation  */

    sprintf(filename, "potential.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(GridPot, sizeof(double), NUMELMNTS, readmag);
    if(numread < NUMELMNTS)
	printf("\nWARNING: not all field file read in. \n\n");
    fclose(readmag);

/*  load magnetic potential from previous calculation  */

    sprintf(filename, "magnetic_potential.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(MagPot, sizeof(double), 3*NUMELMNTS, readmag);
    if(numread < 3*NUMELMNTS)
	printf("\nWARNING: not all field file read in. \n\n");
    fclose(readmag);

/*  create file to hold computed output from this program  */

    fi_pot = fopen("electron_potential.dat", "w");
//    fi_pot = fopen("test_integral.dat", "w");
    if (!fi_pot)
    {
	printf("can't create electric potential data file. \n");
	exit(0);
    }

/* I need to ping pong solution, so start with nothing and see how fast
   it stabilizes  */

    for( i=0; i<MAXSTEPS+1; i++)
    {
	idex = i*(i+1)*(i+2)/6;
	for( j=0; j<=i; j++)
	{
	    jdex = j*(j+1)/2;
	    for(k=0; k<=j; k++)
	    {
		dex = idex + jdex + k;
		ElecPot0[dex] = -GridPot[dex];
	    }
	}
    }

/*  Compute volume element.  Every shell has same volume in terms of
    r^2*dr*dtheta*dphi.  */

    pi4 = M_PI/4.0;          // ubiquitous constant
    delta_volume = pow(MAXWALL/MAXSTEPS, 3.0);
    delta_volume *= pi4*pi4;
    scale = M_PI*12.0*COILRADIUS*delta_volume;

/*  Find maximum magnetic potential amplitude and save it for integral over
    space.
*/

    abmax = 0.0;
    for(i=0; i<NUMELMNTS; i++)
    {
	abval = MagPot[3*i]*MagPot[3*i] +
	    MagPot[3*i+1]*MagPot[3*i+1] + MagPot[3*i+2]*MagPot[3*i+2];
	if(abval > abmax)
	{
	    abmax = abval;
	    abdex = i;
	}
    }
    printf("abmax = %lf\n", abmax);

/*  Compute integral over 1/48 volume to get normalization constant.  
    All constants moved into scale factor.  */

    alpha = DENSITY_A;
    beta = DENSITY_B;
    for(i=0; i<=MAXSTEPS; i++)
    {
	idex = i*(i+1)*(i+2)/6;
	rprime = MAXWALL/MAXSTEPS*i;
	rprime *= rprime;
	for(k=0; k<=i; k++)
	{
	    kdex = k*(k+1)/2 + idex;
	    for(j=0; j<=k; j++)
	    {
		dex = idex + kdex + j;
		magpot_amp = MagPot[3*dex]*MagPot[3*dex];
		magpot_amp += MagPot[3*dex+1]*MagPot[3*dex+1];
		magpot_amp += MagPot[3*dex+2]*MagPot[3*dex+2];
		temp = -(GridPot[dex] + ElecPot0[dex])/ELECTRON_TEMP;
//		temp += (magpot_amp - abmax)*ELECTRON_TEMP/4.0/CONST_CP;
//		energy[dex] = sqrt(magpot_amp)*exp(temp);
		energy[dex] = exp(temp);
		energy[dex] /=(exp(-alpha*rprime) + beta*rprime);
	    }
	}
    }
    sum = integrate_sphere( energy);
    printf("Test integral = %lg\n", sum);
    scale /= sum;          // modify by density integral

/*  compute symmetry angles to every point.  This
    lookup table is used to compute the distance between all points
    during the integration process.  */

    for (i=1; i<MAXSTEPS; i++)
    {
	acp = pi4/i;
	idex = i*(i+1)*(i+2)/6;
	for( k=0; k<=i; k++)
	{
	    anglprime[0] = acp*k;
	    kdex = k*(k+1)/2 + idex;
	    for(j=0; j<=k; j++)
	    {
		dex = kdex + j;
		anglprime[1] = acp*j;
		numangles[dex] = symmetry(anglprime, &angles[dex*96]);

	    }
	}
    }

/*  perform integral by summing over all symmetry points  */

    ElecPot1[0] = ElecPot0[0];
    fwrite( &ElecPot1[0], sizeof(double), 1, fi_pot);
    for( i=1; i<MAXSTEPS; i++)
    {
	printf("i = %d\n", i);
	r = MAXWALL/MAXSTEPS * i ;  // radius from center
	across = pi4/i;
	for(k=0; k<=i; k++)
	{
	    phi = across * k;  //  angle from x axis
	    for( j=0; j<=k; j++)
	    {
		theta = across * j;  //  angle from x,y plane
		edex = 1;
		sum = 0.0;
		for(ip = 1; ip<MAXSTEPS; ip++)
		{
		    rprime = MAXWALL/MAXSTEPS*ip;
		    dtr = r*r + rprime*rprime;
		    dmr = r*rprime;
		    for( kp=0; kp<=ip; kp++)
		    {
			dex = ip*(ip+1)/2;
			for(jp = 0; jp<=kp; jp++)
			{
			    npts = numangles[edex];
			    dsum = 0.0;
			    for(np=0; np<npts; np++)
			    {
				if((i==ip) && (k==kp) && (j==jp))
				    continue;
				thetaprime = angles[96*edex+2*np+1];
				phiprime = angles[96*edex+2*np];
				t1 = cos(theta - thetaprime);
				t2 = cos(phi - phiprime);
				if(((1.0 - t1) < ANGLEPS) && ((1.0 - t2) < ANGLEPS))
				    continue;
				t3 = t1*(1.0 + t2);
				t1 = cos(theta + thetaprime);
				t3 += t1*(t2 - 1.0);
				dist = sqrt(dtr - dmr*t3);
				dsum += 1.0/dist;
			    }
			    dsum *= energy[edex];
			    edex++;
			}
			sum += dsum;
		    }
		}
		ElecPot1[dex] = sum*scale;
		fwrite( &ElecPot1[dex], sizeof(double), 1, fi_pot);
	    }
	}
    }
    fclose(fi_pot);

/*  OK, now let's do it again and see if the new particle distribution and potential match
    with a self consistent description.  If this diverges, the math can't work.
*/

//    for(passnum=2; passnum<5; passnum++)
    {
/*  create file to hold computed output from this program  */

	passnum = 2;   //  it actually converges on first pass
	sprintf(filename, "electron_potential_pass%d.dat", passnum);
	fi_pot = fopen(filename, "w");
	if (!fi_pot)
	{
	    printf("can't create electric potential data file. \n");
	    exit(0);
	}

/* This is start of ping-pong.  I can either change code, use pointers,
   or move data.  For my sanity, I'll move data.  */

	for( i=0; i<MAXSTEPS+1; i++)
	{
	    idex = i*(i+1)*(i+2)/6;
	    for( j=0; j<=i; j++)
	    {
		jdex = j*(j+1)/2;
		for(k=0; k<=j; k++)
		{
		    dex = idex + jdex + k;
		    ElecPot0[dex] = ElecPot1[dex];
		}
	    }
	}

/*  copy everything and recompute the whole 9 yards again.  */

	printf("\n\nStarting Pass %d....\n", passnum);
	scale = M_PI*12.0*COILRADIUS*delta_volume;
	
/*  Compute integral over 1/48 volume to get normalization constant.  
    All constants moved into scale factor.  */

	alpha = DENSITY_A;
	beta = DENSITY_B;
	for(i=0; i<=MAXSTEPS; i++)
	{
	    idex = i*(i+1)*(i+2)/6;
	    rprime = MAXWALL/MAXSTEPS*i;
	    rprime *= rprime;
	    for(k=0; k<=i; k++)
	    {
		kdex = k*(k+1)/2 + idex;
		for(j=0; j<=k; j++)
		{
		    dex = idex + kdex + j;
		    magpot_amp = MagPot[3*dex]*MagPot[3*dex];
		    magpot_amp += MagPot[3*dex+1]*MagPot[3*dex+1];
		    magpot_amp += MagPot[3*dex+2]*MagPot[3*dex+2];
		    temp = -(GridPot[dex] + ElecPot0[dex])/ELECTRON_TEMP;
//		    temp += (magpot_amp - abmax)*ELECTRON_TEMP/4.0/CONST_CP;
//		    energy[dex] = sqrt(magpot_amp)*exp(temp);
		    energy[dex] = exp(temp);
		    energy[dex] /=(exp(-alpha*rprime) + beta*rprime);
		}
	    }
	}
	sum = integrate_sphere( energy);
	printf("Test integral = %lg\n", sum);
	scale /= sum;          // modify by density integral
	
/*  perform integral by summing over all symmetry points  */
	
	ElecPot1[0] = ElecPot0[0];
	fwrite( &ElecPot1[0], sizeof(double), 1, fi_pot);
	for( i=1; i<MAXSTEPS; i++)
	{
	    printf("i = %d\n", i);
	    r = MAXWALL/MAXSTEPS * i ;  // radius from center
	    across = pi4/i;
	    for(k=0; k<=i; k++)
	    {
		phi = across * k;  //  angle from x axis
		for( j=0; j<=k; j++)
		{
		    theta = across * j;  //  angle from x,y plane
		    edex = 1;
		    sum = 0.0;
		    for(ip = 1; ip<MAXSTEPS; ip++)
		    {
			rprime = MAXWALL/MAXSTEPS*ip;
			dtr = r*r + rprime*rprime;
			dmr = r*rprime;
			for( kp=0; kp<=ip; kp++)
			{
			    dex = ip*(ip+1)/2;
			    for(jp = 0; jp<=kp; jp++)
			    {
				npts = numangles[edex];
				dsum = 0.0;
				for(np=0; np<npts; np++)
				{
				    if((i==ip) && (k==kp) && (j==jp))
					continue;
				    thetaprime = angles[96*edex+2*np+1];
				    phiprime = angles[96*edex+2*np];
				    t1 = cos(theta - thetaprime);
				    t2 = cos(phi - phiprime);
				    if(((1.0 - t1) < ANGLEPS) && ((1.0 - t2) < ANGLEPS))
					continue;
				    t3 = t1*(1.0 + t2);
				    t1 = cos(theta + thetaprime);
				    t3 += t1*(t2 - 1.0);
				    dist = sqrt(dtr - dmr*t3);
				    dsum += 1.0/dist;
				}
				dsum *= energy[edex];
				edex++;
			    }
			    sum += dsum;
			}
		    }
		    ElecPot1[dex] = sum*scale;
		    fwrite( &ElecPot1[dex], sizeof(double), 1, fi_pot);
		}
	    }
	}
	fclose(fi_pot);
    }
}