/******************************************************************************** * * * Compute potential 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 potential is scaled using: * * * * psi = V * 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 = Jan. 20, 2008 * * * * This version is from stank on polywell_talk.org. Very much improved, and * * some 50 times faster than the original, it's proof that many heads are * * better than one. * * * *******************************************************************************/ #include #include #include #include #include /* Need one subroutine for each integral to be performed. See notes for derivations. */ #define TOOSMALL (1e-6) #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 ux, uy, uz; double R, k, k2, k2_R; double ux_Minus1_Sqrd; double ux_MinusK2_Sqrd; double ux_Plus1_Sqrd; double ux_PlusK2_Sqrd; double uy_Minus1_Sqrd; double uy_MinusK2_Sqrd; double uy_Plus1_Sqrd; double uy_PlusK2_Sqrd; double uz_Minus1_Sqrd; double uz_MinusK2_Sqrd; double uz_Plus1_Sqrd; double uz_PlusK2_Sqrd; } tIntegralData; /* coil 1 first on positive x axis*/ double e1 (double phi, void *params) { double s1, c1, rx, ry, rz; double t1, t2; tIntegralData *p = (tIntegralData *) params; s1 = sin(phi); c1 = cos(phi); ry = p->uy - p->R * c1; rz = p->uz - p->R * s1; t1 = sqrt(p->ux_Minus1_Sqrd + ry*ry + rz*rz); if (t1 < TOOSMALL) return 0.0; t2 = 1.0/t1; ry = p->uy - p->k2_R * c1; rz = p->uz - p->k2_R * s1; t1 = sqrt(p->ux_MinusK2_Sqrd + ry*ry + rz*rz); return t2 - p->k / t1; } /* coil 2 on negative x axis*/ double e2 (double phi, void *params) { double s1, c1, rx, ry, rz; double t1, t2; tIntegralData *p = (tIntegralData *) params; s1 = sin(phi); c1 = cos(phi); ry = p->ux - p->R * c1; rz = p->uz - p->R * s1; t1 = sqrt(p->ux_Plus1_Sqrd + ry*ry + rz*rz); if (t1 < TOOSMALL) return 0.0; t2 = 1.0/t1; ry = p->uy - p->k2_R * c1; rz = p->uz - p->k2_R * s1; t1 = sqrt(p->ux_PlusK2_Sqrd + ry*ry + rz*rz); return t2 - p->k / t1; } /* coil 3 on positive y axis */ double e3 (double phi, void *params) { double s1, c1, rx, ry, rz; double t1, t2; tIntegralData *p = (tIntegralData *) params; s1 = sin(phi); c1 = cos(phi); rx = p->ux + p->R * c1; rz = p->uz - p->R * s1; t1 = sqrt(rx*rx + p->uy_Minus1_Sqrd + rz*rz); if (t1 < TOOSMALL) return 0.0; t2 = 1.0/t1; rx = p->ux + p->k2_R * c1; rz = p->uz - p->k2_R * s1; t1 = sqrt(rx*rx + p->uy_MinusK2_Sqrd + rz*rz); return t2 - p->k / t1; } /* coil 4 on negative y axis */ double e4 (double phi, void *params) { double s1, c1, rx, ry, rz; double t1, t2; tIntegralData *p = (tIntegralData *) params; s1 = sin(phi); c1 = cos(phi); rx = p->ux + p->R * c1; rz = p->uz - p->R * s1; t1 = sqrt(rx*rx + p->uy_Plus1_Sqrd + rz*rz); if (t1 < TOOSMALL) return 0.0; t2 = 1.0/t1; rx = p->ux + p->k2_R * c1; rz = p->uz - p->k2_R * s1; t1 = sqrt(rx*rx + p->uy_PlusK2_Sqrd + rz*rz); return t2 - p->k / t1; } /* coil 5 on positive z axis */ double e5 (double phi, void *params) { double s1, c1, rx, ry, rz; double t1, t2; tIntegralData *p = (tIntegralData *) params; s1 = sin(phi); c1 = cos(phi); rx = p->ux - p->R * c1; ry = p->uy - p->R * s1; t1 = sqrt(rx*rx + ry*ry + p->uz_Minus1_Sqrd); if (t1 < TOOSMALL) return 0.0; t2 = 1.0/t1; rx = p->ux - p->k2_R * c1; ry = p->uy - p->k2_R * s1; t1 = sqrt(rx*rx + ry*ry + p->uz_MinusK2_Sqrd); return t2 - p->k / t1; } /* coil 6 on negative z axis */ double e6 (double phi, void *params) { double s1, c1, rx, ry, rz; double t1, t2; tIntegralData *p = (tIntegralData *) params; s1 = sin(phi); c1 = cos(phi); rx = p->ux - p->R * c1; ry = p->uy - p->R * s1; t1 = sqrt(rx*rx + ry*ry + p->uz_Plus1_Sqrd); if (t1 < TOOSMALL) return 0.0; t2 = 1.0/t1; rx = p->ux - p->k2_R * c1; ry = p->uy - p->k2_R * s1; t1 = sqrt(rx*rx + ry*ry + p->uz_PlusK2_Sqrd); return t2 - p->k / t1; } main() { int i, j, k; double ex, ey, ez, error, r, theta, phi; double alpha, alpha2, extotal, across, up; tIntegralData iData; FILE *efield; gsl_integration_workspace *w; gsl_function f1; w = gsl_integration_workspace_alloc(1000); alpha = MAXWALL; // radius of sphere in units of L alpha2 = alpha*alpha; iData.R = COILRADIUS; // R = 0.9 iData.k = alpha / sqrt(1.0 + iData.R * iData.R); // image charge effective magnitude iData.k2 = iData.k * iData.k; iData.k2_R = iData.k2 * iData.R; efield = fopen("potential.dat", "w"); if (!efield) { printf("can't create electric 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