/********************************************************************************
*										*
*	Plot Magnetic field data for Bussard fusor.                           	*
*   Add electric field and electron density later.                              *
*                                                                               *
*        Author = Mike Rosing                                                   *
*         date  = Sept. 12, 2007                                                *
*										*
*********************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <png.h>
#include <math.h>

#define UNUMPOINTS 512
#define HNUMPOINTS 512
png_byte magplt[UNUMPOINTS][HNUMPOINTS][3];
png_bytep mptr[UNUMPOINTS];

#define MAXDIM 400
//#define MAXDIM 100

/*  MAXWALL is the ratio of outer containment sphere radius to L (center
    of sphere to center of coil distance)  */

#define MAXWALL  (3.0)

#define MAXSTEPS 100000
//#define NUMELMNTS 401*402*403*3/6
#define NUMELMNTS ((MAXDIM+1)*(MAXDIM+2)*(MAXDIM+3)*3/6)

/*  globals set in main only once  */

double dr;
double field[NUMELMNTS];

#define TSIZE (MAXDIM/2)

double vectors[TSIZE*TSIZE*TSIZE*3];

/*  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;		
}

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;
}


/*  subroutine to return offset into E or B field buffers.
    Enter with x, y, z position vector.
    Returns with offset to start of Fx, Fy, Fz data or -1 if
      input point is out of bounds.
*/

int bufroffset(double x, double y, double z, int *mapping)
{
    unsigned long mask1, mask2, mask3, mask;
    int offset, i, j, k;
    double r, theta, phi, xm, ym, zm;
    double dtheta, dphi;

/*  First map all octants back to zeroth octant  */

    if (x < 0) x = -x;
    if (y < 0) y = -y;
    if (z < 0) z = -z;

/*  find sextant within base octant.  Each bit position is 
    for a particular sextant  */

    if (x < y) mask1 = 0xE;
    else mask1 = 0x31;
    if (x < z) mask2 = 0x1C;
    else mask2 = 0x23;
    if( y < z) mask3 = 0x38;
    else mask3 = 0x7;
    mask = mask1 & mask2 & mask3;
    if( !mask)
    {
	printf("Error in bufroffset: x = %lf  y = %lf  z = %lf\n",
	       x, y, z);
	exit(0);
    }
    if( mask & 1)
    {
	xm = x;
	ym = y;
	zm = z;
	mapping[0] = 0;
	mapping[1] = 1;
	mapping[2] = 2;
    }
    else if (mask & 2)
    {
	xm = y;
	ym = x;
	zm = z;
	mapping[0] = 1;
	mapping[1] = 0;
	mapping[2] = 2;
    }
    else if (mask & 4)
    {
	xm = y;
	ym = z;
	zm = x;
	mapping[0] = 1;
	mapping[1] = 2;
	mapping[2] = 0;
    }
    else if (mask & 8)
    {
	xm = z;
	ym = y;
	zm = x;
	mapping[0] = 2;
	mapping[1] = 1;
	mapping[2] = 0;
    }
    else if (mask & 0x10)
    {
	xm = z;
	ym = x;
	zm = y;
	mapping[0] = 2;
	mapping[1] = 0;
	mapping[2] = 1;
    }
    else if (mask & 0x20)
    {
	xm = x;
	ym = z;
	zm = y;
	mapping[0] = 0;
	mapping[1] = 2;
	mapping[2] = 1;
    }
    r = sqrt(xm*xm + ym*ym + zm*zm);
    theta = atan2(ym, xm);
    phi = atan2(zm,  sqrt(xm*xm + ym*ym));
    i = lrint(r/dr);
    if (i > MAXDIM) return -1;
    dtheta = 4.0*i/M_PI;
    j = lrint(theta*dtheta);
    if (j > MAXDIM) return -1;
    dphi = dtheta;
    k = lrint(phi*dphi);
    if (k > MAXDIM) return -1;
    offset = (i*(i+1)*(i+2)/6 + j*(j+1)/2 + k)*3;
    return offset;
}

main()
{
	FILE  *readmag;
	int numread;
	int test, kmin, kmax, jmin, jmax;
	int i, j, k;
	char filename[64], input[8];
	double uangle, hangle, tlimit;
	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  ix, iy, iz, xb, yb, zb, index, xc, yc, zc;
	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, 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];
	double across, up;

/*  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;
	color[0] = 255;
	color[1] = 127;
	color[2] = 0;

	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(field, sizeof(double), NUMELMNTS, readmag);
    if(numread < NUMELMNTS)
	printf("\nWARNING: not all field file read in. \n\n");
    fclose(readmag);

    dr = MAXWALL/MAXDIM;

/*  test to see if data is where you think it is  

    index = 5;
    while(index > 0)
    {
	printf(" Enter position within volume.\n");
	printf("x: ");
	scanf("%lf", &x);
	printf("y: ");
	scanf("%lf", &y);
	printf("z: ");
	scanf("%lf", &z);
	index = bufroffset(x, y, z, map);
	printf("index = %d  map[0] = %d  map[1] = %d  map[2] = %d\n", index, map[0], map[1], map[2]);
	printf("bx = %lf  by = %lf  bz = %lf\n", field[index+map[0]], field[index+map[1]], 
	       field[index + map[2]]);
    }
    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);
*/
	uangle = 1.0;
	hangle = 1.0;
	f = 0.2;

	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);

	tlimit = TSIZE*dr;

/*  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[7], v[7], color);
	    drawline( h[0], v[0], h[4], v[4], white);
	    drawline( h[0], v[0], h[6], v[6], white);
	    drawline( h[0], v[0], h[5], v[5], white);
	    derr = 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], green);
	    drawline( h[1], v[1], h[5], v[5], white);
	    drawline( h[3], v[3], h[6], v[6], white);
	    drawline( h[3], v[3], h[5], v[5], white);
	    drawline( h[2], v[2], h[4], v[4], green);
	    drawline( h[2], v[2], h[6], v[6], white);
	    drawline( h[7], v[7], h[4], v[4], white);
	    drawline( h[7], v[7], h[5], v[5], white);
	    drawline( h[7], v[7], h[6], v[6], white);

	}

/*  compute vector length for large blocks and map to screen */

	vectormax = 0.0;
	numpoints = 0;
	for(i=0; i<TSIZE; i+=4)
	{
//	    printf("i= %d\n",i);
	    xb = i - TSIZE/2;  //i*8
//	    x = (i - TSIZE/2)*dr;
	    for(j=0; j<TSIZE; j+=4)
	    {
		yb = j - TSIZE/2; //*8
//		y = (j - TSIZE/2)*dr;
		for(k=0; k<TSIZE; k+=4)
		{
		    zb = k - TSIZE/2; //*8
//		    z = (k - TSIZE/2)*dr;
		    bx = 0.0;
		    by = 0.0;
		    bz = 0.0;
		    for(ix = 0; ix<4; ix++)
		    {
			x = (xb + ix)*dr*2.0;
//			xc = xb + ix;
			for(iy = 0; iy<4; iy++)
			{
			    y = (yb + iy)*dr*2.0;
			    for(iz = 0; iz<4; iz++)
			    {
				z = (zb + iz)*dr*2.0;
				index = bufroffset(x, y, z, map);
				if (index < 0) continue;
/*				printf("x = %lf y = %lf z = %lf index = %d  %d %d %d\n",
				       x, y, z, index, map[0], map[1], map[2]);
				printf("f[0] = %lf  f[1] = %lf f[2] = %lf\n", 
				       field[index+map[0]], field[index+map[1]], field[index+map[2]]);
				scanf("%s", input);
				
				zc = zb + iz;
				index = (xc*(xc+1)*(xc+2)/6 + yc*(yc+1)/2 + zc)*3;
				
*/
				numpoints++;
				bx += field[index+map[0]];
				by += field[index+map[1]];
				bz += field[index+map[2]];
			    }
			}
		    }
		    bx /= 4*4*4;
		    by /= 4*4*4;
		    bz /= 4*4*4;
		    vectors[3*(i + TSIZE*(j + TSIZE*k))] = bx;
		    vectors[3*(i + TSIZE*(j + TSIZE*k)) + 1] = by;
		    vectors[3*(i + TSIZE*(j + TSIZE*k)) + 2] = bz;

/*		    index = bufroffset(x, y, z, map);
		    if (index < 0) continue;
		    numpoints++;

		    bx = field[index+map[0]];
		    by = field[index+map[1]];
		    bz = field[index+map[2]];
		    vectors[3*(i + TSIZE*(j + TSIZE*k))] = bx;
		    vectors[3*(i + TSIZE*(j + TSIZE*k)) + 1] = by;
		    vectors[3*(i + TSIZE*(j + TSIZE*k)) + 2] = bz;
*/
		    vtst = sqrt(bx*bx + by*by + bz*bz);
		    if (vtst > vectormax) vectormax = vtst;
		}
	    }
	}
	across = M_PI/4.0;
	up = atan(1.0/sqrt(2.0));
	printf("found %d points\n", numpoints);
	vectormax /= scalevector;
	for(i=0; i<TSIZE; i++)
	{
	    x = (i - TSIZE/2)*dr*2.0;
//	    color[0] = 255 - 255/TSIZE*i;
	    color[0] = 255;

	    for(j=0; j<TSIZE; j++)
	    {
//		if (i) y = across/i * j;  //  angle from x axis
//		else y = 0.0;
		y = (j - TSIZE/2)*dr*2.0;
		color[1] = 255 - 255/TSIZE*j;
		for(k=0; k<TSIZE; k++)
		{
//		    if(j) z = up/i * k;  //  angle from x,y plane
//		    else z = 0.0;
		    z = (k - TSIZE/2)*dr*2.0;
		    color[2] = 255 - 255/TSIZE*k;
		    vectorhead[0] = x + //(x*cos(z)*cos(y)) +
			vectors[3*(i + TSIZE*(j + TSIZE*k))]/vectormax*2.0; 
		    vectortail[0] = x - //(x*cos(z)*cos(y)) - 
			vectors[3*(i + TSIZE*(j + TSIZE*k))]/vectormax*2.0; 
		    vectorhead[1] = y + //(x*cos(z)*sin(y)) +
			vectors[3*(i + TSIZE*(j + TSIZE*k)) + 1]/vectormax*2.0; 
		    vectortail[1] = y - //(x*cos(z)*sin(y)) -
			vectors[3*(i + TSIZE*(j + TSIZE*k)) + 1]/vectormax*2.0; 
		    vectorhead[2] = z + //(x*sin(z)) +
			vectors[3*(i + TSIZE*(j + TSIZE*k)) + 2]/vectormax*2.0; 
		    vectortail[2] =  z - //(x*sin(z)) -
			vectors[3*(i + TSIZE*(j + TSIZE*k)) + 2]/vectormax*2.0; 
		    if (transform3d( vectorhead, w, dptr, params, &h[0], &v[0]) &&
			transform3d( vectortail, w, dptr, params, &h[1], &v[1]))
		    {
//			printf(" p_0 = %d %d  p_1 = %d %d\n", h[0], v[0], h[1], v[1]);
			drawline( h[0], v[0], h[1], v[1], color);
		    }
//		    printf("vechead_x = %lf vechead_y = %lf vechead_z = %lf\n",
//			   vectorhead[0], vectorhead[1], vectorhead[2]);

		}
	    }
	}

	for( j=0; j<UNUMPOINTS; j++) mptr[j] = (png_bytep) magplt + j*3*HNUMPOINTS;

	write_png("eplot.png", (png_bytep*)mptr, UNUMPOINTS);
}