/* plot_pde.c                              */
/* Rotating pde solution demo with plot_pde*/


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

#define bool int /* remove for C++ */
#define false 0
#define true 1

#ifndef M_PI
#define M_PI 3.14159265
#endif

static int  winWidth, winHeight;
static float  angle = 0.0, axis[3], trans[3];
static bool  trackingMouse = false;
static bool  redrawContinue = false;
static bool   plot_pdeMove = false;
static double xmin,xmax,dxw,ymin,ymax,dyw,zmin,zmax,dzw;
static int    nx, ny, nz;
static FILE * infile;
static char   file_name[127];
static char   input_line[255];
static char   cxmin[7],cxmax[8],cymin[7],cymax[8],czmin[7],czmax[8];
static int num_pts;
typedef struct {GLfloat x; GLfloat y; GLfloat z;
                GLfloat nx; GLfloat ny; GLfloat nz;} dpts;
static dpts * data_points; /* malloc'd space for vertices */

/* Draw the pde.*/
GLfloat vertices[][3] = {
  {-1.0,-1.0,-1.0},{1.0,-1.0,-1.0}, {1.0,1.0,-1.0}, {-1.0,1.0,-1.0}, 
  {-1.0,-1.0,1.0}, {1.0,-1.0,1.0}, {1.0,1.0,1.0}, {-1.0,1.0,1.0}
};
GLfloat colors[][3] = {
  {0.0,0.0,0.0},{1.0,0.0,0.0}, {1.0,1.0,0.0}, {0.0,1.0,0.0}, 
  {0.0,0.0,1.0}, {1.0,0.0,1.0}, {1.0,1.0,1.0}, {0.0,1.0,1.0}
};

static writeText(GLfloat x, GLfloat y, GLfloat z, GLfloat size, char * msg)
{                /* floats with figure */
  char * p;

  glPushMatrix();
    glColor3f(0.0,0.0,0.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();
} /* end writeText */

static drawText(GLfloat x, GLfloat y, GLfloat z, GLfloat size, char * msg)
{                /* fixed position */
  char * p;

  glPushMatrix();
    glLoadIdentity();
    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 polygon(int a, int b, int c , int d, int face)
{
  /* draw a polygon via list of vertices */
  glBegin(GL_POLYGON);
    glColor3fv(colors[a]);
    glVertex3fv(vertices[a]);
    glColor3fv(colors[b]);
    glVertex3fv(vertices[b]);
    glColor3fv(colors[c]);
    glVertex3fv(vertices[c]);
    glColor3fv(colors[d]);
    glVertex3fv(vertices[d]);
  glEnd();
}

void colorpde(void)
{
  double size=1.5;

  writeText(xmin-0.1*dxw, ymin-0.1*dyw, zmin-0.1*dzw, size, cxmin); 
  writeText(xmax-0.1*dxw, ymin-0.1*dyw, zmin-0.1*dzw, size, cxmax); 
  writeText(xmin-0.3*dxw, ymin+0.1*dyw, zmin-0.1*dzw, size, cymin); 
  writeText(xmin-0.3*dxw, ymax-0.1*dyw, zmin-0.1*dzw, size, cymax); 
  writeText(xmax+0.1*dxw, ymin-0.1*dyw, zmin+0.1*dzw, size, czmin); 
  writeText(xmax+0.1*dxw, ymin-0.1*dyw, zmax-0.1*dzw, size, czmax); 


  /* draw axes */
  glColor3f(0.0, 0.0, 0.0);
  glLineWidth(2.0);
  glBegin(GL_LINES);
    glVertex3f(xmin, ymin, zmin);
    glVertex3f(xmax, ymin, zmin);
    glVertex3f(xmin, ymin, zmin);
    glVertex3f(xmin, ymax, zmin);
    glVertex3f(xmin, ymin, zmin);
    glVertex3f(xmin, ymin, zmax);
  glEnd();
  glLineWidth(1.0);
  glBegin(GL_LINES);
    glVertex3f(xmin, ymax, zmin);
    glVertex3f(xmax, ymax, zmin);
    glVertex3f(xmin, ymin, zmax);
    glVertex3f(xmax, ymin, zmax);
    glVertex3f(xmin, ymax, zmax);
    glVertex3f(xmax, ymax, zmax);
    glVertex3f(xmax, ymin, zmin);
    glVertex3f(xmax, ymax, zmin);
    glVertex3f(xmin, ymin, zmax);
    glVertex3f(xmin, ymax, zmax);
    glVertex3f(xmax, ymin, zmax);
    glVertex3f(xmax, ymax, zmax);
    glVertex3f(xmax, ymin, zmin);
    glVertex3f(xmax, ymin, zmax);
    glVertex3f(xmin, ymax, zmin);
    glVertex3f(xmin, ymax, zmax);
    glVertex3f(xmax, ymax, zmin);
    glVertex3f(xmax, ymax, zmax);
  glEnd();

  /* map vertices to faces */
  polygon(1,0,3,2,0);
  polygon(3,7,6,2,1);
  polygon(7,3,0,4,2);
  polygon(2,6,5,1,3);
  polygon(4,5,6,7,4);
  polygon(5,4,0,1,5);
}

/* These functions implement a simple plot_pde-like motion control. */
float lastPos[3] = {0.0F, 0.0F, 0.0F};
int curx, cury;
int startX, startY;

void plot_pde_ptov(int x, int y, int width, int height, float v[3])
{
  float d, a;

  /* project x,y onto a hemi-sphere centered within width, height */
  v[0] = (2.0F*x - width) / width;
  v[1] = (height - 2.0F*y) / height;
  d = (float) sqrt(v[0]*v[0] + v[1]*v[1]);
  v[2] = (float) cos((M_PI/2.0F) * ((d < 1.0F) ? d : 1.0F));
  a = 1.0F / (float) sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
  v[0] *= a;
  v[1] *= a;
  v[2] *= a;
}

void mouseMotion(int x, int y)
{
  float curPos[3], dx, dy, dz;

  plot_pde_ptov(x, y, winWidth, winHeight, curPos);
  if(trackingMouse)
  {
    dx = curPos[0] - lastPos[0];
    dy = curPos[1] - lastPos[1];
    dz = curPos[2] - lastPos[2];
    if (dx || dy || dz)
    {
      angle = 90.0F * sqrt(dx*dx + dy*dy + dz*dz);
      axis[0] = lastPos[1]*curPos[2] - lastPos[2]*curPos[1];
      axis[1] = lastPos[2]*curPos[0] - lastPos[0]*curPos[2];
      axis[2] = lastPos[0]*curPos[1] - lastPos[1]*curPos[0];
      lastPos[0] = curPos[0];
      lastPos[1] = curPos[1];
      lastPos[2] = curPos[2];
    }
  } 
  glutPostRedisplay();
}

void startMotion(int x, int y)
{
  trackingMouse = true;
  redrawContinue = false;
  startX = x; startY = y;
  curx = x; cury = y;
  plot_pde_ptov(x, y, winWidth, winHeight, lastPos);
  plot_pdeMove=true;
}

void stopMotion(int x, int y)
{
  trackingMouse = false;
  if (startX != x || startY != y)
  {
    redrawContinue = true;
  }
  else
  {
    angle = 0.0F;
    redrawContinue = false;
    plot_pdeMove = false;
  }
}

void display(void)
{
  double size=2.0;
  glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
  glColor3f(0.0,0.0,0.0);
  drawText(xmin-0.9*dxw, ymin-0.9*dyw, zmin-0.9*dzw, size,
           "drag mouse to rotate, right click to stop");

  /* view transform */
  if (plot_pdeMove)
  {
    glRotatef(angle, axis[0], axis[1], axis[2]);
  }
  else
  {
    glRotatef(-15.0, 1.0, 0.0, 0.0);
    glRotatef( 30.0, 0.0, 1.0, 0.0);
    glRotatef(  0.0, 0.0, 0.0, 1.0);
  }
  colorpde();
  glutSwapBuffers();
}

void mouseButton(int button, int state, int x, int y)
{
  if(button==GLUT_RIGHT_BUTTON)
  {
    trackingMouse = false;
    plot_pdeMove = false;
    angle = 0.0F;
    redrawContinue = false;
    plot_pdeMove = false;
    return;
  }
  if(button==GLUT_LEFT_BUTTON)
  {
    switch(state) 
    {
      case GLUT_DOWN:
        y=winHeight-y;
        startMotion( x,y);
        break;
      case GLUT_UP:
        stopMotion( x,y);
        break;
    }
  }
}

void myReshape(int w, int h)
{
  glViewport(0, 0, w, h);
  winWidth = w;
  winHeight = h;
}

void spinPde(void)
{
  if(redrawContinue) glutPostRedisplay();
}

void init(int argc, char *argv[])
{
  int i, j, k, ii, ipnt;
  double x, y, z, t;
  double ddx, ddy, ddz;

  /* read data file */
  /* get file name of data to render */
  if(argc<2)
  {
    printf("supply a file name of a pde data file\n");
    exit(0);
  }
  else
  {
    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);
    }
    strcpy(file_name, argv[1]);
    fgets(input_line, 254, infile);
    if(!strncmp("pde", input_line, 3))
    {
      sscanf(input_line, "pde%d %d %d", &nx, &ny, &nz);
      fgets(input_line, 254, infile);
      sscanf(input_line,"%lf %lf %lf %lf %lf %lf",
             &xmin, &xmax, &ymin, &ymax, &zmin, &zmax);
    }
    else
    {
      sscanf(input_line, "%d %d %d", &nx, &ny, &nz);
      fgets(input_line, 254, infile);
      sscanf(input_line,"%lf %lf %lf %lf %lf %lf",
             &xmin, &xmax, &ymin, &ymax, &zmin, &zmax);
    }
    printf("nx=%d, ny=%d, nz=%d\n", nx, ny, nz);
    printf("xmin=%f, xmax=%f, ymin=%f, ymax=%f, zmin=%f, zmax=%f\n",
           xmin, xmax, ymin, ymax, zmin, zmax);
    if(nz==0) /* need to handle 4D later */
    {
      ddx = (xmax-xmin)/(double)(nx-1);
      ddy = (ymax-ymin)/(double)(ny-1);
      num_pts = nx*ny;
      data_points = (dpts *)malloc(sizeof(dpts)*num_pts);
      if(data_points == NULL)
      {
        printf("can not allocate enough space for points\n");
        exit(0);
      }
      ipnt = 0;
      for(i=0; i<nx; i++)
      {
        x = xmin + ddx*i;
        for(j=0; j<ny; j++)
	{
          y = ymin + ddy*j;
          fgets(input_line, 254, infile);
          for(ii=0; ii<100; ii++) if(input_line[ii]=='=') break; 
          sscanf(&input_line[ii+1],"%lf", &z);
          printf("i=%d, j=%d, ii=%d, z=%f \n", i, j, ii, z);
          data_points[ipnt].x = x;
          data_points[ipnt].y = y;
          data_points[ipnt].z = z;
          ipnt++;
	}
      }
    } /* end nz==0 */
  }

  dxw=(xmax-xmin)/2.0;
  dyw=(ymax-ymin)/2.0;
  dzw=(zmax-zmin)/2.0;
  sprintf(cxmin,"%4.2f\0",xmin);
  sprintf(cxmax,"%4.2fX\0",xmax);
  sprintf(cymin,"%4.2f\0",ymin);
  sprintf(cymax,"%4.2fY\0",ymax);
  sprintf(czmin,"%4.2f\0",zmin);
  sprintf(czmax,"%4.2fZ\0",zmax);

  glClearColor(1.0, 1.0, 1.0, 1.0);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glOrtho(xmin-dxw, xmax+dxw, ymin-dyw, ymax+dyw, zmin-dzw, zmax+dzw);
  glMatrixMode(GL_MODELVIEW);
}


int main(int argc, char *argv[])
{
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
  glutInitWindowSize(600, 600);
  glutCreateWindow("plot_pde");
  glutReshapeFunc(myReshape);
  glutDisplayFunc(display);
  glutIdleFunc(spinPde);
  glutMouseFunc(mouseButton);
  glutMotionFunc(mouseMotion);
  glEnable(GL_DEPTH_TEST); 
  init(argc, argv);
  glutMainLoop();
  return 0;
} /* end plot_pde.c */