/******************************************************************************** * * * Plot Magnetic field data for Bussard fusor. * * Add electric field and electron density later. * * * * Author = Mike Rosing * * date = Sept. 12, 2007 * * * *********************************************************************************/ #include #include #include #include #define UNUMPOINTS 512 #define HNUMPOINTS 512 png_byte magplt[UNUMPOINTS][HNUMPOINTS][3]; png_bytep mptr[UNUMPOINTS]; /* This subroutine writes a png data structure to disk. Enter with file name, and pointer to data as well as row pointer data. Returns 1 if successful, 0 if errors. */ int write_png( char *filename, png_bytep row_ptr[], int rowdepth) { FILE *pngfile; png_structp png_ptr; png_infop info_ptr; int check, chk; /* follow png write file commands from example.c */ pngfile = fopen( filename, "wb"); if( pngfile == NULL) { printf(" Can't create png file %s.\n", filename); exit(0); } /* create png_struct. I do as little as possible! */ png_ptr = png_create_write_struct( PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); if( png_ptr == NULL) { printf("Can't create png structure. \n"); fclose( pngfile); } /* allocate and initialize info structure */ info_ptr = png_create_info_struct( png_ptr); if( info_ptr == NULL) { fclose( pngfile); png_destroy_write_struct( &png_ptr, (png_infopp) NULL); printf(" Can't create info pointer.\n"); exit(0); } /* haven't a clue what this does .... */ if( setjmp(png_ptr->jmpbuf)) { fclose( pngfile); printf("Failure in setjmp (file reading?)\n"); png_destroy_write_struct( &png_ptr, (png_infopp) NULL); exit(0); } /* set up output control */ png_init_io( png_ptr, pngfile); /* set up image form. 8 bit RGB for now */ png_set_IHDR( png_ptr, info_ptr, rowdepth, rowdepth, 8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE); /* write file header info */ png_write_info( png_ptr, info_ptr); /* write image info */ png_write_image( png_ptr, row_ptr); printf("wrote image data\n"); /* finish file data */ png_write_end( png_ptr, info_ptr); /* eliminate all data allocated on heap */ png_destroy_write_struct( &png_ptr, (png_infopp) NULL); fclose(pngfile); } /* Add color with maximum saturation. Enter with pointer to RGB triple, and red, green, blue new values to sum into it. */ void setcolor( png_byte rgb[3], int red, int green, int blue) { if( rgb[0] + red > 255) rgb[0] = 255; else rgb[0] += (png_byte) red; if( rgb[1] + green > 255) rgb[1] = 255; else rgb[1] += (png_byte)green; if( rgb[2] + blue > 255) rgb[2] = 255; else rgb[2] += (png_byte) blue; } /* show the points of tau actually selected. This uses global data and writes global data. For now, red is imaginary axis marker and green is real axis marker. Find limits of tau, convert values to integers and color each point. */ /* void settau() { FLOAT imin, imax, rmin, rmax; FLOAT ifactr, rfactr, constant, temp; INDEX mi, mr, j, k; initminmax( &imin, &imax); initminmax( &rmin, &rmax); for( j=0; j 255) mr = 255; if( mi > 255) mi = 255; setcolor( tauplt[mi][mr], 63 + k*(192/tlimit) , 255 - j*(192/tlimit), 0); } } for( j=0; j<256; j++) tptr[j] = (png_bytep) tauplt + j*3*256; } */ /* compute a 3x3 determinant. Enter with "a" pointing to a 9 element array of doubles in order 11, 12, 13, 21, 22, 23, 31, 32, 33. 0 1 2 3 4 5 6 7 8 */ 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; } /* Transform a 3D position to a 2D screen position. Input endpoint "w", "dptr" matrix, point to be transformed "p" and parameters qmax, rmax, smax and determinant of dptr "divisor" as params[0], params[1], params[2], params[3], params[4] = tan(uangle), params[5] = tan(hangle). outputs integer h, v referenced to center of screen for horizontal and vertical offset. returns 1 if point can be plotted, 0 if outside view volume. */ int transform3d ( double *plot, double *w, double *dptr, double *params, int *h, int *v) { double b[3], cptr[9], rtst, stst; double q, r, s; int j; b[0] = plot[0] - w[0]; b[1] = plot[1] - w[1]; b[2] = plot[2] - w[2]; for(j=0; j<9; j++) cptr[j] = dptr[j]; cptr[0] = b[0]; cptr[3] = b[1]; cptr[6] = b[2]; q = det(cptr)/params[3]; for(j = 0; j<9; j++) cptr[j] = dptr[j]; cptr[1] = b[0]; cptr[4] = b[1]; cptr[7] = b[2]; r = det(cptr)/params[3]; if (r < 0.0) rtst = -r; else rtst = r; for(j = 0; j<9; j++) cptr[j] = dptr[j]; cptr[2] = b[0]; cptr[5] = b[1]; cptr[8] = b[2]; s = det(cptr)/params[3]; if (s < 0.0) stst = -s; else stst = s; /* printf("q = %lf r = %lf s = %lf\n", q, r, s); printf("params[4] = %lf params[5] = %lf\n", params[4], params[5]); printf("params[1] = %lf params[2] = %lf\n", params[1], params[2]); */ // do testing here and print results!!@! if( q > params[0] || q < 0.0) { // printf("point %d not in volume q = %lf qmax = %lf\n", // i, q[i], qmax); return 0; } if( rtst > params[1] + (params[0] - q)*params[4]) { // printf("point %d not in volume r = %lf rmax = %lf\n", // i, r[i], rmax); return 0; } if( stst > params[2] + (params[0] - q)*params[5]) { // printf("point %d not in volume s = %lf smax = %lf\n", // i, s[i], smax); return 0; } *h = s/params[2] * HNUMPOINTS/2; *v = r/params[1] * UNUMPOINTS/2; return 1; } /* draw a line in magplt buffer. Enter with end points and color. */ int drawline( int h1, int v1, int h2, int v2, int color[3]) { int vcount, hcount, vcheck, hcheck; int checklength, mark; int hdex, vdex, start, end; double vfrac, hfrac, slope, bcept; /* translate input from center reference to screen reference */ h1 = HNUMPOINTS/2 + h1; h2 = HNUMPOINTS/2 + h2; v1 = UNUMPOINTS/2 - v1; v2 = UNUMPOINTS/2 - v2; /* see if vector is long enough to draw */ vcount = v2 - v1; hcount = h2 - h1; // printf("h1 = %d v1 = %d h2 = %d v2 = %d\n", h1, v1, h2, v2); if ((v1 < 0) && (v2 < 0)) return 0; if ((h1 < 0) && (h2 < 0)) return 0; if ((h1 > HNUMPOINTS) && (h2 > HNUMPOINTS)) return 0; if ((v1 > UNUMPOINTS) && (v2 > UNUMPOINTS)) return 0; checklength = vcount*vcount + hcount*hcount; if (checklength < 1) return 0; /* draw a line between head and tail onscreen */ if (vcount < 0) vcheck = -vcount; else vcheck = vcount; if (hcount < 0) hcheck = -hcount; else hcheck = hcount; if (vcheck > UNUMPOINTS) return 0; if (hcheck > HNUMPOINTS) return 0; if (hcheck >= vcheck) { if (h1 < h2) { start = h1; end = h2; } else { start = h2; end = h1; } if(start < 0) start = 0; if(end >= HNUMPOINTS) end = HNUMPOINTS -1; vfrac = v2 - v1; hfrac = h2 - h1; slope = vfrac/hfrac; bcept = v1 - slope*h1; for (hdex=start; hdex<=end; hdex++) { vdex = slope*(double)hdex + bcept; if ((vdex > UNUMPOINTS) || (vdex < 0)) break; /* if(hdex < start+2) { magplt[vdex][hdex][0] = 255; magplt[vdex][hdex][1] = 0; magplt[vdex][hdex][2] = 0; } else */ { magplt[vdex][hdex][0] = color[0]; magplt[vdex][hdex][1] = color[1]; magplt[vdex][hdex][2] = color[2]; } } } else { if (v1 < v2) { start = v1; end = v2; } else { start = v2; end = v1; } vfrac = v2 - v1; hfrac = h2 - h1; slope = hfrac/vfrac; bcept = h1 - slope*v1; if (end >= UNUMPOINTS) end = UNUMPOINTS - 1; if (start < 0) start = 0; for(vdex=start; vdex<=end; vdex++) { hdex = slope*(double)vdex + bcept; if( (hdex > HNUMPOINTS) || (hdex < 0)) break; /* if(vdex < start+4) { magplt[vdex][hdex][0] = 255; magplt[vdex][hdex][1] = 0; magplt[vdex][hdex][2] = 0; } else */ { magplt[vdex][hdex][0] = color[0]; magplt[vdex][hdex][1] = color[1]; magplt[vdex][hdex][2] = color[2]; } } } return 1; } main() { FILE *readmag; int numread, tlimit; int test, kmin, kmax, jmin, jmax; int i, j, k; char filename[64], input[8]; double uangle, hangle; double *magblock, *vectors; double vectorhead[3], vectortail[3], bx, by, bz; double ex, ey, ez, f, cx, cy, cz, ax, ay, az; double wx, wy, wz, ux, uy, uz, hx, hy, hz; double ulen, uscale, eyelength, hscale; int xb, yb, zb, ix, iy, iz; double px, py, pz, eyefrac; double amin[3], divisor, dptr[9], cptr[9]; double cprmx, cprmy, cprmz, uprmx, uprmy, uprmz, hprmx, hprmy, hprmz; double ptest[8][3], q[8], r[8], s[8], qmax, rmax, smax; double vtst; int numpoints, h[8], v[8], red[3], green[3], blue[3], white[3]; double params[6], w[3], vectormax, scalevector; tlimit = 256; magblock = (double *)calloc(256*256*256, 3*sizeof(double)); if( magblock == NULL) { printf("failed to allocate memory space \n"); exit(0); } vectors = (double *)calloc(256*256*256, 3*sizeof(double)); if( vectors == NULL) { printf("failed to allocate vectors memory space \n"); exit(0); } /* set up some simple color references */ red[0] = 255; red[1] = 0; red[2] = 0; green[0] = 0; green[1] = 255; green[2] = 0; blue[0] = 0; blue[1] = 0; blue[2] = 255; white[0] = 255; white[1] = 255; white[2] = 255; printf("name of magnetic file data? "); scanf( "%s", filename); readmag = fopen(filename, "rb"); if( !readmag) { printf("can't read in mag data file %s.\n", filename); exit(0); } /* gather view angle data */ printf("enter eye position point:\n"); printf("x: "); scanf("%lf", &ex); printf("y: "); scanf("%lf", &ey); printf("z: "); scanf("%lf", &ez); printf("enter end of viewing volume point: \n"); printf("x: "); scanf("%lf", &wx); printf("y: "); scanf("%lf", &wy); printf("z: "); scanf("%lf", &wz); w[0] = wx; w[1] = wy; w[2] = wz; printf("theta up: "); scanf("%lf", &uangle); printf("theta horizontal: "); scanf("%lf", &hangle); printf("fractional distance to screen: "); scanf("%lf", &f); printf("Vector length scale factor: "); scanf("%lf", &scalevector); /* compute geometric factors */ ax = ex - wx; ay = ey - wy; az = ez - wz; eyelength = sqrt( ax*ax + ay*ay + az*az); printf("ax = %lf ay = %lf az = %lf\n", ax, ay, az); cx = (1.0 - f)*ax; cy = (1.0 - f)*ay; cz = (1.0 - f)*az; printf("cx = %lf cy = %lf cz = %lf\n", cx, cy, cz); /* determine "up" based on smallest direction cosine */ if (ax < 0) amin[0] = -ax; else amin[0] = ax; if (ay < 0) amin[1] = -ay; else amin[1] = ay; if (az < 0) amin[2] = -az; else amin[2] = az; if ((amin[2] <= amin[0]) && (amin[2] <= amin[1])) { ulen = ax*ax + ay*ay; uscale = f*tan(uangle)/sqrt(ulen); ux = -ax*az*uscale; uy = -ay*az*uscale; uz = ulen*uscale; hscale = f*tan(hangle)/sqrt(ulen)*eyelength; hx = ay*hscale; hy = -ax*hscale; hz = 0.0; } else if ((amin[0] <= amin[1]) && (amin[0] <= amin[2])) { ulen = ay*ay + az*az; uscale = f*tan(uangle)/sqrt(ulen); ux = ulen*uscale; uy = -ax*ay*uscale; uz = - ax*az*uscale; hscale = f*tan(hangle)/sqrt(ulen)*eyelength; hx = 0.0; hy = az*hscale; hz = -ay*hscale; } else if ((amin[1] <= amin[0]) && (amin[1] <= amin[2])) { ulen = ax*ax + az*az; uscale = f*tan(uangle)/sqrt(ulen); ux = -ax*ay*uscale; uy = ulen*uscale; uz = -ay*az*uscale; hscale = f*tan(hangle)/sqrt(ulen)*eyelength; hx = -az*hscale; hy = 0.0; hz = ax*hscale; } else { printf("sorting test failed!!\n"); printf("amin[0] = %lf amin[1] = %lf amin[2] = %lf\n", amin[0], amin[1], amin[2]); exit(0); } /* compute view volume unit vectors. w is origin, c is screen. */ divisor = sqrt(cx*cx + cy*cy + cz*cz); qmax = divisor; params[0] = qmax; cprmx = cx/divisor; cprmy = cy/divisor; cprmz = cz/divisor; printf("cprm x = %lf y = %lf z = %lf\n", cprmx, cprmy, cprmz); divisor = sqrt(ux*ux + uy*uy + uz*uz); rmax = divisor; params[1] = rmax; uprmx = ux/divisor; uprmy = uy/divisor; uprmz = uz/divisor; printf("uprm x = %lf y = %lf z = %lf\n", uprmx, uprmy, uprmz); divisor = sqrt(hx*hx + hy*hy + hz*hz); smax = divisor; params[2] = smax; hprmx = hx/divisor; hprmy = hy/divisor; hprmz = hz/divisor; printf("hprm x = %lf y = %lf z = %lf\n", hprmx, hprmy, hprmz); dptr[0] = cprmx; dptr[1] = uprmx; dptr[2] = hprmx; dptr[3] = cprmy; dptr[4] = uprmy; dptr[5] = hprmy; dptr[6] = cprmz; dptr[7] = uprmz; dptr[8] = hprmz; divisor = det(dptr); params[3] = divisor; if (divisor < 1e-6) { printf("view volume determinant too small!!\n"); printf("divisor = %lf\n\n", divisor); exit(0); } params[4] = tan(uangle); params[5] = tan(hangle); /* see if any portion of object space is in view volume */ ptest[0][0] = 0.0; ptest[0][1] = 0.0; ptest[0][2] = 0.0; ptest[1][0] = tlimit; ptest[1][1] = 0.0; ptest[1][2] = 0.0; ptest[2][0] = 0.0; ptest[2][1] = tlimit; ptest[2][2] = 0.0; ptest[3][0] = 0.0; ptest[3][1] = 0.0; ptest[3][2] = tlimit; ptest[4][0] = tlimit; ptest[4][1] = tlimit; ptest[4][2] = 0.0; ptest[5][0] = tlimit; ptest[5][1] = 0.0; ptest[5][2] = tlimit; ptest[6][0] = 0.0; ptest[6][1] = tlimit; ptest[6][2] = tlimit; ptest[7][0] = tlimit; ptest[7][1] = tlimit; ptest[7][2] = tlimit; numpoints = 0; for(i=0; i<8; i++) { if (transform3d ( ptest[i], w, dptr, params, &h[i], &v[i])) { printf("point %d is in volume\n", i); numpoints++; } } if (numpoints == 8) { drawline( h[0], v[0], h[1], v[1], red); drawline( h[0], v[0], h[2], v[2], green); drawline( h[0], v[0], h[3], v[3], blue); /* drawline( h[1], v[1], h[4], v[4]); drawline( h[1], v[1], h[5], v[5]); drawline( h[3], v[3], h[6], v[6]); drawline( h[3], v[3], h[5], v[5]); drawline( h[2], v[2], h[4], v[4]); drawline( h[2], v[2], h[6], v[6]); drawline( h[7], v[7], h[4], v[4]); drawline( h[7], v[7], h[5], v[5]); drawline( h[7], v[7], h[6], v[6]); */ } /* note first test - bx only */ numread = fread( magblock, 3*sizeof(double), tlimit*tlimit*tlimit , readmag); if( numread < tlimit*tlimit*tlimit) printf("Warning, did not read all of tau file.\n"); /* compute vector length for large blocks and map to screen */ vectormax = 0.0; for(i=0; i vectormax) vectormax = vtst; } } } vectormax /= scalevector; for(i=0; i