/********************************************************** * * * Compute force volume in Vlasov equation for * * a few specific points in 6 coil volume. * * * * Author = Mike Rosing * * date = Oct. 17, 2007 * * * *********************************************************/ #include #include main() { FILE *magdat, *efld, *forcedat; double bx, by, bz, ux, uy, uz, fx, fy, fz; double ex, ey, ez, volts, amps, me, ce, param; int i, j, k, offset; char st[2]; me = 9.10938189e-31; // mass electron ce = 1.602176462e-19; // charge of electron magdat = fopen("dipole.dat", "r"); if(!magdat) { printf("can't open mag data\n"); exit(0); } efld = fopen("e_field.dat", "r"); if(!magdat) { printf("can't open e_field data\n"); exit(0); } st[0] = 'n'; while (st[0] == 'n') { printf("input i offset: "); scanf("%d", &i); printf("input j offset: "); scanf("%d", &j); printf("input k offset: "); scanf("%d", &k); offset = 3*sizeof(double)*(i + 256*(j + 256*k)); // printf("desired set point is %d\n", offset); rewind(magdat); fseek( magdat, offset, SEEK_SET); rewind(efld); fseek( efld, offset, SEEK_SET); // offset = ftell(magdat); // printf("actuall set point is %d\n", offset); fread( (double*)&bx, sizeof(double), 1, magdat); fread( (double*)&by, sizeof(double), 1, magdat); fread( (double*)&bz, sizeof(double), 1, magdat); fread( (double*)&ex, sizeof(double), 1, efld); fread( (double*)&ey, sizeof(double), 1, efld); fread( (double*)&ez, sizeof(double), 1, efld); printf(" mag field values are bx = %lg by = %lg bz = %lg\n", bx, by, bz); printf(" e field values are ex = %lg ey = %lg ez = %lg\n", ex, ey, ez); printf("hit n to retry, y to continue: "); scanf("%s", st); } /* now compute u vector space and u X B */ forcedat = fopen("force.dat", "w"); if (! forcedat) { printf("can't create force file\n"); exit(0); } printf("Voltage on grid: "); scanf ("%lf", &volts); printf("Amp turns in grid: "); scanf("%lf", &s); param = sqrt(2.0*me*volts/ce)/amps/(2.0e-7); printf("containment parameter = %lf\n", param); for(k=0; k<128; k++) { uz = k/32.0 - 2.0; for(j=0; j<128; j++) { uy = j/32.0 - 2.0; for(i=0; i<128; i++) { ux = i/32.0 - 2.0; fx = param*ex + uy*bz - uz*by; fy = param*ey + uz*bx - ux*bz; fz = param*ez + ux*by - uy*bx; fwrite( &fx, sizeof(double), 1, forcedat); fwrite( &fy, sizeof(double), 1, forcedat); fwrite( &fz, sizeof(double), 1, forcedat); } } } fclose(forcedat); }