/********************************************************************************
 *                                                                              *
 *    Compute magnetic 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.  *
 *                                                                              *
 *    Current density and effective length are combined to create unitless      *
 *  magnetic field.  Actual field is scaled using:                              *
 *                                                                              *
 *                 B = (mu_0 * I_0 )/( 4.0 * PI * L) * b                        *
 *                                                                              *
 *  where b is what this program computes, I_0 is the coil current, L is the    *
 *  distance from center of fusor to center of coil and                         *
 *  mu_0 is 4*PI*1e-7 kg m / C^2                                                *
 *                                                                              *
 *                     Author = Mike Rosing                                     *
 *                      date  = Sept. 8, 2007                                   *
 *                                                                              *
 *******************************************************************************/

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

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

main()
{
    int i, j, k;
    double p[4], bx, by, bz, error;
    double bxtotal, bytotal, bztotal;
    FILE *bfield;
    gsl_integration_workspace *w;
    gsl_function f1;
    
    w = gsl_integration_workspace_alloc(1000);
    p[3] = 0.9;  // R = 0.9

    bfield = fopen("dipole.dat", "w");
    if (!bfield)
    {
	printf("can't create mag field data file. \n");
	exit(0);
    }
    for( k=0; k<256; k++)
    {
	printf("%d.\n", k);
	p[2] = 2.0/256.0*k ;   //  uz = 2.0/256 * k
	for( j=0; j<256; j++)
	{
	    p[1] = 2.0/256.0*j;  // uy = 2.0/256 * j
	    for( i=0; i<256; i++)
	    {
		p[0] = 2.0/256.0 * i ;  // ux = 2.0/256 * i
		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);
	    }
	}
    }
    gsl_integration_workspace_free(w);
    fclose(bfield);
}