/* plot4d_gl.c   plot up to 4d data, mouse select choice
 *               reads from stdin, first line  4 nx ny nz nt
 *                                 many lines  x y z t u
 *                                   sorted x first, y next, ...
 *               plot4d_gl <  file.dat > file.out
 *               converted from plot4d.java
 */

#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 i4x(int i, int ii, int iii, int iiii);
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);
void repaint();
void read_data();
 
  /* mouse and display variables */
static  int xm, ym, bm; /* mouse position */
static  int height = 800;
static  int width = 800;
static  double xmin = 0.0; /* scaled for plotting */
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 dxmin, dxmax, dymin, dymax, 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 xw, yh, xoff, yoff;
static  int dim = 4; /* should be read in */
static  int nn[4]; /* nx, ny, nz, nt in var order */
static  double vmin[4]; /* xmin, ymin, zmin, tmin */
static  double vmax[4]; /* xmax, ymax, zmax, tmax */
static  char cmin[4][5] = {"xmin", "ymin", "zmin", "tmin"};
static  char cmax[4][5] = {"xmax", "ymax", "zmax", "tmax"};
static  GLfloat zop[3];
static  int var1 = 0; /* x=0, y=1, z=2, t=3 front axis mouse change */
static  int var2 = 1; /* x=0, y=1, z=2, t=3 side axis */
static  int var3 = 2; /* x=0, y=1, z=2, t=3 scan-run */
static  int var4 = 3; /* x=0, y=1, z=2, t=3 other-run */
static  double v1, v2, v3, v4;
static  int v3step = 1;
static  int v4step = 1;
static  int nx = 5; /* read in */
static  int ny = 5;
static  int nz = 5;
static  int nt = 5;
static  int nxyzt; /* nx*ny*nz*nt */
static  int nv1 = 5; /* nx, ny, nz, nt based on var1, var2, var3, var4 */
static  int nv2 = 5;
static  int nv3 = 5;
static  int nv4 = 5;
static  double * pt; /* double[nxyzt][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("plot4d_gl.c runing\n");
  read_data();
  nn[0] = nx;
  nn[1] = ny;
  nn[2] = nz;
  nn[3] = nt;
  nv1 = nn[0];
  nv2 = nn[1];
  nv3 = nn[2];
  nv4 = nn[3];
  vmin[0] = xmin; /* in case not normalized */
  vmin[1] = ymin;
  vmin[2] = zmin;
  vmin[3] = tmin;
  vmax[0] = xmax;
  vmax[1] = ymax;
  vmax[2] = zmax;
  vmax[3] = tmax;
  v3step = nn[var3]-1;
  v4step = nn[var4]-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;
  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 */
      v4step = nv4; /* optional, could go from current position */
      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))
    {
      /* step2D */
      v4step++;
      if(v4step>=nv4) v4step=0;  
      repaint();
    }
    else if(in_rect(width-110, 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>3) var3 = 0;
      while(var3==var1 || var3==var2)
      {
        var3++;
        if(var3>3) var3 = 0;
      }
      var4 = 0;
      while(var4==var1 || var4==var2 || var4==var3)
      {
        var4++;
      }
      v3step = nn[var3]-1;
      v4step = nn[var4]-1;
      repaint();
    }
    else if(in_rect(width-220, 10, 65, 20))
    {
      /* next side */
      var2++;
      if(var2>3) var2 = 0;
      if(var2==var1) var2++;
      while(var1==var3 || var2==var3)
      {
        var3++;
        if(var3>3) var3 = 0;
      }
      var4 = 0;
      while(var4==var1 || var4==var2 || var4==var3)
      {
        var4++;
      }
      v3step = nn[var3]-1;
      v4step = nn[var4]-1;
      repaint();
    }
    else if(in_rect(width-295, 10, 70, 20))
    {
      /* next front */
      var1++;
      if(var1>3) var1 = 0;
      if(var1==var2) var2++;
      if(var2>3) var2 = 0;
      while(var1==var3 || var2==var3)
      {
        var3++;
        if(var3>3) var3 = 0;
      }
      var4 = 0;
      while(var4==var1 || var4==var2 || var4==var3)
      {
        var4++;
      }
      v3step = nn[var3]-1;
      v4step = nn[var4]-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;
  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, "step2D");

  glColor3f(0.6, 0.6, 0.6);
  drawBox(width-110,70,50,20);
  glColor3f(0.0, 0.0, 0.0);
  printstring(width-105, 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-73), 300, cmin[var3]);
  printstring((int)(width-73), 650, cmax[var3]);
  printstring((int)(width-36), 300, cmin[var4]);
  printstring((int)(width-36), 650, cmax[var4]);

  printstring(250,height-20,"xmin=");
  sprintf(fstring,"%f\0", dxmin);
  printstring(290,height-20,fstring);
  printstring(400,height-20,"xmax=");
  sprintf(fstring,"%f\0", dxmax);
  printstring(440,height-20,fstring);
  printstring(250,height-40,"ymin=");
  sprintf(fstring,"%f\0", dymin);
  printstring(290,height-40,fstring);
  printstring(400,height-40,"ymax=");
  sprintf(fstring,"%f\0", dymax);
  printstring(440,height-40,fstring);
  printstring(250,height-60,"zmin=");
  sprintf(fstring,"%f\0", dzmin);
  printstring(290,height-60,fstring);
  printstring(400,height-60,"zmax=");
  sprintf(fstring,"%f\0", dzmax);
  printstring(440,height-60,fstring);
  printstring(250,height-80,"tmin=");
  sprintf(fstring,"%f\0", dtmin);
  printstring(290,height-80,fstring);
  printstring(400,height-80,"tmax=");
  sprintf(fstring,"%f\0", dtmax);
  printstring(440,height-80,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];
  v1min = vmin[var1];
  v2min = vmin[var2];
  v3min = vmin[var3];
  v4min = vmin[var4];
  v1max = vmax[var1];
  v2max = vmax[var2];
  v3max = vmax[var3];
  v4max = vmax[var4];
  /* 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 color position */
  i4 = i4x(0, 0, 0, v4step);
  v4 = pt[i4*5+var4];
  set_color((int)(ncol*((v4-v4min)/(v4max-v4min))+0.5));
  fillRect(width-31,315+v4step*(300/(nv4-1)),25,30);

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

  /* plot function or data with u value colors */
  glLineWidth(2.0); /* for data */
  iii = v3step;
  iiii = v4step;

  if(debug) printf("paint var1=%d, var2=%d, var3=%d, var4=%d \n",
                   var1, var2, var3, var4);
  for(i=0; i<nv1-1; i++) /* var1 == front */
  {
    for(ii=0; ii<nv2-1; ii++) /* var2 == side */
    {
      i4 = i4x(i, ii, iii, iiii);   /* depends on var1..4 */
      i5 = i4x(i+1, ii, iii, iiii); /* depends on var1..4 */
      i6 = i4x(i, ii+1, iii, iiii); /* depends on var1..4 */
      v1 = pt[i4*5+var1];
      v2 = pt[i4*5+var2];
      v3 = pt[i4*5+var3];
      v4 = pt[i4*5+var4];
      u1 = pt[i4*5+dim];

      if(debug) printf("i=%d, ii=%d, iii=%d, iiii=%d \n", i, ii, iii, iiii);
      if(debug) printf("i4=%d, i5=%d, i6=%d \n", i4, i5, i6);
      if(debug) printf("v1=%f, v2=%f, v3=%f, v4=%f, u1=%f \n",
                         v1, v2, v3, v4, u1);
      u2 = pt[i5*5+dim];
      if(debug) printf("v1p=%f, v2=%f, u=%f \n", pt[i5*5+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*5+var1], v2, u2));  
      glEnd();

      u2 = pt[i6*5+dim];
      if(debug) printf("v1=%f, v2p=%f, u=%f \n", v1, pt[i6*5+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*5+var2], u2));  
      glEnd();
    } /* end ii */
  } /* end i */
  for(i=0; i<nv1-1; i++)
  {
    ii = nv2-1;
    i4 = i4x(i, ii, iii, iiii);
    i5 = i4x(i+1, ii, iii, iiii);
    v1 = pt[i4*5+var1];
    v2 = v2max;
    v3 = pt[i4*5+var3];
    v4 = pt[i4*5+var4];
    u1 = pt[i4*5+dim];
    u2 = pt[i5*5+dim];
    if(debug) printf("i4=%d, i5=%d, v1=%f, v2=%f \n", i4, i5, v1, v2);
    if(debug) printf("v3=%f, v4=%f, u1=%f, u2=%f \n", v3, v4, 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*5+var1], v2, u2));  
    glEnd();
  }
  for(ii=0; ii<nv2-1; ii++)
  {
    i = nv1-1;
    i4 = i4x(i, ii, iii, iiii);
    i6 = i4x(i, ii+1, iii, iiii);
    v2 = pt[i4*5+var2];
    v1 = v1max;
    v3 = pt[i4*5+var3];
    v4 = pt[i4*5+var4];
    u1 = pt[i4*5+dim];
    u2 = pt[i6*5+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, u1=%f, u2=%f \n", v3, v4, 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*5+var2], u2));
      glEnd();
    }
    else
    {
    }
  } /* end ii loop */
  
  glFlush();  
  if(buffer) glutSwapBuffers();
  if(running)
  {
    while(now<next) now = (double)time(NULL);
  }
} /* end paint */

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

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)
{
  return i*ny*nz*nt + ii*nz*nt + iii*nt + iiii;
} /* end S */

void repaint()
{
  glutPostRedisplay();
}

void read_data()
{
  int i, n_input;

  printf("read_data running \n");
  scanf("%d %d %d %d %d", &dim, &nx, &ny, &nz, &nt);
  printf("dim=%d \n", dim);
  printf("nx=%d \n", nx);
  printf("ny=%d \n", ny);
  printf("nz=%d \n", nz);
  printf("nt=%d \n", nt);
  nxyzt = nx*ny*nz*nt;
  printf("about to allocate %d points times %d doubles, %d bytes each for %d bytes \n",
	 nxyzt, dim+1, sizeof(double), nxyzt*(dim+1)*sizeof(double));
  pt = (double *)malloc(nxyzt*(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", &pt[n_input*5+0], &pt[n_input*5+1],
         &pt[n_input*5+2], &pt[n_input*5+3], &pt[n_input*5+4]);
    if(n_input==0)
    {
      dxmin = pt[n_input*5+0];
      dxmax = pt[n_input*5+0];
      dymin = pt[n_input*5+1];
      dymax = pt[n_input*5+1];
      dzmin = pt[n_input*5+2];
      dzmax = pt[n_input*5+2];
      dtmin = pt[n_input*5+3];
      dtmax = pt[n_input*5+3];
      dumin = pt[n_input*5+4];
      dumax = pt[n_input*5+4];
      printf("x=%f, y=%f, z=%f, t=%f, u=%f \n",
              dxmin, dymin, dzmin, dtmin, dumin);
    }
    else
    {
      dxmin = min(dxmin,pt[n_input*5+0]);
      dxmax = max(dxmax,pt[n_input*5+0]);
      dymin = min(dymin,pt[n_input*5+1]);
      dymax = max(dymax,pt[n_input*5+1]);
      dzmin = min(dzmin,pt[n_input*5+2]);
      dzmax = max(dzmax,pt[n_input*5+2]);
      dtmin = min(dtmin,pt[n_input*5+3]);
      dtmax = max(dtmax,pt[n_input*5+3]);
      dumin = min(dumin,pt[n_input*5+4]);
      dumax = max(dumax,pt[n_input*5+4]);
    }
    n_input++;
    if(n_input==nxyzt) break;
  } /* end while */

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

} /* end read_data */
/* end plot4d_gl.c */