/******************************************************************* * * * The purpose of this program is to compute the E and B * * fields for use in the particle tracking code. The fields are * * stored in spherical coordinates. They are unusual coordinates * * with "theta" being in the x,y plane and "phi" going up from * * the x,y plane to the point at z. The zeroth quadrant of this * * spherical system holds all the E and B data, and Indrek's * * vector form is used to find the lookup point within that block.* * * * The B and E field computation code is copied from mag_6.c and * * efield.c, but the dimensions of the data are different. * * * * Author = Mike Rosing * * date = Jan. 1, 2008 * * * ******************************************************************/ #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 b1x (double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux - 1.0; ry = uy - R*c1; rz = uz - R*s1; t1 = c1*ry + s1*rz; // printf("t1=%lf\n", t1); t2 = rx*rx + ry*ry + rz*rz; // printf("phi = %lf ux = %lf uy = %lf uz = %lf R = %lf\n", // phi, ux, uy, uz, R); // scanf("%lf", &R); if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; // printf("f = %lf\n", f); return f; } double b1y ( double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux - 1.0; ry = uy - R*c1; rz = uz - R*s1; t1 = -c1*rx; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } double b1z ( double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux - 1.0; ry = uy - R*c1; rz = uz - R*s1; t1 = -s1*rx; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } /* coil 2 on negative x axis*/ double b2x (double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux + 1.0; ry = uy - R*c1; rz = uz - R*s1; t1 = -(c1*ry + s1*rz); t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } double b2y ( double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux + 1.0; ry = uy - R*c1; rz = uz - R*s1; t1 = c1*rx; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } double b2z ( double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux + 1.0; ry = uy - R*c1; rz = uz - R*s1; t1 = s1*rx; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } /* coil 3 on positive y axis */ double b3x (double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux + R*c1; ry = uy - 1.0; rz = uz - R*s1; t1 = c1*ry; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } double b3y ( double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux + R*c1; ry = uy - 1.0; rz = uz - R*s1; t1 = -c1*rx + s1*rz; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } double b3z ( double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux + R*c1; ry = uy - 1.0; rz = uz - R*s1; t1 = -s1*ry; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } /* coil 4 on negative y axis */ double b4x (double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux + R*c1; ry = uy + 1.0; rz = uz - R*s1; t1 = -c1*ry; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } double b4y ( double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux + R*c1; ry = uy + 1.0; rz = uz - R*s1; t1 = c1*rx - s1*rz; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } double b4z ( double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux + R*c1; ry = uy + 1.0; rz = uz - R*s1; t1 = s1*ry; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } /* coil 5 on positive z axis */ double b5x (double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux - R*c1; ry = uy - R*s1; rz = uz - 1.0; t1 = -c1*rz; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } double b5y ( double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux - R*c1; ry = uy - R*s1; rz = uz - 1.0; t1 = -s1*rz; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } double b5z ( double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux - R*c1; ry = uy - R*s1; rz = uz - 1.0; t1 = s1*ry + c1*rx; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } /* coil 6 on negative z axis */ double b6x (double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux - R*c1; ry = uy - R*s1; rz = uz + 1.0; t1 = c1*rz; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } double b6y ( double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux - R*c1; ry = uy - R*s1; rz = uz + 1.0; t1 = s1*rz; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } double b6z ( double phi, void *params) { double s1, c1, rx, ry, rz, f; double t1, t2, t3, *ptr; double ux, uy, uz, R; ptr = (double *)params; ux = ptr[0]; uy = ptr[1]; uz = ptr[2]; R = ptr[3]; s1 = sin(phi); c1 = cos(phi); rx = ux - R*c1; ry = uy - R*s1; rz = uz + 1.0; t1 = -s1*ry - c1*rx; t2 = rx*rx + ry*ry + rz*rz; if (t2 < TOOSMALL) return 0.0; t3 = sqrt(t2); f = t1/t2/t3; return f; } /* 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; double p[5], ex, ey, ez; double alpha, alpha2, extotal, eytotal, eztotal; FILE *efield; double r, theta, phi, across, up; double bx, by, bz, error; double bxtotal, bytotal, bztotal; FILE *bfield; gsl_integration_workspace *w; gsl_function f1; w = gsl_integration_workspace_alloc(1000); /* generate both magnetic field and electric field. Parameter passing uses first 5 values p[0] -> p[3] for both B and E, but E requires one more: p[0] = ux p[1] = uy p[2] = uz p[3] = coil radius ratio p[4] = image charge effective magnitude */ 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 bfield = fopen("b_field.dat", "w"); if (!bfield) { printf("can't create mag field data file. \n"); exit(0); } efield = fopen("e_field.dat", "w"); if (!efield) { printf("can't create electric field data file. \n"); exit(0); } across = M_PI/4.0; // x, y plane angle up = across; // z up from xy plane angle for( i=0; i