/********************************************************************
 *                                                                  *
 *      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 angles[2*48*NUMELMNTS];
int numangles[NUMELMNTS];
double test_integral[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

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

/*  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 rtlde, rmax, abserr, densparms[2];
    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;
    double delta_volume, 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*4.0*M_PI);
    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 element.  Every shell has same volume in terms of
    r^2*dr*dtheta*dphi.  Also compute scale factor:

            scale = e*N_e / (4*pi*epsilon_0* MaGrid_Voltage * L)

I'm taking N_e ~ 1e17, Voltage ~50kV, L ~ 2 m for scale ~ 4e3
*/

    pi4 = M_PI/4.0;          // ubiquitous constant
    delta_volume = pow(MAXWALL/MAXSTEPS, 3.0);
    delta_volume *= pi4*pi4;
    scale = -4e3 * delta_volume;

/*  Compare gsl integral to 1/48 estimate integral.  */

    alpha = densparms[0];
    beta = densparms[1];
    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 = kdex+j;
		temp = (GridPot[dex] + ElecPot0[dex])/ELECTRON_TEMP;
		test_integral[dex] = exp(-temp)/(exp(-alpha*rprime) + beta*rprime);
	    }
	}
    }
    sum = integrate_sphere( test_integral);
    printf("Test integral = %lf\n", sum);
    scale /= sum;          // modify by density integral

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

    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] = 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]*dens[dex];
			    edex++;
			    dex++;
			}
			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;  // converges on first pass even for high temperature electrons
	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 = -4e3 * delta_volume;
	
/*  Compare gsl integral to 1/48 estimate integral.  */
	
	alpha = densparms[0];
	beta = densparms[1];
	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 = kdex + j;
		    temp = (GridPot[dex] + ElecPot0[dex])/ELECTRON_TEMP;
		    test_integral[dex] = exp(-temp)/(exp(-alpha*rprime) + beta*rprime);
		}
	    }
	}
	sum = integrate_sphere( test_integral);
	printf("Test integral = %lf\n", sum);
	scale /= sum;          // modify by density integral
	
/*  compute energy function table.  This depends on all dimensions.  */

	energy[0] = 1.0;
	for (i=1; i<MAXSTEPS; i++)
	{
	    idex = i*(i+1)*(i+2)/6;
	    for( k=0; k<=i; 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);
		    
		}
	    }
	}
	
/*  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]*dens[dex];
				edex++;
				dex++;
			    }
			    sum += dsum;
			}
		    }
		    ElecPot1[dex] = sum*scale;
		    fwrite( &ElecPot1[dex], sizeof(double), 1, fi_pot);
		}
	    }
	}
	fclose(fi_pot);
    }
}