/********************************************************************************
 *                                                                              *
 *    Compute electric field 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 field is scaled using:                                       *
 *                                                                              *
 *                                E = V/ L * 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  = Sept. 29, 2007                                  *
 *                                                                              *
 *******************************************************************************/

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

/*  coil 1 first  on positive x axis*/

double e1x (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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*rx/t1/sqrt(t1);
    rx = ux - k2;
    ry = uy - k2*R*c1;
    rz = uz - k2*R*s1;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*rx/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

double e1y ( 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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*ry/t1/sqrt(t1);
    rx = ux - k2;
    ry = uy - k2*R*c1;
    rz = uz - k2*R*s1;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*ry/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

double e1z ( 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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*rz/t1/sqrt(t1);
    rx = ux - k2;
    ry = uy - k2*R*c1;
    rz = uz - k2*R*s1;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*rz/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

/*  coil 2   on negative x axis*/
    
double e2x (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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*rx/t1/sqrt(t1);
    rx = ux + k2;
    ry = uy - k2*R*c1;
    rz = uz - k2*R*s1;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*rx/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

double e2y ( 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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*ry/t1/sqrt(t1);
    rx = ux + k2;
    ry = uy - k2*R*c1;
    rz = uz - k2*R*s1;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*ry/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

double e2z ( 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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*rz/t1/sqrt(t1);
    rx = ux + k2;
    ry = uy - k2*R*c1;
    rz = uz - k2*R*s1;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*rz/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

/*  coil 3  on positive y axis */
    
double e3x (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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*rx/t1/sqrt(t1);
    rx = ux + k2*R*c1;
    ry = uy - k2;
    rz = uz - k2*R*s1;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*rx/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

double e3y ( 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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*ry/t1/sqrt(t1);
    rx = ux + k2*R*c1;
    ry = uy - k2;
    rz = uz - k2*R*s1;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*ry/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

double e3z ( 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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*rz/t1/sqrt(t1);
    rx = ux + k2*R*c1;
    ry = uy - k2;
    rz = uz - k2*R*s1;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*rz/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

/*  coil 4  on negative y axis */
    
double e4x (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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*rx/t1/sqrt(t1);
    rx = ux + k2*R*c1;
    ry = uy + k2;
    rz = uz - k2*R*s1;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*rx/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

double e4y ( 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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*ry/t1/sqrt(t1);
    rx = ux + k2*R*c1;
    ry = uy + k2;
    rz = uz - k2*R*s1;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*ry/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

double e4z ( 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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*rz/t1/sqrt(t1);
    rx = ux + k2*R*c1;
    ry = uy + k2;
    rz = uz - k2*R*s1;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*rz/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

/*  coil 5  on positive z axis */
    
double e5x (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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*rx/t1/sqrt(t1);
    rx = ux - k2*R*c1;
    ry = uy - k2*R*s1;
    rz = uz - k2;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*rx/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

double e5y ( 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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*ry/t1/sqrt(t1);
    rx = ux - k2*R*c1;
    ry = uy - k2*R*s1;
    rz = uz - k2;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*ry/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

double e5z ( 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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*rz/t1/sqrt(t1);
    rx = ux - k2*R*c1;
    ry = uy - k2*R*s1;
    rz = uz - k2;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*rz/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

/*  coil 6  on negative z axis */
    
double e6x (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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*rx/t1/sqrt(t1);
    rx = ux - k2*R*c1;
    ry = uy - k2*R*s1;
    rz = uz + k2;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*rx/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

double e6y ( 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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*ry/t1/sqrt(t1);
    rx = ux - k2*R*c1;
    ry = uy - k2*R*s1;
    rz = uz + k2;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*ry/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

double e6z ( 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 = rx*rx + ry*ry + rz*rz;
    if (t1 < TOOSMALL) return 0.0;
    t2 = R*rz/t1/sqrt(t1);
    rx = ux - k2*R*c1;
    ry = uy - k2*R*s1;
    rz = uz + k2;
    t1 = rx*rx + ry*ry + rz*rz;
    t3 = R*k*rz/t1/sqrt(t1);
    f = t2 - t3;
    return f;
}

main()
{
    int i, j, k, kdex, jdex, kup, kdown, jup, jdown;
    int test1, test2;
    double p[5], ex, ey, ez, error, zlen, ylen, xlen;
    double alpha, alpha2, extotal, eytotal, eztotal;
    FILE *efield;
    gsl_integration_workspace *w;
    gsl_function f1;
    
    w = gsl_integration_workspace_alloc(1000);
    alpha = 2.0; // radius of sphere in units of L 
    alpha2 = alpha*alpha;
    p[3] = 0.9;  // R = 0.9
    p[4] = alpha / sqrt(1.0 + p[3]*p[3]);  //  image charge effective magnitude

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

/*  compute potential from geometry of coils  */

    for( k=0; k<MAXSTEPS; k++)
    {
	printf("%d.\n", k);
	p[2] = 2.0/256.0*k;   //  uz = 2.0/256 * k
	zlen = p[2] * p[2];
	for( j=0; j<MAXSTEPS; j++)
	{
	    p[1] = 2.0/256.0*j;  // uy = 2.0/256 * j
	    ylen = p[1] * p[1];
	    for( i=0; i<MAXSTEPS; i++)
	    {
		p[0] = 2.0/256.0 * i;  // ux = 2.0/256 * i
		xlen = p[0] * p[0];
		if (xlen + ylen + zlen > alpha2) 
		{
		    extotal = 0.0;  // point outside sphere
		    eytotal = 0.0;
		    eztotal = 0.0;
		}
		else
		{
		    f1.params = p;
		    f1.function = & e1x;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ex, &error);
		    f1.function = & e1y;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ey, &error);
		    f1.function = & e1z;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ez, &error);
		    extotal = ex;
		    eytotal = ey;
		    eztotal = ez;

		    f1.function = & e2x;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ex, &error);
		    f1.function = & e2y;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ey, &error);
		    f1.function = & e2z;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ez, &error);
		    extotal += ex;
		    eytotal += ey;
		    eztotal += ez;
		    
		    f1.function = & e3x;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ex, &error);
		    f1.function = & e3y;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ey, &error);
		    f1.function = & e3z;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ez, &error);
		    
		    extotal += ex;
		    eytotal += ey;
		    eztotal += ez;
		    
		    f1.function = & e4x;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ex, &error);
		    f1.function = & e4y;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ey, &error);
		    f1.function = & e4z;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ez, &error);
		    extotal += ex;
		    eytotal += ey;
		    eztotal += ez;
		    
		    f1.function = & e5x;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ex, &error);
		    f1.function = & e5y;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ey, &error);
		    f1.function = & e5z;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ez, &error);
		    
		    extotal += ex;
		    eytotal += ey;
		    eztotal += ez;
		    
		    f1.function = & e6x;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ex, &error);
		    f1.function = & e6y;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ey, &error);
		    f1.function = & e6z;
		    gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
					 GSL_INTEG_GAUSS41, w, &ez, &error);
		    extotal += ex;
		    eytotal += ey;
		    eztotal += ez;
		    
		}
		fwrite( &extotal, sizeof(double), 1, efield);
		fwrite( &eytotal, sizeof(double), 1, efield);
		fwrite( &eztotal, sizeof(double), 1, efield);
	    }
	}
    }
    gsl_integration_workspace_free(w);
    fclose(efield);
}