/********************************************************************************
*										*
*	Plot Magnetic or electric field data for Bussard fusor.               	*
*    Modified from plotfield.c to include 3D arrows for *.obj type files.       *
*                                                                               *
*        Author = Mike Rosing                                                   *
*         date  = July. 12, 2008                                                *
*										*
*********************************************************************************/

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

//#define MAXDIM 400
#define MAXDIM 25

/*  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 ((MAXDIM+1)*(MAXDIM+2)*(MAXDIM+3)/6)

typedef struct
{
    double x, y, z;
} VECTOR;

/*  globals set in main only once  */

double dr;
VECTOR mag_field[NUMELMNTS], elc_field[NUMELMNTS];
double magmag[NUMELMNTS], elcmag[NUMELMNTS];

static VECTOR arrow_def[7];

void init_arrow_def()
{
  arrow_def[0].x = 3.24 - 1.62;
  arrow_def[0].y = 0.0;
  arrow_def[0].z = 0.0;
  arrow_def[1].x = 3.24 - 0.81 - 1.62;
  arrow_def[1].y = 0.309;
  arrow_def[1].z = 0.309;
  arrow_def[2].x = 3.24 - 0.81 - 1.62;
  arrow_def[2].y = -0.309;
  arrow_def[2].z = 0.309;
  arrow_def[3].x = 3.24 - 0.81 - 1.62;
  arrow_def[3].y = 0.309;
  arrow_def[3].z = -0.309;
  arrow_def[4].x = 3.24 - 0.81 - 1.62;
  arrow_def[4].y = -0.309;
  arrow_def[4].z = -0.309;
  arrow_def[5].x = 0.0 - 1.62;
  arrow_def[5].y = 0.309;
  arrow_def[5].z = 0.0;
  arrow_def[6].x = 0.0 - 1.62;
  arrow_def[6].y = -0.309;
  arrow_def[6].z = 0.0;
}

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, int elcmag)
{
  static int arrowdex=1;
  VECTOR arrow, dir[3];
  int i, j, k, ldex;
  double len, x, y, z, scale;
  int path;

  if(color < 20) scale = color/20.0;
  else scale = 1.0;

  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;

  for(i=0; i<3; i++)
  {
    dir[i].x *= scale;
    dir[i].y *= scale;
    dir[i].z *= scale;
  }
/*  compute rotation of arrow  */

  if (color > 127) color = 127;
  if(color < 0) color = 0;
  if(elcmag)
    fprintf(arrowfile, "usemtl elc%03d\n", color);
  else
    fprintf(arrowfile, "usemtl mag%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, *arrowfile;
  int numread;
  int i, j, k;
  char filename[64];
  int  ix, iy, iz, xb, yb, zb, index, xc, yc, zc;
  double vtst, x, y, z, bx, by, bz;
  double r, theta, phi;
  int numpoints;
  double elcvectormax, magvectormax, scalevector;
  VECTOR color_map, plot_dir, plot_loc;
  double across, up;

/*  generate different colors for magnetic and electric fields  */

  arrowfile = fopen("electromag_color", "w");
  if(!arrowfile)
  {
    printf("can't create electromag_color file\n");
    exit(0);
  }
  fprintf(arrowfile, "#Electromagnetic color data\n");
  for(i=0; i<64; i++)
  {
    color_map.x = 0.0;
    color_map.y = i/64.0;
    color_map.z = 1.0;
    fprintf(arrowfile, "newmtl elc%03d\n", i);
    fprintf(arrowfile, "Ka %f  %f  %f\n", color_map.x, color_map.y, color_map.z);
    fprintf(arrowfile, "Kd %f  %f  %f\n", 1.0-color_map.x, 
	    1.0-color_map.y, 1.0-color_map.z);
    fprintf(arrowfile, "illum 1\n");
    color_map.x = 1.0;
    color_map.z = 1.0 - i/64.0;
    fprintf(arrowfile, "newmtl mag%03d\n", i);
    fprintf(arrowfile, "Ka %f  %f  %f\n", color_map.x, color_map.y, color_map.z);
    fprintf(arrowfile, "Kd %f  %f  %f\n", 1.0-color_map.x, 
	    1.0-color_map.y, 1.0-color_map.z);
    fprintf(arrowfile, "illum 1\n");
  }
  for(i=64; i<128; i++)
  {
    color_map.x = 0.0;
    color_map.y = 1.0;
    color_map.z = 1.0 - (i - 64)/64.0;
    fprintf(arrowfile, "newmtl elc%03d\n", i);
    fprintf(arrowfile, "Ka %f  %f  %f\n", color_map.x, color_map.y, color_map.z);
    fprintf(arrowfile, "Kd %f  %f  %f\n", 1.0-color_map.x, 
	    1.0-color_map.y, 1.0-color_map.z);
    fprintf(arrowfile, "illum 1\n");
    color_map.x = 1.0;
    color_map.y = color_map.z;
    color_map.z = 0.0;
    fprintf(arrowfile, "newmtl mag%03d\n", i);
    fprintf(arrowfile, "Ka %f  %f  %f\n", color_map.x, color_map.y, color_map.z);
    fprintf(arrowfile, "Kd %f  %f  %f\n", 1.0-color_map.x, 
	    1.0-color_map.y, 1.0-color_map.z);
    fprintf(arrowfile, "illum 1\n");
  }
  fclose(arrowfile);

/*  read in magnetic field data  */

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

/*  read in electric field data  */

  printf("name of electric field file data? ");
  scanf( "%s", filename);
  readmag = fopen(filename, "rb");
  if( !readmag)
  {
    printf("can't read in electric data file %s.\n", filename);
    exit(0);
  }
  printf("opening data files and reading them in.\n");
  numread = fread(elc_field, sizeof(VECTOR), NUMELMNTS, readmag);
  if(numread < NUMELMNTS)
    printf("\nWARNING: not all field file read in. \n\n");
  fclose(readmag);

  init_arrow_def();

//  dr = MAXWALL/MAXDIM;

/*  Plot each vector as an arrow. */
/*  first find maximum vector length for electric and magnetic field vectors  */

  elcvectormax = 0.0;
  magvectormax = 0.0;
  numpoints = 0;
  for(i=0; i<MAXDIM; i++)
  {
    ix = i*(i+1)*(i+2)/6;
    for(j=0; j<=i; j++)
    {
      iy = j*(j+1)/2;
      for(k=0; k<=j; k++)
      {
	iz = ix + iy + k;
	bx = mag_field[iz].x;
	by = mag_field[iz].y;
	bz = mag_field[iz].z;
	vtst = sqrt(bx*bx + by*by + bz*bz);
	magmag[iz] = vtst;
	if (vtst > magvectormax) magvectormax = vtst;
	bx = elc_field[iz].x;
	by = elc_field[iz].y;
	bz = elc_field[iz].z;
	vtst = sqrt(bx*bx + by*by + bz*bz);
	elcmag[iz] = vtst;
	if (vtst > elcvectormax) elcvectormax = vtst;
	numpoints++;
      }
    }
  }
  across = M_PI/4.0;
  up = atan(1.0/sqrt(2.0));
  scalevector = 6.0;
  printf("found %d points\n", numpoints);
  printf("name of picture file: ");
  scanf( "%s", filename);
  arrowfile = fopen(filename, "w");
  if( !arrowfile)
  {
    printf("can't create in picture data file %s.\n", filename);
    exit(0);
  }
  fprintf(arrowfile, "mtllib  electromag_color\n");
  for(i=0; i<MAXDIM; i++)
  {
    r = i*scalevector;
    ix = i*(i+1)*(i+2)/6;
    for(j=0; j<=i; j++)
    {
      iy = j*(j+1)/2;
      if(i) phi = across/i * j;
      else phi = 0.0;
      for(k=0; k<=j; k++)
      {
	iz = ix + iy + k;
	if(i) theta = up/i * k;
	else theta = 0.0;
	plot_loc.x = r*cos(theta)*cos(phi);  //  ux
	plot_loc.y = r*sin(theta)*cos(phi);  //  uy
	plot_loc.z = r*sin(phi);             //  uz
	vtst = magmag[iz];
	plot_dir.x = mag_field[iz].x/vtst;
	plot_dir.y = mag_field[iz].y/vtst;
	plot_dir.z = mag_field[iz].z/vtst;
	index = vtst/magvectormax*128;
	draw_arrow( arrowfile, plot_dir, plot_loc, index, 0);
	vtst = elcmag[iz];
	plot_dir.x = elc_field[iz].x/vtst;
	plot_dir.y = elc_field[iz].y/vtst;
	plot_dir.z = elc_field[iz].z/vtst;
	index = vtst/elcvectormax*128;
	draw_arrow( arrowfile, plot_dir, plot_loc, index, 1);
      }
    }
  }
  fclose(arrowfile);
  printf("all done!\n");
}