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