/******************************************************************************** * * * Plot Magnetic field data for Bussard fusor. * * Modified from plotfield.c to plot tombo_field.dat. * * Modified from plot_tombo_png.c to create 3D arrow plots for Blender. * * * * Author = Mike Rosing * * date = July 8, 2008 * * * *********************************************************************************/ #include #include #include #define MAXDIM 256 //#define MAXDIM 16 /* MAXWALL is the ratio of outer containment sphere radius to L (center of sphere to center of coil distance) */ #define MAXWALL (2.0) #define NUMELMNTS (MAXDIM*MAXDIM*MAXDIM) /* globals set in main only once */ #define TSIZE (MAXDIM/16) //#define TSIZE MAXDIM typedef struct { double x, y, z; } VECTOR; VECTOR Bfield[MAXDIM][MAXDIM][MAXDIM]; VECTOR vectors[TSIZE][TSIZE][TSIZE]; double vector_color[TSIZE][TSIZE][TSIZE]; static VECTOR arrow_def[7]; void init_arrow_def() { arrow_def[0].x = 3.24; arrow_def[0].y = 0.0; arrow_def[0].z = 0.0; arrow_def[1].x = 3.24 - 0.81; arrow_def[1].y = 0.309; arrow_def[1].z = 0.309; arrow_def[2].x = 3.24 - 0.81; arrow_def[2].y = -0.309; arrow_def[2].z = 0.309; arrow_def[3].x = 3.24 - 0.81; arrow_def[3].y = 0.309; arrow_def[3].z = -0.309; arrow_def[4].x = 3.24 - 0.81; arrow_def[4].y = -0.309; arrow_def[4].z = -0.309; arrow_def[5].x = 0.0; arrow_def[5].y = 0.309; arrow_def[5].z = 0.0; arrow_def[6].x = 0.0; arrow_def[6].y = -0.309; arrow_def[6].z = 0.0; } double det( double *a) { double t1, t2, t3, t4, t5, t6, t7; t1 = a[0] * a[4] * a[8]; t2 = a[1] * a[5] * a[6]; t3 = a[2] * a[3] * a[7]; t4 = a[2] * a[4] * a[6]; t5 = a[0] * a[5] * a[7]; t6 = a[1] * a[3] * a[8]; t7 = t1 + t2 + t3 - t4 - t5 - t6; return t7; } double absx(double x) { if (x < 0.0) return -x; return x; } /* draw an arrow in 3D. Enter with direction and location within volume. Also include color as magnitude reference. Returns 1 on success, 0 on failure (direction is zero length) */ int draw_arrow( FILE *arrowfile, VECTOR direction, VECTOR location, int color) { static int arrowdex=1; VECTOR arrow, dir[3]; int i, j, k, ldex; double len, x, y, z, scale; int path; dir[0].x = direction.x; dir[0].y = direction.y; dir[0].z = direction.z; if (absx(dir[0].x) > 1e-4) { dir[1].y = dir[0].y; dir[1].z = dir[0].z; dir[1].x = -(dir[0].y*dir[0].y + dir[0].z*dir[0].z)/dir[0].x; path = 0; if(absx(dir[1].x) < 1e-4) { dir[1].y = dir[0].x; dir[1].x = dir[0].y; dir[1].z = dir[0].z; path = 1; } } else if (absx(dir[0].y) > 1e-4) { dir[1].x = dir[0].x; dir[1].z = dir[0].z; dir[1].y = -(dir[0].x*dir[0].x + dir[0].z*dir[0].z)/dir[0].y; path = 2; if(absx(dir[1].y) < 1e-4) { dir[1].y = dir[0].z; dir[1].x = dir[0].x; dir[1].z = dir[0].y; path = 3; } } else if (absx(dir[0].z) > 1e-4) { dir[1].x = dir[0].x; dir[1].y = dir[0].y; dir[1].z = -(dir[0].x*dir[0].x + dir[0].y*dir[0].y)/dir[0].z; path = 4; if(absx(dir[1].z) < 1e-4) { dir[1].z = dir[0].x; dir[1].y = dir[0].y; dir[1].x = dir[0].z; path = 5; } } else return 0; // Nothing to plot since there is no direction /* renormalize direction vectors */ len = sqrt(dir[1].x*dir[1].x + dir[1].y*dir[1].y + dir[1].z*dir[1].z); dir[1].x /= len; dir[1].y /= len; dir[1].z /= len; dir[2].x = dir[0].y*dir[1].z - dir[0].z*dir[1].y; dir[2].y = dir[0].z*dir[1].x - dir[0].x*dir[1].z; dir[2].z = dir[0].x*dir[1].y - dir[0].y*dir[1].x; if((absx(dir[2].x < 1e-4) && (absx(dir[2].y) < 1e-4) && (absx(dir[2].z) < 1e-4))) return 0; /* compute rotation of arrow */ if (color > 255) color = 255; if(color < 0) color = 0; fprintf(arrowfile, "usemtl color%03d\n", color); for( ldex=0; ldex<7; ldex++) { arrow.x = dir[0].x*arrow_def[ldex].x + dir[1].x*arrow_def[ldex].y + dir[2].x*arrow_def[ldex].z; arrow.y = dir[0].y*arrow_def[ldex].x + dir[1].y*arrow_def[ldex].y + dir[2].y*arrow_def[ldex].z; arrow.z = dir[0].z*arrow_def[ldex].x + dir[1].z*arrow_def[ldex].y + dir[2].z*arrow_def[ldex].z; arrow.x += location.x; arrow.y += location.y; arrow.z += location.z; if((absx(arrow.x) > 180.0) || (absx(arrow.y) > 180.0) || (absx(arrow.z) > 180.0)) { printf("path = %d ldex = %d location = %f %f %f\n", path, ldex, location.x, location.y, location.z); printf("direction = %f %f %f\n", dir[0].x, dir[0].y, dir[0].z); printf("direction = %f %f %f\n", dir[1].x, dir[1].y, dir[1].z); printf("direction = %f %f %f\n", dir[2].x, dir[2].y, dir[2].z); printf("\n"); } /* output complete arrow to file */ fprintf( arrowfile, "v %f %f %f\n", arrow.x, arrow.y, arrow.z); } fprintf(arrowfile, "f %d %d %d\n", arrowdex, arrowdex+1, arrowdex+2); fprintf(arrowfile, "f %d %d %d\n", arrowdex, arrowdex+3, arrowdex+4); fprintf(arrowfile, "f %d %d %d\n", arrowdex, arrowdex+5, arrowdex+6); arrowdex += 7; return 1; } main() { FILE *readmag, *debug; int numread; int test, kmin, kmax, jmin, jmax; int i, j, k; char filename[64], input[8]; double uangle, hangle, tlimit, dr; double vectorhead[3], vectortail[3], bx, by, bz; int ix, iy, iz, xb, yb, zb, index, xc, yc, zc; double vtst, x, y, z; int numpoints, h[8], v[8], red[3], green[3], blue[3], white[3]; double params[6], w[3], vectormax, scalevector; int map[3], derr, color[3]; VECTOR color_map[256], location; /* create color arrow map */ for(i=0; i<43; i++) { color_map[i].x = 1.0; color_map[i].y = i/43.0; color_map[i].z = 0.0; } for(i=43; i<85; i++) { color_map[i].x = 1.0 - (i-43)/42.0; color_map[i].y = 1.0; color_map[i].z = 0.0; } for(i=85; i<128; i++) { color_map[i].x = 0.0; color_map[i].y = 1.0; color_map[i].z = (i - 85)/43.0; } for(i=128; i<170; i++) { color_map[i].x = 0.0; color_map[i].y = 1.0 - (i - 128)/42.0; color_map[i].z = 1.0; } for(i=170; i<213; i++) { color_map[i].x = (i - 170)/43.0; color_map[i].y = 0.0; color_map[i].z = 1.0; } for(i=213; i<256; i++) { color_map[i].x = 1.0; color_map[i].y = 0.0; color_map[i].z = 1.0 - (i - 213)/43.0; } /* output color table for arrows */ sprintf(filename, "arrow_color_map"); readmag = fopen(filename, "w"); if(!readmag) { printf("can't create arrow color map file\n"); exit(0); } for(i=0; i<256; i++) { fprintf(readmag, "newmtl color%03d\n", i); fprintf(readmag, "Ka %f %f %f\n", color_map[i].x, color_map[i].y, color_map[i].z); fprintf(readmag, "Kd %f %f %f\n", 1.0-color_map[i].x, 1.0-color_map[i].y, 1.0-color_map[i].z); fprintf(readmag, "illum1\n"); } fclose(readmag); init_arrow_def(); /* read in data from magnetic field file */ printf("name of field file data? "); scanf( "%s", filename); readmag = fopen(filename, "rb"); if( !readmag) { printf("can't read in mag data file %s.\n", filename); exit(0); } printf("opening data files and reading them in.\n"); numread = fread(Bfield, sizeof(VECTOR), NUMELMNTS, readmag); if(numread < NUMELMNTS) printf("\nWARNING: not all field file read in. \n\n"); fclose(readmag); dr = 2.0*MAXWALL/MAXDIM; /* compute vector length for large blocks and compute color of associated arrow */ vectormax = 0.0; numpoints = 0; for(i=0; i> 4; for(j=0; j> 4; for(k=0; k> 4; bx = 0.0; by = 0.0; bz = 0.0; for(ix = 0; ix<16; ix++) { for(iy = 0; iy<16; iy++) { for(iz = 0; iz<16; iz++) { numpoints++; bx += Bfield[i+ix][j+iy][k+iz].x; by += Bfield[i+ix][j+iy][k+iz].y; bz += Bfield[i+ix][j+iy][k+iz].z; } } } if(isnan(bx) || isnan(by) || isnan(bz)) printf("%d %d %d\n", xb, yb, zb); // bx /= 16*16*16; // by /= 16*16*16; // bz /= 16*16*16; vtst = sqrt(bx*bx + by*by + bz*bz); if(vtst > 1e-5) { vectors[xb][yb][zb].x = bx/vtst; vectors[xb][yb][zb].y = by/vtst; vectors[xb][yb][zb].z = bz/vtst; if (vtst > vectormax) vectormax = vtst; vector_color[xb][yb][zb] = vtst; } else { vectors[xb][yb][zb].x = 0.0; vectors[xb][yb][zb].y = 0.0; vectors[xb][yb][zb].z = 0.0; vector_color[xb][yb][zb] = 0.0; } } } } printf("found %d points\n", numpoints); printf("vectormax = %f\n", vectormax); printf("name of picture file: "); scanf( "%s", filename); readmag = fopen(filename, "w"); if( !readmag) { printf("can't create in picture data file %s.\n", filename); exit(0); } fprintf(readmag, "mtllib arrow_color_map\n"); /* debug = fopen("debug.dat", "w"); if( !debug) { printf("can't create debug file\n"); exit(0); }*/ for(i=0; i