/* plot_ucd.c  Rotating, wire frame or with lighting */
/*             UCD format   3D boundary file  *.inp  */
/*             link with opengl(mesa) and glut       */

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

static FILE * infile;    /* input file name on command line */
static FILE * outfile;   /* optional, second, output file on command line */
static int num_pts;
static int num_pys;
typedef struct {GLfloat x; GLfloat y; GLfloat z;} dpts;
typedef struct {GLfloat r; GLfloat g; GLfloat b; GLfloat a;} rgba;
typedef struct {int n1; int n2; 
                GLfloat r; GLfloat g; GLfloat b; GLfloat a;} wcol;
static dpts * data_points;
static rgba * color_points;
static int    num_pts_col;
static int  * data_verts;
static rgba * color_verts;
static int    num_pys_col;
static wcol * color_wire;
static int    num_wire_col;
static int debug = 0;
static GLfloat size = 1.0;     /* auto size data */
static GLfloat scale = 1.0;    /* user scale for plot */
static GLfloat roll = 0.0;     /* rotate about Z axis 'r' */
static GLfloat pitch = 0.0;    /* rotate about X axis 'p' */
static GLfloat heading = 0.0;  /* rotate about Y axis 'h' */
static GLfloat pos[3] = {0.0, 0.0, 0.0};
static int background = 1;     /* background color 'b' */
static int wire = 1;           /* draw wireframe 'w' */
static int vert = 0;           /* draw vertices 'v' */
static int colorz = 0;         /* color based on z 'c' */
static GLfloat zavg = 0.0;
static GLfloat zmin = 0.0;
static GLfloat zmax = 0.0;

static GLfloat amb[4] = { 0.329412, 0.223529, 0.027451, 1.0}; /* brass */
static GLfloat dif[4] = { 0.780392, 0.568627, 0.113725, 1.0}; /* material */
static GLfloat spe[4] = { 0.992157, 0.941176, 0.807843, 1.0};
static GLfloat shiney = 0.21794872;

                                                  /* light source    */
static GLfloat ambient[] = {0.7, 0.0, 0.0, 1.0};  /* red   ambient    */
static GLfloat diffuse[] = {0.5, 0.4, 0.4, 1.0};  /* mid    diffuse  */
static GLfloat specular[] = {0.4, 0.4, 0.5, 1.0}; /* mid    specular */
static GLfloat position[] = {0.0, 0.0, 3.0, 0.0};

static GLdouble m[16];   /* 4 by 4 transformation matrix */
static int win_w, win_h;

static void draw_solid_dat(void)
{
  int i, j, k, pts, pt;
  GLfloat v[3][3];
  int kk[3];
  GLfloat ax, ay, az, bx, by, bz, nx, ny, nz, s;
  
  glLightfv(GL_LIGHT0, GL_AMBIENT, ambient);
  glLightfv(GL_LIGHT0, GL_DIFFUSE, diffuse);
  glLightfv(GL_LIGHT0, GL_POSITION, position);

  glFrontFace(GL_CCW);
  glEnable(GL_LIGHTING);
  glEnable(GL_LIGHT0);
  glEnable(GL_NORMALIZE);
  glEnable(GL_AUTO_NORMAL);
  glEnable(GL_DEPTH_TEST); 
  glLightfv(GL_LIGHT0, GL_POSITION, position);
  glPushMatrix();
    glLoadIdentity();
    glRotatef(pitch, 1.0, 0.0, 0.0);
    glRotatef(heading, 0.0, 1.0, 0.0);
    glRotatef(roll, 0.0, 0.0, 1.0);
    glTranslatef(pos[0], pos[1], pos[2]);
    glGetDoublev(GL_MODELVIEW_MATRIX, m); /* save transformation */

    glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, amb);
    /*
    glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, dif);
    glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, spe);
    glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, shiney * 64.0);
    */

    /* use pys to do surfaces */
    k = 0;
    for(i=0; i<num_pys; i++)
    {
      pts = data_verts[k++];
      if(pts==1) {k++; continue;} /* must have at lease 3 points */
      if(pts==2) {k += 2; continue;}
      if(debug>2) printf("pts=%d\n", pts);
      if(num_pys_col>0) 
      {
        amb[0] = color_verts[i].r;
        amb[1] = color_verts[i].g;
	amb[2] = color_verts[i].b;
        amb[3] = color_verts[i].a;
        glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, amb);
      }
      glBegin(GL_POLYGON);
        for(j=0; j<3; j++)
        {
          pt = data_verts[k++];
          if(debug>2) printf("j=%d, pt=%d\n", j, pt);
          kk[j] = pt;
          v[j][0] = data_points[pt-1].x;
          v[j][1] = data_points[pt-1].y;
          v[j][2] = data_points[pt-1].z;
        }
        /* compute normals, assuming planar, normalize */
        ax = v[2][0] - v[1][0];
        ay = v[2][1] - v[1][1];
        az = v[2][2] - v[1][2];
        bx = v[1][0] - v[0][0];
        by = v[1][1] - v[0][1];
        bz = v[1][2] - v[0][2];
        nx = ay*bz-az*by;
        ny = az*bx-ax*bz;
        nz = ax*by-ay*bx;
        s = sqrt(nx*nx+ny*ny+nz*nz);
        nx = nx / s;
        ny = ny / s;
        nz = nz / s;
        if(debug>1)
        {
          printf("\nv0x=%f, v0y=%f, v0z=%f\n", v[0][0], v[0][1], v[0][2]);
          printf("v1x=%f, v1y=%f, v1z=%f\n", v[1][0], v[1][1], v[1][2]);
          printf("v2x=%f, v2y=%f, v2z=%f\n", v[2][0], v[2][1], v[2][2]);
          printf("nx=%f, ny=%f, nz=%f\n", nx, ny, nz); 
        }
        for(j=0; j<3; j++)
        {
         glNormal3f(nx, ny, nz);
         glVertex3f(scale*data_points[kk[j]-1].x,
                    scale*data_points[kk[j]-1].y,
                    scale*data_points[kk[j]-1].z);
        }
        for(j=3; j<pts; j++)
        {
          pt = data_verts[k++];
          if(debug>2) printf("j=%d, pt=%d\n", j, pt);
          glNormal3f(nx, ny, nz);
          glVertex3f(scale*data_points[pt-1].x,
                     scale*data_points[pt-1].y,
                     scale*data_points[pt-1].z);
        }
      glEnd();
    } /* end loop on num_pys */
  glPopMatrix();
} /* end draw_solid_dat */

static void draw_wire_dat(void)
{
  int i, j, k, pts, pt, pt0;

  glDisable(GL_LIGHTING);
  glDisable(GL_LIGHT0);
  glDisable(GL_NORMALIZE);
  glDisable(GL_AUTO_NORMAL);
  if(background)  glColor3f(0.0, 0.0, 0.0);
  if(!background) glColor3f(1.0, 1.0, 1.0);
  glPointSize(3.0);
  glPushMatrix();
    glLoadIdentity();
    glRotatef(pitch, 1.0, 0.0, 0.0);
    glRotatef(heading, 0.0, 1.0, 0.0);
    glRotatef(roll, 0.0, 0.0, 1.0);
    glTranslatef(pos[0], pos[1], pos[2]);
    glGetDoublev(GL_MODELVIEW_MATRIX, m); /* save transformation */

    /* use pys to do edges */
    k = 0;
    for(i=0; i<num_pys; i++)
    {
      pts = data_verts[k++];
      if(debug>2) printf("pts=%d\n", pts);
      if(num_pys_col>0) glColor4f(color_verts[i].r,color_verts[i].g,
				  color_verts[i].b,color_verts[i].a);
      glBegin(GL_LINE_STRIP);
        for(j=0; j<pts; j++)
        {
          pt = data_verts[k++]; /* packed data */
          if(j==0) pt0 = pt; /* close polygon */
          if(debug>2) printf("j=%d, pt=%d\n", j, pt);
          glVertex3f(scale*data_points[pt-1].x,
                     scale*data_points[pt-1].y,
                     scale*data_points[pt-1].z);
        }
        glVertex3f(scale*data_points[pt0-1].x,
                   scale*data_points[pt0-1].y,
                   scale*data_points[pt0-1].z);
      glEnd();
    } /* end loop on num_pys */
  glPopMatrix();
} /* end draw_wire_dat */

static void draw_vert_dat(void)
{
  int i, j, k, pts, pt;

  glDisable(GL_LIGHTING);
  glDisable(GL_LIGHT0);
  glDisable(GL_NORMALIZE);
  glDisable(GL_AUTO_NORMAL);
  if(background)  glColor3f(0.0, 0.0, 0.0);
  if(!background) glColor3f(1.0, 1.0, 1.0);
  glPointSize(3.0);
  glPushMatrix();
    glLoadIdentity();
    glRotatef(pitch, 1.0, 0.0, 0.0);
    glRotatef(heading, 0.0, 1.0, 0.0);
    glRotatef(roll, 0.0, 0.0, 1.0);
    glTranslatef(pos[0], pos[1], pos[2]);
    glGetDoublev(GL_MODELVIEW_MATRIX, m); /* save transformation */

    /* use pys to do surfaces */
    k = 0;
    for(i=0; i<num_pys; i++)
    {
      pts = data_verts[k++];
      if(debug>2) printf("pts=%d\n", pts);
      glBegin(GL_POINTS);
        for(j=0; j<pts; j++)
        {
          pt = data_verts[k++];
          if(debug>2) printf("j=%d, pt=%d\n", j, pt);
          glVertex3f(scale*data_points[pt-1].x,
                     scale*data_points[pt-1].y,
                     scale*data_points[pt-1].z);
        }
      glEnd();
    } /* end loop on num_pys */
  glPopMatrix();
} /* end draw_vert_dat */

static writeText(GLfloat x, GLfloat y, GLfloat z, GLfloat size, char * msg)
{            /* uses background and wire frame */
  char * p;
  GLfloat wht[] = {1.0, 1.0, 1.0, 1.0};
  GLfloat blk[] = {0.0, 0.0, 0.0, 0.0};

  glPushMatrix();
    glLoadIdentity ();
    if(background)
    {
      if(!wire)
          { 
            glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, blk); 
            glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, blk);
          }
          else
          {
            glColor3f(0.0, 0.0, 0.0);
          }
    }
    else
    {
      if(!wire)
          { 
        glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, wht); 
        glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, wht);
      }
          else
          {
            glColor3f(1.0, 1.0, 1.0);
          }
    } 
    glEnable(GL_LINE_SMOOTH);
    glEnable(GL_BLEND);
    glTranslatef(x, y, z);
    glScalef(size/1000.0, size/1000.0, 0.0);
    for(p=msg; *p; p++)
      glutStrokeCharacter(GLUT_STROKE_ROMAN, *p);
  glPopMatrix();
}

void display(void)
{
  if(background)  glClearColor(1.0, 1.0, 1.0, 1.0); 
  if(!background) glClearColor(0.0, 0.0, 0.0, 1.0); 
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

  if(vert)      draw_vert_dat();
  else if(wire) draw_wire_dat();
  else          draw_solid_dat();
  writeText(-size*0.95, -size*0.95+size/16.0, 2.0*size-1.0, size/2.5,
     "Key r,R roll, p,P pitch, h,H heading, l,L large  b toggles background"); 
  writeText(-size*0.95, -size*0.95, 2.0*size-1.0, size/2.5,
            "Key x,X y,Y z,Z to position, w for wireframe, v for vertices");  
  glFlush();
  glutSwapBuffers();
} /* end display */


static GLdouble dist(GLdouble a1x, GLdouble a1y, GLdouble a1z, /* line p1 */
                     GLdouble a2x, GLdouble a2y, GLdouble a2z, /* line p2 */
                     GLdouble a3x, GLdouble a3y, GLdouble a3z) /* point   */
{
  GLdouble d, ax, ay, az, bx, by, bz, cx, cy, cz;
  
  ax = a2x-a1x;      ay = a2y-a1y;      az = a2z-a1z;
  bx = a1x-a3x;      by = a1y-a3y;      bz = a1z-a3z;
  cx = ay*bz-az*by;  cy = ax*bz-az*bx;  cz = ax*by-ay*bx;
  d =  sqrt(cx*cx+cy*cy+cz*cz)/sqrt(ax*ax+ay*ay+az*az);
  return d;
}

void mouse(int button, int state, int x, int y) 
{
   GLint viewport[4];
   GLdouble mvmatrix[16], projmatrix[16];
   GLint realy;  /*  OpenGL y coordinate position  */
   GLdouble wx0, wy0, wz0;  /*  returned world x, y, z coords at near */
   GLdouble wx1, wy1, wz1;  /*  returned world x, y, z coords at far */
   GLdouble d, dmin;
   int i, imin;

   if(button == GLUT_LEFT_BUTTON && state == GLUT_DOWN)
   {
     glGetIntegerv (GL_VIEWPORT, viewport);
     glGetDoublev(GL_PROJECTION_MATRIX, projmatrix);
     /*  note viewport[3] is height of window in pixels  */
     realy = viewport[3] - (GLint) y - 1;
     gluUnProject((GLdouble) x, (GLdouble) realy, 0.0, /* near */
                  m, projmatrix, viewport, &wx0, &wy0, &wz0); 
     gluUnProject((GLdouble) x, (GLdouble) realy, 1.0, /* far */
                  m, projmatrix, viewport, &wx1, &wy1, &wz1); 
     imin = 0;
     for(i=0; i<num_pts; i++)
     {
       d = dist(wx0, wy0, wz0, wx1, wy1, wz1,
                data_points[i].x, data_points[i].y, data_points[i].z);
       if(i==0) dmin = d;
       if(debug) printf("distance, d to point %d is %g\n", i, d);
       if(d < dmin) { dmin = d; imin = i; }
     }
     printf("%g close to data_points[%d]\n", dmin, imin);
  }
  if(button==GLUT_MIDDLE_BUTTON && state == GLUT_DOWN) ; /* commit */
  if(button==GLUT_RIGHT_BUTTON && state == GLUT_DOWN) ;  /* drag */
}

void myReshape(int w, int h)
{
  win_w = w;
  win_h = h;
  glViewport(0, 0, w, h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  if(w <= h)
      glOrtho(-size, size, -size * (GLfloat)h / (GLfloat)w,
              size * (GLfloat)h / (GLfloat)w, -5.0*size, 5.0*size);
  else
      glOrtho(-size * (GLfloat)w / (GLfloat)h,
              size * (GLfloat)w / (GLfloat)h, -size, size, -5.0*size, 5.0*size);
  glMatrixMode(GL_MODELVIEW);
}

/* x,X,y,Y,z and Z keys for moving position
   r,R,p,P,h and H keys for rool, pitch and heading
   w toggle wire frame, v toggle vertices
   b toggle background, c toggle color in z */
void keyboard (unsigned char key, int x, int y)
{
  if(key == 'x') pos[0] -= 0.125;
  if(key == 'X') pos[0] += 0.125;
  if(key == 'y') pos[1] -= 0.125;
  if(key == 'Y') pos[1] += 0.125;
  if(key == 'z') pos[2] -= 0.125;
  if(key == 'Z') pos[2] += 0.125;
  if(key == 'r') roll    -= 5.0;
  if(key == 'R') roll    += 5.0;
  if(key == 'p') pitch   -= 5.0;
  if(key == 'P') pitch   += 5.0;
  if(key == 'h') heading -= 5.0;
  if(key == 'H') heading += 5.0;
  if(key == 'l') scale *= 0.7071;
  if(key == 'L') scale *= 1.414;
  if(key == 'b' || key == 'B') background = 1-background;
  if(key == 'w' || key == 'W') wire = 1-wire;
  if(key == 'v' || key == 'V') vert = 1-vert;
  if(key == 'c' || key == 'C') colorz = 1-colorz;
  glutPostRedisplay();
} /* end keyboard */

#undef abs
#define abs(x) (((x)<0.0)?(-(x)):(x))

static void ucd_read(int argc, char *argv[])
{
  int i, j, k, pts, pt;
  int kk, len;
  char line[255];
  char *stat = NULL;
  float x,y,z;
  int icell, matl;
  char typ[8];
  /*  cell type name, integer number, number of vertices */
  int cvert[8];
  char *ctyps[9] = {"none", "pt",  "line",  "tri", "quad",
                   "tet",  "pyr", "prism", "hex"};
  int  ctypi[9] = { 0,      1,     2,       3,     4,
                    5,      6,     7,       8};
  int  ctypn[9] = { 0,      1,     2,       3,     4,
                    4,      5,     6,       8};
  int nodecomp, noden, cellcomp, celln, modelcomp, modeln;
  int p1, p2, p3, p4, p5, p6, p7, p8;
  int tnum_pys; /* may be more than one polygon per cell */

  /* get file name of data to render */
  if(argc<2)
  {
    printf("supply a file name of a UCD .inp file\n");
    exit(0);
  }
  printf("open %s for reading\n", argv[1]);
  infile = fopen(argv[1], "r");
  if(infile == NULL)
  {
    printf("could not open %s for reading\n", argv[1]);
    exit(0);
  }
  /* skip comments */
  while(1) /* read comments from UCD file */
  {  
    stat = fgets(line, 254, infile);
    if(stat==NULL) 
    {
      printf("plot_ucd, premature EOF in %s\n", argv[1]);
      exit(1);
    }
    printf("input %s",line); /* optional */
    if(line[0]!='#') break;
  } /* end while to get past comments */

  sscanf(line, "%d %d %d %d %d",
         &num_pts, &tnum_pys, &num_pts_col, &num_pys_col, &num_wire_col);

  printf("nodes=%d, cells=%d \n", num_pts, tnum_pys);
  data_points = (dpts *)malloc(sizeof(dpts)*(num_pts+2));
  data_verts  = (int *)malloc(sizeof(int)*16*(tnum_pys+2));
                       /* only room 4 * pts<=4 */
  if(data_points == NULL || data_verts == NULL)
  {
    printf("can not allocate enough space for data\n");
    exit(0);
  }
  /* used for display of vertices */
  if(num_pts_col>0) color_points = (rgba *)malloc(sizeof(rgba)*(num_pts+2));

  /* used for solid display, surface color, per "cell" */
  if(num_pys_col>0)  color_verts  = (rgba *)malloc(sizeof(rgba)*(num_pys+2));

  /* used for wire frame color */
  if(num_wire_col>0) color_wire =(wcol *)malloc(sizeof(wcol)*(num_wire_col+2));

  for(i=0; i<num_pts; i++) /* nodes required */
  {
    fscanf(infile, "%d %f %f %f\n", &k, &data_points[i].x,
           &data_points[i].y, &data_points[i].z);
    if(debug) printf("k=%d, x=%f, y=%f, z=%f\n", k, data_points[i].x,
                     data_points[i].y, data_points[i].z);
    if(abs(data_points[i].x)>size) size = abs(data_points[i].x);
    if(abs(data_points[i].y)>size) size = abs(data_points[i].y);
    if(abs(data_points[i].z)>size) size = abs(data_points[i].z);
  }
  printf("scaling to size=%g\n", size);
  size = size * 1.2; /* for spinning */
  position[2] = 1.2*size;

  /* use pys to read in cells = boundary surfaces */
  kk = 0;
  num_pys = 0;
  for(i=0; i<tnum_pys; i++) /* cells required */
  {
    stat = fgets(line, 254, infile);
    if(debug) printf("%s", line);
    sscanf(line,"%d %d %s", &icell, &matl, typ);
    /* ignored, could test matl==0 for boundary */
    k = 0;
    for(j=1; j<9; j++)
    {
      if(strcmp(ctyps[j],typ)==0)
      {
        k = ctypi[j];
        break;
      }
    }
    switch(k)
    {
    case 1: /* pt */
      num_pys += 1;
      data_verts[kk++] = 1;
      sscanf(line,"%d %d %s %d", &icell, &matl, typ, &data_verts[kk++]);
      break;
    case 2: /* line */
      num_pys += 1;
      data_verts[kk++] = 2;
      sscanf(line,"%d %d %s %d %d", &icell, &matl, typ, &data_verts[kk++],
             &data_verts[kk++]);
      break;
    case 3: /* tri */
      num_pys += 1;
      data_verts[kk++] = 3;
      sscanf(line,"%d %d %s %d %d %d", &icell, &matl, typ, &data_verts[kk++],
             &data_verts[kk++], &data_verts[kk++]);
      break;
    case 4: /* quad */
      num_pys += 1;
      data_verts[kk++] = 4;
      sscanf(line,"%d %d %s %d %d %d %d", &icell, &matl, typ, &data_verts[kk++],
             &data_verts[kk++], &data_verts[kk++], &data_verts[kk++]);
      break;
    case 5: /* tet 4 tri */
      num_pys += 4;
      sscanf(line,"%d %d %s %d %d %d %d", &icell, &matl, typ,
             &p1, &p2, &p3, &p4);
      data_verts[kk++] = 3;
      data_verts[kk++] = p1;
      data_verts[kk++] = p2;
      data_verts[kk++] = p3;
      data_verts[kk++] = 3;
      data_verts[kk++] = p1;
      data_verts[kk++] = p2;
      data_verts[kk++] = p4;
      data_verts[kk++] = 3;
      data_verts[kk++] = p2;
      data_verts[kk++] = p3;
      data_verts[kk++] = p4;
      data_verts[kk++] = 3;
      data_verts[kk++] = p3;
      data_verts[kk++] = p1;
      data_verts[kk++] = p4;
      break;
    case 6: /* pyr quad plus 4 tri*/
      num_pys += 5;
      sscanf(line,"%d %d %s %d %d %d %d %d", &icell, &matl, typ,
             &p1, &p2, &p3, &p4, &p5);
      data_verts[kk++] = 4;
      data_verts[kk++] = p1;
      data_verts[kk++] = p2;
      data_verts[kk++] = p3;
      data_verts[kk++] = p4;
      data_verts[kk++] = 3;
      data_verts[kk++] = p1;
      data_verts[kk++] = p2;
      data_verts[kk++] = p5;
      data_verts[kk++] = 3;
      data_verts[kk++] = p2;
      data_verts[kk++] = p3;
      data_verts[kk++] = p5;
      data_verts[kk++] = 3;
      data_verts[kk++] = p3;
      data_verts[kk++] = p4;
      data_verts[kk++] = p5;
      data_verts[kk++] = 3;
      data_verts[kk++] = p4;
      data_verts[kk++] = p1;
      data_verts[kk++] = p5;
      break;
    case 7: /* prism 2 tri plus 3 quad*/
      num_pys += 5;
      sscanf(line,"%d %d %s %d %d %d %d %d %d", &icell, &matl, typ,
             &p1, &p2, &p3, &p4, &p5, &p6);
      data_verts[kk++] = 3;
      data_verts[kk++] = p1;
      data_verts[kk++] = p2;
      data_verts[kk++] = p3;
      data_verts[kk++] = 3;
      data_verts[kk++] = p4;
      data_verts[kk++] = p5;
      data_verts[kk++] = p6;
      data_verts[kk++] = 4;
      data_verts[kk++] = p1;
      data_verts[kk++] = p2;
      data_verts[kk++] = p5;
      data_verts[kk++] = p4;
      data_verts[kk++] = 4;
      data_verts[kk++] = p2;
      data_verts[kk++] = p3;
      data_verts[kk++] = p6;
      data_verts[kk++] = p5;
      data_verts[kk++] = 4;
      data_verts[kk++] = p3;
      data_verts[kk++] = p1;
      data_verts[kk++] = p4;
      data_verts[kk++] = p6;
      break;
    case 8: /* hex 6 quad */
      num_pys += 6;
      sscanf(line,"%d %d %s %d %d %d %d %d %d %d %d", &icell, &matl, typ,
             &p1, &p2, &p3, &p4, &p5, &p6, &p7, &p8);
      data_verts[kk++] = 4;
      data_verts[kk++] = p1;
      data_verts[kk++] = p2;
      data_verts[kk++] = p3;
      data_verts[kk++] = p4;
      data_verts[kk++] = 4;
      data_verts[kk++] = p5;
      data_verts[kk++] = p6;
      data_verts[kk++] = p7;
      data_verts[kk++] = p8;
      data_verts[kk++] = 4;
      data_verts[kk++] = p1;
      data_verts[kk++] = p2;
      data_verts[kk++] = p6;
      data_verts[kk++] = p5;
      data_verts[kk++] = 4;
      data_verts[kk++] = p2;
      data_verts[kk++] = p3;
      data_verts[kk++] = p7;
      data_verts[kk++] = p6;
      data_verts[kk++] = 4;
      data_verts[kk++] = p3;
      data_verts[kk++] = p4;
      data_verts[kk++] = p8;
      data_verts[kk++] = p7;
      data_verts[kk++] = 4;
      data_verts[kk++] = p4;
      data_verts[kk++] = p1;
      data_verts[kk++] = p5;
      data_verts[kk++] = p8;
      break;
    default:
      printf("bad cell type %s\n", line);
    } /* end switch */
  } /* end cells */

  /* used for color display of vertices, optional */
  if(num_pts_col>0) /* color_points = (rgba *)  n_ndata*/
  {
    stat = fgets(line, 254, infile);
    if(debug) printf("%s", line);
    sscanf(line,"%d %d", &nodecomp, &noden); /* expect 1 3 or 1 4 */
    if(nodecomp!=1 || noden<3 || noden>4)
    {
      printf("UDF format not correct for color plot, ignoring\n");
      num_pts_col = 0;
      num_pys_col = 0;
      num_wire_col = 0;
      goto done;
    }
    stat = fgets(line, 254, infile); /* text */
    for(i=0; i<num_pts; i++) /* nodes required */
    {
      stat = fgets(line, 254, infile);
      if(debug) printf("%s", line);
      if(noden==3)
      {
        sscanf(line, "%d %f %f %f", &k, &color_points[i].r,
               &color_points[i].g, &color_points[i].b);
        if(debug) printf("k=%d, r=%f, g=%f, b=%f\n", k, color_points[i].r,
                         color_points[i].g, color_points[i].b);
        color_points[i].a = 1.0;
      }
      else if(noden==4)
      {
        sscanf(line, "%d %f %f %f %f", &k, &color_points[i].r,
               &color_points[i].g, &color_points[i].b, &color_points[i].a);
        if(debug) printf("k=%d, r=%f, g=%f, b=%f\n", k, color_points[i].r,
                         color_points[i].g, color_points[i].b,
                         color_points[i].a);
      }
    } /* end loop on num_pts */
  } /* end num_pts>0 */

  /* used for solid display, surface color, per "cell", optional */
  if(num_pys_col>0) /* color_verts  = (rgba *)   n_cdata*/
  {
    stat = fgets(line, 254, infile);
    if(debug) printf("%s", line);
    sscanf(line,"%d %d", &cellcomp, &celln); /* expect 1 3 or 1 4 */
    if(cellcomp!=1 || celln<3 || celln>4)
    {
      printf("UDF format not correct for color plot, ignoring\n");
      num_pys_col = 0;
      num_wire_col = 0;
      goto done;
    }
    stat = fgets(line, 254, infile); /* text */
    for(i=0; i<num_pys; i++) /* color optional */
    {
      stat = fgets(line, 254, infile);
      if(debug) printf("%s", line);
      if(celln==3)
      {
        sscanf(line, "%d %f %f %f", &k, &color_verts[i].r,
               &color_verts[i].g, &color_verts[i].b);
        if(debug) printf("k=%d, r=%f, g=%f, b=%f\n", k, color_verts[i].r,
                         color_verts[i].g, color_verts[i].b);
        color_verts[i].a = 1.0;
      }
      else if(celln==4)
      {
        sscanf(line, "%d %f %f %f %f", &k, &color_verts[i].r,
               &color_verts[i].g, &color_verts[i].b, &color_verts[i].a);
        if(debug) printf("k=%d, r=%f, g=%f, b=%f\n", k, color_verts[i].r,
                         color_verts[i].g, color_verts[i].b,
                         color_verts[i].a);
      }
    } /* end loop on num_pts */
  } /* end num_pys_col>0 */

  /* used for wire frame color, optional   n_mdata*/
  if(num_wire_col>0) /* color_wire =(wcol *) n1, n2, r, g, b, a */
  {
  }

  done: fclose(infile);
}  /* end ucdread */

static void ucd_write(int argc, char *argv[])
{
  
  /* get file name to write */
  if(argc<3)
  {
    printf("no output of UCD .inp file\n");
    return;
  }
  printf("open %s for writing\n", argv[2]);
  outfile = fopen(argv[2], "w");
  if(outfile == NULL)
  {
    printf("could not open %s for writing\n", argv[2]);
    return;
  }
  /* fprintf(outfile, "%d %d \n", num_pts, num_pys); */
  printf("not implemented yet.\n");
  /* fflush(outfile); */
  /* fclose(outfile); */
}  /* end ucd_write */

int main(int argc, char *argv[])
{
  glutInit(&argc, argv);
  /* need both double buffering and z buffer */
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
  glutInitWindowSize(600, 600);
  glutCreateWindow(argv[0]);
  glutReshapeFunc(myReshape);
  glutDisplayFunc(display);
  glutMouseFunc(mouse);       /* "mouse" setup    */
  glutKeyboardFunc(keyboard); /* "keyboard" setup */
  glEnable(GL_DEPTH_TEST);    /* Enable hidden--surface--removal */
  ucd_read(argc, argv);       /* read a UCD .inp file */
  /* ucd_write(argc, argv);*/ /* write a .inp file when second file name */
  glutMainLoop();
  return 0;
} /* end plot_ucd.c */