/*******************************************************************
 *                                                                 *
 *      The purpose of this program is to compute the E and B      *
 *  fields for use in the particle tracking code.  The fields are  *
 *  stored in spherical coordinates. They are unusual coordinates  *
 *  with "theta" being in the x,y plane and "phi" going up from    *
 *  the x,y plane to the point at z.  The zeroth quadrant of this  *
 *  spherical system holds all the E and B data, and Indrek's      *
 *  vector form is used to find the lookup point within that block.*
 *                                                                 *
 *  The B and E field computation code is copied from mag_6.c and  *
 *  efield.c, but the dimensions of the data are different.        *
 *                                                                 *
 *                        Author = Mike Rosing                     *
 *                         date = Jan. 1, 2008                     *
 *                                                                 *
 ******************************************************************/

#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 b1x (double phi, void *params)
{

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

    ptr = (double *)params;
    ux = ptr[0];
    uy = ptr[1];
    uz = ptr[2];
    R = ptr[3];
    s1 = sin(phi);
    c1 = cos(phi);
    rx = ux - 1.0;
    ry = uy - R*c1;
    rz = uz - R*s1;
    t1 = c1*ry + s1*rz;
//    printf("t1=%lf\n", t1);
    t2 = rx*rx + ry*ry + rz*rz;
//    printf("phi = %lf ux = %lf uy = %lf uz = %lf R = %lf\n", 
//	   phi, ux, uy, uz, R);
//    scanf("%lf", &R);
    if (t2 < TOOSMALL) return 0.0;
    t3 = sqrt(t2);
    f = t1/t2/t3;
//    printf("f = %lf\n", f);
    return f;
}

double b1y ( double phi, void *params)
{

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

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

double b1z ( double phi, void *params)
{

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

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

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

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

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

double b2y ( double phi, void *params)
{

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

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

double b2z ( double phi, void *params)
{

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

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

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

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

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

double b3y ( double phi, void *params)
{

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

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

double b3z ( double phi, void *params)
{

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

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

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

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

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

double b4y ( double phi, void *params)
{

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

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

double b4z ( double phi, void *params)
{

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

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

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

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

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

double b5y ( double phi, void *params)
{

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

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

double b5z ( double phi, void *params)
{

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

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

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

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

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

double b6y ( double phi, void *params)
{

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

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

double b6z ( double phi, void *params)
{

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

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


/*  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;
    double p[5], ex, ey, ez;
    double alpha, alpha2, extotal, eytotal, eztotal;
    FILE   *efield;
    double r, theta, phi, across, up;
    double bx, by, bz, error;
    double bxtotal, bytotal, bztotal;
    FILE   *bfield;

    gsl_integration_workspace *w;
    gsl_function f1;
    
    w = gsl_integration_workspace_alloc(1000);

/*  generate both magnetic field and electric field.  Parameter passing uses
    first 5 values p[0] -> p[3] for both B and E, but E requires one more:
    p[0] = ux
    p[1] = uy
    p[2] = uz
    p[3] = coil radius ratio
    p[4] = image charge effective magnitude
*/

    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

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

    efield = fopen("e_field.dat", "w");
    if (!efield)
    {
	printf("can't create electric field data file. \n");
	exit(0);
    }
    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 = & b1x;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &bx, &error);
		f1.function = & b1y;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &by, &error);
		f1.function = & b1z;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &bz, &error);
		bxtotal = bx;
		bytotal = by;
		bztotal = bz;

		f1.function = & b2x;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &bx, &error);
		f1.function = & b2y;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &by, &error);
		f1.function = & b2z;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &bz, &error);
		bxtotal += bx;
		bytotal += by;
		bztotal += bz;

		f1.function = & b3x;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &bx, &error);
		f1.function = & b3y;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &by, &error);
		f1.function = & b3z;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &bz, &error);
		bxtotal += bx;
		bytotal += by;
		bztotal += bz;

		f1.function = & b4x;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &bx, &error);
		f1.function = & b4y;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &by, &error);
		f1.function = & b4z;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &bz, &error);
		bxtotal += bx;
		bytotal += by;
		bztotal += bz;

		f1.function = & b5x;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &bx, &error);
		f1.function = & b5y;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &by, &error);
		f1.function = & b5z;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &bz, &error);
		bxtotal += bx;
		bytotal += by;
		bztotal += bz;

		f1.function = & b6x;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &bx, &error);
		f1.function = & b6y;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &by, &error);
		f1.function = & b6z;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &bz, &error);
		bxtotal += bx;
		bytotal += by;
		bztotal += bz;
		bxtotal *= p[3];  // NOTE: forgot R * dphi in integral
		bytotal *= p[3];  // but constant comes outside ok
		bztotal *= p[3];
		fwrite( &bxtotal, sizeof(double), 1, bfield);
		fwrite( &bytotal, sizeof(double), 1, bfield);
		fwrite( &bztotal, sizeof(double), 1, bfield);

		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);
    fclose(bfield);
}