/******************************************************************************** * * * Compute magnetic 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. * * * * Current density and effective length are combined to create unitless * * magnetic field. Actual field is scaled using: * * * * B = (mu_0 * I_0 )/( 4.0 * PI * L) * b * * * * where b is what this program computes, I_0 is the coil current, L is the * * distance from center of fusor to center of coil and * * mu_0 is 4*PI*1e-7 kg m / C^2 * * * * Author = Mike Rosing * * date = Sept. 8, 2007 * * * *******************************************************************************/ #include #include #include #include /* Need one subroutine for each integral to be performed. See notes for derivations. */ #define TOOSMALL (1e-6) /* 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; } main() { int i, j, k; double p[4], bx, by, bz, error; double bxtotal, bytotal, bztotal; FILE *bfield; gsl_integration_workspace *w; gsl_function f1; w = gsl_integration_workspace_alloc(1000); p[3] = 0.9; // R = 0.9 bfield = fopen("dipole.dat", "w"); if (!bfield) { printf("can't create mag field data file. \n"); exit(0); } for( k=0; k<256; k++) { printf("%d.\n", k); p[2] = 2.0/256.0*k ; // uz = 2.0/256 * k for( j=0; j<256; j++) { p[1] = 2.0/256.0*j; // uy = 2.0/256 * j for( i=0; i<256; i++) { p[0] = 2.0/256.0 * i ; // ux = 2.0/256 * i f1.params = p; f1.function = & b1x; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &bx, &error); f1.function = & b1y; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &by, &error); f1.function = & b1z; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &bz, &error); bxtotal = bx; bytotal = by; bztotal = bz; f1.function = & b2x; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &bx, &error); f1.function = & b2y; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &by, &error); f1.function = & b2z; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &bz, &error); bxtotal += bx; bytotal += by; bztotal += bz; f1.function = & b3x; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &bx, &error); f1.function = & b3y; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &by, &error); f1.function = & b3z; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &bz, &error); bxtotal += bx; bytotal += by; bztotal += bz; f1.function = & b4x; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &bx, &error); f1.function = & b4y; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &by, &error); f1.function = & b4z; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &bz, &error); bxtotal += bx; bytotal += by; bztotal += bz; f1.function = & b5x; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &bx, &error); f1.function = & b5y; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &by, &error); f1.function = & b5z; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &bz, &error); bxtotal += bx; bytotal += by; bztotal += bz; f1.function = & b6x; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &bx, &error); f1.function = & b6y; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &by, &error); f1.function = & b6z; gsl_integration_qag( &f1, 0.0, 2.0*M_PI, 1e-6, 1e-5, 1000, GSL_INTEG_GAUSS41, w, &bz, &error); bxtotal += bx; bytotal += by; bztotal += bz; bxtotal *= p[3]; // NOTE: forgot R * dphi in integral bytotal *= p[3]; // but constant comes outside ok bztotal *= p[3]; fwrite( &bxtotal, sizeof(double), 1, bfield); fwrite( &bytotal, sizeof(double), 1, bfield); fwrite( &bztotal, sizeof(double), 1, bfield); } } } gsl_integration_workspace_free(w); fclose(bfield); }