/******************************************************************************** * * * Compute electric 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. * * * * Potential and effective length are combined to create unitless electric * * field. Actual field is scaled using: * * * * E = V/ L * 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 = Sept. 29, 2007 * * * *******************************************************************************/ #include #include #include #include #include /* Need one subroutine for each integral to be performed. See notes for derivations. */ #define TOOSMALL (1e-6) #define MAXSTEPS 256 /* coil 1 first on positive x axis*/ double e1x (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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*rx/t1/sqrt(t1); rx = ux - k2; ry = uy - k2*R*c1; rz = uz - k2*R*s1; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*rx/t1/sqrt(t1); f = t2 - t3; return f; } double e1y ( 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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*ry/t1/sqrt(t1); rx = ux - k2; ry = uy - k2*R*c1; rz = uz - k2*R*s1; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*ry/t1/sqrt(t1); f = t2 - t3; return f; } double e1z ( 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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*rz/t1/sqrt(t1); rx = ux - k2; ry = uy - k2*R*c1; rz = uz - k2*R*s1; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*rz/t1/sqrt(t1); f = t2 - t3; return f; } /* coil 2 on negative x axis*/ double e2x (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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*rx/t1/sqrt(t1); rx = ux + k2; ry = uy - k2*R*c1; rz = uz - k2*R*s1; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*rx/t1/sqrt(t1); f = t2 - t3; return f; } double e2y ( 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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*ry/t1/sqrt(t1); rx = ux + k2; ry = uy - k2*R*c1; rz = uz - k2*R*s1; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*ry/t1/sqrt(t1); f = t2 - t3; return f; } double e2z ( 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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*rz/t1/sqrt(t1); rx = ux + k2; ry = uy - k2*R*c1; rz = uz - k2*R*s1; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*rz/t1/sqrt(t1); f = t2 - t3; return f; } /* coil 3 on positive y axis */ double e3x (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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*rx/t1/sqrt(t1); rx = ux + k2*R*c1; ry = uy - k2; rz = uz - k2*R*s1; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*rx/t1/sqrt(t1); f = t2 - t3; return f; } double e3y ( 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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*ry/t1/sqrt(t1); rx = ux + k2*R*c1; ry = uy - k2; rz = uz - k2*R*s1; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*ry/t1/sqrt(t1); f = t2 - t3; return f; } double e3z ( 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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*rz/t1/sqrt(t1); rx = ux + k2*R*c1; ry = uy - k2; rz = uz - k2*R*s1; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*rz/t1/sqrt(t1); f = t2 - t3; return f; } /* coil 4 on negative y axis */ double e4x (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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*rx/t1/sqrt(t1); rx = ux + k2*R*c1; ry = uy + k2; rz = uz - k2*R*s1; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*rx/t1/sqrt(t1); f = t2 - t3; return f; } double e4y ( 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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*ry/t1/sqrt(t1); rx = ux + k2*R*c1; ry = uy + k2; rz = uz - k2*R*s1; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*ry/t1/sqrt(t1); f = t2 - t3; return f; } double e4z ( 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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*rz/t1/sqrt(t1); rx = ux + k2*R*c1; ry = uy + k2; rz = uz - k2*R*s1; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*rz/t1/sqrt(t1); f = t2 - t3; return f; } /* coil 5 on positive z axis */ double e5x (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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*rx/t1/sqrt(t1); rx = ux - k2*R*c1; ry = uy - k2*R*s1; rz = uz - k2; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*rx/t1/sqrt(t1); f = t2 - t3; return f; } double e5y ( 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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*ry/t1/sqrt(t1); rx = ux - k2*R*c1; ry = uy - k2*R*s1; rz = uz - k2; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*ry/t1/sqrt(t1); f = t2 - t3; return f; } double e5z ( 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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*rz/t1/sqrt(t1); rx = ux - k2*R*c1; ry = uy - k2*R*s1; rz = uz - k2; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*rz/t1/sqrt(t1); f = t2 - t3; return f; } /* coil 6 on negative z axis */ double e6x (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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*rx/t1/sqrt(t1); rx = ux - k2*R*c1; ry = uy - k2*R*s1; rz = uz + k2; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*rx/t1/sqrt(t1); f = t2 - t3; return f; } double e6y ( 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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*ry/t1/sqrt(t1); rx = ux - k2*R*c1; ry = uy - k2*R*s1; rz = uz + k2; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*ry/t1/sqrt(t1); f = t2 - t3; return f; } double e6z ( 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 = rx*rx + ry*ry + rz*rz; if (t1 < TOOSMALL) return 0.0; t2 = R*rz/t1/sqrt(t1); rx = ux - k2*R*c1; ry = uy - k2*R*s1; rz = uz + k2; t1 = rx*rx + ry*ry + rz*rz; t3 = R*k*rz/t1/sqrt(t1); f = t2 - t3; return f; } main() { int i, j, k, kdex, jdex, kup, kdown, jup, jdown; int test1, test2; double p[5], ex, ey, ez, error, zlen, ylen, xlen; double alpha, alpha2, extotal, eytotal, eztotal; FILE *efield; gsl_integration_workspace *w; gsl_function f1; w = gsl_integration_workspace_alloc(1000); alpha = 2.0; // radius of sphere in units of L alpha2 = alpha*alpha; p[3] = 0.9; // R = 0.9 p[4] = alpha / sqrt(1.0 + p[3]*p[3]); // image charge effective magnitude efield = fopen("e_field.dat", "w"); if (!efield) { printf("can't create electric field data file. \n"); exit(0); } /* compute potential from geometry of coils */ for( k=0; k alpha2) { extotal = 0.0; // point outside sphere eytotal = 0.0; eztotal = 0.0; } else { f1.params = p; f1.function = & e1x; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ex, &error); f1.function = & e1y; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ey, &error); f1.function = & e1z; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ez, &error); extotal = ex; eytotal = ey; eztotal = ez; f1.function = & e2x; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ex, &error); f1.function = & e2y; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ey, &error); f1.function = & e2z; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ez, &error); extotal += ex; eytotal += ey; eztotal += ez; f1.function = & e3x; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ex, &error); f1.function = & e3y; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ey, &error); f1.function = & e3z; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ez, &error); extotal += ex; eytotal += ey; eztotal += ez; f1.function = & e4x; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ex, &error); f1.function = & e4y; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ey, &error); f1.function = & e4z; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ez, &error); extotal += ex; eytotal += ey; eztotal += ez; f1.function = & e5x; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ex, &error); f1.function = & e5y; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ey, &error); f1.function = & e5z; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ez, &error); extotal += ex; eytotal += ey; eztotal += ez; f1.function = & e6x; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ex, &error); f1.function = & e6y; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ey, &error); f1.function = & e6z; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &ez, &error); extotal += ex; eytotal += ey; eztotal += ez; } fwrite( &extotal, sizeof(double), 1, efield); fwrite( &eytotal, sizeof(double), 1, efield); fwrite( &eztotal, sizeof(double), 1, efield); } } } gsl_integration_workspace_free(w); fclose(efield); }