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

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

//#define MAXSTEPS 100
#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 dens[MAXSTEPS*(MAXSTEPS+1)/2], energy[NUMELMNTS];
double delta_volume[MAXSTEPS], 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.1

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

/*  params[0] = a in density distribution, 
    params[1] = b
*/

double density( double p, void *params)
{
    double p2, *b;
    
    b = (double*) params;
    p2 = p*p;
    return (p2/(exp(-b[0]*p2) + b[1]*p2));
}

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

main()
{
    FILE  *readmag, *fi_pot, *denscheck;
    int numread;
    char filename[64], input[8];
    double rtlde, rmax, abserr, densparms[2];
    double across, r, phi, theta;
    int err, numval, i, j, k, idex, jdex, dex;
    double pi4, alpha, beta;
    double rprime, thetaprime, phiprime;
    double anglprime[2], sympts[48], t1, t2, t3, dist, dtr, dmr;
    int ip, jp, kp, npts, np, edex, kdex;
    double dsum, sum, dr, temp, acp;

/* debug symmetry function  */

/*    ip = 3;
    pi4 = M_PI/4.0;          // ubiquitous constant
    {
	acp = pi4/ip;
	for( kp=0; kp<=ip; kp++)
	{
	    anglprime[0] = acp*kp;
	    for(jp = 0; jp<=kp; jp++)
	    {
		anglprime[1] = acp*jp;
		npts = symmetry(anglprime, sympts);
		for(np=0; np<npts; np++)
		    printf("%lf  %lf\n", sympts[2*np], sympts[2*np+1]);
		printf("numpoints is %d\n\n", npts);
	    }
	}
	exit(0);
    }
*/	
    gsl_function F;

/*  compute normalization coefficient for density distribution chosen  */

    densparms[0] = 9.0;
    densparms[1] = 10.0;
    F.function = &density;
    F.params = densparms;
    rmax = MAXWALL;
    err = gsl_integration_qng( &F, 0.0, rmax, 
			       1e-6, 1e-6, &rtlde, &abserr, (size_t*)&numval);
    printf("Density integral result = %lf\n", rtlde);
    printf("absolute error = %lf\n", abserr);

/*  load grid potential from previous calculation  */

    sprintf(filename, "potential.dat");
    readmag = fopen(filename, "rb");
    if( !readmag)
    {
	printf("can't read in mag data file %s.\n", filename);
	exit(0);
    }
    printf("opening 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);

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

    fi_pot = fopen("electron_potential.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 elements.  Each radius is different but all angles
    are the same.  Multiply each shell by its corresponding volume element */

    pi4 = M_PI/4.0;          // ubiquitous constant
    delta_volume[0] = 1.0;
    dr = MAXWALL/MAXSTEPS*pow(pi4, 2.0);
    for(i=1; i<MAXSTEPS; i++)
	delta_volume[i] = dr/i/i;

/*  compute density distribution including radial integration factors.
    Create a lookup table over r, theta and assume uniform in phi.
    This complicates things a touch later since all other tables loop
    over phi before theta.
*/

    alpha = densparms[0];
    beta = densparms[1];
    dens[0] = 0.0;
    dex = 1;
    for(i=1; i<MAXSTEPS; i++)
    {
	rprime = MAXWALL/MAXSTEPS*i;
	rprime *= rprime;
	across = pi4/i;
	for( j=0; j<=i; j++)
	{
	    thetaprime = across*j;
	    dens[dex] = rprime*cos(thetaprime);
	    dens[dex] /= exp(-alpha*rprime) + beta*rprime;
	    dex++;
	}
    }

/*  compute energy function table.  This depends on all dimensions.
    At the same time, compute symmetry angles to every point.  This
    lookup table is used to compute the distance between all points
    during the integration process.  */

    energy[0] = 1.0;
    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;
		temp = (GridPot[dex] + ElecPot0[dex])/ELECTRON_TEMP;
		energy[dex] = exp(-temp);
		anglprime[1] = acp*j;
		numangles[dex] = symmetry(anglprime, &angles[dex*96]);

	    }
	}
    }

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

    ElecPot1[0] = 0.0;
    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);
/*				if(dtr-dmr*t3 < 1e-7)
				{
				    printf("edex = %d  dex = %d np=%d\n", edex, dex, np);
				    printf("i=%d k=%d j=%d ip=%d kp=%d jp=%d\n", 
					   i, k, j, ip, kp, jp);
				    printf("theta = %lf  phi = %lf  thetaprime = %lf phiprime=%lf\n",
					   theta, phi,thetaprime,  phiprime);
				    printf("dtr = %lf  dmr= %lf  t3 = %lf\n", dtr, dmr, t3);
				} 
*/
				dist = sqrt(dtr - dmr*t3);
				dsum += 1.0/dist;
			    }
			    dsum *= energy[edex]*dens[dex];
			    edex++;
			    dex++;
			}
			sum += dsum;
		    }
		    sum *= delta_volume[ip];
		}
		ElecPot1[dex] = sum;
		fwrite( &ElecPot1[dex], sizeof(double), 1, fi_pot);
	    }
	}
    }
/*    for( i=0; i<MAXSTEPS; 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;
		fwrite( &ElecPot1[dex], sizeof(double), 1, fi_pot);
	    }
	}
	}*/
    fclose(fi_pot);
}