/******************************************************************************** * * * 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 * * * *******************************************************************************/ #include #include #include #include #include /* 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 e1 (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 = sqrt(rx*rx + ry*ry + rz*rz); if (t1 < TOOSMALL) return 0.0; t2 = 1.0/t1; rx = ux - k2; ry = uy - k2*R*c1; rz = uz - k2*R*s1; t1 = sqrt(rx*rx + ry*ry + rz*rz); t3 = k/t1; f = t2 - t3; return f; } /* coil 2 on negative x axis*/ double e2 (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 = sqrt(rx*rx + ry*ry + rz*rz); if (t1 < TOOSMALL) return 0.0; t2 = 1.0/t1; rx = ux + k2; ry = uy - k2*R*c1; rz = uz - k2*R*s1; t1 = sqrt(rx*rx + ry*ry + rz*rz); t3 = k/t1; f = t2 - t3; return f; } /* coil 3 on positive y axis */ double e3 (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 = sqrt(rx*rx + ry*ry + rz*rz); if (t1 < TOOSMALL) return 0.0; t2 = 1.0/t1; rx = ux + k2*R*c1; ry = uy - k2; rz = uz - k2*R*s1; t1 = sqrt(rx*rx + ry*ry + rz*rz); t3 = k/t1; f = t2 - t3; return f; } /* coil 4 on negative y axis */ double e4 (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 = sqrt(rx*rx + ry*ry + rz*rz); if (t1 < TOOSMALL) return 0.0; t2 = 1.0/t1; rx = ux + k2*R*c1; ry = uy + k2; rz = uz - k2*R*s1; t1 = sqrt(rx*rx + ry*ry + rz*rz); t3 = k/t1; f = t2 - t3; return f; } /* coil 5 on positive z axis */ double e5 (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 = sqrt(rx*rx + ry*ry + rz*rz); if (t1 < TOOSMALL) return 0.0; t2 = 1.0/t1; rx = ux - k2*R*c1; ry = uy - k2*R*s1; rz = uz - k2; t1 = sqrt(rx*rx + ry*ry + rz*rz); t3 = k/t1; f = t2 - t3; return f; } /* coil 6 on negative z axis */ double e6 (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 = sqrt(rx*rx + ry*ry + rz*rz); if (t1 < TOOSMALL) return 0.0; t2 = 1.0/t1; rx = ux - k2*R*c1; ry = uy - k2*R*s1; rz = uz + k2; t1 = sqrt(rx*rx + ry*ry + rz*rz); t3 = k/t1; f = t2 - t3; return f; } main() { int i, j, k; double p[5], ex, ey, ez, error, r, theta, phi; double alpha, alpha2, extotal, across, up; 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; p[3] = COILRADIUS; // R = 0.9 p[4] = alpha / sqrt(1.0 + p[3]*p[3]); // image charge effective magnitude efield = fopen("potential.dat", "w"); if (!efield) { printf("can't create electric potential data file. \n"); exit(0); } /* compute potential from geometry of coils */ across = M_PI/4.0; // x, y plane angle up = across; // z up from xy plane angle for( i=0; i