/****************************************************************************
 *                                                                          *
 *    Compute magnetic potential for Polywell MaGrid.  Each point in space  *
 *  is part of a 1/48 block in r, theta, phi coordinates where r is radius  *
 *  from center, theta is angle up from x,y plane and phi is angle from x   *
 *  axis.  Actual potential is computed in cartesian space, so vectors are  *
 *  referenced to x, y and z axies.  This makes computation easier, but I   *
 *  may revisit that later.                                                 *
 *                                                                          *
 *  Thanks to stank from polywell-talk.org for the fundamental integration  *
 *  parameter passing method.  No way I could have figured it out myself!   *
 *                                                                          *
 *                   Author = Mike Rosing                                   *
 *                    date = Feb 16, 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-20)
#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	theta, phi;
    double      k1, k2, c1;
    double      cz1, cz2;
} Integrand;

/*  coil 1 first  on positive x axis*/

double a1y (double d, void *params)
{	
    double cp1, sd, cd, t1;
    double ct1, ct2, tprm, pprm;
    Integrand	*p = (Integrand *) params;
	
    sd = sin(d);
    cd = cos(d);
    tprm = asin(p->c1*sd);
    pprm = atan(COILRADIUS * cd);
    ct1 = cos(p->theta - tprm);
    ct2 = cos(p->theta + tprm);
    cp1 = cos(p->phi - pprm);
    t1 = p->k1 - p->k2*(ct1*(1.0 + cp1) + ct2*(cp1 - 1.0));
    if (t1 < TOOSMALL) return 0.0;
    return sd / sqrt(t1);
}

double a1z (double d, void *params)
{	
    double cp1, sd, cd, t1;
    double ct1, ct2, tprm, pprm;
    Integrand	*p = (Integrand *) params;
	
    sd = sin(d);
    cd = cos(d);
    tprm = asin(p->c1*sd);
    pprm = atan(COILRADIUS * cd);
    ct1 = cos(p->theta - tprm);
    ct2 = cos(p->theta + tprm);
    cp1 = cos(p->phi - pprm);
    t1 = p->k1 - p->k2*(ct1*(1.0 + cp1) + ct2*(cp1 - 1.0));
    if (t1 < TOOSMALL) return 0.0;
    return -cd / sqrt(t1);
}

/*  coil 2 next,  on negative x axis*/

double a2y (double d, void *params)
{	
    double cp1, sd, cd, t1;
    double ct1, ct2, tprm, pprm;
    Integrand	*p = (Integrand *) params;
	
    sd = sin(d);
    cd = cos(d);
    tprm = asin(p->c1*sd);
    pprm = atan(COILRADIUS * cd);
    ct1 = cos(p->theta - tprm);
    ct2 = cos(p->theta + tprm);
    cp1 = -cos(p->phi - pprm);
    t1 = p->k1 - p->k2*(ct1*(1.0 + cp1) + ct2*(cp1 - 1.0));
    if (t1 < TOOSMALL) return 0.0;
    return -sd / sqrt(t1);
}

double a2z (double d, void *params)
{	
    double cp1, sd, cd, t1;
    double ct1, ct2, tprm, pprm;
    Integrand	*p = (Integrand *) params;
	
    sd = sin(d);
    cd = cos(d);
    tprm = asin(p->c1*sd);
    pprm = atan(COILRADIUS * cd);
    ct1 = cos(p->theta - tprm);
    ct2 = cos(p->theta + tprm);
    cp1 = -cos(p->phi - pprm);
    t1 = p->k1 - p->k2*(ct1*(1.0 + cp1) + ct2*(cp1 - 1.0));
    if (t1 < TOOSMALL) return 0.0;
    return cd / sqrt(t1);
}

/*  coil 3 on positive y axis*/

double a3x (double d, void *params)
{	
    double cp1, sd, cd, t1;
    double ct1, ct2, tprm, pprm;
    Integrand	*p = (Integrand *) params;
	
    sd = sin(d);
    cd = cos(d);
    tprm = asin(p->c1*sd);
    pprm = atan(COILRADIUS * cd);
    ct1 = cos(p->theta - tprm);
    ct2 = cos(p->theta + tprm);
    cp1 = sin(p->phi - pprm);
    t1 = p->k1 - p->k2*(ct1*(1.0 + cp1) + ct2*(cp1 - 1.0));
    if (t1 < TOOSMALL) return 0.0;
    return sd / sqrt(t1);
}

double a3z (double d, void *params)
{	
    double cp1, sd, cd, t1;
    double ct1, ct2, tprm, pprm;
    Integrand	*p = (Integrand *) params;
	
    sd = sin(d);
    cd = cos(d);
    tprm = asin(p->c1*sd);
    pprm = atan(COILRADIUS * cd);
    ct1 = cos(p->theta - tprm);
    ct2 = cos(p->theta + tprm);
    cp1 = sin(p->phi - pprm);
    t1 = p->k1 - p->k2*(ct1*(1.0 + cp1) + ct2*(cp1 - 1.0));
    if (t1 < TOOSMALL) return 0.0;
    return -cd / sqrt(t1);
}

/*  coil 4 on negative y axis*/

double a4x (double d, void *params)
{	
    double cp1, sd, cd, t1;
    double ct1, ct2, tprm, pprm;
    Integrand	*p = (Integrand *) params;
	
    sd = sin(d);
    cd = cos(d);
    tprm = asin(p->c1*sd);
    pprm = atan(COILRADIUS * cd);
    ct1 = cos(p->theta - tprm);
    ct2 = cos(p->theta + tprm);
    cp1 = -sin(p->phi - pprm);
    t1 = p->k1 - p->k2*(ct1*(1.0 + cp1) + ct2*(cp1 - 1.0));
    if (t1 < TOOSMALL) return 0.0;
    return -sd / sqrt(t1);
}

double a4z (double d, void *params)
{	
    double cp1, sd, cd, t1;
    double ct1, ct2, tprm, pprm;
    Integrand	*p = (Integrand *) params;
	
    sd = sin(d);
    cd = cos(d);
    tprm = asin(p->c1*sd);
    pprm = atan(COILRADIUS * cd);
    ct1 = cos(p->theta - tprm);
    ct2 = cos(p->theta + tprm);
    cp1 = -sin(p->phi - pprm);
    t1 = p->k1 - p->k2*(ct1*(1.0 + cp1) + ct2*(cp1 - 1.0));
    if (t1 < TOOSMALL) return 0.0;
    return cd / sqrt(t1);
}

/*  coil 5 on positive z axis, use cz1 and cz2 instead of theta  */

double a5x (double phi, void *params)
{
    double cp1, sd, cd, t1;
    double ct1, ct2, tprm, pprm;
    Integrand	*p = (Integrand *) params;
	
    sd = sin(phi);

    cp1 = cos(p->phi - phi);
    t1 = p->k1 - p->k2*(p->cz1*(1.0 + cp1) + p->cz2*(cp1 - 1.0));
    if (t1 < TOOSMALL) return 0.0;
    return sd / sqrt(t1);
}

double a5y (double phi, void *params)
{
    double cp1, sd, cd, t1;
    double ct1, ct2, tprm, pprm;
    Integrand	*p = (Integrand *) params;
	
    cd = cos(phi);

    cp1 = cos(p->phi - phi);
    t1 = p->k1 - p->k2*(p->cz1*(1.0 + cp1) + p->cz2*(cp1 - 1.0));
    if (t1 < TOOSMALL) return 0.0;
    return -cd / sqrt(t1);
}

/*  coil 6 on negative z axis, use cz1 and cz2 instead of theta 
    don't forget to change theta' !! */

double a6x (double phi, void *params)
{
    double cp1, sd, t1;
    Integrand	*p = (Integrand *) params;
	
    sd = sin(phi);

    cp1 = cos(p->phi - phi);
    t1 = p->k1 - p->k2*(p->cz1*(1.0 + cp1) + p->cz2*(cp1 - 1.0));
    if (t1 < TOOSMALL) return 0.0;
    return -sd / sqrt(t1);
}

double a6y (double phi, void *params)
{
    double cp1, cd, t1;
    Integrand	*p = (Integrand *) params;
	
    cd = cos(phi);

    cp1 = cos(p->phi - phi);
    t1 = p->k1 - p->k2*(p->cz1*(1.0 + cp1) + p->cz2*(cp1 - 1.0));
    if (t1 < TOOSMALL) return 0.0;
    return cd / sqrt(t1);
}

main()
{
    int i, j, k;
    double ax, ay, az, error, r, theta, phi, tzprm;
    double Rpls1, sqrRpls1, atotal[3], across, up;
    Integrand	iData;
    FILE *afield;
    gsl_integration_workspace *w;
    gsl_function f1;
    
    w = gsl_integration_workspace_alloc(1000);
	
    Rpls1 = COILRADIUS*COILRADIUS + 1.0;
    sqrRpls1 = sqrt(Rpls1);
    tzprm = atan(1/COILRADIUS);
    iData.c1 = COILRADIUS / sqrRpls1;
	
    afield = fopen("magnetic_potential.dat", "w");
    if (!afield)
    {
		printf("can't create magnetic 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<3; i++) atotal[i] = 0.0;
    fwrite( &atotal, sizeof(double), 3, afield);
    for( i=1; i<MAXSTEPS+1; i++)
    {
	printf("i = %d\n", i);
	r = MAXWALL/MAXSTEPS * i ;  // radius from center
	for( k=0; k<=i; k++)
	{
	    if (i) phi = up / i * k;  //  angle from x axis
	    else phi = 0.0;
	    iData.phi = phi;
	    for(j=0; j<=k; j++)
	    {
		if(i) theta = across / i * j;  //  angle from x,y plane
		else theta = 0.0;
		iData.theta = theta;

		/*  compute integration constants for each point  */
		
		iData.k1 = r*r + Rpls1;  //  r^2 + 1 + R^2
		iData.k2 = r*sqrRpls1;  //  r*sqrt(1 + R^2)
		iData.c1 = COILRADIUS/sqrRpls1;  //  R/sqrt(1 + R^2)
		iData.cz1 = cos(theta - tzprm);
		iData.cz2 = cos(theta + tzprm);
		
		f1.params = &iData;
		f1.function = & a1y;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &ay, &error);
		atotal[1] = ay;
		
		f1.function = & a1z;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &az, &error);
		atotal[2] = az;

		f1.function = & a2y;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &ay, &error);
		atotal[1] += ay;
		
		f1.function = & a2z;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &az, &error);
		atotal[2] += az;

		f1.function = & a3x;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &ax, &error);
		atotal[0] = ax;
		
		f1.function = & a3z;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &az, &error);
		atotal[2] += az;
		
		f1.function = & a4x;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &ax, &error);
		atotal[0] += ax;
		
		f1.function = & a4z;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &az, &error);
		atotal[2] += az;

		f1.function = & a5x;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &ax, &error);
		atotal[0] += ax;
		
		f1.function = & a5y;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &ay, &error);
		atotal[1] += ay;

		f1.function = & a6x;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &ax, &error);
		atotal[0] += ax;
		
		f1.function = & a6y;
		gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000,
				     GSL_INTEG_GAUSS41, w, &ay, &error);
		atotal[1] += ay;

		atotal[0] *= COILRADIUS;  // multiply by R outside all integrals
		atotal[1] *= COILRADIUS;  // multiply by R outside all integrals
		atotal[2] *= COILRADIUS;  // multiply by R outside all integrals
		fwrite( &atotal, sizeof(double), 3, afield);
	    }
	}
    }
    gsl_integration_workspace_free(w);
    fclose(afield);
}