/* plot6d_gl.c   plot up to 6d data, mouse select choice
 *               reads from stdin, first line  6 nv nw nx ny nz nt
 *                                 many lines  v w x y z t u
 *                                 sorted v first, w next, ...
 *               plot6d_gl <  file.dat > file.out
 *               converted from plot5d_gl.c
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <GL/glut.h>

#undef  max
#define max(a,b) ((a>b)?(a):(b))
#undef  min
#define min(a,b) ((a<b)?(a):(b))

void reshape(int w, int h);
void do_scan();
void mouse(int button, int state, int x, int y);
int in_rect(int x, int y, int w, int h);
void paint();
int i6x(int i, int ii, int iii, int iiii, int iiiii, int iiiiii);
void set_color(int i);
void fillRect(int x, int y, int w, int h);
void drawBox(int x, int y, int w, int h);
void printstring(int x, int y, char *string);
int S(int i, int ii, int iii, int iiii, int iiiii, int iiiiii);
void repaint();
void read_data();
 
  /* mouse and display variables */
static  int xm, ym, bm; /* mouse position */
static  int height = 800;
static  int width = 900;
static  double vmin = 0.0; /* scaled for plotting */
static  double vmax = 1.0;
static  double wmin = 0.0;
static  double wmax = 1.0;
static  double xmin = 0.0;
static  double xmax = 1.0;
static  double ymin = 0.0;
static  double ymax = 1.0;
static  double zmin = 0.0;
static  double zmax = 1.0;
static  double tmin = 0.0;
static  double tmax = 1.0;
static  double umin = 0.0;
static  double umax = 1.0;
static  double scale = 400.0;
static  double xoff = 30.0;
static  double yoff = 30.0;
static  double dvmin, dvmax, dwmin, dwmax, dxmin, dxmax, dymin, dymax;
static  double dzmin, dzmax, dtmin, dtmax, dumin, dumax, duh;
static  double v1min = 0.0; /* computed based on mouse click */
static  double v1max = 1.0;
static  double v2min = 0.0;
static  double v2max = 1.0;
static  double v3min = 0.0;
static  double v3max = 1.0;
static  double v4min = 0.0;
static  double v4max = 1.0;
static  double v5min = 0.0;
static  double v5max = 1.0;
static  double v6min = 0.0;
static  double v6max = 1.0;
static  double xw, yh, xoff, yoff;
static  int dim = 6; /* should be read in */
static  int nn[6]; /* nw, nx, ny, nz, nt in var order */
static  double vamin[6]; /* wmin, xmin, ymin, zmin, tmin */
static  double vamax[6]; /* wmax, xmax, ymax, zmax, tmax */
static  char cmin[6][5] = {"vmin", "wmin", "xmin", "ymin", "zmin", "tmin"};
static  char cmax[6][5] = {"vmax", "wmax", "xmax", "ymax", "zmax", "tmax"};
static  GLfloat zop[3];
static  int var1 = 0; /* v=0,w=0,x=1,y=2,z=3,t=4 front axis mouse change */
static  int var2 = 1; /* v=0,w=0,x=1,y=2,z=2,t=4 side axis */
static  int var3 = 2; /* v=0,w=0,x=1,y=2,z=3,t=4 scan1-run */
static  int var4 = 3; /* v=0,w=0,x=1,y=2,z=3,t=4 scan2-run */
static  int var5 = 4; /* v=0,w=0,x=1,y=2,z=3,t=4 scan3-run */
static  int var6 = 5; /* v=0,w=0,x=1,y=2,z=3,t=4 scan4-run */
static  double v1, v2, v3, v4, v5, v6;
static  int v3step = 1; /* start of scan variable */
static  int v4step = 1;
static  int v5step = 1;
static  int v6step = 1;
static  int nv = 7; /* read in */
static  int nw = 7;
static  int nx = 7;
static  int ny = 7;
static  int nz = 7;
static  int nt = 7;
static  int nvwxyzt; /* nv*nw*nx*ny*nz*nt */
static  int nv1 = 7; /* nv, nw, nx, ny, nz, nt based on */
                     /* var1, var2, var3, var4, var5, var6 */
static  int nv2 = 7;
static  int nv3 = 7;
static  int nv4 = 7;
static  int nv5 = 7;
static  int nv6 = 7;
static  double * pt; /* double[nvwxyzt][dim+1]; */
static  int debug = 0;
static  int running = 0;  /* "run" dynamic */
static  int stop = 0; /* quit "run" */
static  int wireframe = 1;
static  double speed = 0.05;
static  int solid = 0; /* 0 is wireframe */
static  int buffer = 1; /* 0 is no double buffer */

int main(int argc, char *argv[])
{
  printf("plot6d_gl.c runing\n");
  read_data();
  nn[0] = nv;
  nn[1] = nw;
  nn[2] = nx;
  nn[3] = ny;
  nn[4] = nz;
  nn[5] = nt;
  nv1 = nn[0];
  nv2 = nn[1];
  nv3 = nn[2];
  nv4 = nn[3];
  nv5 = nn[4];
  nv6 = nn[5];
  vamin[0] = vmin; /* in case not normalized */
  vamin[1] = wmin;
  vamin[2] = zmin;
  vamin[3] = ymin;
  vamin[4] = zmin;
  vamin[5] = tmin;
  vamax[0] = vmax;
  vamax[1] = wmax;
  vamax[2] = xmax;
  vamax[3] = ymax;
  vamax[4] = zmax;
  vamax[5] = tmax;
  v3step = nn[var3]-1;
  v4step = nn[var4]-1;
  v5step = nn[var5]-1;
  v6step = nn[var6]-1;
  scale = 0.5*height;

  glutInit(&argc, argv);
  if(argc>1)
  {
    if(argv[1][0]=='s') solid = 1;
    if(argv[1][0]=='b') buffer = 1;
    if(argv[1][0]=='a') { solid = 1; buffer = 1;}
  }
  if(solid && buffer)
     glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH );
  if(solid && !buffer)
     glutInitDisplayMode(GLUT_RGB | GLUT_SINGLE | GLUT_DEPTH );
  if(!solid & buffer) glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE );
  if(!solid & !buffer) glutInitDisplayMode(GLUT_RGB | GLUT_SINGLE );
  glutInitWindowSize(width, height);
  glutInitWindowPosition(10,10); 
  glutCreateWindow(argv[0]);
  glutReshapeFunc(reshape);
  glutDisplayFunc(paint);
  glutMouseFunc(mouse);
  glutIdleFunc(do_scan);
  if(solid) glEnable(GL_DEPTH_TEST);
  glViewport(0, 0, width, height);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glOrtho(0.0, (GLfloat)width, 0.0, (GLfloat)height,
          -0.1, (GLfloat)height+0.1);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
  glClearColor(1.0, 1.0, 1.0, 0.0);
  glutMainLoop();
  return 0;
} /* end main */

void reshape(int w, int h)
{
  glViewport(0, 0, w, h);
  width = w;
  height = h;
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glOrtho(0.0, (GLfloat)width, 0.0, (GLfloat)height,
	  -0.1, (GLfloat)height+0.1);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
  glutPostRedisplay();
}

void do_scan()
{
  if(running==0) return;
  if(stop==1) return;
  v6step++;
  if(v6step>=nv6)
  {
    v6step = 0;
    v5step++;
    if(v5step>=nv5)
    {
      v5step=0;
      v4step++;
      if(v4step>=nv4)
      {
        v4step = 0;
        v3step++;
        if(v3step>=nv3)
	{
          v3step=0;
	}
      }
    }
  }
  repaint();
} /* end do_scan */

void mouse(int button, int state, int x, int y)
{
  xm = x;
  ym = height-y;
  bm = 1;  
  if(state==GLUT_DOWN)
  {
    /* printf("mouse down x=%d, y=%d, ym=%d \n", x, y, ym); */

    if(in_rect(width-75, 10, 30, 20))
    {
      /* dynamic display do_scan */
      v6step = nv6; /* optional, could go from current position */
      v5step = nv5;
      v4step = nv4;
      v3step = nv3;
      running = 1;
      stop = 0;
    }
    else if(in_rect(width-75, 40, 30, 20))
    {
      /* stop */
      stop = 1;
      running = 0;
      printf("check stop \n");  
    }
    else if(in_rect(width-120, 40, 50, 20))
    {
      /* faster */
      speed -= 0.01;  
      printf("faster speed = %f \n",speed);
      if(speed<0) speed = 0.0;
    }
    else if(in_rect(width-180, 40, 50, 20))
    {
      /* slower */
      speed += 0.01;
      printf("slower speed = %f \n",speed);
      if(speed>250) speed = 0.3;  
    }
    else if(in_rect(width-55, 70, 50, 20))
    {
      /* step4D */
      v6step++;
      if(v6step>=nv6) v6step=0;  
      repaint();
    }
    else if(in_rect(width-110, 70, 50, 20))
    {
      /* step3D */
      v5step++;
      if(v5step>=nv5) v5step=0;  
      repaint();
    }
    else if(in_rect(width-165, 70, 50, 20))
    {
      /* step2D */
      v4step++;
      if(v4step>=nv4) v4step=0;  
      repaint();
    }
    else if(in_rect(width-220, 70, 50, 20))
    {
      /* step1D */
      v3step++;
      if(v3step>=nv3) v3step=0;  
      repaint();
    }
    else if(in_rect(width-140, 10, 60, 20))
    {
      /* next scan */
      var3++;
      if(var3>dim-1) var3 = 0;
      while(var3==var1 || var3==var2)
      {
        var3++;
        if(var3>dim-1) var3 = 0;
      }
      var4 = 0;
      while(var4==var1 || var4==var2 || var4==var3)
      {
        var4++;
      }
      var5 = 0;
      while(var5==var1 || var5==var2 || var5==var3 || var5==var4)
      {
        var5++;
      }
      var6 = 0;
      while(var6==var1 || var6==var2 || var6==var3 || var6==var4 || var6==var5)
      {
        var6++;
      }
      v3step = nn[var3]-1;
      v4step = nn[var4]-1;
      v5step = nn[var5]-1;
      v6step = nn[var6]-1;
      repaint();
    }
    else if(in_rect(width-220, 10, 65, 20))
    {
      /* next side */
      var2++;
      if(var2>dim-1) var2 = 0;
      if(var2==var1) var2++;
      while(var1==var3 || var2==var3)
      {
        var3++;
        if(var3>dim-1) var3 = 0;
      }
      var4 = 0;
      while(var4==var1 || var4==var2 || var4==var3)
      {
        var4++;
      }
      var5 = 0;
      while(var5==var1 || var5==var2 || var5==var3 || var5==var4)
      {
        var5++;
      }
      var6 = 0;
      while(var6==var1 || var6==var2 || var6==var3 || var6==var4 || var6==var5)
      {
        var6++;
      }
      v3step = nn[var3]-1;
      v4step = nn[var4]-1;
      v5step = nn[var5]-1;
      v6step = nn[var6]-1;
      repaint();
    }
    else if(in_rect(width-295, 10, 70, 20))
    {
      /* next front */
      var1++;
      if(var1>dim-1) var1 = 0;
      if(var1==var2) var2++;
      if(var2>dim-1) var2 = 0;
      while(var1==var3 || var2==var3)
      {
        var3++;
        if(var3>dim-1) var3 = 0;
      }
      var4 = 0;
      while(var4==var1 || var4==var2 || var4==var3)
      {
        var4++;
      }
      var5 = 0;
      while(var5==var1 || var5==var2 || var5==var3 || var5==var4)
      {
        var5++;
      }
      var6 = 0;
      while(var6==var1 || var6==var2 || var6==var3 || var6==var4 || var6==var5)
      {
        var6++;
      }
      v3step = nn[var3]-1;
      v4step = nn[var4]-1;
      v5step = nn[var5]-1;
      v6step = nn[var6]-1;
      repaint();
    }
    if(in_rect(width-35, 10, 30, 20))
    {
      /* quit  run or job */
      if(running) stop=1;
      else exit(0);
    }
  } /* end if down */
} /* end mousePressed */ 

GLfloat * Op(double x, double y, double u) /* Ortho plot */
{
  zop[0] = xoff+x*scale+0.5*y*scale;
  if(zop[0]>0.9*width) zop[0] = 0.9*width;
  zop[1] = yoff+u*scale+0.5*y*scale;
  if(zop[0]>0.9*height) zop[0] = 0.9*height;
  zop[2] = y*scale;
  if(!solid) zop[2] = 0.0;
  if(debug) printf("Op returns %f, %f, %f \n", zop[0], zop[1], zop[2]);
  return zop;  
} /* end Op */

int in_rect(int x, int y, int w, int h)
{
  return xm>x && xm<x+w && ym>y && ym<y+h;
} /* end in_rect */

void paint()
{
  double z, t, u; 
  double u1, u2, hu, dhu;
  double v1, v2, v3, v4;
  int i4, i5, i6;
  int i, ii, iii, iiii, iiiii, iiiiii;
  int hgt = 210;
  int ncol = 10;
  double now, next;
  char fstring[60];

  dhu = (dumax-dumin)/(double)(ncol-1);
  now = (double)time(NULL);
  next = now+speed;
  if(solid)  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  if(!solid) glClear(GL_COLOR_BUFFER_BIT );
  glLoadIdentity();

  /* glRotatef(30.0, 0.0, 1.0, 0.0); */
  glLineWidth(5.0); /* for box */

  /* buttons */
  glColor3f(1.0, 0.0, 0.0);
  drawBox(width-35,10,30,20);
  glColor3f(0.0, 0.0, 0.0);
  printstring(width-30, 15, "quit");

  glColor3f(0.0, 1.0, 0.0);
  drawBox(width-75,10,30,20);
  glColor3f(0.0, 0.0, 0.0);
  printstring(width-70, 15, "run");

  glColor3f(1.0, 0.0, 0.0);
  drawBox(width-75,40,30,20);
  glColor3f(0.0, 0.0, 0.0);
  printstring(width-70, 45, "stop");

  glColor3f(0.0, 1.0, 0.0);
  drawBox(width-130,40,45,20);
  glColor3f(0.0, 0.0, 0.0);
  printstring(width-120, 45, "faster");

  glColor3f(0.0, 0.0, 1.0);
  drawBox(width-190,40,50,20);
  glColor3f(0.0, 0.0, 0.0);
  printstring(width-180, 45, "slower");

  glColor3f(0.5, 0.5, 0.5);
  drawBox(width-55,70,50,20);
  glColor3f(0.0, 0.0, 0.0);
  printstring(width-50, 75, "step4D");

  glColor3f(0.5, 0.5, 0.5);
  drawBox(width-110,70,50,20);
  glColor3f(0.0, 0.0, 0.0);
  printstring(width-105, 75, "step3D");

  glColor3f(0.5, 0.5, 0.5);
  drawBox(width-165,70,50,20);
  glColor3f(0.0, 0.0, 0.0);
  printstring(width-160, 75, "step2D");

  glColor3f(0.6, 0.6, 0.6);
  drawBox(width-220,70,50,20);
  glColor3f(0.0, 0.0, 0.0);
  printstring(width-215, 75, "step1D");

  glColor3f(0.0, 0.0, 1.0);
  drawBox(width-150,10,65,20);
  glColor3f(0.0, 0.0, 0.0);
  printstring(width-145, 15, "next scan");

  glColor3f(1.0, 1.0, 0.0);
  drawBox(width-225,10,65,20);
  glColor3f(0.0, 0.0, 0.0);
  printstring(width-220, 15, "next side");

  glColor3f(1.0, 0.0, 1.0);
  drawBox(width-305,10,70,20);
  glColor3f(0.0, 0.0, 0.0);
  printstring(width-295, 15, "next front");

  /* axis labels */
  printstring(xoff-15, 15, cmin[var1]);
  printstring((int)(xoff+scale-15), 15, cmax[var1]);
  printstring((int)(xoff+scale+10), (int)yoff, cmin[var2]);
  printstring((int)(xoff+1.5*scale), (int)(yoff+0.5*scale-15), cmax[var2]);
  printstring((int)(width-201), 300, cmin[var3]);
  printstring((int)(width-201), 650, cmax[var3]);
  printstring((int)(width-146), 300, cmin[var4]);
  printstring((int)(width-146), 650, cmax[var4]);
  printstring((int)(width-91), 300, cmin[var5]);
  printstring((int)(width-91), 650, cmax[var5]);
  printstring((int)(width-36), 300, cmin[var6]);
  printstring((int)(width-36), 650, cmax[var6]);

  printstring(250,height-20,"vmin=");
  sprintf(fstring,"%f\0", dvmin);
  printstring(290,height-20,fstring);
  printstring(400,height-20,"vmax=");
  sprintf(fstring,"%f\0", dvmax);
  printstring(440,height-20,fstring);
  printstring(250,height-40,"wmin=");
  sprintf(fstring,"%f\0", dwmin);
  printstring(290,height-40,fstring);
  printstring(400,height-40,"wmax=");
  sprintf(fstring,"%f\0", dwmax);
  printstring(440,height-40,fstring);
  printstring(250,height-60,"xmin=");
  sprintf(fstring,"%f\0", dxmin);
  printstring(290,height-60,fstring);
  printstring(400,height-60,"xmax=");
  sprintf(fstring,"%f\0", dxmax);
  printstring(440,height-60,fstring);
  printstring(250,height-80,"ymin=");
  sprintf(fstring,"%f\0", dymin);
  printstring(290,height-80,fstring);
  printstring(400,height-80,"ymax=");
  sprintf(fstring,"%f\0", dymax);
  printstring(440,height-80,fstring);
  printstring(250,height-100,"zmin=");
  sprintf(fstring,"%f\0", dzmin);
  printstring(290,height-100,fstring);
  printstring(400,height-100,"zmax=");
  sprintf(fstring,"%f\0", dzmax);
  printstring(440,height-100,fstring);
  printstring(250,height-120,"tmin=");
  sprintf(fstring,"%f\0", dtmin);
  printstring(290,height-120,fstring);
  printstring(400,height-120,"tmax=");
  sprintf(fstring,"%f\0", dtmax);
  printstring(440,height-120,fstring);

  printstring((int)(xoff+1.5*scale+10), (int)(yoff+0.5*scale+10), "umin");
  printstring((int)(xoff+1.5*scale+10), (int)(yoff+1.5*scale-15), "umax");

  /* height color scale */
  for(i=0; i<ncol; i++)
  {
    set_color(i);
    fillRect(20,height-220+20*i,50,15);
  }
  glColor3f(0.0, 0.0, 0.0);
  printstring(80, height-20, "u value");
  sprintf(fstring,"%f\0", dumax);
  printstring(80, height-40, fstring);
  sprintf(fstring,"%f\0", dumin);
  printstring(80, height-220, fstring);
  for(i=1; i<ncol-1; i++)
  {
    sprintf(fstring,"%f\0", dumin+dhu*(double)(ncol-1-i));
    printstring(80, height-40-i*20, fstring);
  }



  nv1 = nn[var1];
  nv2 = nn[var2];
  nv3 = nn[var3];
  nv4 = nn[var4];
  nv5 = nn[var5];
  nv6 = nn[var6];
  v1min = vamin[var1];
  v2min = vamin[var2];
  v3min = vamin[var3];
  v4min = vamin[var4];
  v5min = vamin[var5];
  v6min = vamin[var6];
  v1max = vamax[var1];
  v2max = vamax[var2];
  v3max = vamax[var3];
  v4max = vamax[var4];
  v5max = vamax[var5];
  v6max = vamax[var6];
  /* draw box, bottom */
  set_color(0);
  glBegin(GL_LINE_LOOP);
    glVertex3fv(Op(v1min, v2min, umin));  
    glVertex3fv(Op(v1max, v2min, umin));  
    glVertex3fv(Op(v1max, v2max, umin));  
    glVertex3fv(Op(v1min, v2max, umin));    
  glEnd();

  set_color(ncol-1); /* top */
  glBegin(GL_LINE_LOOP); 
    glVertex3fv(Op(v1min, v2min, umax));  
    glVertex3fv(Op(v1max, v2min, umax));
    glVertex3fv(Op(v1max, v2max, umax));  
    glVertex3fv(Op(v1min, v2max, umax));  
    glVertex3fv(Op(v1min, v2min, umax));  
  glEnd();

  /* draw box edges with z value colors */
  u1 = umin;
  hu = (umax-umin)/(ncol);
  for(i=0; i<ncol; i++)
  {
    u2 = u1+hu;
    set_color(i);
    glBegin(GL_LINES); 
      glVertex3fv(Op(v1min, v2min, u1));  
      glVertex3fv(Op(v1min, v2min, u2));  

      glVertex3fv(Op(v1min, v2max, u1));  
      glVertex3fv(Op(v1min, v2max, u2));

      glVertex3fv(Op(v1max, v2max, u1));  
      glVertex3fv(Op(v1max, v2max, u2));  

      glVertex3fv(Op(v1max, v2min, u1));  
      glVertex3fv(Op(v1max, v2min, u2));  
    glEnd();
    u1 = u2;
  } /* end i loop */

  /* var3 and var4 and var5 and var6 color position */
  i4 = i6x(0, 0, 0, 0, 0, v6step);
  v6 = pt[i4*(dim+1)+var6];
  set_color((int)(ncol*((v6-v6min)/(v6max-v6min))+0.5));
  fillRect(width-36,315+v6step*(300/(nv6-1)),25,30);

  i4 = i6x(0, 0, 0, 0, v5step, 0);
  v5 = pt[i4*(dim+1)+var5];
  set_color((int)(ncol*((v5-v5min)/(v5max-v5min))+0.5));
  fillRect(width-91,315+v5step*(300/(nv5-1)),25,30);

  i4 = i6x(0, 0, 0, v4step, 0, 0);
  v4 = pt[i4*(dim+1)+var4];
  set_color((int)(ncol*((v4-v4min)/(v4max-v4min))+0.5));
  fillRect(width-146,315+v4step*(300/(nv4-1)),25,30);

  i4 = i6x(0, 0, v3step, 0, 0, 0);
  v3 = pt[i4*(dim+1)+var3];
  set_color((int)(ncol*((v3-v3min)/(v3max-v3min))+0.5));
  fillRect(width-201,315+v3step*(300/(nv3-1)),25,30);

  /* plot function or data with u value colors */
  glLineWidth(2.0); /* for data */
  iii = v3step;
  iiii = v4step;
  iiiii = v5step;
  iiiiii = v6step;
  if(debug) printf("paint var1=%d, var2=%d, var3=%d, var4=%d, var5=%d, var6=%d \n",
                   var1, var2, var3, var4, var5, var6);
  for(i=0; i<nv1-1; i++) /* var1 == front */
  {
    for(ii=0; ii<nv2-1; ii++) /* var2 == side */
    {
      i4 = i6x(i, ii, iii, iiii, iiiii, iiiiii);   /* depends on var1..4 */
      i5 = i6x(i+1, ii, iii, iiii, iiiii, iiiiii); /* depends on var1..4 */
      i6 = i6x(i, ii+1, iii, iiii, iiiii, iiiiii); /* depends on var1..4 */
      v1 = pt[i4*(dim+1)+var1];
      v2 = pt[i4*(dim+1)+var2];
      v3 = pt[i4*(dim+1)+var3];
      v4 = pt[i4*(dim+1)+var4];
      v5 = pt[i4*(dim+1)+var5];
      v6 = pt[i4*(dim+1)+var6];
      u1 = pt[i4*(dim+1)+dim];

      if(debug) printf("i=%d, ii=%d, iii=%d, iiii=%d, iiiii=%d, iiiiii=%d\n",
                       i, ii, iii, iiii, iiiii, iiiiii);
      if(debug) printf("i4=%d, i5=%d, i6=%d \n", i4, i5, i6);
      if(debug) printf("v1=%f, v2=%f, v3=%f, v4=%f, v5=%f, v6=%f, u1=%f \n",
		       v1, v2, v3, v4, v5, v6, u1);
      u2 = pt[i5*(dim+1)+dim];
      if(debug) printf("v1p=%f, v2=%f, u=%f \n", pt[i5*(dim+1)+var1], v2, u2);
      set_color((int)((ncol*(((u1+u2)/2.0)-umin)/(umax-umin))+0.5));
      glBegin(GL_LINES); 
        glVertex3fv(Op(v1, v2, u1));
        glVertex3fv(Op(pt[i5*(dim+1)+var1], v2, u2));  
      glEnd();

      u2 = pt[i6*(dim+1)+dim];
      if(debug) printf("v1=%f, v2p=%f, u=%f \n", v1, pt[i6*(dim+1)+var2], u2);
      set_color((int)((ncol*(((u1+u2)/2.0)-umin)/(umax-umin))+0.5));
      glBegin(GL_LINES); 
        glVertex3fv(Op(v1, v2, u1));
        glVertex3fv(Op(v1, pt[i6*(dim+1)+var2], u2));  
      glEnd();
    } /* end ii */
  } /* end i */
  for(i=0; i<nv1-1; i++)
  {
    ii = nv2-1;
    i4 = i6x(i, ii, iii, iiii, iiiii, iiiiii);
    i5 = i6x(i+1, ii, iii, iiii, iiiii, iiiiii);
    v1 = pt[i4*(dim+1)+var1];
    v2 = pt[i4*(dim+1)+var2];
    v3 = pt[i4*(dim+1)+var3];
    v4 = pt[i4*(dim+1)+var4];
    v5 = pt[i4*(dim+1)+var5];
    v6 = pt[i4*(dim+1)+var6];
    u1 = pt[i4*(dim+1)+dim];
    u2 = pt[i5*(dim+1)+dim];
    if(debug) printf("i4=%d, i5=%d, v1=%f, v2=%f \n", i4, i5, v1, v2);
    if(debug) printf("v3=%f, v4=%f, v5=%f, v6=%f, u1=%f, u2=%f \n",
                     v3, v4, v5, v6, u1, u2);
    set_color((int)((ncol*(((u1+u2)/2.0)-umin)/(umax-umin))+0.5));
    glBegin(GL_LINES); 
      glVertex3fv(Op(v1, v2, u1));  
      glVertex3fv(Op(pt[i5*(dim+1)+var1], v2, u2));  
    glEnd();
  }
  for(ii=0; ii<nv2-1; ii++)
  {
    i = nv1-1;
    i4 = i6x(i, ii, iii, iiii, iiiii, iiiiii);
    i6 = i6x(i, ii+1, iii, iiii, iiiii, iiiiii);
    v1 = v1max;
    v1 = pt[i4*(dim+1)+var1];
    v2 = pt[i4*(dim+1)+var2];
    v3 = pt[i4*(dim+1)+var3];
    v4 = pt[i4*(dim+1)+var4];
    v5 = pt[i4*(dim+1)+var5];
    v6 = pt[i4*(dim+1)+var6];
    u1 = pt[i4*(dim+1)+dim];
    u2 = pt[i6*(dim+1)+dim];
    /* color = average */
    if(debug) printf("i4=%d, i6=%d, v1=%f, v2=%f \n", i4, i6, v1, v2);
    if(debug) printf("v3=%f, v4=%f, v5=%f, u1=%f, u2=%f \n",
                     v3, v4, v5, u1, u2);
    set_color((int)((ncol*(((u1+u2)/2.0)-umin)/(umax-umin))+0.5));
    if(wireframe)
    {
      glBegin(GL_LINES); 
        glVertex3fv(Op(v1, v2, u1));
        glVertex3fv(Op(v1, pt[i6*(dim+1)+var2], u2));
      glEnd();
    }
    else
    {
    }
  } /* end ii loop */
  
  glFlush();  
  if(buffer) glutSwapBuffers();
  if(running)
  {
    while(now<next) now = (double)time(NULL);
  }
} /* end paint */

int i6x(int i, int ii, int iii, int iiii, int iiiii, int iiiiii)
        /* i is var1, ii is var2, iii is var3, iiii is var4, */
        /* iiiii is var5, iiiiii is var6 */
{
  int lin;
  int offset[] = { nw*nx*ny*nz*nt, nx*ny*nz*nt, ny*nz*nt, nz*nt, nt, 1 };
  lin = i*offset[var1] + ii*offset[var2] + iii*offset[var3] + 
	iiii*offset[var4] + iiiii*offset[var5] + iiiiii*offset[var6];
  if(lin>=nvwxyzt)
  {
    printf("index error, lin=%d >= %d \n", lin, nvwxyzt);
    lin = nvwxyzt-1;
  }
  return lin; 
} /* end i6x */

void set_color(int i)
{
  if(i<0) {glColor3f(0.0, 0.0, 0.0); return;}
  if(i>9) {glColor3f(0.7, 0.7, 0.7); return;}
  switch(i)
  {
  case 0: glColor3f(0.0, 0.0, 0.0); break;
  case 1: glColor3f(1.0, 0.5, 0.5); break; /* should be brown */
  case 2: glColor3f(1.0, 0.0, 0.0); break;
  case 3: glColor3f(1.0, 1.0, 0.0); break;
  case 4: glColor3f(0.0, 1.0, 1.0); break;
  case 5: glColor3f(0.0, 1.0, 0.0); break;
  case 6: glColor3f(0.0, 0.0, 1.0); break;
  case 7: glColor3f(0.5, 0.5, 1.0); break;
  case 8: glColor3f(1.0, 0.0, 1.0); break;
  case 9: glColor3f(0.7, 0.7, 0.7); break;
  }
} /* end set_color */

void fillRect(int x, int y, int w, int h)
{
  glBegin(GL_POLYGON);
    glVertex2f(x, y);
    glVertex2f(x, y+h);
    glVertex2f(x+w, y+h);
    glVertex2f(x+w, y);
  glEnd();
} /* end fillrect */

void drawBox(int x, int y, int w, int h)
{
  glBegin(GL_LINE_LOOP);
    glVertex2f(x-2, y-2);
    glVertex2f(x-2, y+h+2);
    glVertex2f(x+w+2, y+h+2);
    glVertex2f(x+w+2, y-2);
    glVertex2f(x-2, y-2);
  glEnd();
} /* end drawBox */

void printstring(int x, int y, char *string)
{
   int len, i;

   glRasterPos2i(x, y);
   len = (int) strlen(string);
   for (i = 0; i < len; i++)
      glutBitmapCharacter(GLUT_BITMAP_HELVETICA_12, string[i]);
} /* end printstring */

int S(int i, int ii, int iii, int iiii, int iiiii, int iiiiii)
{
  return i*nw*nx*ny*nz*nt + ii*nx*ny*nz*nt + iii*ny*nz*nt +
         iiii*nz*nt + iiiii*nt + iiiiii;
} /* end S */

void repaint()
{
  glutPostRedisplay();
}

void read_data()
{
  int i, n_input;

  printf("read_data running \n");
  scanf("%d %d %d %d %d %d %d", &dim, &nv, &nw, &nx, &ny, &nz, &nt);
  printf("dim=%d \n", dim);
  printf("nv=%d \n", nv);
  printf("nw=%d \n", nw);
  printf("nx=%d \n", nx);
  printf("ny=%d \n", ny);
  printf("nz=%d \n", nz);
  printf("nt=%d \n", nt);
  nvwxyzt = nv*nw*nx*ny*nz*nt;
  printf("about to allocate %d points times %d doubles, %d bytes each for %d bytes \n",
	 nvwxyzt, dim+1, sizeof(double), nvwxyzt*(dim+1)*sizeof(double));
  pt = (double *)malloc(nvwxyzt*(dim+1)*sizeof(double));
  printf("allocation sucessful \n");
  n_input = 0;

  /* read data points, last index varies by 1 */
  while(1)
    {
    scanf("%lf %lf %lf %lf %lf %lf %lf", &pt[n_input*(dim+1)+0],
          &pt[n_input*(dim+1)+1], &pt[n_input*(dim+1)+2],
          &pt[n_input*(dim+1)+3], &pt[n_input*(dim+1)+4],
          &pt[n_input*(dim+1)+5], &pt[n_input*(dim+1)+6]);
    if(n_input==0)
    {
      dvmin = pt[n_input*(dim+1)+0];
      dvmax = pt[n_input*(dim+1)+0];
      dwmin = pt[n_input*(dim+1)+1];
      dwmax = pt[n_input*(dim+1)+1];
      dxmin = pt[n_input*(dim+1)+2];
      dxmax = pt[n_input*(dim+1)+2];
      dymin = pt[n_input*(dim+1)+3];
      dymax = pt[n_input*(dim+1)+3];
      dzmin = pt[n_input*(dim+1)+4];
      dzmax = pt[n_input*(dim+1)+4];
      dtmin = pt[n_input*(dim+1)+5];
      dtmax = pt[n_input*(dim+1)+5];
      dumin = pt[n_input*(dim+1)+6];
      dumax = pt[n_input*(dim+1)+6];
      printf("v=%f w=%f x=%f, y=%f, z=%f, t=%f, u=%f \n",
	     dvmin, dwmin, dxmin, dymin, dzmin, dtmin, dumin);
    }
    else
    {
      dvmin = min(dvmin,pt[n_input*(dim+1)+0]);
      dvmax = max(dvmax,pt[n_input*(dim+1)+0]);
      dwmin = min(dwmin,pt[n_input*(dim+1)+1]);
      dwmax = max(dwmax,pt[n_input*(dim+1)+1]);
      dxmin = min(dxmin,pt[n_input*(dim+1)+2]);
      dxmax = max(dxmax,pt[n_input*(dim+1)+2]);
      dymin = min(dymin,pt[n_input*(dim+1)+3]);
      dymax = max(dymax,pt[n_input*(dim+1)+3]);
      dzmin = min(dzmin,pt[n_input*(dim+1)+4]);
      dzmax = max(dzmax,pt[n_input*(dim+1)+4]);
      dtmin = min(dtmin,pt[n_input*(dim+1)+5]);
      dtmax = max(dtmax,pt[n_input*(dim+1)+5]);
      dumin = min(dumin,pt[n_input*(dim+1)+6]);
      dumax = max(dumax,pt[n_input*(dim+1)+6]);
    }
    n_input++;
    if(n_input==nvwxyzt) break;
  } /* end while */

  printf("read_data finished, n_input=%d \n", n_input);
  printf("vmin=%f, vmax=%f, last=%f \n",
         dvmin, dvmax, pt[(nvwxyzt-1)*(dim+1)+0]);
  printf("wmin=%f, wmax=%f, last=%f \n",
         dwmin, dwmax, pt[(nvwxyzt-1)*(dim+1)+1]);
  printf("xmin=%f, xmax=%f, last=%f \n",
         dxmin, dxmax, pt[(nvwxyzt-1)*(dim+1)+2]);
  printf("ymin=%f, ymax=%f, last=%f \n",
         dymin, dymax, pt[(nvwxyzt-1)*(dim+1)+3]);
  printf("zmin=%f, zmax=%f, last=%f \n",
         dzmin, dzmax, pt[(nvwxyzt-1)*(dim+1)+4]);
  printf("tmin=%f, tmax=%f, last=%f \n",
         dtmin, dtmax, pt[(nvwxyzt-1)*(dim+1)+5]);
  printf("umin=%f, umax=%f, last=%f \n",
         dumin, dumax, pt[(nvwxyzt-1)*(dim+1)+6]);
    
  for(i=0; i<nvwxyzt; i++) /* normalize all in range 0.0 to 1.0 */
  {
    pt[i*(dim+1)+0] = (pt[i*(dim+1)+0]-dvmin)/(dvmax-dvmin);
    pt[i*(dim+1)+1] = (pt[i*(dim+1)+1]-dwmin)/(dwmax-dwmin);
    pt[i*(dim+1)+2] = (pt[i*(dim+1)+2]-dxmin)/(dxmax-dxmin);
    pt[i*(dim+1)+3] = (pt[i*(dim+1)+3]-dymin)/(dymax-dymin);
    pt[i*(dim+1)+4] = (pt[i*(dim+1)+4]-dzmin)/(dzmax-dzmin);
    pt[i*(dim+1)+5] = (pt[i*(dim+1)+5]-dtmin)/(dtmax-dtmin);
    pt[i*(dim+1)+6] = (pt[i*(dim+1)+6]-dumin)/(dumax-dumin);
  }
  /* all data in range 0.0 .. 1.0 for convenient plotting */
  vmin = 0.0;
  wmin = 0.0;
  xmin = 0.0;
  ymin = 0.0;
  zmin = 0.0;
  tmin = 0.0;
  umin = 0.0;
  vmax = 1.0;
  wmax = 1.0;
  xmax = 1.0;
  ymax = 1.0;
  zmax = 1.0;
  tmax = 1.0;
  umax = 1.0;

} /* end read_data */
/* end plot6d_gl.c */