/********************************************************************************
 *                                                                              *
 *    Compute potential for 6 coils.       Output is file of unitless form      *
 *  with ux, uy, uz being relative to coil distance from center.  All coils are *
 *  1 unit from center.  Coil diameter is parameter R.  Maximum physical R is   *
 *  1.0 and coils are assumed uniform in size and current.                      *
 *  Only positive octant is computed.  All other regions are mirrored from it.  *
 *                                                                              *
 *    Potential and effective length are combined to create unitless electric   *
 *  field.  Actual potential is scaled using:                                   *
 *                                                                              *
 *                                psi = V * b                                   *
 *                                                                              *
 *  where b is what this program computes, V is the coil potential and L is     *
 *  the distance from center of fusor to center of coil.  Parameter "alpha"     *
 *  is radius of grounded sphere around center of fusor (in p[4])               *
 *                                                                              *
 *                     Author = Mike Rosing                                     *
 *                      date  = Jan. 20, 2008                                   *
 *                                                                              *
 *  This version is from stank on polywell_talk.org.  Very much improved, and   *
 *  some 50 times faster than the original, it's proof that many heads are      *
 *  better than one.                                                            *
 *                                                                              *
 *******************************************************************************/

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

/*  Need one subroutine for each integral to be performed.  See notes for
 derivations.  */

#define TOOSMALL (1e-6)
#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)

typedef struct integral {
	double	ux, uy, uz;
	double	R, k, k2, k2_R;
	double	ux_Minus1_Sqrd;
	double	ux_MinusK2_Sqrd;
	double	ux_Plus1_Sqrd;
	double	ux_PlusK2_Sqrd;
	double	uy_Minus1_Sqrd;
	double	uy_MinusK2_Sqrd;
	double	uy_Plus1_Sqrd;
	double	uy_PlusK2_Sqrd;
	double	uz_Minus1_Sqrd;
	double	uz_MinusK2_Sqrd;
	double	uz_Plus1_Sqrd;
	double	uz_PlusK2_Sqrd;
} tIntegralData;

/*  coil 1 first  on positive x axis*/

double e1 (double phi, void *params)
{	
    double s1, c1, rx, ry, rz;
    double t1, t2;
	tIntegralData	*p = (tIntegralData *) params;
	
	s1 = sin(phi);
    c1 = cos(phi);
    ry = p->uy - p->R * c1;
    rz = p->uz - p->R * s1;
    t1 = sqrt(p->ux_Minus1_Sqrd + ry*ry + rz*rz);
    if (t1 < TOOSMALL) return 0.0;
    t2 = 1.0/t1;
    ry = p->uy - p->k2_R * c1;
    rz = p->uz - p->k2_R * s1;
    t1 = sqrt(p->ux_MinusK2_Sqrd + ry*ry + rz*rz);
    return t2 - p->k / t1;
}


/*  coil 2   on negative x axis*/

double e2 (double phi, void *params)
{	
    double s1, c1, rx, ry, rz;
    double t1, t2;
	tIntegralData	*p = (tIntegralData *) params;

    s1 = sin(phi);
    c1 = cos(phi);
    ry = p->ux - p->R * c1;
    rz = p->uz - p->R * s1;
    t1 = sqrt(p->ux_Plus1_Sqrd + ry*ry + rz*rz);
    if (t1 < TOOSMALL) return 0.0;
    t2 = 1.0/t1;
    ry = p->uy - p->k2_R * c1;
    rz = p->uz - p->k2_R * s1;
    t1 = sqrt(p->ux_PlusK2_Sqrd + ry*ry + rz*rz);
    return t2 - p->k / t1;
}

/*  coil 3  on positive y axis */

double e3 (double phi, void *params)
{
    double s1, c1, rx, ry, rz;
    double t1, t2;
	tIntegralData	*p = (tIntegralData *) params;
	 
    s1 = sin(phi);
    c1 = cos(phi);
    rx =  p->ux + p->R * c1;
    rz =  p->uz - p->R * s1;
    t1 = sqrt(rx*rx + p->uy_Minus1_Sqrd + rz*rz);
    if (t1 < TOOSMALL) return 0.0;
    t2 = 1.0/t1;
    rx =  p->ux +  p->k2_R * c1;
    rz =  p->uz -  p->k2_R * s1;
    t1 = sqrt(rx*rx + p->uy_MinusK2_Sqrd + rz*rz);
	return t2 - p->k / t1;
}

/*  coil 4  on negative y axis */

double e4 (double phi, void *params)
{
	
    double s1, c1, rx, ry, rz;
    double t1, t2;
	tIntegralData	*p = (tIntegralData *) params;
	
    s1 = sin(phi);
    c1 = cos(phi);
    rx = p->ux + p->R * c1;
    rz = p->uz - p->R * s1;
    t1 = sqrt(rx*rx + p->uy_Plus1_Sqrd + rz*rz);
    if (t1 < TOOSMALL) return 0.0;
    t2 = 1.0/t1;
    rx = p->ux + p->k2_R * c1;
    rz = p->uz - p->k2_R * s1;
    t1 = sqrt(rx*rx + p->uy_PlusK2_Sqrd + rz*rz);
    return t2 - p->k / t1;
}

/*  coil 5  on positive z axis */

double e5 (double phi, void *params)
{	
    double s1, c1, rx, ry, rz;
    double t1, t2;	
	tIntegralData	*p = (tIntegralData *) params;

    s1 = sin(phi);
    c1 = cos(phi);
    rx = p->ux - p->R * c1;
    ry = p->uy - p->R * s1;
    t1 = sqrt(rx*rx + ry*ry + p->uz_Minus1_Sqrd);
    if (t1 < TOOSMALL) return 0.0;
    t2 = 1.0/t1;
    rx = p->ux - p->k2_R * c1;
    ry = p->uy - p->k2_R * s1;
    t1 = sqrt(rx*rx + ry*ry + p->uz_MinusK2_Sqrd);
    return t2 - p->k / t1;
}

/*  coil 6  on negative z axis */

double e6 (double phi, void *params)
{	
    double s1, c1, rx, ry, rz;
    double t1, t2;
	tIntegralData	*p = (tIntegralData *) params;
	
    s1 = sin(phi);
    c1 = cos(phi);
    rx = p->ux - p->R * c1;
    ry = p->uy - p->R * s1;
    t1 = sqrt(rx*rx + ry*ry + p->uz_Plus1_Sqrd);
    if (t1 < TOOSMALL) return 0.0;
    t2 = 1.0/t1;
    rx = p->ux - p->k2_R * c1;
    ry = p->uy - p->k2_R * s1;
    t1 = sqrt(rx*rx + ry*ry + p->uz_PlusK2_Sqrd);
	return t2 - p->k / t1;
}

main()
{
    int i, j, k;
    double ex, ey, ez, error, r, theta, phi;
    double alpha, alpha2, extotal, across, up;
	tIntegralData	iData;
    FILE *efield;
    gsl_integration_workspace *w;
    gsl_function f1;
    
    w = gsl_integration_workspace_alloc(1000);
    alpha = MAXWALL; // radius of sphere in units of L 
    alpha2 = alpha*alpha;
	
    iData.R = COILRADIUS;  // R = 0.9
    iData.k = alpha / sqrt(1.0 + iData.R * iData.R);  //  image charge effective magnitude
	iData.k2 = iData.k * iData.k;
	iData.k2_R = iData.k2 * iData.R;
	
    efield = fopen("potential.dat", "w");
    if (!efield)
    {
		printf("can't create electric potential data file. \n");
		exit(0);
    }
	
	/*  compute potential from geometry of coils  */
	
    up = across = M_PI/4.0;        //  x, y plane angle
	
    for( i=0; i<MAXSTEPS+1; i++)
    {
		printf("i = %d\n", i);
		r = MAXWALL/MAXSTEPS * i ;  // radius from center
		for( j=0; j<=i; j++)
		{
			if (i) theta = across / i * j;  //  angle from x axis
			else theta = 0.0;
			for(k=0; k<=j; k++)
			{
				if(j) phi = up / i * k;  //  angle from x,y plane
				else phi = 0.0;
				
				/*  compute x, y, z from r, theta, phi  */
				
				iData.ux = r*cos(theta)*cos(phi);  //  ux
				iData.uy = r*sin(theta)*cos(phi);  //  uy
				iData.uz = r*sin(phi);             //  uz
				iData.ux_Minus1_Sqrd = (iData.ux - 1.0) * (iData.ux - 1.0);
				iData.ux_MinusK2_Sqrd = (iData.ux - iData.k2) * (iData.ux - iData.k2);
				iData.ux_Plus1_Sqrd = (iData.ux + 1.0) * (iData.ux + 1.0);
				iData.ux_PlusK2_Sqrd = (iData.ux + iData.k2) * (iData.ux + iData.k2);
				iData.uy_Minus1_Sqrd = (iData.uy - 1.0) * (iData.uy - 1.0);
				iData.uy_MinusK2_Sqrd = (iData.uy - iData.k2) * (iData.uy - iData.k2);
				iData.uy_Plus1_Sqrd = (iData.uy + 1.0) * (iData.uy + 1.0);
				iData.uy_PlusK2_Sqrd = (iData.uy + iData.k2) * (iData.uy + iData.k2);
				iData.uz_Minus1_Sqrd = (iData.uz - 1.0) * (iData.uz - 1.0);
				iData.uz_MinusK2_Sqrd = (iData.uz - iData.k2) * (iData.uz - iData.k2);
				iData.uz_Plus1_Sqrd = (iData.uz + 1.0) * (iData.uz + 1.0);
				iData.uz_PlusK2_Sqrd = (iData.uz + iData.k2) * (iData.uz + iData.k2);
				
				f1.params = &iData;
				f1.function = & e1;
				gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
									GSL_INTEG_GAUSS41, w, &ex, &error);
				extotal = ex;
				
				f1.function = & e2;
				gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
									GSL_INTEG_GAUSS41, w, &ex, &error);
				extotal += ex;
				
				f1.function = & e3;
				gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
									GSL_INTEG_GAUSS41, w, &ex, &error);
				extotal += ex;
				
				f1.function = & e4;
				gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
									GSL_INTEG_GAUSS41, w, &ex, &error);
				extotal += ex;
				
				f1.function = & e5;
				gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
									GSL_INTEG_GAUSS41, w, &ex, &error);
				extotal += ex;
				
				f1.function = & e6;
				gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
									GSL_INTEG_GAUSS41, w, &ex, &error);
				extotal += ex;
				extotal *= iData.R;  // multiply by R outside all integrals
				fwrite( &extotal, sizeof(double), 1, efield);
			}
		}
    }
    gsl_integration_workspace_free(w);
    fclose(efield);
}