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

#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 400

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

/*  coil 1 first  on positive x axis*/

double e1 (double phi, void *params)
{

    double s1, c1, rx, ry, rz, f;
    double t1, t2, t3, *ptr;
    double ux, uy, uz, R, k, k2;

    ptr = (double *)params;
    ux = ptr[0];
    uy = ptr[1];
    uz = ptr[2];
    R = ptr[3];
    k = ptr[4];
    k2 = k*k;
    s1 = sin(phi);
    c1 = cos(phi);
    rx = ux - 1.0;
    ry = uy - R*c1;
    rz = uz - R*s1;
    t1 = sqrt(rx*rx + ry*ry + rz*rz);
    if (t1 < TOOSMALL) return 0.0;
    t2 = 1.0/t1;
    rx = ux - k2;
    ry = uy - k2*R*c1;
    rz = uz - k2*R*s1;
    t1 = sqrt(rx*rx + ry*ry + rz*rz);
    t3 = k/t1;
    f = t2 - t3;
    return f;
}


/*  coil 2   on negative x axis*/
    
double e2 (double phi, void *params)
{

    double s1, c1, rx, ry, rz, f;
    double t1, t2, t3, *ptr;
    double ux, uy, uz, R, k, k2;

    ptr = (double *)params;
    ux = ptr[0];
    uy = ptr[1];
    uz = ptr[2];
    R = ptr[3];
    k = ptr[4];
    k2 = k*k;
    s1 = sin(phi);
    c1 = cos(phi);
    rx = ux + 1.0;
    ry = uy - R*c1;
    rz = uz - R*s1;
    t1 = sqrt(rx*rx + ry*ry + rz*rz);
    if (t1 < TOOSMALL) return 0.0;
    t2 = 1.0/t1;
    rx = ux + k2;
    ry = uy - k2*R*c1;
    rz = uz - k2*R*s1;
    t1 = sqrt(rx*rx + ry*ry + rz*rz);
    t3 = k/t1;
    f = t2 - t3;
    return f;
}

/*  coil 3  on positive y axis */
    
double e3 (double phi, void *params)
{

    double s1, c1, rx, ry, rz, f;
    double t1, t2, t3, *ptr;
    double ux, uy, uz, R, k, k2;

    ptr = (double *)params;
    ux = ptr[0];
    uy = ptr[1];
    uz = ptr[2];
    R = ptr[3];
    k = ptr[4];
    k2 = k*k;
    s1 = sin(phi);
    c1 = cos(phi);
    rx = ux + R*c1;
    ry = uy - 1.0;
    rz = uz - R*s1;
    t1 = sqrt(rx*rx + ry*ry + rz*rz);
    if (t1 < TOOSMALL) return 0.0;
    t2 = 1.0/t1;
    rx = ux + k2*R*c1;
    ry = uy - k2;
    rz = uz - k2*R*s1;
    t1 = sqrt(rx*rx + ry*ry + rz*rz);
    t3 = k/t1;
    f = t2 - t3;
    return f;
}

/*  coil 4  on negative y axis */
    
double e4 (double phi, void *params)
{

    double s1, c1, rx, ry, rz, f;
    double t1, t2, t3, *ptr;
    double ux, uy, uz, R, k, k2;

    ptr = (double *)params;
    ux = ptr[0];
    uy = ptr[1];
    uz = ptr[2];
    R = ptr[3];
    k = ptr[4];
    k2 = k*k;
    s1 = sin(phi);
    c1 = cos(phi);
    rx = ux + R*c1;
    ry = uy + 1.0;
    rz = uz - R*s1;
    t1 = sqrt(rx*rx + ry*ry + rz*rz);
    if (t1 < TOOSMALL) return 0.0;
    t2 = 1.0/t1;
    rx = ux + k2*R*c1;
    ry = uy + k2;
    rz = uz - k2*R*s1;
    t1 = sqrt(rx*rx + ry*ry + rz*rz);
    t3 = k/t1;
    f = t2 - t3;
    return f;
}

/*  coil 5  on positive z axis */
    
double e5 (double phi, void *params)
{

    double s1, c1, rx, ry, rz, f;
    double t1, t2, t3, *ptr;
    double ux, uy, uz, R, k, k2;

    ptr = (double *)params;
    ux = ptr[0];
    uy = ptr[1];
    uz = ptr[2];
    R = ptr[3];
    k = ptr[4];
    k2 = k*k;
    s1 = sin(phi);
    c1 = cos(phi);
    rx = ux - R*c1;
    ry = uy - R*s1;
    rz = uz - 1.0;
    t1 = sqrt(rx*rx + ry*ry + rz*rz);
    if (t1 < TOOSMALL) return 0.0;
    t2 = 1.0/t1;
    rx = ux - k2*R*c1;
    ry = uy - k2*R*s1;
    rz = uz - k2;
    t1 = sqrt(rx*rx + ry*ry + rz*rz);
    t3 = k/t1;
    f = t2 - t3;
    return f;
}

/*  coil 6  on negative z axis */
    
double e6 (double phi, void *params)
{

    double s1, c1, rx, ry, rz, f;
    double t1, t2, t3, *ptr;
    double ux, uy, uz, R, k, k2;

    ptr = (double *)params;
    ux = ptr[0];
    uy = ptr[1];
    uz = ptr[2];
    R = ptr[3];
    k = ptr[4];
    k2 = k*k;
    s1 = sin(phi);
    c1 = cos(phi);
    rx = ux - R*c1;
    ry = uy - R*s1;
    rz = uz + 1.0;
    t1 = sqrt(rx*rx + ry*ry + rz*rz);
    if (t1 < TOOSMALL) return 0.0;
    t2 = 1.0/t1;
    rx = ux - k2*R*c1;
    ry = uy - k2*R*s1;
    rz = uz + k2;
    t1 = sqrt(rx*rx + ry*ry + rz*rz);
    t3 = k/t1;
    f = t2 - t3;
    return f;
}

main()
{
    int i, j, k;
    double p[5], ex, ey, ez, error, r, theta, phi;
    double alpha, alpha2, extotal, across, up;
    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;
    p[3] = COILRADIUS;  // R = 0.9
    p[4] = alpha / sqrt(1.0 + p[3]*p[3]);  //  image charge effective magnitude

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

/*  compute potential from geometry of coils  */

    across = M_PI/4.0;        //  x, y plane angle
    up = across;  //  z up from xy 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  */

		p[0] = r*cos(theta)*cos(phi);  //  ux
		p[1] = r*sin(theta)*cos(phi);  //  uy
		p[2] = r*sin(phi);             //  uz
		f1.params = p;
		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 *= p[3];  // multiply by R outside all integrals
		fwrite( &extotal, sizeof(double), 1, efield);
	    }
	}
    }
    gsl_integration_workspace_free(w);
    fclose(efield);
}